begin
integer ii,jj,T;
ii:= 30; jj:= 36; T:= 10;

begin real array BLS(1:jj,1:ii), g,h(0:T);

real procedure W(n,m);
value n,m; integer n,m;
begin integer k,p,q,r,s; real t;
    r:= 2*n-1; p:= n+m; q:= n-m; s:= 0;
    t:= if m=0 the 1 else 2;
    for k:= 2*n step -1 until 1 do
    begin
        k:= t*r; r:= r-s; s:=2-s;
        if p>1 then begin t:= t/p; p:= p-1 end else
        if q>1 then begin t:= t/q; q:= q-1 end;
    end;
    W:= sqrt(t);
end;

real procedure P(n,m,theta); value n,m,theta;
integer n,m; real theta;
begin integer i,nm; real s,t,c,c2;
    c:= cos(theta); c2:= c**2, nm:= n-m; s:= t:= 1;
    for i:= 2 step 2 until nm do
    begin
        t:= t*(i-nm-2)*(nm-i+1)/(2*n-i+1)/i;
        s:= s*c2+t;
    end;
    if nm mod 2 = 1 then s:= s*c;
    P:= s*W(n,m)*sin(theta)**m;
 end;

procedure solve(K,m,T,a,b); value m,T;
real array K,a,b; integer m,T;
begin real div,ratio; integer i,j,h;
    for i:= m step 1 until T do
    begin
        div:= K(i,i);
        for j:= m step 1 until T do
            K(i,j):= K(i,j)/div;
        a(i):= a(i)/div; b(i):= b(i)/div;
        for j:= m step 1 until T do
        if i<>j then
        begin
            ratio:= K(j,i);
            if ratio<>0 then
            begin
                for h:= m step 1 until T do
                    K(j,h):= K(j,h)-ratio*K(i,h);
                a(j):= a(j)-ratio*a(i);
                b(j):= b(j)-ratio*b(i);
            end;
        end;
    end i;
end;

procedure getK(K,m,T,rs); value m,T,rs;
real array K; integer m,T; real rs;
begin integer i,n; real m2; m2:= m**2;
    for n:= m step 1 until T do
    begin
        for i:= m step 1 until T do K(i,n):= 0;
        K(n,n):= ((n+1)*(n+2)+m2)/(2*n+3)
               + ((n-1)*n+m2)/(2*n-1)/rs**(2*n+1);
    end;
    for n:= m+2 step 1 until T do
        K(n,n-2):= -sqrt(((n-1)**2-m2)*(n**2-m2))/(2*n-1);
    for n:= T-2 step -1 until m do
        K(n.n+2):= -sqrt(((n+1)**2-m2)*((n+2)**2-m2))
                         /(2*n+3)/rs**(2*n+5);
end;

procedure gettheta(cth,sth); real array cth,sth;
begin real th; integer i;
    for i:= 1 step 1 until ii do
    begin
        th:= 1.57079637 - arcsin((ii+1-2*i)/ii);
        cth(i):= cos(th); sth(i):= sin(th);
    end;
end;

procedure getmphi(m,cmphi,smphi); value m;
integer m; real array cmphi,smphi;
begin real s,mphi; integer j;
    s:= 2*3.14159265*m/jj;
    for j:= 1 step 1 until jj do
    begin
        mphi:= (jj-j+0.5)*s;
        cmphi(j):= cos(mphi); smphi(j):= sin(mphi);
    end;
end;

procedure getx(x,n,m,sm); value n,m;
real array x,sm; integer n,m;
begin integer i,nm,k; real c,c2,s,u;
    for i:= 1 step 1 until ii do
    begin
        c:= cth(i); c2:= c**2; nm:= n-m;
        s:= u:= 1;
        for k:= 2 step 2 until nm do
        begin
            u:= u*(k-nm-2)*(nm-k+1)/(2*n+1-k)/k;
            s:= s*c2+u;
        end;
        if nm mod 2 = 1 then s:= s*c;
        x(i):= sm(i)*s*W(n,m);
    end i;
end;

procedure getab(BLS,m,T,a,b); value m,T;
real array BLS,a,b; integer m,T;
begin integer i,j,n; real sa,sb,si;
real array x,sm(1:ii), c,s(1:jj);
    for i:= 1 step 1 until ii do
        sm(i):= sth(i)**(m+1);
    getmphi(m,c,s);
    for n:= m step 1 until T do
    begin
        sa:= sb:= 0; getx(n,m,sm);
        for j:= 1 step 1 until jj do
        begin
            si:= 0;
            for i:=1 step 1 until ii do
                si:= si+BLS(j,i)*x(i);
        end;
        a(n):= sa*(2*n+1)/(ii*jj);
        b(n):= sb*(2*n+1)/(ii*jj);
    end n;
end;

procedure getgh(BLS,rs,T,g,h); value rs,T;
real array BLS,g,h; real rs; integer T;
begin integer m;
    for m:= 0 step 1 until T do
    begin integer n;
    real array a,b(m:T),K(m:T,m:T);
        getab(BLS,m,T,a,b);
        getK(K,m,T,rs);
        solve(K,m,T,a,b);
        for n:= m step 1 until T do
        begin
            g(n,m):= a(n);
            h(n,m):= b(n);
        end;
    end m;
end;

comment get BLS somehow;
getgh(BLS) rs:(2.5) order:(T) coeffs:(g,h);
comment print out g and h somehow;

end BLS and g,h;
end;