getglint.f

Go to the documentation of this file.

C Incident ray
C       i = ( sin(PHI)       , 0              , cos(PHI) )
C
C Reflected ray
C       o = ( sin(MU)cos(NU) , sin(MU)sin(NU) , cos(MU)  )
C
C Reflection angle
C       i . o = cos(2 OMEGA) = sin(PHI)sin(MU)cos(NU)+cos(PHI)cos(MU)
C
C Normal to plane of reflection
C       n = ( -sin(BETA)cos(ALPHA) , -sin(BETA)sin(ALPHA) , cos(BETA) )
C
C Reflection plane orientation
C       n = i+o
C
C       n = ( sin(PHI)+sin(MU)cos(NU) , sin(MU)sin(NU) , cos(PHI)+cos(MU) )/
C            (2 cos(OMEGA))
C
C  Tilt of plane
C       cos(BETA) = (cos(PHI)+cos(MU))/(2 cos(OMEGA))
C
C  Rotation of plane
C       (i+n)^2 = 4 cos(OMEGA/2)
C
C       cos(OMEGA) = cos(BETA)cos(PHI)-sin(BETA)sin(PHI)cos(ALPHA)
C
C   See papers by Cox and Munk (1954a,b; 1956) for derivation of equations.
 
 
        subroutine getglint(X1,X2,X3,X4,X5,X6)
C
C  Calculate sun gliter coefficient
C
C   X1  Angle MU       (sensor zenith angle)
C   X2  Angle PHI      (solar zenith angle)
C   X3  Angle NU       (sensor-sun azimuth)
C   X4  Wind speed     (m/s)
C   X5  Wind direction (radians)
C   X6  Radiance       (ignoring the atmosphere)
C
C #include "sensor_cmn.fin"
        save
        real pi
        parameter(pi=3.1415926)
        real X2, X3, X4, X5, acoss, RHO
        real X1, Y1, RAD, X, Y2, OMEGA, BETA, ALPHA, SIGC, SIGU, Y3
        real CHI, ALPHAP, SWIG, ETA, PROB, X6, EXPON
        real Y4
 
        rad(x) = x*pi/180.
 
        y4 = max(x4,0.001)
 
        y1 = rad(x1)
        if (y1.eq.0.) y1 = 1.e-7
        y2 = rad(x2)
        if (y2.eq.0.) y2 = 1.e-7
        y3 = rad(x3)
        omega  = acoss(cos(y1)*cos(y2)-sin(y1)*sin(y2)*cos(y3))/2.
        if (omega.eq.0.) omega = 1.e-7
        beta   = acoss((cos(y1)+cos(y2))/(2.*cos(omega)))
        if (beta.eq.0.) beta = 1.e-7
        alpha  = acoss((cos(beta)*cos(y2)-cos(omega))/(sin(beta)*sin(y2)))
        if (sin(y3).lt.0.) alpha = -alpha
 
c Isotropic wind
c        if (glintOn .ge. 2) then
c           ! from Ebuchi & Kizu
c            SIGC = SQRT(0.0048 + 0.00152*Y4)
c            SIGU = SQRT(0.0053 + 0.000671*Y4)
c        else
c           ! from Cox & Munk
            sigc   = .04964*sqrt(y4)
            sigu   = .04964*sqrt(y4)
c        endif
        chi    = x5
        alphap = alpha+chi
        swig   = sin(alphap)*tan(beta)/sigc
        eta    = cos(alphap)*tan(beta)/sigu
        expon  = -(swig**2+eta**2)/2. 
        if (expon.lt.-30.) expon = -30.         ! trap underflow
        if (expon.gt.+30.) expon = +30.         ! trap overflow
        prob   = exp(expon)/(2.*pi*sigu*sigc)
        call reflec(omega,rho)
 
c Normal distribution
        x6     = rho*prob/(4.*cos(y1)*cos(beta)**4)
 
        return
        end
 
 
        subroutine getglint_iqu(X1,X2,X3,X4,X5,X6,X7,X8)
 
C  Calculate sun gliter coefficient, including polarization
C
C       X1(real) - Angle MU       (sensor zenith angle)
C       X2(real) - Angle PHI      (solar zenith angle)
C       X3(real) - Angle NU       (sensor-sun azimuth)
C       X4(real) - Wind speed     (m/s)
C       X5(real) - Wind direction
C       X6(real) - Sun glitter coefficient
C       X7(real) - Q/I for glitter
C       X8(real) - U/I for glitter
C
C #include "sensor_cmn.fin"
        save
        real pi
        parameter(pi=3.1415926)
        real X2, X3, X4, X5, acoss, RHO_PLUS
        real X1, Y1, RAD, X, Y2, OMEGA, BETA, ALPHA, SIGC, SIGU
        real CHI, ALPHAP, SWIG, ETA, PROB, X6, EXPON
        real X7, X8, RHO_MINUS, SR, CR, ROT_ANG, C2R, S2R, ASINN
        real Y4
 
        rad(x) = x*pi/180.
 
        y4 = max(x4,0.001)
 
        y1 = rad(x1)
        if (y1.eq.0.) y1 = 1.e-7
        y2 = rad(x2)
        if (y2.eq.0.) y2 = 1.e-7
        y3 = rad(x3)
        omega  = acoss(cos(y1)*cos(y2)-sin(y1)*sin(y2)*cos(y3))/2.
        if (omega.eq.0.) omega = 1.e-7
        beta   = acoss((cos(y1)+cos(y2))/(2.*cos(omega)))
        if (beta.eq.0.) beta = 1.e-7
        alpha  = acoss((cos(beta)*cos(y2)-cos(omega))/(sin(beta)*sin(y2)))
        if (sin(y3).lt.0.) alpha = -alpha
 
c Isotropic wind
c        if (glintOn .ge. 2) then
c           ! from Ebuchi & Kizu
c            SIGC = SQRT(0.0048 + 0.00152*Y4)
c            SIGU = SQRT(0.0053 + 0.000671*Y4)
c        else
c           ! from Cox & Munk
            sigc   = .04964*sqrt(y4)
            sigu   = .04964*sqrt(y4)
c        endif
        chi    = x5
        alphap = alpha+chi
        swig   = sin(alphap)*tan(beta)/sigc
        eta    = cos(alphap)*tan(beta)/sigu
        expon  = -(swig**2+eta**2)/2.
        if (expon.lt.-30.) expon = -30.         ! trap underflow
        if (expon.gt.+30.) expon = +30.         ! trap overflow
        prob   = exp(expon)/(2.*pi*sigu*sigc)
        call reflec_both(omega,rho_plus,rho_minus)
 
c Normal distribution
        x6     = rho_plus*prob/(4.*cos(y1)*cos(beta)**4)
 
c Polarization components
        if (omega .gt. .0001) then
          cr = (cos(y2) - cos(2.*omega)*cos(y1))/(sin(2.*omega)*sin(y2))
          sr = sin(y2)*sin(pi-y3) / sin(2.*omega)
          rot_ang = sign(1., cr)*asinn(sr)
        else
          rot_ang = pi/2.
        endif
 
        c2r = cos(2.*rot_ang)
        s2r = sin(2.*rot_ang)
 
        x7 =  c2r * rho_minus / rho_plus ! q_ov_i
        x8 = -s2r * rho_minus / rho_plus ! u_ov_i
 
        return
        end
 
 
        subroutine reflec (X1,X3)
C
C  X1  Incident angle (radians)
C  X3  Reflectance
C
C       n1 sin(x1) = n2 sin(x2)
C
C                    tan(x1-x2)**2
C       Refl(par ) = -------------
C                    tan(x1+x2)**2
C
C                    sin(x1-x2)**2
C       Refl(perp) = -------------
C                    sin(x1+x2)**2
C
C  Where:
C       x1  Incident angle
C       n1  Index refraction of Air
C       x2  Refracted angle
C       n2  Index refraction of Water
C
      real X1, X2, X3, ref
*                                       ! Index refraction of sea water
        ref = 4./3.
        if (x1.lt..00001) then
          x3 = .0204078
        else
          x2 = asin(sin(x1)/ref)
          x3 = (sin(x1-x2)/sin(x1+x2))**2+(tan(x1-x2)/tan(x1+x2))**2
          x3 = x3/2.
        end if
        return
        end
 
 
        subroutine reflec_both (X1,X3,X4)
C
C  X1  Incident angle (radians)
C  X3  Reflectance sum
C  X4  Reflectance difference
C
C       n1 sin(x1) = n2 sin(x2)
C
C                    tan(x1-x2)**2
C       Refl(par ) = -------------
C                    tan(x1+x2)**2
C
C                    sin(x1-x2)**2
C       Refl(perp) = -------------
C                    sin(x1+x2)**2
C
C  Where:
C       x1  Incident angle
C       n1  Index refraction of Air
C       x2  Refracted angle
C       n2  Index refraction of Water
C
      real X1, X2, X3, X4, REF
      real PERP, PAR
c                                       ! Index refraction of sea water
        ref = 4./3.
        if (x1.lt..00001) then
          x3 = .0204078
          x4 = 0.
        else
          x2 = asin(sin(x1)/ref)
          perp = (sin(x1-x2)/sin(x1+x2))**2
          par  = (tan(x1-x2)/tan(x1+x2))**2
          x3 =  perp+par
          x3 = x3/2.
          x4 = -perp+par
          x4 = x4/2.
        end if
        return
        end
 
 
        function acoss (X1)
        real    acoss, x1
C
C   This calculates ACOS(X) when X is near + or - 1
C
        if (x1.ge.1.) then
          acoss = 0.
        else if (x1.le.-1.) then
          acoss = 3.141592654
        else
          acoss = acos(x1)
        end if
        return
        end
 
 
        function asinn (X1)
        real    asinn, x1
C
C   This calculates ASIN(X) when X is near + or - 1
C
        if (x1.ge.1.) then
          asinn = 3.141592654/2.0
        else if (x1.le.-1.) then
          asinn =-3.141592654/2.0
        else
          asinn = asin(x1)
        end if
        return
        end