Improve global B finding in Lgm_Trace

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; ii<Info2->nPnts; 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,29 +310,34 @@ 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);
                 }
         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 );
 

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

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

tests/check_McIlwain_L_01.expected

@@ -1,4 +1,4 @@
-4.000000000000042
-0.000000000000092
+4.000000000000001
+0.000000000000000
 468.232303048585436
 29966.867395109486097

tests/check_McIlwain_L_02.expected

@@ -1,4 +1,4 @@
-11.532670707544700
-12.825003323634284
+11.535395826950113
+12.832099449687490
 117.957459036319008
 29966.867395109486097

tests/check_McIlwain_L_03.expected

@@ -1,4 +1,4 @@
-3.176196186252976
-1.325574405043580
+3.176190063029837
+1.325559723369823
 1698.990212267689230
 30006.443502889760566

tests/check_McIlwain_L_04.expected

@@ -1,4 +1,4 @@
-57.762544970857483
-139.866472042890592
+57.759416444353988
+139.857851735183516
 132.908118809885536
 29966.867395109486097

tests/check_McIlwain_L_05.expected

@@ -1,4 +1,4 @@
-7.000000000000013
-0.000000000000030
+7.000000154304486
+0.000000342779522
 87.366960335596147
 29966.867395109486097

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

tests/check_McIlwain_L_07.expected

@@ -1,4 +1,4 @@
-4.515135188966586
-0.691623196271838
+4.515151500328265
+0.691660560957078
 401.967905765166449
 29966.867395109486097

tests/check_McIlwain_L_08.expected

@@ -1,4 +1,4 @@
-5.029342624138683
-0.911360044797447
+5.029303895910284
+0.911270903147744
 302.289018570471740
 29937.183805846714677