function [cosTeta_new,sinTeta_new,cosFi_new,sinFi_new]=uklad_new(cosTeta0,sinTeta0,...
    cosFi0,sinFi0,cosTeta,sinTeta,cosFi,sinFi);
% input:
 % cosTeta0,sinTeta0,cosFi0,sinFi0, old photon direction
 % cosTeta,sinTeta,cosFi,sinFi % components of scatering direction
% output: 
 % cosTeta_new,sinTeta_new,cosFinew,sinFinew  , new photon direction 
cosTeta_new=cosTeta0*cosTeta+sinTeta0*sinTeta*cosFi;
sinTeta_new=sqrt((1.0-cosTeta_new)*(1.0+cosTeta_new));

if sinTeta_new~=0.0
    if sinTeta0~=0.0
        sinDFi=sinTeta*sinFi/sinTeta_new;
        cosDFi=(cosTeta-cosTeta0*cosTeta_new)/(sinTeta0*sinTeta_new);
        cosFi_new=cosFi0*cosDFi+sinFi0*sinDFi;
        sinFi_new=sinFi0*cosDFi-cosFi0*sinDFi;
    else
        if sinFi0>=0.0
            cosFi_new=cosFi;
            sinFi_new=sinFi;
        else
            cosFi_new=-cosFi;
            sinFi_new=-sinFi;
        end 
    end
else
    cosFi_new=cosFi;
    sinFi_new=sinFi;
end    
        