function rog=sunglint(wspd,nr,ni,azw,ts,tv,fi)
% input parameters:   wspd=speed of the wind (in m/s)
%                     nr=index of refraction of the sea water
%                     ni=extinction coefficient of the sea water
%                     azw=azim. of the sun minus(-) azim. of the wind (in deg.) ------- relative azimuth angle
%                     ts=solar zenith angle (in deg.)
%                     tv=view zenith angle (in deg.)
%                     fi=relative azimuth (sun-satellite)
% output parameters:  rog=reflectance of the sun glint

            
 fac=pi/180;
 phw=azw*fac;
 cs=cos(ts*fac);
 cv=cos(tv*fac);
 ss=sin(ts*fac);
 sv=sin(tv*fac);
 phi=fi*fac;
 Zx=-sv*sin(phi)/(cs+cv);
 Zy=(ss+sv*cos(phi))/(cs+cv);
 tantilt=sqrt(Zx*Zx+Zy*Zy);
 tilt=atan(tantilt);
 
%Anisotropic Gaussian distribution
%    phw=phi_sun-phi_wind
 SigmaC=0.003+0.00192*wspd;
 SigmaU=0.00316*wspd;
 C21=0.01-0.0086*wspd;
 C03=0.04-0.033*wspd;
 C40=0.40;
 C22=0.12;
 C04=0.23;
 xe=(cos(phw)*Zx+sin(phw)*Zy)/sqrt(SigmaC);
 xn=(-sin(phw)*Zx+cos(phw)*Zy)/sqrt(SigmaU);
 xe2=xe*xe;
 xn2=xn*xn;
 coef=1-C21/2.*(xe2-1)*xn-C03/6.*(xn2-3)*xn;
 coef=coef+C40/24.*(xe2*xe2-6*xe2+3);
 coef=coef+C04/24.*(xn2*xn2-6*xn2+3);
 coef=coef+C22/4.*(xe2-1)*(xn2-1);
 proba=coef/2./pi/sqrt(SigmaU)/sqrt(SigmaC)*exp(-(xe2+xn2)/2.);
% Compute Fresnel's coefficient R1
 cos2chi=cv*cs+sv*ss*cos(phi);
 if (cos2chi>1.0)
    cos2chi=0.99999999999;
 end   
 if (cos2chi<-1.0)
    cos2chi=-0.99999999999;
 end   
 coschi=sqrt(0.5*(1+cos2chi));
 sinchi=sqrt(0.5*(1-cos2chi));
 R1=fresn(nr,ni,coschi,sinchi);
      
% Compute Reflectance of the sun glint
  rog=pi*R1*proba/4./cs/cv/(cos(tilt)^4);
