PMM_inc_coef.m

function [int_P1_Q1, int_P1_Q2] = PMM_inc_coef(La, b_x1, b_x2, alpha_ref, beta_ref,...
    N_basis_x, N_basis_y, N_intervals_x, N_intervals_y, Nx, Ny, nx, ny)


n_points_int = 1000;
ksi1 = linspace(-1,1,n_points_int);
ksi2 = linspace(-1,1,n_points_int);

Nmax_x = max(N_basis_x);   %maximum number of Gegenbauer polynomial on all intervals
Nmax_y = max(N_basis_y);   %maximum number of Gegenbauer polynomial on all intervals
Nmax = max(Nmax_x, Nmax_y);
C = zeros(Nmax, n_points_int);

%for m=1:Nmax
for m=1:2
    for i = 1:n_points_int
        %C(m,i) = mfun('G',m-1,La,ksi1(i));
        C(m,i) = gegenbauerC(m-1,La,ksi1(i));
    end
end


%that's where we submit matched coordinates,
%because x(k) and x(k+1) must be constant

%{
lx1=zeros(N_intervals_x,1);
for k=1:N_intervals_x
    lx1(k) = (b_x1(k+1) - b_x1(k))/2;

end
lx2=zeros(N_intervals_y,1);
for k=1:N_intervals_y
    lx2(k) = (b_x2(k+1) - b_x2(k))/2;
end
%}

N_total_x = sum(N_basis_x);  %total number of basis functions
N_total_y = sum(N_basis_y);  %total number of basis functions
N_total_x3 = N_total_x - N_intervals_x;
N_total_y3 = N_total_y - N_intervals_y;

int_inc = zeros(N_total_x3, N_total_y3);

%{
x1 = zeros(N_intervals_x,1);
x2 = zeros(N_intervals_y,1);
x = zeros(N_intervals_x,1);
y = zeros(N_intervals_y,1);
%}

%for k1 = 1:N_intervals_x   %for each interval
%for k2 = 1:N_intervals_y 
 %   for j1 = (Nx(k1)+1):(Nx(k1)+nx(k1))
  %  for j2 = (Ny(k2)+1):(Ny(k2)+ny(k2))
for k1 = 1:1   %for each interval
for k2 = 1:1 
    for j1 = (Nx(k1)+1):(Nx(k1)+2)
    for j2 = (Ny(k2)+1):(Ny(k2)+2) 
        for i1 = 1:(n_points_int-1)   %points of integration on the interval
        for i2 = 1:(n_points_int-1)
            x1 = ( (b_x1(k1+1)-b_x1(k1))*ksi1(i1)+(b_x1(k1+1)+b_x1(k1)) )/2;
            x2 = ( (b_x2(k2+1)-b_x2(k2))*ksi2(i2)+(b_x2(k2+1)+b_x2(k2)) )/2;
            %x(k1) = PMM_matched_coordinates_ellipse(x1(k1), x2(k2), P1, P2, R1, R2)
            %y(k2) = PMM_matched_coordinates_ellipse(x1(k1), x2(k2), P1, P2, R1, R2)
            x = x1;
            y = x2;
        
            num1 = j1 - Nx(k1) ;   %true number of Gegenbauer polynomial
            num2 = j2 - Ny(k2) ;            
            int_inc(j1, j2) = int_inc(j1,j2) + C(num1,i1)*C(num2,i2)*...
                exp(1j*alpha_ref*x)*exp(1j*beta_ref*y)*...
                (1-ksi1(i1)^2)^(La-0.5)*(1-ksi2(i2)^2)^(La-0.5)*...
                (ksi1(i1+1)-ksi1(i1))*(ksi2(i2+1)-ksi2(i2));          
        end
        end
    end
    end
end
end

   
%to define h=<Pi,Pj> we should start from here

p = zeros(Nmax,1);     
norm = zeros(Nmax,1);
for i=0:(Nmax-1)
    p(i+1) = gamma(i+2*La)/(gamma(2*La)*gamma(i+1));
    %p(i)=Ci i-th Gegenbauer polynomial at 1
    norm(i+1) = pi^0.5*p(i+1)*gamma(La+0.5)/(gamma(La)*(i+La));
    %<Cn,Cm> = delta(n,m)*norm(n)
end

%define h=<Pi,Pj>, integral over (1,-1), without lx, ly!
hx = zeros(N_total_x3);
hy = zeros(N_total_y3);
for k=1:N_intervals_x
    for i=(Nx(k)+1):(Nx(k)+nx(k))
    hx(i,i) = norm(i-Nx(k));
    end
end
for k=1:N_intervals_y
    for i=(Ny(k)+1):(Ny(k)+ny(k))
    hy(i,i) = norm(i-Ny(k));
    end
end

coeff = zeros(N_total_x3*N_total_y3,1);
%{
for k1 = 1:N_intervals_x   %for each interval
for k2 = 1:N_intervals_y 
        for j1 = (Nx(k1)+1):(Nx(k1)+nx(k1))
        for j2 = (Ny(k2)+1):(Ny(k2)+ny(k2))
%}
for k1 = 1:1   %for each interval
for k2 = 1:1
        for j1 = (Nx(k1)+1):(Nx(k1)+2)
        for j2 = (Ny(k2)+1):(Ny(k2)+2)
            row = j2 + (j1-1)*N_total_y3;
            coeff(row) = int_inc(j1,j2)/(hx(j1,j1)*hy(j2,j2));
        end
        end
end
end

int_P1_Q1 = coeff(1);
int_P1_Q2 = coeff(2);