clear all
%close all

rdx = 1;
for n=[20,50,100],
delta=0.01;
a=1;
dy=a/n;
y=zeros(1,n+1);
r=zeros(1,n+1);
integrand=zeros(1,n+1);

I=zeros(1,n+1);
C=zeros(1,n+1);
Cexp=zeros(1,n+1);

ijdx=0;

noisemax=delta;
for idx=1:n+1
  y(idx)=(idx-1)*dy; %Shifts the value
  r(idx)=y(idx);
  Cr = 0.2; %2g/L %constant
  I(idx) = 2*Cr*sqrt(a.*a - y(idx).*y(idx));%+rand*noisemax; %constant
  %I(idx)=a*sqrt(a^2-y(idx)*y(idx))+y(idx)*y(idx)*log((a+sqrt(a^2-y(idx)*y(idx)))/y(idx))+rand*noisemax;%linear
  %I(idx)=1/2*sqrt(pi)*exp(-4*y(idx)^2)*erf(2*sqrt(1-y(idx)^2));%+rand*noisemax; %guassian intensity
end

%trapezoid rule
for jdx=1:n
    ijdx=0;
    for idx=jdx+2:n
        ijdx=ijdx+dy*y(idx)*I(idx)/sqrt(y(idx)^2-r(jdx)^2);
    end
    %end points for trapezoid rule
    ijdx=ijdx+(1/2)*dy*y(n+1)*I(n+1)/sqrt(y(n+1)^2-r(jdx)^2); 
    ijdx=ijdx+(1/2)*dy*y(jdx+1)*I(jdx+1)/sqrt(y(jdx+1)^2-r(jdx)^2);
    %account for region r,r+eps
    %if n < n/2,
        ijdx=ijdx+((I(jdx)+I(jdx+1))/2)*sqrt(-r(jdx)^2+(r(jdx)+dy)^2);
    %else
    %   ijdx=ijdx+((4*I(jdx)+I(jdx+1))/5)*sqrt(-r(jdx)^2+(r(jdx)+dy)^2);
    %end
    %ijdx=ijdx+I(jdx)*sqrt(-r(jdx)^2+(r(jdx)+dy)^2);
    integrand(jdx)=ijdx;
end

for idx=2:n-1
    
    %Cexp(idx) = r(idx); %linear
    %Cexp(idx)=exp(-4*r(idx)^2); %gaussian
    C(idx)=(-1/(pi*r(idx)))*((integrand(idx+1)-integrand(idx-1))/(2*dy));
end
figure(2),
hold on
plot(y,C)
rdx = rdx +1;
end

Cexp = 0.2*ones(length(C));
plot(y,Cexp)
hold off

figure(2),
xlabel('radius'), ylabel('concentration')
legend('n=20','n=50','n=100','analytic concentration')
%ylim([0 1.2])