test_Lstar failures

[jtniehof]

test_Lstar fails quite badly; I've got five failures of which three are off by a little under 10% and two are way off (10-20%). I don't feel reasonable just loosening these checks...L* of 5.7 vs. 4.2 is pretty big.

[drsteve]

Addressing just one of the test failures (I skipped the Hilton approximation tests on this run):

FAIL: testLgm_B_OP77_1 (__main__.Lstar_Tests)
This is a regression functional test for LstarVersusPA (regression)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "test_Lstar.py", line 56, in testLgm_B_OP77_1
    self.assertAlmostEqual(4.86203516, ans[90]['LMcIlwain'], places=4)
AssertionError: 4.86203516 != 4.550315708950278 within 4 places

The test is looking for a McIlwain L of 4.86203516 and an L* of 4.571340083084.
With a simple example in C I get very close to the expected values:

./Lstar 


                        Computing L* for, Date: 20101012, UT: 0
                =========================================================================
                Date (yyyymmdd):                         20101012
                UTC (hours):                             0
                Pitch Angle (deg.):                      90
                Initial Position, vin (Re):              < -4.31276, 1, 0.200202 >
                Mirror Mag. Field Strength, Bm (nT):     247.134

                L*, Dipole Approximation.
                  Magnetic Flux:                         -41207.336411971838970
                  L*:                                    4.566722246561874
                  L* (Using McIlwain M):                 4.752002242905313

                L*, Full Field.
                  Magnetic Flux:                         -40841.016242890211288
                  L*:                                    4.607683089836693
                  L* (Using McIlwain M):                 4.794624940018057

                Drift Orbit Type: CLOSED



L* = 4.60768
Pitch Angle: 90    McIlwain L  = 4.862654033   ( I, Bm, M = 0.704204 317.752 29937.3 )

So I believe there's something that's no longer working properly in the Python. I'll paste the C code in in the next comment so that the point of this comment doesn't get too lost. Basically, small changes are expected as there have been substantial library changes since the tests were last touched. The Python doesn't seem to be reproducing the C result though, and that's a problem.

(EDIT: Missed pasting the last line of output from the C code...)

[drsteve]

C code for the Lstar program in my previous comment:

#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <Lgm_CTrans.h>
#include <Lgm_LstarInfo.h>


int main( int argc, char *argv[] ){

    double           UTC, alpha, a;
    double           L, I, Bm, M;
    long int         Date;
    int              nAlpha, Kp;
    int              ZeroOneTwo = 1, ThreeFourFive=0, SixSevenEight=0;
    Lgm_Vector       Psm, P, v1, v2, v3;
    Lgm_LstarInfo *LstarInfo = InitLstarInfo(0);

    // Date and UTC
    Date       = 20101012;
    UTC        = 0.0;
    Lgm_Set_Coord_Transforms( Date, UTC, LstarInfo->mInfo->c );

    // Position in SM
    Psm.x = -4.2; Psm.y = 1.0; Psm.z = 1.0;
    Lgm_Convert_Coords( &Psm, &P, SM_TO_GSM, LstarInfo->mInfo->c );

    // pitch angles to compute
    LstarInfo->PitchAngle = 90.0;

    Lgm_SetLstarTolerances( 1, 24, LstarInfo );

    LstarInfo->VerbosityLevel = 1;
    LstarInfo->mInfo->VerbosityLevel = 0;
    LstarInfo->mInfo->Bfield        = Lgm_B_OP77;
    LstarInfo->mInfo->InternalModel = LGM_IGRF;
    LstarInfo->mInfo->Kp = 4;

    /*
     * Compute L*s, Is, Bms, Footprints, etc...
     * These quantities are stored in the MagEphemInfo Structure
     */
    LstarInfo->ShabanskyHandling = LGM_SHABANSKY_IGNORE;
    Lstar( &P, LstarInfo );
    printf("L* = %g\n", LstarInfo->LS);
    L = Lgm_McIlwain_L( Date, UTC, &P, LstarInfo->PitchAngle, 0, &I, &Bm, &M, LstarInfo->mInfo );
    printf("Pitch Angle: %g    McIlwain L  = %.10g   ( I, Bm, M = %g %g %g )\n", LstarInfo->PitchAngle, L, I, Bm, M);


    FreeLstarInfo( LstarInfo );

    return(0);

}

[jtniehof]

Thanks Steve. That does suggest a Python problem and gives me something to work with.

I was hoping to start de-duping e.g. the various field model wrappers.

This may have bearing on #17...most of those are small changes that I'd sort of expect from library updates, but there was one that was fairly large, so unless you know that related to a known substantial change, it might be suspect.

I did try going back to see if earlier revs worked, but that put me into problems with the build and it didn't seem worth that deep a dive.

[jtniehof]

The Python code, incidentally, is about twice as long as your C code. Which is bizarre. It is trying to do more stuff.

[drsteve]

I would expect some extra code for ctypes wrapping and user-friendly outputs, but yeah, this needs some rationalization. I likely won't have any time before the 21st, but if it's still an issue then I can take a look.

[jtniehof]

I think I've got a handle on this.

The Python code first does a simple closed field line trace:
trace, northern, southern, minB, Lsimple = Closed_Field.Closed_Field(MagEphemInfo, extended_out=True)

Then in calculating the L*:

MagEphemInfo.LstarInfo.contents.PitchAngle = pa
Ls_vec = Lgm_Vector.Lgm_Vector(*minB)
LS_Flag = Lgm_Lstar( pointer(Ls_vec), MagEphemInfo.LstarInfo)

So L* is calculated for a position at the minimum B along the field line, but with the pitch angle at this point equal to the pitch angle given at the observation point.

I'm not sure why the L* is calculated from that point rather than just handing it off to Lgm, but I imagine the PA has to be corrected. @balarsen?

[drsteve]

The Lstar calculation starts with a trace of the field line, so a lot of the work done in the Python in redundant. I also changed things a while back so Lstar should actually calculate Bm internally. In short, a lot of that setup can be removed now.

More importantly, yes, the L* will be wrong if it's being calculated from Pmin but using alpha from the initial location. The position and pitch angle should both be local (or both converted to (Pmin, apha_eq), but that's extra work and may increase the amount of numerical error).

[jtniehof]

Okay, I'll roll up my sleeves and start ripping stuff out.

[balarsen]

I suspect that the C library has changed a lot and gotten better as well as understanding of how to wrap it. Feel free to wholesale rip.

…

On Fri, Apr 12, 2019 at 12:59 PM Jon Niehof ***@***.***> wrote: Okay, I'll roll up my sleeves and start ripping stuff out. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#16 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAfB7U5JGL-dRngXyK90d_C3VpCBjQtWks5vgNeCgaJpZM4aiBPN> .

[jtniehof]

Changing either the pitch angle to the "fixed" one or the location to the Pgsm dramatically improves the match on Lstar. Changing the location is a better match, although still only good to two decimals. PA has to be changed to match the Hilton/McIlwain L.

But it does seem like pulling as much as possible out of the LstarInfo as populated by the C is better. We'll give that a shot.

[jtniehof]

Going from real position instead of minB, and using the same LStar tolerances, makes the L* match.

The tests fail on the Hilton and McIlwain because they use I0 as calculated from the initial L* trace, which appears always to be 0. It looks like something in that trace does a PA90 trace from the minimum point, since checking LStarInfo->I0 in the C always returns 0.

[jtniehof]

Working entirely in C, I can't even get the McIlwain values as in Steve's example to match McIlwain values for the same inputs in check_McIlwain_L. The L* calculation seems to be doing something to the CTrans or MagModelInfo structures; if I call Lgm_McIlwain_L before Lstar I get different output than calling it after.

[jtniehof]

https://github.com/jtniehof/LANLGeoMag/tree/lgmpy_tests

[jtniehof]

Tracked down part of the remaining problem: I need to call McIlwain_L before Lstar and before SetLstarTolerances, as SetLstarTolerances affects quite a bit in the MagModelInfo. Hopefully setting that up before calling the non-L* test will let those match, and then I can potentially track down what the actual Lstar calculation is doing to things.