function [i1,i2,CosTETA]=angles(Teta,Fi,Teta_p,Fi_p)
% calculate angles for polarized transfer equation based on Liou book
% p. 324
m=cos(Teta);
mp=cos(Teta_p);
CosTETA=m*mp+sqrt((1-m^2)*(1-mp^2))*cos(Fi-Fi_p);  % cos of scattering angle
if CosTETA==1    % particular cases
    i1=0;
    i2=pi;
else  
  cosi1=(m*sqrt(1-mp^2)-mp*sqrt(1-m^2)*cos(Fi-Fi_p))/sqrt(1-CosTETA^2);
  cosi2=(mp*sqrt(1-m^2)-m*sqrt(1-mp^2)*cos(Fi-Fi_p))/sqrt(1-CosTETA^2);
 
 if (cosi1>1)       % check if cos is out of cosine domain
     cosi1=1;
 elseif (cosi1<-1)
     cosi1=-1;
 end
 if (cosi2>1)
     cosi2=1;
 elseif (cosi2<-1)
     cosi2=-1;
 end
 i1=acos(cosi1);
 i2=acos(cosi2);
end 
