From 673ef8b80c96000ceb9da50198f93deddc98a113 Mon Sep 17 00:00:00 2001 From: "Michael G. Henderson" Date: Wed, 15 Nov 2023 10:11:45 -0700 Subject: [PATCH 1/2] Added some error messages. Fixed some global B finding steps Include merge of prior SKM/GC changes with cherry-picked MGH updates --- libLanlGeoMag/Lgm_Trace.c | 76 ++++++++++++++++++++++++++++++++++----- 1 file changed, 68 insertions(+), 8 deletions(-) diff --git a/libLanlGeoMag/Lgm_Trace.c b/libLanlGeoMag/Lgm_Trace.c index 784a8b58..3fe9f59f 100644 --- a/libLanlGeoMag/Lgm_Trace.c +++ b/libLanlGeoMag/Lgm_Trace.c @@ -114,6 +114,7 @@ int Lgm_Trace( Lgm_Vector *u, Lgm_Vector *v1, Lgm_Vector *v2, Lgm_Vector *v3, do Lgm_Vector u_scale, P, gpp; +u_scale.x = u_scale.y = u_scale.z = 1.0; /* * Determine our initial geocentric radius in km. (u is assumed to be in @@ -135,6 +136,7 @@ int Lgm_Trace( Lgm_Vector *u, Lgm_Vector *v1, Lgm_Vector *v2, Lgm_Vector *v3, do // must be inside the Earth, which is no good -- bail with // LGM_INSIDE_EARTH error code + if ( Info->VerbosityLevel > 2 ) { printf("Initial point is inside the Earth\n"); } return( LGM_INSIDE_EARTH ); } else if ( Rinitial < WGS84_A ) { @@ -148,6 +150,7 @@ int Lgm_Trace( Lgm_Vector *u, Lgm_Vector *v1, Lgm_Vector *v2, Lgm_Vector *v3, do if ( GeodHeight < 0.0 ) { // inside the Earth, which is no good -- bail with error + if ( Info->VerbosityLevel > 2 ) { printf("Initial point is inside the Earth\n"); } return( LGM_INSIDE_EARTH ); } @@ -228,6 +231,9 @@ Info->Hmax = 0.10; flag2 = Lgm_TraceToEarth( u, v2, Height, -sgn, TOL1, Info ); +if (Info->VerbosityLevel == -100){ +printf("u = %g %g %g v2, Height, sgn, TOL1 = %g %g %g, %g, %g, %g\n", u->x, u->y, u->z, v2->x, v2->y, v2->z, Height, sgn, TOL1); +} Info->Snorth = Info->Trace_s; // save distance from u to northern footpoint location. double MIKEA = Info->Trace_s; Info->v2_final = *v2; @@ -247,6 +253,55 @@ double MIKEB = Info->Trace_s; if ( flag1 && flag2 ) { + /* + * Pre-trace the FL to find the true global Bmin. This is necessary for FLs + * that have multiple minima on the -- probably wasteful otherwise. + * We could also add code to detect number of minima here. + * + * 1. Copy the Info structure so we dont mess anything up. + * 2. Start from the southern footpoint (v1) and tracev up to north (we know the distance.) + * 3. save the various info as well as how far we are in s. + */ + double tSS, tBmin, s_approx; + Lgm_Vector u_approx; + int ii, iBmin; + Lgm_MagModelInfo *Info2 = Lgm_CopyMagInfo( Info ); + +//printf("Info2->s[%d] = %g\n", iBmin, Info2->s[iBmin]); + tSS = Info->Snorth + Info->Ssouth; + Lgm_TraceLine3( v1, tSS, 500, 1.0, 1e-7, FALSE, Info2 ); + tBmin = 9e99; + iBmin = -1; + for (ii=0; iinPnts; ii++){ + if (Info2->Bmag[ii] <= tBmin) { + tBmin = Info2->Bmag[ii]; + iBmin = ii; + } + } + if ( (iBmin >= 0) && (iBmin < Info2->nPnts) ){ + u_approx.x = Info2->Px[iBmin]; + u_approx.y = Info2->Py[iBmin]; + u_approx.z = Info2->Pz[iBmin]; + s_approx = Info2->s[iBmin]; // distancev along FL from v1 to approximate global Bmin. + } else { + u_approx = *u; + s_approx = 0.0; + } +//printf("u_approx = %g %g %g\n", u_approx.x, u_approx.y, u_approx.z); +//printf("Info2->s[%d] = %g\n", iBmin, Info2->s[iBmin]); +//printf("Info2->Bmag[%d] = %g\n", iBmin-1, Info2->Bmag[iBmin-1]); +//printf("Info2->Bmag[%d] = %g\n", iBmin, Info2->Bmag[iBmin]); +//printf("Info2->Bmag[%d] = %g\n", iBmin+1, Info2->Bmag[iBmin+1]); + Lgm_FreeMagInfo( Info2 ); + + // We should probably put in a test here to see if the final point from + // TraceLine3() is different than what the other routines gives. + // This could be the basis for dynamically setting the tolerances. + + + + + /* * Closed FL -- attempt to trace to Eq Plane. */ @@ -255,7 +310,7 @@ double MIKEB = Info->Trace_s; /* * Check the result of Lgm_TraceToMinBSurf to make sure it returns a closed field line */ - flag3 = Lgm_TraceToMinBSurf( u, v3, 0.1, TOL2, Info ); + flag3 = Lgm_TraceToMinBSurf( &u_approx, v3, 0.1, TOL2, Info ); if (flag3 != LGM_CLOSED) { printf("Major problem in Lgm_Trace(): northern and southern footpoints found, but TraceToMinBSurf does not return LGM_CLOSED, rather returns %d. Returning %d.\n", flag3, flag3); return(flag3); @@ -263,21 +318,26 @@ double MIKEB = Info->Trace_s; Info->v3_final = *v3; Info->Pmin = *v3; //Info->Smin = Info->Trace_s; // save location of Bmin. NOTE: Smin is measured from the southern footpoint. +//printf("About to call Info->Bfield( v3, &Bvec, Info );\n"); Info->Bfield( v3, &Bvec, Info ); +//printf("Done calling Info->Bfield( v3, &Bvec, Info );\n"); Info->Bvecmin = Bvec; Info->Bmin = Lgm_Magnitude( &Bvec ); - //printf("Bmin = %.15lf\n", Info->Bmin ); +// printf("Bmin = %.15lf\n", Info->Bmin ); +// printf("Info->Trace_s = %.15lf\n", Info->Trace_s ); +//printf("s_approx - Info->Trace_s = %g\n", s_approx - Info->Trace_s); +//exit(0); /* * Various FL arc lengths... * Snorth - distance from S/C to northern footpoint (set above) * Ssouth - distance from S/C to southern footpoint (set above) * Stotal - Total FL length ( Snorth + Ssouth ) (only compute for closed FL) - * Smin - distance from southern footpoint to S/C (Ssouth - what we - * got from Lgm_TraceToMinBSurf() because we started at S/C) + * Smin - Distance from southern footpoint to Pmin along the FL. (in Re). */ - Info->Stotal = Info->Snorth + Info->Ssouth; // Total FL length - Info->Smin = Info->Ssouth - Info->Trace_s; // length from south foot to S/C + Info->Stotal = Info->Snorth + Info->Ssouth; // Total FL length + //Info->Smin = Info->Ssouth - Info->Trace_s; // length from south foot to Pmin +Info->Smin = s_approx - Info->Trace_s; // length from south foot to Pmin Info->Trace_s = Info->Stotal; //printf("Info->Ssouth, Info->Trace_s = %g %g\n", Info->Ssouth, Info->Trace_s); @@ -363,8 +423,8 @@ double MIKEB = Info->Trace_s; */ Info->v1_final = *v1; Info->v2_final = *v2; - Info->Stotal = Info->Snorth + Info->Ssouth; // Total FL length within bounding box - Info->Trace_s = Info->Stotal; + Info->Stotal = Info->Snorth + Info->Ssouth; // Total FL length within bounding box + Info->Trace_s = Info->Stotal; return( LGM_OPEN_IMF ); From 95c5f8f3d7d24c639c1a3ec6db03e0fba9dbd34d Mon Sep 17 00:00:00 2001 From: Steve Morley Date: Mon, 20 May 2024 10:47:50 -0600 Subject: [PATCH 2/2] Update test results (ClosedField, McIlwain_L, Lstar) for expected change Minor channges in tracing result in different output. Changes are minor and expected. Regression tests updated. --- tests/check_ClosedField_01.expected | 2 +- tests/check_Lstar.expected | 10 +++++----- tests/check_McIlwain_L_01.expected | 4 ++-- tests/check_McIlwain_L_02.expected | 4 ++-- tests/check_McIlwain_L_03.expected | 4 ++-- tests/check_McIlwain_L_04.expected | 4 ++-- tests/check_McIlwain_L_05.expected | 4 ++-- tests/check_McIlwain_L_06.expected | 6 +++--- tests/check_McIlwain_L_07.expected | 4 ++-- tests/check_McIlwain_L_08.expected | 4 ++-- 10 files changed, 23 insertions(+), 23 deletions(-) diff --git a/tests/check_ClosedField_01.expected b/tests/check_ClosedField_01.expected index b62ea7ef..7f51ee02 100644 --- a/tests/check_ClosedField_01.expected +++ b/tests/check_ClosedField_01.expected @@ -43,4 +43,4 @@ # "units": "km" # } #} -2010-12-12T00:00:00.000000 T89 IGRF 2.000000 1.000000 2.000000 2.000000 LGM_CLOSED -0.062044 0.394560 0.930652 0.766626 0.335260 -0.571420 2.279626 2.809833 1.124498 +2010-12-12T00:00:00.000000 T89 IGRF 2.000000 1.000000 2.000000 2.000000 LGM_CLOSED -0.062044 0.394560 0.930652 0.766626 0.335260 -0.571420 2.279594 2.809784 1.124482 diff --git a/tests/check_Lstar.expected b/tests/check_Lstar.expected index 980e435d..5b582ec0 100644 --- a/tests/check_Lstar.expected +++ b/tests/check_Lstar.expected @@ -40,8 +40,8 @@ # "units": "km" # } #} -2010-10-12T00:00:00.000000 T89 IGRF 4.000000 -4.200000 1.000000 1.000000 90.000000 5.030699 5.030821 4.449974 -2010-10-12T00:00:00.000000 T89 IGRF 5.000000 -4.200000 1.000000 1.000000 90.000000 5.181867 5.182040 4.351110 -2010-10-12T00:00:00.000000 OP77 IGRF 4.000000 -4.200000 1.000000 1.000000 90.000000 4.862530 4.862611 4.564142 -2010-10-12T00:00:00.000000 OP77 IGRF 5.000000 -4.200000 1.000000 1.000000 90.000000 4.862530 4.862611 4.564142 -2010-10-12T00:00:00.000000 NULL CDIP 2.000000 -3.000000 0.000000 0.000000 90.000000 3.041069 3.041069 2.981678 +2010-10-12T00:00:00.000000 T89 IGRF 4.000000 -4.200000 1.000000 1.000000 90.000000 5.030914 5.031036 4.450113 +2010-10-12T00:00:00.000000 T89 IGRF 5.000000 -4.200000 1.000000 1.000000 90.000000 5.182028 5.182201 4.351208 +2010-10-12T00:00:00.000000 OP77 IGRF 4.000000 -4.200000 1.000000 1.000000 90.000000 4.862773 4.862854 4.564350 +2010-10-12T00:00:00.000000 OP77 IGRF 5.000000 -4.200000 1.000000 1.000000 90.000000 4.862773 4.862854 4.564350 +2010-10-12T00:00:00.000000 NULL CDIP 2.000000 -3.000000 0.000000 0.000000 90.000000 3.041360 3.041360 2.981987 diff --git a/tests/check_McIlwain_L_01.expected b/tests/check_McIlwain_L_01.expected index 61cbca89..4197f1e5 100644 --- a/tests/check_McIlwain_L_01.expected +++ b/tests/check_McIlwain_L_01.expected @@ -1,4 +1,4 @@ -4.000000000000042 -0.000000000000092 +4.000000000000001 +0.000000000000000 468.232303048585436 29966.867395109486097 diff --git a/tests/check_McIlwain_L_02.expected b/tests/check_McIlwain_L_02.expected index 09b816e0..5aaf503b 100644 --- a/tests/check_McIlwain_L_02.expected +++ b/tests/check_McIlwain_L_02.expected @@ -1,4 +1,4 @@ -11.532670707544700 -12.825003323634284 +11.535395826950113 +12.832099449687490 117.957459036319008 29966.867395109486097 diff --git a/tests/check_McIlwain_L_03.expected b/tests/check_McIlwain_L_03.expected index b222c2f3..441b8562 100644 --- a/tests/check_McIlwain_L_03.expected +++ b/tests/check_McIlwain_L_03.expected @@ -1,4 +1,4 @@ -3.176196186252976 -1.325574405043580 +3.176190063029837 +1.325559723369823 1698.990212267689230 30006.443502889760566 diff --git a/tests/check_McIlwain_L_04.expected b/tests/check_McIlwain_L_04.expected index a016897c..a24658ff 100644 --- a/tests/check_McIlwain_L_04.expected +++ b/tests/check_McIlwain_L_04.expected @@ -1,4 +1,4 @@ -57.762544970857483 -139.866472042890592 +57.759416444353988 +139.857851735183516 132.908118809885536 29966.867395109486097 diff --git a/tests/check_McIlwain_L_05.expected b/tests/check_McIlwain_L_05.expected index e77c240f..add85656 100644 --- a/tests/check_McIlwain_L_05.expected +++ b/tests/check_McIlwain_L_05.expected @@ -1,4 +1,4 @@ -7.000000000000013 -0.000000000000030 +7.000000154304486 +0.000000342779522 87.366960335596147 29966.867395109486097 diff --git a/tests/check_McIlwain_L_06.expected b/tests/check_McIlwain_L_06.expected index 0e6c6b6d..c592deb9 100644 --- a/tests/check_McIlwain_L_06.expected +++ b/tests/check_McIlwain_L_06.expected @@ -1,6 +1,6 @@ -7.967070727674994 -10.931271506834396 +7.968043287070341 +10.933854073716107 628.079558197794540 29966.867395109486097 628.079558197794540 -21.678768796737319 +21.664215975438278 diff --git a/tests/check_McIlwain_L_07.expected b/tests/check_McIlwain_L_07.expected index b40c717d..146a1733 100644 --- a/tests/check_McIlwain_L_07.expected +++ b/tests/check_McIlwain_L_07.expected @@ -1,4 +1,4 @@ -4.515135188966586 -0.691623196271838 +4.515151500328265 +0.691660560957078 401.967905765166449 29966.867395109486097 diff --git a/tests/check_McIlwain_L_08.expected b/tests/check_McIlwain_L_08.expected index 20afe780..27d18efb 100644 --- a/tests/check_McIlwain_L_08.expected +++ b/tests/check_McIlwain_L_08.expected @@ -1,4 +1,4 @@ -5.029342624138683 -0.911360044797447 +5.029303895910284 +0.911270903147744 302.289018570471740 29937.183805846714677