LANLGeoMag version 2.0.0-alpha

.github/workflows/lgm-ci.yml

@@ -27,7 +27,7 @@ jobs:
     - name: compile_test
       run: |
         autoreconf -i
-        ./configure --disable-perl
+        ./configure
         make
         sudo make install
     - name: unit_tests

ChangeLog

@@ -1,34 +1,68 @@
+2024-02-29  Steve Morley <smorley@lanl.gov>
+    * Release: v2.0.0-alpha
+
+    * Continuous integration: Migrated to GitHub actions
+
+    * Installation: Perl extension disabled by default (--disable-perl argument removed from configure). To enable perl extension, use --enable-perl
+
+    * IGRF: Updated to DGRF12/IGRF13. Update extrapolation to use secular variation.
+
+    * Lgm_QinDenton: Extra checks in Lgm_QinDenton to reduce unnecessary output and check for inputs with standard QD naming convention ([1min, 5min, hour], in that priority order)
+
+    * New functions: Lgm_Terminator, FluxTubeVolume, Lgm_ReplaceSubString2, Lgm_B_FromScatteredData6 added to library
+
+    * Experimental Gaussian Process Regression functionality added to library
+
+    * Updates to enable second-order evaluation of fist adiabatic invariant
+
+    * magstep: Bulirsch-Stoer relative tolerance set to zero by default (only use absolute error)
+
+    * API change: Add RelaxTolerance argument to FindShellLine
+
+    * API change: Add ErrorStatus to all function variants of ComputeI_FromMltMlat
+
+    * TA16: Tsyganenko-Andreeva 2016 model added. Currently requires TA coeficient input files.
+
+    * TS04: Trap added for clean return if evaluating deep inside Earth
+
+    * JensonCain1960: New internal field model added
+
+    * Bug fixes: Correct altitude for lambda integrand in L* calculation, Correct angle in Lgm_EarthEclipse, Correct intersection calculation for Lgm_Ellipse
+
+    * Other: Preliminary version of field-aligned (PQB) coordinates, expanded KdTree functionality, updates to various examples, tests updated for expected changes, numerous minor documentation updates, etc.
+
+
 2018-03-09  Steve Morley <smorley@lanl.gov>
-        * Release: v1.5.16
+    * Release: v1.5.16
 
-        * tests: revised and expanded test suite, github repo set to use Travis
-        CI
+    * tests: revised and expanded test suite, github repo set to use Travis
+    CI
 
-        * Lgm_QinDenton: updated defaults, changed handling of missing data
-        files
+    * Lgm_QinDenton: updated defaults, changed handling of missing data
+    files
 
-        * Lgm_Ellipsoid: added shapes, intersections, etc.
+    * Lgm_Ellipsoid: added shapes, intersections, etc.
 
-        * Tools: Command line utils for Magnetic Ephemeris files and Last 
-        Closed Drift Shell calculation
+    * Tools: Command line utils for Magnetic Ephemeris files and Last 
+    Closed Drift Shell calculation
 
-        * ViewDriftShell: Updates to improve portability and UI
+    * ViewDriftShell: Updates to improve portability and UI
 
-        * ComputeLstar: Behavior in Shabansky orbits now set by flag in 
-        LstarInfo struct
+    * ComputeLstar: Behavior in Shabansky orbits now set by flag in 
+    LstarInfo struct
 
-        * Lgm_CTrans: Added JPL DE421 to obtain Sun vector (now have DE, high 
-        accuracy and low accuracy)
+    * Lgm_CTrans: Added JPL DE421 to obtain Sun vector (now have DE, high 
+    accuracy and low accuracy)
 
-        * TS07: TS07 magnetic field model fully enabled. Requires coefficient
-        files in same location/format as IRBEMlib.
+    * TS07: TS07 magnetic field model fully enabled. Requires coefficient
+    files in same location/format as IRBEMlib.
 
-        * Other: More changes since this file was last updated than can be 
-        encapsulated here. OSS release and GitHub hosting means that commit 
-        history is now publicly available.
+    * Other: More changes since this file was last updated than can be 
+    encapsulated here. OSS release and GitHub hosting means that commit 
+    history is now publicly available.
 
 2011-07-08  Mike Henderson <mghenderson@lanl.gov>
-        * configure.ac: added LGM_INDEX_DATA_DIR to locate where to find index data (e.g. Kp, Dst, AE). 
+    * configure.ac: added LGM_INDEX_DATA_DIR to locate where to find index data (e.g. Kp, Dst, AE). 
 
 2011-04-25  Mark Galassi  <markgalassi@lanl.gov>
 	* configure.ac, lanlgeomag.spec.in: started maintaining a .spec.in

Examples/Bfield/GradBvec.c

@@ -0,0 +1,106 @@
+#include <Lgm_MagModelInfo.h>
+
+/* BAL 04-Jan-2011 modified for more cases */
+
+int main(){
+
+
+    long int            Date;
+    double              UTC;
+    double              GradBvec[3][3];
+    Lgm_Vector          u, B, CurlB, GradB;
+    Lgm_MagModelInfo    *mInfo;
+    int                 i, j;
+
+
+    Date = 20050831;
+    UTC  = 9.0;
+
+    mInfo = Lgm_InitMagInfo( );
+    Lgm_Set_Coord_Transforms( Date, UTC, mInfo->c );
+
+    //mInfo->Bfield = Lgm_B_T89;
+    //mInfo->Bfield = Lgm_B_TS04;
+    //mInfo->Bfield = Lgm_B_OP77;
+    //mInfo->Bfield = Lgm_B_T89;
+    mInfo->Bfield = Lgm_B_TS04;
+    mInfo->P      = 4.1011111111111118;
+    mInfo->Dst    = 7.7777777777777777;
+    mInfo->By     = 3.7244444444444444;
+    mInfo->Bz     = -0.12666666666666665;
+    mInfo->W[0]   = 0.12244444444444445;
+    mInfo->W[1]   = 0.2514;
+    mInfo->W[2]   = 0.089266666666666661;
+    mInfo->W[3]   = 0.047866666666666668;
+    mInfo->W[4]   = 0.22586666666666666;
+    mInfo->W[5]   = 1.0461333333333334;
+
+
+    printf("%5s", "Kp");
+    printf(" %10s", "Ux (Re)");
+    printf(" %10s", "Uy (Re)");
+    printf(" %10s", "Uz (Re)");
+    printf(" %10s", "Bx (nT)");
+    printf(" %10s", "By (nT)");
+    printf(" %10s", "Bz (nT)");
+    printf(" %10s", "Bmag (nT)");
+    printf("   %15s", "(\u2207B)_x");
+    printf("   %15s", "(\u2207B)_y");
+    printf("   %15s", "(\u2207B)_z\n");
+
+
+    for (i=0; i<=5; i++) {
+        mInfo->Kp = i;
+        u.x = -6.6; u.y =  0.0;  u.z =  0.0;
+        mInfo->Bfield( &u, &B, mInfo );
+        Lgm_CurlB( &u, &CurlB, LGM_DERIV_SIX_POINT, 1e-3, mInfo );
+        Lgm_GradB( &u, &GradB, LGM_DERIV_SIX_POINT, 1e-3, mInfo );
+        Lgm_GradBvec( &u, GradBvec, LGM_DERIV_SIX_POINT, 1e-3, mInfo );
+
+        printf("\n\n");
+        printf( "\t Kp        = %5d\n", mInfo->Kp);
+        printf( "\t u         = %10g %10g %10g\n", u.x, u.y, u.z );
+        printf( "\t B         = %10g %10g %10g\n", B.x, B.y, B.z );
+        printf( "\t |B|       = %10g\n", Lgm_Magnitude( &B ) );
+        printf( "\t CurlB     = %15g %15g %15g\n", CurlB.x, CurlB.y, CurlB.z );
+        printf( "\t GradB     = %15g %15g %15g\n", GradB.x, GradB.y, GradB.z );
+        printf( "\t             (%8g %8g %8g)\n", GradBvec[0][0], GradBvec[0][1], GradBvec[0][2] );
+        printf( "\t Grad Bvec = (%8g %8g %8g)\n", GradBvec[1][0], GradBvec[1][1], GradBvec[1][2] );
+        printf( "\t             (%8g %8g %8g)\n", GradBvec[2][0], GradBvec[2][1], GradBvec[2][2] );
+        double Bmag = Lgm_Magnitude( &B );
+        printf("\n\t Show that GradB can be derived from GradBvec...\n");
+        printf("\t GradB_X: %g\n", (B.x*GradBvec[0][0] + B.y*GradBvec[1][0] + B.z*GradBvec[2][0])/Bmag );
+        printf("\t GradB_Y: %g\n", (B.x*GradBvec[0][1] + B.y*GradBvec[1][1] + B.z*GradBvec[2][1])/Bmag );
+        printf("\t GradB_Z: %g\n", (B.x*GradBvec[0][2] + B.y*GradBvec[1][2] + B.z*GradBvec[2][2])/Bmag );
+        printf("\n\t Show that CurlB can be derived from GradBvec...\n");
+        printf("\t CurlB_X: %g\n", GradBvec[2][1] - GradBvec[1][2] );
+        printf("\t CurlB_Y: %g\n", GradBvec[0][2] - GradBvec[2][0] );
+        printf("\t CurlB_Z: %g\n", GradBvec[1][0] - GradBvec[0][1] );
+    }
+
+/*
+    for (j=0; j<100; j++){
+        mInfo->Kp = 3;
+        for (i=0; i<13; i++) {
+            u.x = -1. - (double)i * 0.5;
+            u.y =  0.0;  u.z =  0.0;
+            mInfo->Bfield( &u, &B, mInfo );
+            Lgm_GradB( &u, &GradB, LGM_DERIV_SIX_POINT, 1e-3, mInfo );
+            printf( "%5d", mInfo->Kp);
+            printf( " %10g %10g %10g", u.x, u.y, u.z );
+            printf( " %10g %10g %10g", B.x, B.y, B.z );
+            printf( " %10g", Lgm_Magnitude( &B ) );
+            printf( " %15g %15g %15g \n", GradB.x, GradB.y, GradB.z );
+        }
+    }
+*/
+
+
+
+
+
+    Lgm_FreeMagInfo( mInfo );
+
+
+    exit(0);
+}

Examples/Bfield/GradBvec2.c

@@ -0,0 +1,115 @@
+#include <Lgm_MagModelInfo.h>
+#include <Lgm_DynamicMemory.h>
+
+/* BAL 04-Jan-2011 modified for more cases */
+
+int main(){
+
+
+    long int            Date;
+    double              UTC;
+    double              GradBvec[4][4];
+    Lgm_Vector          u, B, CurlB, GradB;
+    Lgm_MagModelInfo    **mInfo;
+    int                 i, j, n;
+
+    j = 167;
+    n = 8640;
+    LGM_ARRAY_1D( mInfo, n, Lgm_MagModelInfo * );
+
+    /*
+     * Initialize an array of mInfo[i] structures...
+     */
+    Date = 20050831;
+    for ( i=0; i<8640; ++i ) {
+        UTC  = (double)i*10.0/3600.0;
+        mInfo[i] = Lgm_InitMagInfo( );
+        mInfo[i]->Kp = i%5;
+        Lgm_Set_Coord_Transforms( Date, UTC, mInfo[i]->c );
+        //printf("%8ldT%02d%02d%02d %g\n", mInfo[i]->c->UTC.Date, mInfo[i]->c->UTC.Hour, mInfo[i]->c->UTC.Minute, (int)mInfo[i]->c->UTC.Second, mInfo[i]->c->psi*DegPerRad );
+        Lgm_MagModelInfo_Set_MagModel( LGM_IGRF, LGM_EXTMODEL_T89, mInfo[i] );
+        mInfo[i]->P      = 4.1011111111111118;
+        mInfo[i]->Dst    = 7.7777777777777777;
+        mInfo[i]->By     = 3.7244444444444444;
+        mInfo[i]->Bz     = -0.12666666666666665;
+        mInfo[i]->W[0]   = 0.12244444444444445;
+        mInfo[i]->W[1]   = 0.2514;
+        mInfo[i]->W[2]   = 0.089266666666666661;
+        mInfo[i]->W[3]   = 0.047866666666666668;
+        mInfo[i]->W[4]   = 0.22586666666666666;
+        mInfo[i]->W[5]   = 1.0461333333333334;
+    }
+
+
+    printf("%5s", "Kp");
+    printf(" %10s", "Ux (Re)");
+    printf(" %10s", "Uy (Re)");
+    printf(" %10s", "Uz (Re)");
+    printf(" %10s", "Bx (nT)");
+    printf(" %10s", "By (nT)");
+    printf(" %10s", "Bz (nT)");
+    printf(" %10s", "Bmag (nT)");
+    printf("   %15s", "(\u2207B)_x");
+    printf("   %15s", "(\u2207B)_y");
+    printf("   %15s", "(\u2207B)_z\n");
+
+
+    for (i=0; i<=5; i++) {
+        u.x = -6.6; u.y =  0.3;  u.z =  -2.0;
+        mInfo[j]->Bfield( &u, &B, mInfo[j] );
+        Lgm_CurlB( &u, &CurlB, LGM_DERIV_SIX_POINT, 1e-3, mInfo[j] );
+        Lgm_GradB( &u, &GradB, LGM_DERIV_SIX_POINT, 1e-3, mInfo[j] );
+        Lgm_GradBvec2( j, &u, GradBvec, LGM_DERIV_SIX_POINT, 1e-3, n, 10.0, &mInfo[0] );
+
+        printf("\n\n");
+        printf( "\t Kp        = %5d\n", mInfo[i]->Kp);
+        printf( "\t u         = %10g %10g %10g\n", u.x, u.y, u.z );
+        printf( "\t B         = %10g %10g %10g\n", B.x, B.y, B.z );
+        printf( "\t |B|       = %10g\n", Lgm_Magnitude( &B ) );
+        printf( "\t CurlB     = %15g %15g %15g\n", CurlB.x, CurlB.y, CurlB.z );
+        printf( "\t GradB     = %15g %15g %15g\n", GradB.x, GradB.y, GradB.z );
+        printf( "\t             (%8g %8g %8g %8g)\n", GradBvec[0][0], GradBvec[0][1], GradBvec[0][2], GradBvec[0][3] );
+        printf( "\t Grad Bvec = (%8g %8g %8g %8g)\n", GradBvec[1][0], GradBvec[1][1], GradBvec[1][2], GradBvec[1][3] );
+        printf( "\t             (%8g %8g %8g %8g)\n", GradBvec[2][0], GradBvec[2][1], GradBvec[2][2], GradBvec[2][3] );
+        printf( "\t             (%8g %8g %8g %8g)\n", GradBvec[3][0], GradBvec[3][1], GradBvec[3][2], GradBvec[3][3] );
+        double Bmag = Lgm_Magnitude( &B );
+        printf("\n\t Show that GradB can be derived from GradBvec...\n");
+        printf("\t GradB_T: %g\n", (B.x*GradBvec[1][0] + B.y*GradBvec[2][0] + B.z*GradBvec[3][0])/Bmag );
+        printf("\t GradB_X: %g\n", (B.x*GradBvec[1][1] + B.y*GradBvec[2][1] + B.z*GradBvec[3][1])/Bmag );
+        printf("\t GradB_Y: %g\n", (B.x*GradBvec[1][2] + B.y*GradBvec[2][2] + B.z*GradBvec[3][2])/Bmag );
+        printf("\t GradB_Z: %g\n", (B.x*GradBvec[1][3] + B.y*GradBvec[2][3] + B.z*GradBvec[3][3])/Bmag );
+        printf("\n\t Show that CurlB can be derived from GradBvec...\n");
+        printf("\t CurlB_X: %g\n", GradBvec[3][2] - GradBvec[2][3] );
+        printf("\t CurlB_Y: %g\n", GradBvec[1][3] - GradBvec[3][1] );
+        printf("\t CurlB_Z: %g\n", GradBvec[2][1] - GradBvec[1][2] );
+    }
+
+/*
+    for (j=0; j<100; j++){
+        mInfo[i]->Kp = 3;
+        for (i=0; i<13; i++) {
+            u.x = -1. - (double)i * 0.5;
+            u.y =  0.0;  u.z =  0.0;
+            mInfo[i]->Bfield( &u, &B, mInfo[i] );
+            Lgm_GradB( &u, &GradB, LGM_DERIV_SIX_POINT, 1e-3, mInfo[i] );
+            printf( "%5d", mInfo[i]->Kp);
+            printf( " %10g %10g %10g", u.x, u.y, u.z );
+            printf( " %10g %10g %10g", B.x, B.y, B.z );
+            printf( " %10g", Lgm_Magnitude( &B ) );
+            printf( " %15g %15g %15g \n", GradB.x, GradB.y, GradB.z );
+        }
+    }
+*/
+
+
+
+
+
+    for ( i=0; i<8640; ++i ) {
+        Lgm_FreeMagInfo( mInfo[i] );
+    }
+    LGM_ARRAY_1D_FREE( mInfo );
+
+
+    exit(0);
+}

Examples/Bfield/Makefile

@@ -2,7 +2,7 @@ HDF5FLAGS = `pkg-config hdf5 --cflags --libs 2>/dev/null`
 LGMFLAGS = `pkg-config lgm --cflags --libs`
 
 
-all: B Bcompare B2 GradB GradB2 B_Cross_GradB_Over_B CurlB CurlB2 DivB
+all: B Bcompare B2 GradBvec GradBvec2 GradB GradB2 B_Cross_GradB_Over_B CurlB CurlB2 DivB
 
 
 B: B.c
@@ -14,6 +14,12 @@ Bcompare: Bcompare.c
 B2: B2.c
 	gcc B2.c $(LGMFLAGS) $(HDF5FLAGS) -o B2 -Wall
 
+GradBvec: GradBvec.c
+	gcc GradBvec.c $(LGMFLAGS) $(HDF5FLAGS) -o GradBvec -Wall
+
+GradBvec2: GradBvec2.c
+	gcc GradBvec2.c $(LGMFLAGS) $(HDF5FLAGS) -o GradBvec2 -Wall
+
 GradB: GradB.c
 	gcc GradB.c $(LGMFLAGS) $(HDF5FLAGS) -o GradB -Wall
 
@@ -34,4 +40,4 @@ DivB: DivB.c
 
 
 clean:
-	rm *~ *.o B Bcompare B2 GradB GradB2 B_Cross_GradB_Over_B CurlB CurlB2 DivB
+	rm *~ *.o B Bcompare B2 GradBvec GradBvec2 GradB GradB2 B_Cross_GradB_Over_B CurlB CurlB2 DivB

Examples/Geodetic/Geodetic2.c

@@ -0,0 +1,39 @@
+// Geodetic.c
+#include <Lgm_CTrans.h>
+int main(){
+
+    double      GeodLat, GeodLong, GeodHeight, r;
+    double      GeoLat, GeoLong, GeoHeight;
+    Lgm_Vector  v;
+
+    GeodLat    = 35.0 + 53.0/60.0 + 17.0/3600.0;   // degrees
+    GeodLong   = -106.0 - 18.0/60.0 - 23.0/3600.0; // degrees
+    GeodHeight = 20.0; // km
+GeoLat    = 56.02;   // degrees
+GeoLong   = -141.42; // degrees
+GeoHeight = 850.0; // km
+v.x = (1.0+GeoHeight/WGS84_A)*cos( GeoLong * RadPerDeg )*cos( GeoLat*RadPerDeg );
+v.y = (1.0+GeoHeight/WGS84_A)*sin( GeoLong * RadPerDeg )*cos( GeoLat*RadPerDeg );
+v.z = (1.0+GeoHeight/WGS84_A)*sin( GeoLat*RadPerDeg );
+    Lgm_WGS84_to_GEOD( &v, &GeodLat, &GeodLong, &GeodHeight );
+
+//    Lgm_GEOD_to_WGS84( GeodLat, GeodLong, GeodHeight, &v );
+
+    printf( "            Geodetic Position\n");
+    printf( "   ===================================\n");
+    printf( "      GeodLat    : %lf degrees\n", GeodLat );
+    printf( "      GeodLong   : %lf degrees\n", GeodLong );
+    printf( "      GeodHeight : %lf km\n\n", GeodHeight );
+
+    printf( "            Geocentric Position\n");
+    printf( "   ===================================\n");
+    printf( "      X         : %-15lf meters\n", v.x*WGS84_A*1000.0);   // meters
+    printf( "      Y         : %-15lf meters\n", v.y*WGS84_A*1000.0);   // meters
+    printf( "      Z         : %-15lf meters\n", v.z*WGS84_A*1000.0); // meters
+    r = sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
+    printf( "      Radius    : %lf Re (%lf km)\n", r, r*WGS84_A );
+    printf( "      Latitude  : %lf degrees\n", DegPerRad*asin( v.z/r ) );
+    printf( "      Longitude : %lf degrees\n\n", DegPerRad*atan2( v.y, v.x ) );
+
+    exit(0);
+}

Examples/Geodetic/Makefile

@@ -4,5 +4,8 @@
 Geodetic: Geodetic.c
 	gcc Geodetic.c `pkg-config --cflags --libs lgm` -o Geodetic
 
+Geodetic2: Geodetic2.c
+	gcc Geodetic2.c `pkg-config --cflags --libs lgm` -o Geodetic2
+
 clean:
-	rm Geodetic
+	rm Geodetic Geodetic2