%Applies Shepelevich/Paramonov procedure to obtain AD solutions for
%spheroids by using a polydispersion of spheres.
%E. Boss, UMaine, 2006

epsilon=0.33; %ratio of a to b, a 0.5 of equatorial axis, b 0.5 of rotation axis.
n=1.05;
np=0.01;
lambda=0.66*0.75; %wavelength
k=2*pi/lambda; %wavenumber

a=1;
b=a/epsilon;
V=4/3*pi*b*a^2;
if epsilon < 1 %prolate
    e=sqrt(1-(a^2/b^2));
    S=2*pi*a*(a+b*asin(e)/e);
else  %oblate
    e=sqrt(1-(b^2/a^2));
    S=pi*(2*a^2+b^2/e*log((1+e)/(1-e)));
end

alpha=k*b;

alpha_epsilon=alpha*epsilon;

if epsilon>1
QS = quad(@(x)area(x,epsilon,alpha,k),alpha,alpha_epsilon*0.999999999);       
QV = quad(@(x)volume(x,epsilon,alpha,k),alpha,alpha_epsilon*0.999999999);       
Qext=quad(@(x)C_ext(x,epsilon,alpha,k,n,np),alpha,alpha_epsilon*0.999999999)/(S/4);
Qabs=quad(@(x)C_abs(x,epsilon,alpha,k,np),alpha,alpha_epsilon*0.999999999)/(S/4);
else
QS = quad(@(x)area(x,epsilon,alpha,k),alpha,alpha_epsilon*1.000000000001);       
QV = quad(@(x)volume(x,epsilon,alpha,k),alpha,alpha_epsilon*1.00000000001);       
Qext=quad(@(x)C_ext(x,epsilon,alpha,k,n,np),alpha,alpha_epsilon*1.00000000001)/(S/4);
Qabs=quad(@(x)C_abs(x,epsilon,alpha,k,np),alpha,alpha_epsilon*1.00000000001)/(S/4);
end

QS/S
QV/V
    


