; This file will calculate the emissivity of the ocean surface
; in both V and H polarizations at a view angle of 53.1 degrees.
;
;   It contains four subroutines:
;   eps_sea_wentz  : calculates dielectric constant of ocean surface
;   fresnel        : calculates specular emissivity of ocean surface.
;   wind_roughen   : computes wind-roughened emissivity of ocean surface.
;   emiss          : calculates full emissivity of ocean surface in V & H,
;                    taking into account both whitecaps (foam) and wind
;                    roughening.

;   NOTE: You only need to call "emiss"!  It in turn will call the other 
;         functions as needed.
;
;   Calling syntax:
;
;   emiss, freq, Tsurf, wind, ev, eh
;   
;   freq : frequency in GHz
;   Tsurf: Surface Temperature in K
;   wind : Surface wind speed in m/s
;   ev, eh: Ocean surface emissivity for V,H polarization


function eps_sea_wentz, f, T 
; This function returns the complex dielectric constant of sea water,
; using the model of Wentz (1996; AMSR Ocean Algorithm ATBD, RSS Tech.
; Report 120296)

; Inputs
;  real, intent (in) :: f     ; [GHz]    Frequency (valid from 0 to 1000 GHz)
;  real, intent (in) :: T     ; [°C]     Temperature

; Output
;  complex :: eps             ;          Dielectric constant

  clight = 2.998e10  ; [cm/s]
  pi = 3.14159265
  salinity = 32.54

; Eq. (36) -- static dielectric constant of distilled water
  eps_s0 = 87.90 * exp( - 0.004585 * T)

; Eq. (38) -- relaxation wavelength of distilled water
  lambda_r0 = 3.30 * exp( - 0.0346 * T + 0.00017 * T^2)

; Eq. (41)
  C = 0.5536 * salinity

; Eq. (42)
  delta = 25.0 - T

; Eq. (40)
  zeta = 2.03e-2 + 1.27e-4*delta + 2.46e-6*delta^2  $
       - C * (3.34e-5 - 4.60e-7*delta + 4.60e-8*delta^2)

; Eq. (39)
  sigma = 3.39e9 * (C^0.892) * exp( - delta * zeta)

; Eq. (43)
  eps_s = eps_s0 * exp(- 3.45e-3*salinity + 4.69e-6*salinity^2 + 1.36e-5*salinity*T)

; Eq. (44)
  lambda_r = lambda_r0 - 6.54e-3*(1.0 - 3.06e-2*T + 2.0e-4*T^2)*salinity

; radiation wavelength (cm)
  lambda = clight / (f * 1.0e9)

; from p. 20
  eta = 0.012
  eps_inf = 4.44

; Eq. (35)
  term1 = complex(0.0, - lambda_r/lambda)
  term1 = term1 ^ (1.0 - eta)
  term2 = complex(0.0, - 2.0 * sigma * lambda / clight)
  eps = eps_inf + (eps_s - eps_inf)/(1.0 + term1) - term2

  return , eps
end

PRO fresnel, XPERM, ANGLE, ev, eh
; Calculates the vertical (RV) and Horizontal (RH) reflection
; coefficients from a dielectric surface given the
; complex electrical permittivity of the dielectric.
; This uses the standard Fresnel relations at a dielectric interface.
;
;  INPUTS:
;   XPERM : (Complex) electrical permittivity
;   ANGLE : Viewing angle from Nadir (degrees)
;  
; OUTPUTS:
;     RV : Reflectivity for vertical (V) polarization
;     RH : Reflectivity for horizontal (H) polarization
; 
; NOTE: Emissivity = 1 - Reflectivity

       PS2 = sin(angle*!pi/180.)^2
       PCC = cos(angle*!pi/180.)
       PERM1 = SQRT(XPERM - PS2)
       PERM2  = XPERM*PCC
       RHTH = (   PCC - PERM1)/(   PCC + PERM1)
       RVTH = (PERM2 - PERM1)/(PERM2 + PERM1)
       
       RH =(real_part(RHTH))^2+(imaginary(RHTH))^2
       RV =(real_part(RVTH))^2+(imaginary(RVTH))^2
       
       ev = 1. - rv
       eh = 1. - rh       
END

PRO wind_roughen, wind, freq, ev, eh
; This routines calculates the wind-roughened ocean
; surface emissivity at a view angle of 53.1 degrees, 
; given the specular ocean surface emissivity, 
; surface wind speed in m/s, and the frequency of observation in GHz.
;
; It is based on the FASTEM-2 algorithm of Deblonde and English 
; (UK Met Office, 2000).
C = [ -9.4674e-03, -1.5567e-03,  3.6254e-06, -4.6501e-02, -1.8435e-03,  6.3099e-06,  2.9533e-04,  4.4369e-06, -1.4016e-08, $
       5.3421e-02,  3.8867e-04, -3.5054e-06,  9.3315e-02,  1.3192e-03, -5.3577e-06,  2.3784e-05,  8.6950e-07,  2.8249e-09 ]
C = reform(C, 3, 3, 2, /over)

xcorr = fltarr(2)
z =fltarr(3)
for jp = 0, 1 do begin ; 2 polarizations
   z[*] = c[0,*,jp] + c[1,*,jp]*freq + c[2,*,jp]*freq^2
   xcorr[jp] = z[0] + z[1] * wind + z[2] * wind^2
endfor
xcorr = xcorr/100.
ev = ev + xcorr[0]
eh = eh + xcorr[1]

END

PRO emiss, freq, Tsurf, Wind, ev, eh
; Calculates the ocean surface emissivity in V and H polarizations 
; at a view angle of 53.1 degrees.  It takes into account both
; wind roughening of the ocean surface, and foam coverage as parameterized
; by Monahan and Muircheartaigh, J. Phys. Oceanogr., Vol 10, pp 2094-2099, 1980.
; A microwave foam emissivity of 1.0 is assumed.

; 1 - Get the dielectric constant at these frequencies
; and this surface temperature.
eps = eps_sea_wentz(freq, Tsurf-273.15)

; 2 - Get the specular emissivity in V and H.
fresnel, eps, 53.1, ev, eh

; 3 - Compute the wind-roughened emissivities.
wind_roughen, wind, freq, ev, eh

; 4 - Compute the foam fraction and adjust emissivities accordingly.
frac = (2.95d-6 * wind^3.52) < 0.6
ev = ev*(1.-frac) + frac
eh = eh*(1.-frac) + frac

END