function [deff_prime,deff,ave_area,d_mean]=d_prime_to_eff_table(ice_distribution);
%[d_eff_prime,d_eff,ave_area,d_mean]=d_prime_to_eff_table(ice_distribution);
% Sucessive values in the output arrays correspond to increasing values
% of the parameter (1/b) in the gamma size distribution defined in
% "ice_distribution" using power law descriptions of crystal volume and
% area as a function of particle size.
% size distribution of ice crystals
%   N(D) ~ D^alpha * exp(-b*D^g_ice)
%   where D=diameter of smallest sphere enclosing ice crystal
%   and where we shall assume that D and Dr are expressed in microns.
% Projected area of ice crystals
%   for D < Dr;    A=sigma_a*Dr^(2-delta_a1)*D^delta_a1   
%       D>= Dr;    A=sigma_a*Dr^(2-delta_a2)*D^delta_a2
% Volume of ice in crystal
%       D < Dr;    V=sigma_v*Dr^(3-delta_v1)*D^delta_v1
%       D>= Dr;    V=sigma_v*Dr^(3-delta_v2)*D^delta_v2
% ave_area= ave area of ice particle corresponding to d_eff_prime (microns^2).
% sample call:
%s=struct('alpha_ice',2,'g_ice',1,'sigma_a',1,'delta_a1',2,'delta_a2',2,'sigma_v',1,'delta_v1',3,'delta_v2',3,'h20_depol_threshold',.1,'Dr',100);[dep,de]=d_prime_to_eff_table(s);
%structure ice_distribution
%   .alpha_ice
%   .g_ice    
%   .sigma_a
%   .delta_a1  
%   .delta_a2 
%   .sigma_v
%   .delta_v1  
%   .delta_v2 
%   .h20_depol_threshold  %range [0,1]
%   .Dr       


alpha_ice=ice_distribution.alpha_ice;
g_ice=ice_distribution.g_ice;
sigma_a=ice_distribution.sigma_a;
delta_a1=ice_distribution.delta_a1;
delta_a2=ice_distribution.delta_a2;
sigma_v=ice_distribution.sigma_v;
delta_v1=ice_distribution.delta_v1;
delta_v2=ice_distribution.delta_v2;
h20_depol_threshold=ice_distribution.h20_depol_threshold;
Dr=ice_distribution.Dr;


%find maximum value of x needed to cover particles up to 1 cm.
xmax=100;
deff_prime_max=0;
while deff_prime_max<10000
[deff_prime_max,ave_area]=compute_deff_prime(alpha_ice,g_ice,sigma_a,delta_a1,delta_a2 ...
				,sigma_v,delta_v1,delta_v2,Dr,xmax);
xmax=xmax*1.2;
end
xmax

% make vector of x=(1/b)*Dr^gamma
x=logspace(1, log(xmax)/log(10), 200);

%clear x
%xmax=max(3,log10(1000^g_ice))
% x=logspace(-2, 2*xmax, 400);


[deff_prime,ave_area]=compute_deff_prime(alpha_ice,g_ice,sigma_a,delta_a1,delta_a2 ...
				,sigma_v,delta_v1,delta_v2,Dr,x);

[deff,d_mean]=compute_deff(alpha_ice,g_ice,sigma_a,delta_a1,delta_a2 ...
				,sigma_v,delta_v1,delta_v2,Dr,x);


if 0
  figure(2)
  hold off
  plot(deff_prime_raw,'+b')
  hold on
  plot(deff_raw*1.05,'+r')
  grid
end 



%deff=interp1(deff_prime_raw,deff_raw,1:500,'linear',nan);  


if 0
  figure(3)
  plot(deff_prime_raw,deff_raw,'+',1:500,deff,'r')
  ax=axis;
  axis([1 500 0 deff(500)])
  grid
end

if 0
figure(1)
max(deff_prime)
max(deff)
max(ave_area)
plot(deff_prime,deff./deff_prime,deff_prime,ave_area./((pi/4)*deff.^2))
ylabel('d_e_f_f / d_e_f_f prime')
xlabel('d_e_f_f prime')
ax=axis;
axis([1 500 0 ax(4)*1.1])
h=line([Dr Dr], [0 ax(4)*1.1]);
set(h,'color','r')
grid
end

%b_table=deff;
return
%***************************************************************




function [deff_prime,ave_area]=compute_deff_prime(alpha_ice,g_ice,sigma_a,delta_a1,delta_a2...
		,sigma_v,delta_v1,delta_v2,Dr,x);


alpha_a1_g_ice=(alpha_ice+delta_a1+1)/g_ice;
alpha_a2_g_ice=(alpha_ice+delta_a2+1)/g_ice;

t_k=(1./x)*Dr^g_ice; %coordinate transform for limit of integration

gamma_fun_a1=gamma(alpha_a1_g_ice);
inc_gamma_fun_a1=gammainc(t_k,alpha_a1_g_ice);
gamma_fun_a2=gamma(alpha_a2_g_ice);
inc_gamma_fun_a2=gammainc(t_k,alpha_a2_g_ice);

alpha_v1_g_ice=(alpha_ice+2*delta_v1+1)/g_ice;
alpha_v2_g_ice=(alpha_ice+2*delta_v2+1)/g_ice;



gamma_fun_v1=gamma(alpha_v1_g_ice);
inc_gamma_fun_v1=gammainc(t_k,alpha_v1_g_ice);
gamma_fun_v2=gamma(alpha_v2_g_ice);
inc_gamma_fun_v2=gammainc(t_k,alpha_v2_g_ice);


A=Dr^(-2*delta_v1)*inc_gamma_fun_v1.*gamma_fun_v1;
B=Dr^(-2*delta_v2)*(1-inc_gamma_fun_v2).*gamma_fun_v2;
C=Dr^(-delta_a1)*inc_gamma_fun_a1.*gamma_fun_a1;
G=Dr^(-delta_a2)*(1-inc_gamma_fun_a2).*gamma_fun_a2;

area_factor=(C.*x.^(delta_a1/g_ice)+G.*x.^(delta_a2/g_ice));

deff_prime=(A.*x.^(2*delta_v1/g_ice)+B.*x.^(2*delta_v2/g_ice))...
           ./area_factor;

deff_prime=Dr*deff_prime.^0.25*(sigma_v^2/sigma_a)^0.25;


%compute average area of particle in square microns.
ave_area=sigma_a*(pi/4)*Dr^2*area_factor./gamma((alpha_ice+1)/g_ice);



return


%*************************************************************

function [deff,d_mean]=compute_deff(alpha_ice,g_ice,sigma_a,delta_a1, delta_a2...
		,sigma_v,delta_v1,delta_v2,Dr,x);

t_k=(1./x)*Dr^g_ice;

gamma_fun_a1=gamma((alpha_ice+delta_a1+1)/g_ice);
inc_gamma_fun_a1=gammainc(t_k,(alpha_ice+delta_a1+1)/g_ice);
gamma_fun_a2=gamma((alpha_ice+delta_a2+1)/g_ice);
inc_gamma_fun_a2=gammainc(t_k,(alpha_ice+delta_a2+1)/g_ice);

gamma_fun_v1=gamma((alpha_ice+delta_v1+1)/g_ice);
inc_gamma_fun_v1=gammainc(t_k,(alpha_ice+delta_v1+1)/g_ice);
gamma_fun_v2=gamma((alpha_ice+delta_v2+1)/g_ice);
inc_gamma_fun_v2=gammainc(t_k,(alpha_ice+delta_v2+1)/g_ice);

A=Dr^(-delta_v1)*inc_gamma_fun_v1.*gamma_fun_v1;
B=Dr^(-delta_v2)*(1-inc_gamma_fun_v2).*gamma_fun_v2;
C=Dr^(-delta_a1)*inc_gamma_fun_a1*gamma_fun_a1;
G=Dr^(-delta_a2)*(1-inc_gamma_fun_a2).*gamma_fun_a2;



deff=(A.*x.^(delta_v1/g_ice)+B.*x.^(delta_v2/g_ice))...
           ./(C.*x.^(delta_a1/g_ice)+G.*x.^(delta_a2/g_ice));
deff=deff*Dr*sigma_v/sigma_a;

%for ice regions
%dmean=x.^g_ice*gamma((alpha_ice+2)/g_ice)/gamma((alpha_ice+1)/g_ice);
%corrected 28-Nov-07 ewe
d_mean=x.^(1/g_ice)*gamma((alpha_ice+2)/g_ice)/gamma((alpha_ice+1)/g_ice);


return