Feature expansion and some cleanup

Examples/Bfield/B.c

@@ -35,13 +35,14 @@ int main(){
      * Or there is a setter routine ...
      */
     //Lgm_MagModelInfo_Set_MagModel( LGM_EDIP, LGM_EXTMODEL_T89, mInfo );
+    Lgm_MagModelInfo_Set_MagModel( LGM_EDIP, LGM_EXTMODEL_TS04, mInfo );
 
 
 
     /*
      * For TS07, the coeffs need to be initialized for each new time...
      */
-    Lgm_SetCoeffs_TS07( Date, UTC, &mInfo->TS07_Info );
+    //Lgm_SetCoeffs_TS07( Date, UTC, &mInfo->TS07_Info );
 
 
 
@@ -83,7 +84,7 @@ int main(){
 
     for (i=0; i<=5; i++) {
         mInfo->Kp = i;
-        u.x = -6.6; u.y =  0.0;  u.z =  0.0;
+        u.x = -5.8; u.y =  0.0;  u.z =  0.0;
         //Lgm_Convert_Coords( &u, &ugsm, GEO_TO_GSM, mInfo->c );
         Lgm_Convert_Coords( &u, &ugsm, SM_TO_GSM, mInfo->c );
         mInfo->Bfield( &ugsm, &B, mInfo );
@@ -112,6 +113,35 @@ int main(){
 
 
 
+    double Tilt, Phi, Theta, MLT, R;
+    R = 6.6;
+    FILE *fp1 = fopen( "Tilt1.txt", "w" );
+    FILE *fp2 = fopen( "Tilt2.txt", "w" );
+    for (MLT = 18.0; MLT <= 30.0; MLT += 0.1 ) {
+
+        Phi   = 180.0 + (MLT-18.0)*15.0;
+        Theta = 0.0;
+
+        u.x = R * cos( Phi*RadPerDeg ) * cos( Theta *RadPerDeg );
+        u.y = R * sin( Phi*RadPerDeg ) * cos( Theta *RadPerDeg );
+        u.z = R * sin( Theta *RadPerDeg );
+printf("u = %g %g %g\n", u.x, u.y, u.z);
+
+        Lgm_Convert_Coords( &u, &ugsm, SM_TO_GSM, mInfo->c );
+
+        mInfo->Bfield( &ugsm, &B, mInfo );
+
+        Tilt = DegPerRad * atan2( B.z, sqrt(B.x*B.x + B.y*B.y) );
+        fprintf( fp1, "%g %g\n", MLT, Tilt );
+
+
+        Tilt = DegPerRad * asin( B.z/sqrt(B.x*B.x + B.y*B.y + B.z*B.z) );
+        fprintf( fp2, "%g %g\n", MLT, Tilt );
+
+    }
+    fclose(fp1);
+    fclose(fp2);
+
 
 
 

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 B2 GradB GradB2 B_Cross_GradB_Over_B CurlB CurlB2 DivB
+all: B B2 GradBvec GradBvec2 GradB GradB2 B_Cross_GradB_Over_B CurlB CurlB2 DivB
 
 
 B: B.c
@@ -11,6 +11,12 @@ B: B.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
 
@@ -31,4 +37,4 @@ DivB: DivB.c
 
 
 clean:
-	rm *~ *.o B B2 GradB GradB2 B_Cross_GradB_Over_B CurlB CurlB2 DivB
+	rm *~ *.o B 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);
+}