function [eff_diameter,num_particles,LWC,mean_diameter]=d_eff_from_d_eff_prime(eff_diameter_prime...
          ,beta_a_backscat,depol,beta_a,pparm);
%[eff_diameter,num_particles,LWC,mean_diameter]=d_eff_from_d_eff_prime(eff_diameter_prime...
%          ,beta_a_backscat,beta_a,depol,pparm);
%beta_a_backscat     = lidar backscatter cross section, 1/(m sr)
%beta_a              = lidar scattering cross section, 1/m
%depol               = lidar measured circular depolarization. 
%pparm               = structure containing the following elements(shown
%                          with sample values) 
%              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: 0.1000
%                     Dr: 100
%            alpha_water: 2
%                g_water: 1
%               p180_ice: 0.0350
%   where:
%h20_depol_threshold = linear depolarization threshold, linear_depol<threshold identified as liquid. 
%alpha_water         = user supplied size distribution shape parameter alpha for liquid
%alpha_ice           = user supplied size
%g_water             = user supplled size distribution parameter g for liquid
%g_ice               = user supplied size distribution shape parameter g for ice
%p180_ice            = p(180)/4pi = backscatter phase function for ice crystals
%                      this can be left empty and beta_a will be used
%                      directly rather than beta_a_backscat/p180_ice.
%sigma_a, delta_a    = projected area of ice crystal = sigma_a*(pi/4)*D^delta_a
%                      where D is maximum dimension of ice crystal
%                      see Donovan and van Lammeren, JGR, Vol 106,
%                      pp 27433, no D21, Nov 16 2001.
%                      for spheres, sigma_a=1, delta_a=2
%sigma_v, delta_v   = volume of ice crystal = sigma_v*(pi/6)D^delta_a
%                      for spheres sigma_v=1, delta_v=3
%eff_diameter_prime = lidar-radar effective diameter (microns)
%                   = (9*<volume^2>/(pi*<projected area>))^.25
%LWC                = liquid water content gr/m^3
%mean_diameter      = mean value of particle diameter(ie mean of max dimension)
%eff_diameter       = effective diameter (microns) 
%                   = 2*<volume>/(3*<area>)
%                     computed from lidar and radar backscatter cross sections
%                     assuming and a particle sizes given by a modified gamma
%                     distribution as presented by Deirmendjian,
%                     'Electromagnetic Scattering on Spherical
%                     Polydispersions', Elsevier, NY, 1969:
%
%                     n(D) = a * D^alpha * exp(-b*D^g)      (equation 1) 
%    
%                      Where:
%                      D    = Maximum dimension of particle
%                      a    = Num_particles*g*b^((alpha+1)/g)*gamma((alpha+1)/g)
%                      n    = number of particles per unit volume
%                      b    = parameter computed from lidar-radar signal ratio
%                      gamma= the gamma function
%     
% Assumptions:
%      - Size distribution given by equation 1
%      - cloud at a given data point is either all water or all ice
%      - when depolarization is < h20_depol_threshold, cloud is water
%      - backscatter phase function for water droplets, p(180)/4pi = 0.05
%      - particles are large compared to the lidar wavelength such that
%              the optical scattering cross section is twice the
%              projected area of the particle.
%      - particles are small compared to the radar wavelength
%      - radar attenuation is negligible
%      - the volume of ice in a particle is related to its max dimension, D, by:
%                   
%                 Volume = sigma_v * pi/6 *Dr^(3-delta_v)* D^delta_v    
%                    for water, sigma_v =1, delta_v=3 thus Volume= pi/6 * D^3
%
%      - the projected area of an ice particle is related to D by:
%
%                 Area = sigma_a *pi/4 *Dr^(2-delta_a)* D^delta_a        
%                    for water, sigma_a =1, delta_a=2, thus Area=(pi/4) * D^2
%      - because the power laws for volume and area are often different
%        for small and large particles two values can be specified for
%        delta_a and delta_v:
%           delta_a=delta_a1, delta_v=delta_v1   for D < Dr    (microns)
%           delta_a=delta_a2, delta_v=delta_v2   for D>= Dr   
%          
%      - sigma_a=the area fill fraction at D=Dr, it is the projected area of an
%        ice crystal of size D=Dr divided by (pi/4)*Dr^2.
%      - sigma_v=volume fill fraction at D=Dr, it is the volume of ice
%        within an ice particle of size D=Dr divided by (pi/6)*Dr^3.



h20_depol_threshold=pparm.h20_depol_threshold;

%convert depol threshold to work with circular depolarization
h20_depol_threshold=2*h20_depol_threshold/(1-h20_depol_threshold);

p180_ice=pparm.p180_ice;
alpha_water=pparm.alpha_water;
g_water=pparm.g_water;

density_ice=0.92; %g/cm^3
warning off MATLAB:divideByZero





%this returns vectors of d_prime and the corresponding values of d_eff
%which can be used as a lookup table and covert from microns to meters.
[d_prime_table,d_eff_table,ave_area_table,mean_diameter_table]=d_prime_to_eff_table(pparm); 
%figure(1000);
% plot(d_prime_table,mean_diameter_table);grid;xlabel('d_e_f_f prime');ylabel('d_e_f_f')
% pause
%convert from microns to meters
d_prime_table=d_prime_table/1e6;
d_eff_table=d_eff_table/1e6;
ave_area_table=ave_area_table/1e12;
mean_diameter_table=mean_diameter_table/1e6;

eff_diameter=zeros(size(eff_diameter_prime)); %create empty array

%convert eff_diamter_prime from microns to meters
eff_diameter_prime=eff_diameter_prime/1e6;

%convert d_eff_prime in those portions of the cloud dominated by ice 
%to d_eff by interpolating lookup table
[ntimes,nalts]=size(eff_diameter_prime);
reshape(eff_diameter_prime,1,ntimes*nalts);
eff_diameter=zeros(size(eff_diameter_prime)); %create empty array

eff_diameter=interp1(d_prime_table,d_eff_table,eff_diameter_prime);


reshape(eff_diameter,ntimes,nalts);
reshape(eff_diameter_prime,ntimes,nalts);


  if ~isempty(p180_ice)
    beta_a=beta_a_backscat/p180_ice;
  end
  

%compute beta_a in water clouds using a fixed value of p180.
    beta_a(depol< h20_depol_threshold)=beta_a_backscat(depol<h20_depol_threshold)/.05;
  

%set size distribution factor in water regions
c_factor=gamma((alpha_water+4)/g_water)/gamma((alpha_water+3)/g_water);
c_factor=c_factor*(gamma((alpha_water+3)/g_water)/gamma((alpha_water+7)/g_water))^.25;
eff_diameter(depol<h20_depol_threshold)=...
    c_factor*eff_diameter_prime(depol<h20_depol_threshold);
 


    %assume particles large compared to lidar lidar wavelength
    %beta_a = num_particles * 2 * <area per particle>;

    %in ice regions
     area_per_particle=zeros(size(eff_diameter_prime));
     area_per_particle(depol>=h20_depol_threshold)=...
	 interp1(d_prime_table,ave_area_table,eff_diameter_prime(depol>=h20_depol_threshold));
     
     num_particles=beta_a./(2*area_per_particle); %1/m^3
    
      
     %eff_diameter = 3*<volume per particle>/(2*<area per particle>)
     %LWC   = num_particles * <volume per particle> * density of water
     %multiply by 10^6 to convert m^3 to cm^3 yielding gr/m^3
     
     LWC = (2/3)*density_ice*num_particles .* eff_diameter .*area_per_particle *1e6;
     
     %in water regions
     area_per_particle(depol<h20_depol_threshold)=...
	 (pi/4)*(eff_diameter(depol<h20_depol_threshold)*gamma((alpha_water+3)/g_water)...
	  /(gamma((alpha_water+4)/g_water))).^2 ...
	 *gamma((alpha_water+3)/g_water)/gamma((alpha_water+1)/g_water);
     
     num_particles(depol<h20_depol_threshold)=...
	 beta_a_backscat(depol<h20_depol_threshold)...
	 ./(.05*2*area_per_particle(depol<h20_depol_threshold)); %1/m^3
     
     
     
     %eff_diameter = 3*<volume per particle>/(2*<area per particle>)
     %LWC   = num_particles * <volume per particle> * density of water
     %multiply by 10^6 to convert m^3 to cm^3 yielding gr/m^3

    
     LWC(depol<h20_depol_threshold)=(2/3)*num_particles(depol<h20_depol_threshold)... 
	 .* eff_diameter(depol<h20_depol_threshold)...
	 .*area_per_particle(depol<h20_depol_threshold) *1e6;
	 

     
     %compute mean diameter of particles at each point in image
     %ice regions of cloud
     mean_diameter=zeros(size(eff_diameter_prime));
     mean_diameter=interp1(d_prime_table,mean_diameter_table,eff_diameter_prime);
 
    
     %for water regions of cloud
     mean_dia_const=gamma((alpha_water+3)/g_water)/gamma((alpha_water+4)/g_water)...
			   *gamma((alpha_water+2)/g_water)/gamma((alpha_water+1)/g_water)
               
     mean_diameter(depol<h20_depol_threshold)=mean_dia_const*eff_diameter(depol<h20_depol_threshold);
     mean_diameter(~isreal(mean_diameter))=nan;
     mean_diameter=mean_diameter*1e6; %change units to microns

     h20_depol_threshold
    
     
     %eff_diameter(~isreal(eff_diameter))=nan;      
     eff_diameter(~isreal(eff_diameter))=-inf;
     eff_diameter=eff_diameter*1e6; %change units to microns
     num_particles=num_particles/1000; %change to number particles/litter.