Rewrite RoC calculation

[Ian Culverwell]

The calculation of radius of curvature in ropp_utils/coordinates/curvature.f90 is a bit dodgy: Ak2 is indeterminate (0/0) at lat = pi/2, and Ak1 includes such huge quantities (like ((SQRT(Re4))3) that it causes a floating point overflow if evaluated at single precision!

Suggest rewriting in terms of the principal radii R_NS and R_EW.

(Prompted by discussions with Josep Aparicio.)

[Ian Culverwell]

Replace

! 1.3 Calculate radius of curvature
! ---------------------------------
     
  Ak1 = Rp*Re**4 / (Sqrt(Re**4 + (Rp**2 - Re**2)*Sum(x(1:2)**2)))**3
  Ak2 = Cos(rlat) / Sqrt(Sum(X(1:2)**2))
  roc  = 1.0_wp / Abs(Ak1*(Cos(Theta))**2 + Ak2*(Sin(Theta))**2)

which gives, at lat=0.14123822136156511942E+01 , azi=0.40621548728487049118E+01,

RoC = 0.63624451366924149916E+07,

by

! 1.3 Calculate radius of curvature
! ---------------------------------

  tphi = (Re/Rp)**2 * x(3) / SQRT(x(1)**2 + x(2)**2)
  cphi = 1.0_wp / SQRT(1.0_wp + tphi**2)
  sphi = tphi*cphi
  r_NS = Re * (Rp/Re)**2 / (SQRT(cphi**2 + (sphi*Rp/Re)**2))**3
  r_EW = Re              /  SQRT(cphi**2 + (sphi*Rp/Re)**2)
  roc  = 1.0_wp / ((COS(theta)**2/r_NS) + (SIN(theta)**2/r_EW))

which gives, at lat=0.14123822136156511942E+01 , azi=0.40621548728487049118E+01,

RoC = 0.63624451366924159229E+07

ie the same within 1 nanometre, which is probably good enough.

[Ian Culverwell]

Close ticket.

[Ian Culverwell]

Actually it's done a lot neater by writing (see r4172)

  tphi2 = (Re/Rp)**4 * x(3)**2 / (x(1)**2 + x(2)**2)
  cphi2 = 1.0_wp / (1.0_wp + tphi2)
  sphi2 = tphi2*cphi2
  r_NS = Re * (Rp/Re)**2 / (SQRT(cphi2 + sphi2*(Rp/Re)**2))**3
  r_EW = Re              /  SQRT(cphi2 + sphi2*(Rp/Re)**2)
  roc  = 1.0_wp / ((COS(theta)**2/r_NS) + (SIN(theta)**2/r_EW))

Makes no difference to RoC, so reclose ticket.

[Ian Culverwell]

Actually, this still goes wrong at the pole, because tphi2 is infinite there.

What we want is:

lambda = (Re/Rp)^2

cphi2 = (x^2 + y^2)/(x^2 + y^2 + (lambda*z)^2) = 0,1 at pole,eq - OK

sphi2 = ((lambda*z)^2)/(x^2 + y^2 + (lambda*z)^2) = 1,0 at pole,eq - OK

r_EW = 1/sqrt(cphi2 + sphi2/lambda) = Re/Rp,1 at pole,eq - OK

r_NS = r_EW^3 / lambda = Re/Rp,(Rp/Re)^2 at pole,eq - OK

roc = Re/(cos(theta)^2/r_NS + sin(theta)^2/r_EW) = Re^2/Rp,Re-->Rp^2/Re at pole,eq - OK

Reopening for ROPP9.0.

[Ian Culverwell]

Trying out the above, we get

At lat=pi/2:
------------
 **** PREVIOUS CALCULATION ****
roc =                            NaN

 **** LATEST CALCULATION ****
roc =            0.63995936257584933192E+07
cf
Re/(Rp/Re) =     0.63995936257584933192E+07

At lat=0, azi=0:
----------------
 **** PREVIOUS CALCULATION ****
roc =                 0.63354393272928195074E+07
cf
1/(1/Re*(Rp/Re)**2) = 0.63354393272928195074E+07

 **** LATEST CALCULATION ****
roc =                 0.63354393272928213701E+07
cf
Re/(1/1/(Re/Rp)**2) = 0.63354393272928213701E+07

(2 nm difference.)

At lat=0, azi=pi/2:
-------------------
 **** PREVIOUS CALCULATION ****
roc =       0.63781370000000000000E+07
cf
1/(1/Re) =  0.63781370000000000000E+07

 **** LATEST CALCULATION ****
roc =     0.63781370000000000000E+07
cf
Re =      0.63781370000000000000E+07

(No difference.)

At lat=0.14123822136156511942E+01, 
---------------------------------
azi=0.40621548728487049118E+01:
-------------------------------
 **** PREVIOUS CALCULATION ****
roc =     0.63624451366924159229E+07
 **** LATEST CALCULATION ****
roc =     0.63624451366924187168E+07

(3 nm difference.)

Good enough, so commit at r4746.

[Ian Culverwell]

Closing ticket (again).