% x(t)=a*sin(t)
% z(t)=c*cos(t)  parametric equations for curve C
clear

N=101; % number of points on C
NN=30; % number of basis vectors for flow
c=1; a=3.0; % when changing c and a, change line below and @dtds
L=quad('sqrt((1)^2*sin(x).^2+(3.0)^2*cos(x).^2)',0,pi) %length of C
%a=1; b=1; % when changing a and b, change line below and @dtds
%L=quad('sqrt(1^2*sin(x).^2+1^2*cos(x).^2)',0,pi) %length of C
[s,t]=ode45(@dtds,[0,L],0)  %solve for t as a function of arclength s
T=spline(s,t,linspace(0,L,101)) %parameter t for equidistant points along C
plot(a*sin(T),c*cos(T),'.') % how does it look?
x=a*sin(T);
z=c*cos(T); %Cartesian coordinates for equidistant points
r=sqrt(x.^2+z.^2);
th=acos(z./r);  % spherical polar coordinates for equidistant points
den=sqrt(c^2*sin(T).^2+a^2*cos(T).^2) % a frequently occurring denominator
for i=1:NN
    [phi1a,phi1b,phi2,phi3]=eigbas(0,i,cos(th));
    %[phi1a,phi1b,phi2,phi3]=eigbas(0,i,cos(th));
    Mat1(:,i)=[phi1a(1,:)./r.^i+phi1b(1,:)./r.^(i+2)  phi1a(2,:)./r.^i+phi1b(2,:)./r.^(i+2)]';
    Mat2(:,i)=[phi2(1,:)./r.^(i+2)  phi2(2,:)./r.^(i+2)]';
end
A=[Mat1 Mat2]\[cos(th)';-sin(th)']
%c=[Mat1 Mat2]\[cos(th)'; -sin(th)']
hold off
plot([Mat1 Mat2]*A)
hold on
plot([cos(th) -sin(th)],'r')
%plot([cos(th) -sin(th)],'r')
%plot(rMat')
%plot(thMat')
%c=[rMat' thMat']\
A(1)/2/sqrt(pi)
R=(8/3)/(2*c/(a^2-c^2)+(2*a^2-4*c^2)/(a^2-c^2)^(3/2)*atan(sqrt(a^2-c^2)/c))