Possible error in linearisation in ropp_pp_geometric_optics

[Ian Culverwell]

User João Francisco Galera Monico has uncovered a potential error in the linearisation of eqns 2.18-2.21, ROPP user guide III, in ropp_geometric_optics.f90.

He writes:

Hello Kristian...

I have been studying several topics realted with ROPP. In the subroutine ropp_pp_bending_angle_go, which call ropp_pp_geometric_optics, I found that maybe we have a problem.

In the linearization we have:

A(1,1:3) = -v_leo/(c_light - DOT_PRODUCT(v_leo,U_leo))

but from my derivation it is

A(1,1:3) = -v_leo/(c_light - DOT_PRODUCT(v_gns,U_gns)).

Could you verify that?

Thank you

J F Galera Monico

I think he's right, but need to investigate further.

[Ian Culverwell]

This comes from the expression for the non-relativistic fractional Doppler shift = (c - v_L.u_L)/(c - vG.uG).

This and 5 other expressions are solved simultaneously to give (u_L, u_G) = unit vectors in dirn of motion of L & G. In the code this is DF as a function of (DZ): DF(DZ)=0. This is solved iteratively, by Newton Raphson, so at each step we have DDZ = -(dDF/dDZ)^{-1 DF, and then DZ --> DZ + DDZ. It's the matrix A = dDF/dDZ that appears to be wrong (and maybe in OCC too).}

However, as long as it converges to DDZ=0 then it shouldn't matter too much if dDF/dDZ is a bit wrong. So I think this may be less serious than it looks.

I am concerned that correcting it might affect the degree of convergence to DDZ=0, since the code uses a fixed number of iterations (=5). This needs to be checked.

[Ian Culverwell]

Details_on_the_subroutine_ropp.docx

[Ian Culverwell]

See attached explanation from João Francisco.

I (idculv) agree about A(1,1:3) being wrong.

The same error appears in occ_refraction.f90, so someone will have to speak to Michael Gorbunov about it.

Initially, there are apparently small differences in refractivity profiles - evidence of insensitivity of NR iteration to the estimate of df/dx? (Eg when finding root of f(x)=x*x-2=0, get convergence with f'(x)=2x, 5x, (1.5+r)x where r~U(0,1), ...)

Need more evidence that correcting it won't mess up things before making a change.

[Ian Culverwell]

Further theoretical evidence comes from the (belated) realisation that u_L and u_G are unit vectors; hence |V_L.u_L | </= | |V_L| | << c. So the difference in A(1,1:3) between the two methods of calculation would be negligible. (Michael Gorbunov agrees.) Indeed, we could probably just set A(1,1:3) = -V_L/c and A(1, 4:6) = V_G/c, which at least has the virtue of symmetry.

[Ian Culverwell]

These changes have been made to ropp_pp_geometric_optics.f90 and ropp_pp_geometric_optics_adj.f90 in the ROPP5.1_prototype.

How to test? These routines are called by ropp_pp_radiooptic_analysis.f90 and ropp_pp_DCT.f90.

The former is called from ropp_pp_spectra_tool.f90, which is tested as part of the general ropp_pp testing. So we should see a difference there. From all that has been said above, it should be small. Is it?

Yes - very small. The difference in the ROanalysis_??_L?.dat files, which ropp_pp_spectra_tool generates, is less than the 5 or 6 significant figures that are printed out by default. In fact, you have to artificially increase the output precision to 10 significant figures before you see occasional (~0.2% of the time) differences in the ROanalysis_ep_L1/2.dat files.

So: it's correct, and it won't stop anything passing the standard ropp_pp testing. Put it in the code.

[Ian Culverwell]

Footnote: although this does indeed have a small impact on the results (see above), it's large enough to cause ropp5.1 to fail some of the it-pp-?? tests in the ropp test folder, if the results are judged against uncorrected OCC results. I suggest this is evidence that the pass criterion in these tests (<0.1% frac diff) is too stringent.