Introduce wave optics propagator tool (WOPT)

[cburrows]

​This branch includes the developments of a wave optics propagator that is intended to become part of the ropp_pp module. This takes a vertical refractivity profile as input and simulates the phase and amplitude values that such an atmosphere would produce for a fixed GNSS and a LEO in a circular orbit.

[Ian Culverwell]

This has been ROPPified. In doing so, I noticed a crazy sensitivity, which is worth recording before things are overwritten.

The command

ropp_pp_wopt_tool CASE_07_ref1d.nc -o ../data/CASE_07_ref1d_out2.nc

generates this excess phase

and this SNR:

Both of these look reasonably sensible.

But simply changing

r_gps(:,2) = - SQRT(2.66E7**2-r_gps(:,1)**2)     ! radius of GPS = 26600 km

to

r_gps(:,2) = - SQRT(2.66E7_wp**2-r_gps(:,1)**2)     ! radius of GPS = 26600 km

in ropp_pp_wopt_tool.f90 (ie just promoting the GNSS altitude 2.66e7 m from single precision to double precision) generates excess phases that differ from the above by this much:

This is crazy: about 500 m difference when the phase is 10 000 m, ie about 5% difference in excess phase.

The SNRs differ from the single precision result by this much:

At the end of the (setting) occultation, this is around 0.001 in 0.01, ie ~ 10%.

This seems a very large 'structural uncertainty'. The differences grow rapidly until dim_lev1a ~ 14500, towards the end of the (setting) occultation, so perhaps it occurs when s_geom = |r_GPS - r_LEO| exceeds some threshold.

[Ian Culverwell]

CASE_07_excess_phase.png

[Ian Culverwell]

CASE_07_SNR.png

[Ian Culverwell]

CASE_07_excess_phase_diff.png

[Ian Culverwell]

CASE_07_SNR_diff.png

[Ian Culverwell]

Fortunately the differences in refractivity

and bending angle

resulting from this change, are pretty small.

(These are obtained by running, for example,

../tools/ropp_pp_occ_tool -w CASE_07_ref1d_out.nc -o occ_out.nc -c ../config/wopt_pp.cf -d

)

[Ian Culverwell]

CASE_07_refrac_diff.png

[Ian Culverwell]

CASE_07_bangle_diff.png

[Ian Culverwell]

I meant the differences grow rapidly after dim_lev1a ~ 14500. They're very small (1e-5 m in phase, 1e-7 V/V in SMR) until then.

[Ian Culverwell]

Oddly, for CASE 12 (a more difficult one), we don't see anything like the same difference in excess phase

although the difference in SNR is similar to CASE_07:

[Ian Culverwell]

CASE_12_excess_phase_diff.png

[Ian Culverwell]

CASE_12_SNR_diff.png

[Ian Culverwell]

I mean and

Weirder and weirder.

[Ian Culverwell]

Increasing nsmooth (phase modelling window in ropp_wopt_phase_and_amplitude_LEO) from 4 to 8 reduces the difference in excess phase to ~15 m, and with negligible impact on retrieved refractivities.

[Ian Culverwell]

it_wopt_02_CASE_07_nsmooth=8.jpg

[Ian Culverwell]

CASE_07_excess_phase_diff_nsmooth=8.png

[Ian Culverwell]

CASE_07_SNR_diff_nsmooth=8.png

[Ian Culverwell]

A bit more sensitivity came to light when testing CASE12 (a difficult one). When I use

IF (ABS(U(i)) > 1.0E-5*Amp_max)

in ropp_pp_wopt_propagate_to_leo.f90 and propagate through WOPT and then OCC, the resulting refracs look bas cf: the input ones: When I reduce this threshold, which allows more points to contribute to the signal at the LEO, to

IF (ABS(U(i)) > 1.0E-6*Amp_max)

then the fit is closer: This looks better, and is closer to wnat I got with Sean's original code (with 1.e-5, however):

[Ian Culverwell]

ropp_pp_comp_CASE_12_idc_1.0E-5.jpg

[Ian Culverwell]

ropp_pp_comp_CASE_12_idc_1.0E-6.jpg

[Ian Culverwell]

ropp_pp_comp_CASE_12_sbh.jpg

[Ian Culverwell]

ropp_pp_wopt_tool_out_CASE_12_occ_idc.out-ropp_pp_wopt_tool_out_CASE_12_occ_sbh.out.png

[Ian Culverwell]

I think the reduction in threshold is justifiable, but I don't understand the different behaviour of both codes with 1e-5. What have I done to change this? One clue: the output from the OCC processing is very different:

[Ian Culverwell]

Note in particular the failure of the Global MSIS search for Sean's data. (Of course, iono corrn should have no effect as there's no ionosphere in this problem: phase_L1 = phase_L2.)

[Ian Culverwell]

Note that if I revert my code to use Sean's parameters, namely:

real(wp) :: f_L1, c_light, pi1, radius_earth
f_L1 = 1.57542e9_wp
c_light = 3.0e8_wp
pi1 =  3.141592653589793238_wp
radius_earth = 6.378135e6_wp  
radius_gps = 2.66e7
theta_min = 0.451953697_wp 
theta_dot = 1.036e-3_wp 
  REAL(wp), PARAMETER      :: amplitude_cutoff = 1.0E-5_wp 

then I get very similar results to Sean:

And the occ processing is also very similar:

So I evidently haven't introduced any errors in my translation of the code: the differences in these parameters are deliberate (mainly done to match the values used in the ROPP OCC code).

The difference in the OCC processing probably accounts for the difference in the retrieved refracs, and I think we can just say that in this case Sean's happened to give a closer fit to the input refracs than mine.

[Ian Culverwell]

ropp_pp_comp_CASE_12_idc_1.0E-5_SBH_params.jpg

[Ian Culverwell]

ropp_pp_wopt_tool_out_CASE_12_occ_idc_SBH_params.out-ropp_pp_wopt_tool_out_CASE_12_occ_sbh.out.png

[Ian Culverwell]

Summary

• The new code matches the original pretty closely - no apparent errors
• This is true even on very challenging cases such as CASE 12
• There are good grounds for setting amplitude_cutoff = 1.0E-6_wp

[Ian Culverwell]

Commit various aspects of these changes (to https://trac.romsaf.org/ropp/browser/ropp_src/branches/dev/Share/cb_wopt_test):

• Tool and man page at r4975.
• Support routines at r4976 and r4981.
• Test routines at r4977.
• Test datasets at r4978 and r4984.
• Configuration namelists at r4979.
• Reference manual routines at r4980.
• Top level ropp_pp routines at r4982.
• Change to datum_hmsl in coordinates/earth.f90 at r4983.

[Ian Culverwell]

The tolerance on |difference in phase_L2_1| in ropp_io_fields_compare.f90 has been loosened (to match that of phase_L1) at r4990.

A fix to prevent calculation of ATAN2(0, 0) in ropp_pp_wopt_phase_and_amplitude.f90 has been made at r4991.

A cosmetic change ymax --> ylimit has been made to the WOPT tool at r4992.

[Ian Culverwell]

Another bug in ropp_io_fields_compare.f90 has been fixed at r4995.

[Ian Culverwell]

Remove ymax-dependency, and add extra diagnostics, at r4996.

[Ian Culverwell]

Chris and Sean squashed some of the bugs in this while I was away.

• Update ropp_io/data/eum_test_l.bfr_ref (as a result of IP precision (r4987) and software_version changes (r4997) at r5023.

• Include higher precision Fresnel integral estimates at r5024.

• Zero the 'edge' contribution to the 'linear' part of the integral giving the signal at the LEO at r5026 and r5028. (Not doing so is equivalent to assuming a linear variation in phase and amplitude out to 'infinity' (1E20 in practice). Omitting the edge contribution is equivalent to using a flat extrapolation. Keeping it in reveals a strange and nasty compiler sensitivity.)

[Ian Culverwell]

ROPP9.0_beta, which contains these latest changes, has been tagged at r5039.

[Ian Culverwell]

Thought for future developments.

• Why not do a cubic spline fit over each of the mini-screens in ropp_pp_wopt_propagate_to_leo.f90? The relevant integrals involve exponentials and Fresnel integrals, so the hard work has already been done. This might reduce the sensitivity to nsample (and the compiler!) because the 1st and 2nd derivatives of phase and amplitude would then be guaranteed to be continuous across the 'sample' boundaries.

[Ian Culverwell]

Actually a cubic spline might be a step too far: we would end up with integrals involving exp(i(1+x+x²+x³)), which is a bit involved. But quadratic splines (C1 piecewise quadratics) might be possible. They would involve integrals of quadratic*exp(i*quadratic), and they're not too hard.

But in fact a better method might be some sort of least squares fit which was nonetheless cts across the 'sample' boundaries. Can we do some sort of constrained least squares optimisation?

[Ian Culverwell]

CASE_05_ref1d.nc and CASE_51_ref1d.nc have been added to the repository at r5055. ROPP now includes the following examples:

CASE_51 = cat 1
CASE_07 = cat 2
CASE_12 = cat 3
CASE_05 = cat 4

For the definitions of the various categories, see the attached wop_studyv5.pdf. Basically, the higher the cat, the higher the maximum refractivity gradient.

Sean has passed all 55 profiles through his version of the WOPT; the results are plotted in the attached wopt55.pdf.

[Ian Culverwell]

wop_studyv5.pdf

[Ian Culverwell]

wopt55.pdf

[Ian Culverwell]

Sean's 'quick FSI' code has been plumbed in to ROPP at r5057. It is invoked by calling ropp_pp_wopt_tool with the -f option.

[Ian Culverwell]

The ROPP PP UG has been updated to describe this option at r5058 and r5059.

[cburrows]

A new test has been produced in the Test Folder for this new functionality: IT-PP-06 (note that the previous test 'IT-PP-06' has been renamed 'IT-APPS-01' to reflect the addition of the new 'ropp_apps' module).

The test runs the WOPT on a series of 1D refractivity profiles from the 55-profile dataset to produce simulated phases/amplitudes. The new FSI tool then retrieves bending angles from these, and the existing 'invert' tool produces refractivities which are compared against the original profiles. The comparisons are within acceptable bounds.

Closing ticket.

[Ian Culverwell]

I wondered if a continuous piecewise linear fit of the signal on the final screen might generate better results. It doesn't. See the attached Continuous_LSQfit.pdf for details. So abandon.

For the record, in case the idea proves useful elsewhere, the code to do this is contained within the attached ropp_pp/wopt/ropp_pp_wopt_propagate_to_leo.f90.

[Ian Culverwell]

ropp_pp_wopt_propagate_to_leo.f90

[Ian Culverwell]

Continuous_LSQfit.pdf

[Ian Culverwell]

Reclosing ticket.