function [tau,omega,P11,G]=MIXoptics(lambda,AOT,SSA,cosTETA,g,p0)
% Computes single scattering properties for mixture aerosol and air molecule
%------------------------------------------------------------------------

tau_r=rayleigh(lambda,p0);     % Rayleigh optical depth
tau=AOT+tau_r;                 % Total aerosol plus molecular optical depth
omega=(1*tau_r+SSA*AOT)/tau;   % Mean single scattering albedo
G=g*AOT*SSA/(tau_r+AOT*SSA);   % Mean assymetry parameter
% Computation on mean scattering phase function
Snum=100;  
Teta=acos(cosTETA);   
if 1>2    % MIE camputation
  x=4.0;  % size parameter (2*pi*r/lambda);
  n=1.45+0.001i;      % aerosol refractive index
 % Run MIE code for aerosol particles 
 [Paer,ext,SSA]=MieOptics(lambda,n,0.5,2,Teta);
else % Henyey-Greenstein phase function.
    Paer=HGfun(Teta*180/pi,g);
end    
%Pa=interp1(Teta,Paer,TETA_s);
  
Pm=3/4*(1+cosTETA.^2);   % Rayleigh phase function
Pa=Paer;
P11=(Pm*1*tau_r+Pa*AOT*SSA)/tau;    % Phase function for mixture

