function Tnew = updateTNonSquareGrid(Fij, h1, h2, T1, T2)

global HUGE

if(T1 == HUGE & T2 == HUGE)
    Tnew = HUGE;
    return;
end
h1sq = h1*h1;
h2sq = h2*h2;
a = (h1sq + h2sq);
Tdiff = T1 - T2;
Discrim = Fij*Fij*a - Tdiff*Tdiff;
if Discrim > 0.0
    % both i and j directions approximate the gradient - solve quadratic equation
    Tnew = (T1*h2sq + T2*h1sq + h1*h2*sqrt(Discrim))/a;
else
    % only i or j direction approximates the gradient - just add appropriate one to F*h
    Tnew = min(Fij*h1 + T1, Fij*h2 + T2);
end         


