% Temperature of water in Celsius
T=5;

% Salinity in psu
S=35;

% Frequency in hertz
freq=1413e6;

% Observation angle in radians
th=60*pi/180;

e0=(87.134-1.949e-1*T-1.276e-2*T^2+2.491e-4*T^3)*(1.+1.613e-5*T*S-3.656e-3*S+3.21e-5*S^2-4.232e-7*S^3);

einf=4.9;

relax=(1.1109e-10-3.824e-12*T+6.398e-14*T^2-5.096e-16*T^3)*(1.+2.282e-5*T*S-7.638e-4*S-7.760e-6*S^2+1.105e-8*S^3);

sigma=S*(.18252-1.4619e-3*S+2.093e-5*S^2-1.282e-7*S^3)*exp((T-25)*(2.033e-2+1.266e-4*(25-T)+2.464e-6*(25-T)^2-S*(1.849e-5-2.551e-7*(25-T)+2.551e-8*(25-T)^2))); 

[r c]=size(freq);

epsr=(e0-einf)./(1+(relax*freq).^2)+einf

epsi=(relax*freq.*(e0-einf))./(1+(relax*freq).^2)+(sigma/2./pi)./freq/8.854e-12

eps=epsr+j*epsi;

rh=(cos(th)-sqrt(eps-sin(th)^2))./(cos(th)+sqrt(eps-sin(th)^2));
rv=(eps*cos(th)-sqrt(eps-sin(th)^2))./(eps*cos(th)+sqrt(eps-sin(th)^2));

th0=(1-abs(rh).^2)*(273.16+T);
tv0=(1-abs(rv).^2)*(273.16+T);


[th0 tv0]