function eps_swd_l91sd(f, T) result(eps) 

! This function returns the complex dielectric constant of suspended 
! (pure) water droplets as a function of frequency f (GHz) and 
! temperature T (deg. C), using the single Debye model of Liebe, Hufford, 
! and Manabe (Int. J. IR & mm Waves V.12 (7) July 1991).

  implicit none

  real, intent (in) :: f  ! [GHz]  Frequency (valid from 0 to 1000 GHz)
  real, intent (in) :: T  ! [°C]   Temperature
  complex :: eps          !        Dielectric constant
  real :: theta           !        Inverse temperature variable 
  real :: eps0            !        Static dielectric constant
  real :: epsinf          !        High frequency dielectric constant
  real :: gamma           ! [GHz]  Relaxation frequency

  complex, parameter :: i = (0.0, 1.0)

  theta = 1.0 - 300.0/(273.15 + T)
  eps0 = 77.66 - 103.3*theta
  epsinf = 0.066 * eps0
  gamma = 20.27 + 146.5 * theta + 314.0 * theta**2

  eps = (eps0 - epsinf)/(1.0 - i * f/gamma) + epsinf

  return
end function eps_swd_l91sd


!-----------------------------------------------------------------------
function eps_swd_l91dd(f, T) result(eps) 

! This function returns the complex dielectric constant of suspended 
! (pure) water droplets as a function of frequency f (GHz) and 
! temperature T (deg. C), using the double Debye model of Liebe, Hufford, 
! and Manabe (Int. J. IR & mm Waves V.12 (7) July 1991).

  implicit none

  real, intent (in) :: f  ! [GHz]  Frequency (valid from 0 to 1000 GHz)
  real, intent (in) :: T  ! [°C]   Temperature
  complex :: eps          !        Dielectric constant
  real :: theta           !        Inverse temperature variable 
  real :: eps0            !        Static dielectric constant
  real :: eps1            !        First high-frequency constant
  real :: eps2            !        Second high-frequency constant
  real :: gamma1          ! [GHz]  Primary relaxation frequency
  real :: gamma2          ! [GHz]  Secondary relaxation frequency

  theta = 1.0 - 300.0/(273.15 + T)
  eps0 = 77.66 - 103.3*theta

  eps1 = 0.0671 * eps0 
  gamma1 = 20.20 + 146.4*theta + 316.0*theta*theta 
  eps2 = 3.52 + 7.52*theta 
  gamma2 = 39.8*gamma1 

  eps = eps0 - f * ((eps0 - eps1)/cmplx(f,gamma1)  &
                  + (eps1 - eps2)/cmplx(f,gamma2))

  return
end function eps_swd_l91dd