
% Interactive program to compute aerosol optical properties from OPAC
% Require aermod.mat data file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
load aermod
dent=[2.0 1.8 1.0 2.2 2.2 2.6 2.6 2.6 2.6 1.7]; % aerosol dry (RH=0%) density
% higroscopic growth factor -- k=r/r_dry for 7 relative humidity
%                               0.0 50  70  80  90  95  99%
% for sulfate  nr 10            1.0 1.4 1.6 1.7 1.9 2.3 3.3
% for organic carbon nr 2       1.0 1.2 1.3 1.4 1.6 1.9 2.5
% for black carbon-soot nr 3    1.0 1.0 1.0 1.2 1.4 1.5 1.9
% for sea salt nr 4 and 5       1.0 1.6 1.8 2.0 2.4 2.9 4.8
% for dust growth factor is 1
disp('Aerosol optical propertice')
%wav=input('Select wavelength in [nm]: ');
RH=input('Select relative humidity: ');
if RH>99
   RH=99;
end
 x=1;
 k=0;
 disp ('___________________________________')
 disp('AEROSOL MODELS')
 disp('1  *** insoluble');
 disp('2  *** watersoluble');
 disp('3  *** soot');
 disp('4  *** sea salt accululation mod');
 disp('5  *** sea salt coarse mod');
 disp('6  *** mineral nucleation mod');
 disp('7  *** accululation mod');
 disp('8  *** mineral coarse mod');
 disp('9  *** mineral transported');
 disp('10 *** sulfate droplets');
 disp('____________________________________')

while x~=0
 x0=input('Select aerosol model: ');
 if x0==0
  disp(['You select ',num2str(k),' aerosol models'])
  break
 end
 k=k+1;
 x(k)=x0;
 N(k)=input('Select number of particles in cm^3: ');
end
A=zeros(10,8);
A(1,1)=1;
A(2,:)=[2:9];
A(3,1)=14;
A(4,:)=[15:22];
A(5,:)=[23:30];
A(6,1)=12;
A(7,1)=10;
A(8,1)=11;
A(9,1)=13;
A(10,:)=[31:38];
B=A(x,:);
nmodels=length(x);
rh=[0,50,70,80,90,95,98,99];
for i=1:nmodels
   a=find(B(i,:)>0);
    b=B(i,a);
   if length(a)==1
      datext(i,:)=ext(b,:);
      datsingle(i,:)=single(b,:);
      datg(i,:)=g(b,:);
   else
      datext(i,:)=interp1(rh,ext(b,:),RH);
      datsingle(i,:)=interp1(rh,single(b,:),RH);
      datg(i,:)=interp1(rh,g(b,:),RH);
   end
end
extnew=datext'*N';
singlenew=(datext.*datsingle)'*N'./extnew;
gnew=(datext.*datsingle.*datg)'*N'./(extnew.*singlenew);

disp(['extincion for 550 nm ',num2str(extnew(7))])
disp(['single scatter albedo for 550 nm ',num2str(singlenew(7))])