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

N=101; % number of points on C
NN=20; % number of rings of Stokeslets
c=1; a=6.0; % when changing c and a, change line below and @dtds
L=quad('sqrt((1)^2*sin(x).^2+(6.0)^2*cos(x).^2)',0,pi) %length of C
rho=linspace(0,a,NN+1); %radii of rings of Stokeslets
rho(NN+1)=[]; % don't want rho=a, so discard
%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,N)) %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

% add up contribution of each ring of Stokeslets at each point of C
phi=linspace(0,2*pi,101);  % variable of integration
cp=cos(phi);
sp=sin(phi);
for i=1:NN
    for j=1:N
        r=(sqrt(x(j)^2+rho(i)^2+z(j)^2-2*x(j)*rho(i)*cp));
        r3=r.^3;
        r5=r.^5;
        M(j,i)=simpson(z(j)*(x(j)-rho(i)*cp)./r3,0,2*pi);
        M(j,NN+i)=simpson(-3*z(j)*(x(j)-rho(i)*cp)./r5,0,2*pi);
        M(N+j,i)=simpson(z(j)^2./r3+1./r,0,2*pi);
        M(N+j,NN+i)=simpson(((x(j)-rho(i)*cp).^2+(-rho(i)*sp).^2-2*z(j)^2)./r5,0,2*pi);
    end
end
        
rhs=[zeros(N,1); ones(N,1)];
A=M\rhs;
plot(M*A)
Rfit=8*pi*sum(A(1:NN))/3
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))