function liquid_refractive_index_mw_turner, freq_in, temp_in

  ;This function calculates the complex permittivity as a function of temperature and frequency.
  ; To do so, we use the two-term Debye function given in Turner Kneifel Cadeddu. (2016) in JTech.  
  ;
  ; Inputs: freq = frequency in GHz
  ;         temp = temperature in Kelvins
  ;
  ; Outputs: a = real part of the complex permittivity
  ;          b = imaginary part of the complex permittivity
  ; Note: The permittivity (epsilon) is a complex number in the form of a - bi.
  ;
freq=freq_in*1.d9
temp=double(temp_in)-273.15d

  S0=8.79d1
  S1=-4.0440d-1
  S2=9.5873d-4
  S3=-1.328d-6
  
  a1=8.111d1
  b1=4.434d-3
  c1=1.302d-13
  d1=6.627d2
  a2=2.025d
  b2=1.073d-2
  c2=1.012d-14
  d2=6.089d2
  tc=1.342d2
  
  c=29979245800.d ; cm/s
  rho_l=1.
  
  
  eps_s=S0+(S1*temp)+(S2*temp^2)+(S3*temp^3)
  
  tau1=c1*exp(d1/(temp+tc))
  tau2=c2*exp(d2/(temp+tc))
  
  delta1=a1*exp(-b1*temp)
  delta2=a2*exp(-b2*temp)
  
  aa1=((tau1^2)*delta1)/((1.d)+(((2.d)*!dpi*freq*tau1)^2))
  aa2=((tau2^2)*delta2)/((1.d)+(((2.d)*!dpi*freq*tau2)^2))
  
  bb1=(tau1*delta1)/((1.d)+(((2.d)*!dpi*freq*tau1)^2))
  bb2=(tau2*delta2)/((1.d)+(((2.d)*!dpi*freq*tau2)^2))
  
  eps=eps_s-((((2.d)*!dpi*freq)^2)*(aa1+aa2))
  
  eps_prime=((((2.d)*!dpi*freq))*(bb1+bb2))
  
  eps_imag=complex(eps,eps_prime)

  a_ref=sqrt(( sqrt(eps^2+eps_prime^2)+eps   )/2)
  b_ref=sqrt(( sqrt(eps^2+eps_prime^2)-eps   )/2)
  x=complex(a_ref,b_ref)
  k_sqrd_liq=abs((((x^2)-1.)/((x^2)+2.))^2)
  
  alpha_l=((6.d)*!dpi*freq/(rho_l*c))*imaginary((eps_imag-(1.d))/(eps_imag+(2.d)))  ; mass absorption coefficient in cm^2/g

  ;print, eps, eps_prime, a_ref, b_ref, k_sqrd_liq, alpha_l


  return, x
end