% EXAMPLE SCRIPT TO CALCULATE REFLECTIVITY OF ICE PARTICLES FROM A MEASURED PSD

% DEFINE WAVELENGTH OF INTEREST (example here is 35GHz or 8.6mm)
wavelength=8.6e-3 % IN METRES

% THIS BIT IS THE PSD DATA [SWAP THIS FOR YOUR OWN PSD DATA / MODEL] - IN M^(-4)
PSD=[131000000 61400000 18600000 54500000 17800000 61100000 17100000 11100000 27100000 42200000 10300000 48600000 3620000 7010000 2510000 1360000 334000 265000 123000 24300 7980 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% PARTICLE SIZE BINS [METRES]
D=1e-6*[75.0   125.0   175.0   225.0   275.0   325.0   375.0   437.5   512.5   587.5   662.5   750.0   850.0   950.0  1100.0  1300.0  1500.0  1700.0  2000.0  2400.0  2800.0  3200.0  3600.0  4000.0  4400.0  4800.0  5500.0  6500.0  7500.0  8500.0  9500.0 11000.0 13000.0 15000.0 17000.0 19000.0 22500.0 27500.0];
% WIDTH OF THE SIZE BINS [METRES]
dD=1e-3*[0.0250 0.0250 0.0250 0.0250 0.0250 0.0250 0.0250 0.0375 0.0375 0.0375 0.0375 0.0500 0.0500 0.0500 0.1000 0.1000 0.1000 0.1000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.5000 0.5000 0.5000 0.5000 0.5000 1.0000 1.0000 1.0000 1.0000 1.0000 2.5000 2.5000];

% THIS loops through the particle size bins and adds up the contributions to the total reflectivity
Z_total=0;
for jD=1:length(D) % 
	mass=0.0185*D(jD)^1.9; % Brown and Francis mass-size relationship
	Z_particle=radar_agg_westbrook(mass,D(jD),wavelength); % calculate reflectivity in mm^6 of an aggregate this size
	Z_contribution=Z_particle*PSD(jD)*dD(jD); % multiply by the number of particles of this size in the cloud
	Z_total=Z_total+Z_contribution; % add on to the total
end
dBZ=10*log10(Z_total);

[' The answer is...' num2str(dBZ) 'dBZ']

