function [ext,omega,g,angs,lambda]=aerfun(f,num,N)
% Aerosol optical properties  from OPAC data dabe
% Function require aermod.mat data file
% input : 
%        f - relative humidity
%        num - vector include number of aeroosl model see OPAC description
%        N - particle number density [par/cm^3]
% output : 
%         ext - aerosol extinction coefficient [1/km]
%         omega - aerosol single scatterin albedo
%         g - asymmetry parameter
%         angs - Angstrom exponent 500/1000 nm
%         lambda - wavelength in um

load aermod
RH=f;
if RH>99
   RH=99;
end   
nmodels=length(num);
x=num;
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,:);
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,:)=omega(b,:);
      datg(i,:)=g(b,:);
   else
      datext(i,:)=interp1(rh,ext(b,:),RH);  
      datsingle(i,:)=interp1(rh,omega(b,:),RH);  
      datg(i,:)=interp1(rh,g(b,:),RH);  
   end
   extmod(i)=datext(i,7)*N(i);   
end   
ext=datext'*N';
omega=(datext.*datsingle)'*N'./extnew;
g=(datext.*datsingle.*datg)'*N'./(extnew.*singlenew);
angs=log(extnew(6)/extnew(14))/log(1000/500);
