function th=eosben07_th(p,r,s)
%EOSBEN07_TH Potential temperature.
%   EOSBEN07_TH(P, R, TH) computes the potential temperature (referenced
%   to 0 dbar) [degC] at pressure, P [dbar], density, R [kg m-3], and
%   salinity, S [g kg-1].

eosben07_const

a=a14-a24*r;
b=a12-a22*r+(a15-a25*r).*s+(b12-b22*r).*p;
c=a11-a21*r+(b11-b21*r).*p+(a13-a23*r+(a16-a26*r).*s+(b13-b23*r).*p).*s;

th=(-b-sqrt(b.*b-4*a.*c))./(2*a);
