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])