DRAFT: Shabansky fix

Examples/CoordQuickStart/CoordQuickStart.c

@@ -9,15 +9,28 @@ int main( ) {
     long int    Date;
     double      UTC, Lat, Lon, r;
 
-    Date = 20170101;    // Jan 1, 2000
+/*
+Date: 2014-08-29 15:32:13.811285
+GMST: 3.67943418927
+ITRF: [-13710.99888512  -3465.32964842  10000.        ]
+TEME: [ 10000.  10000.  10000.]
+*/
+
+
+    //Date = 20170101;    // Jan 1, 2000
+    Date = 20140829;    // Jan 1, 2000
     //Date = 20000101;    // Jan 1, 2000
-    UTC  = 12.0+0.0/60.0;         // Universal Time Coordinated (in decimal hours)
+    //UTC  = 12.0+0.0/60.0;         // Universal Time Coordinated (in decimal hours)
+    UTC  = 15.0+32.0/60.0+13.811285/3600.0;
     Ugsm.x = -6.6; Ugsm.y = 3.4; Ugsm.z = -2.3; // Set a vector in GSM coordinates
 
     // Set up all the necessary variables to do transformations for this Date and UTC
     /* Options for setting Sun/Moon pos are:  */
-    //Lgm_Set_CTrans_Options(LGM_EPH_HIGH_ACCURACY, LGM_PN_IAU76, c); /* Uses LGM high accuracy analytic solution */
-    Lgm_Set_CTrans_Options(LGM_EPH_DE, LGM_PN_IAU76, c); /* Uses JPL Development Ephemeris */
+    Lgm_Set_CTrans_Options(LGM_EPH_HIGH_ACCURACY, LGM_PN_IAU76, c); /* Uses LGM high accuracy analytic solution */
+    //Lgm_Set_CTrans_Options(LGM_EPH_DE, LGM_PN_IAU76, c); /* Uses JPL Development Ephemeris */
+    //Lgm_Set_CTrans_Options(LGM_EPH_LOW_ACCURACY, LGM_PN_IAU76, c); /* Uses Basic Ephemeris Calcs*/
+
+
     //Lgm_Set_CTrans_Options(LGM_EPH_LOW_ACCURACY, LGM_PN_IAU76, c); /* Same as NOT calling Set_CTrans_Options */
     Lgm_Set_Coord_Transforms( Date, UTC, c );
 

Examples/EopExample/QuickStartEop.c

@@ -10,6 +10,9 @@ int main( ) {
 
     Date = 20040812;                  // August 12, 2004
     UTC  = 12.34567;                  // Universal Time Coordinated (in decimal hours)
+Date = 20140829;    // Jan 1, 2000
+UTC  = 15.0+32.0/60.0+13.811285/3600.0;
+
     JD = Lgm_Date_to_JD( Date, UTC, c ); // Compute JD
     Ugsm.x = -6.6; Ugsm.y = 3.4; Ugsm.z = -2.3; // Set a vector in GSM coordinates
 

Examples/FluxTubeVolume/FluxTubeVolume.c

@@ -0,0 +1,132 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <Lgm_CTrans.h>
+#include <Lgm_QinDenton.h>
+#include <Lgm_ElapsedTime.h>
+#include <Lgm_FluxToPsd.h>
+#include <omp.h>
+
+
+int main(){
+
+    double              x, y, GeodLat, GeodLong, GeodHeight ;
+    double              UTC, JD;
+    long int            Date;
+    int                 Flag;
+    double              V, r, mlat, MLT, cl, sl, Phi, Sb;
+
+
+    double              Hdid, Hnext, s, Inc;
+    int                 reset=0, Ni;
+    Lgm_Vector          u_scale;
+    int                 EnhancedFlag;
+    Lgm_QinDentonOne    p;
+    Lgm_Vector          u, u_sm, v1, v2, v3, Bvec;
+    Lgm_MagModelInfo    *mInfo2, *mInfo = Lgm_InitMagInfo();
+    char    Filename[1024];
+    int     i, ii, nDivs, Status;
+    FILE    *fpout;
+    FILE    *fpout2 = fopen( "FL.txt", "w" );
+
+
+    Date       = 20130324;
+    UTC        = 23.462;
+    UTC        = 23.662;
+    JD = Lgm_Date_to_JD( Date, UTC, mInfo->c );
+    Lgm_Set_Coord_Transforms( Date, UTC, mInfo->c );
+    printf("Tilt = %g\n", mInfo->c->psi*DegPerRad );
+    //exit(0);
+
+    Lgm_get_QinDenton_at_JD( JD, &p, 1, 1 );
+    Lgm_set_QinDenton( &p, mInfo );
+
+    Lgm_MagModelInfo_Set_MagModel( LGM_IGRF, LGM_EXTMODEL_TS04, mInfo );
+    Lgm_MagModelInfo_Set_MagModel( LGM_CDIP, LGM_EXTMODEL_T87, mInfo );
+    Lgm_MagModelInfo_Set_MagModel( LGM_IGRF, LGM_EXTMODEL_T96, mInfo );
+    Lgm_MagModelInfo_Set_MagModel( LGM_IGRF, LGM_EXTMODEL_T89, mInfo );
+    Lgm_MagModelInfo_Set_MagModel( LGM_IGRF, LGM_EXTMODEL_TS04, mInfo );
+    Lgm_MagModelInfo_Set_MagModel( LGM_IGRF, LGM_EXTMODEL_OP77, mInfo );
+    Lgm_MagModelInfo_Set_MagModel( LGM_IGRF, LGM_EXTMODEL_T01S, mInfo );
+
+    char *s1, *s2, *s3, *s4;
+    Lgm_Get_ExtMagModelStrings( &s1, &s2, &s3, &s4, mInfo );
+    printf( "s1 = %s\n", s1 );
+    printf( "s2 = %s\n", s2 );
+    printf( "s3 = %s\n", s3 );
+    printf( "s4 = %s\n", s4 );
+    sprintf( Filename, "V_%s.txt", mInfo->ExtMagModelStr1 );
+    fpout = fopen( Filename, "w" );
+
+
+    //mInfo->fp = fpout2;
+    mInfo->SavePoints = FALSE;
+    MLT = 0.0;
+    r = 1.0 + 200.0/Re;
+    GeodHeight = 200.0;
+    Inc = 0.1;
+    Ni = (80.0-60.0)/Inc;
+mInfo->Lgm_MagStep_BS_atol = 1e-6;
+mInfo->Lgm_MagStep_BS_rtol = 0.0;
+mInfo->VerbosityLevel = 0;
+
+    #pragma omp parallel private(Status,mlat,Phi,cl,sl,u_sm,u,v1,v2,v3,mInfo2,Flag,nDivs,V,Sb)
+    #pragma omp for schedule(dynamic, 20)
+    for ( i=0; i<=Ni; i++ ) {
+
+        mlat = 60.0 + Inc*i;
+
+        mInfo2 = Lgm_CopyMagInfo( mInfo );
+
+        /*
+         * Trace from given SM position
+         */
+        printf("mlat, MLT = %g %g     ", mlat, MLT);
+        Phi = 15.0*(MLT-12.0)*RadPerDeg;
+        cl = cos( mlat * RadPerDeg ); sl = sin( mlat * RadPerDeg );
+        u_sm.x = r*cl*cos(Phi); u_sm.y = r*cl*sin(Phi); u_sm.z = r*sl;
+        Lgm_Convert_Coords( &u_sm, &u, SM_TO_GSM, mInfo2->c );
+        //printf("u_sm  = %g %g %g\n", u_sm.x, u_sm.y, u_sm.z );
+        //printf("u     = %g %g %g\n", u.x, u.y, u.z );
+
+        Flag = Lgm_Trace( &u, &v1, &v2, &v3, GeodHeight, 1e-7, 1e-7, mInfo2 );
+mInfo2->Bm = mInfo2->Ellipsoid_Footprint_Bn;
+
+
+        if ( Flag == LGM_CLOSED ) {
+            nDivs = mInfo2->Stotal/0.1;
+            if ( nDivs < 200 ) nDivs = 200;
+            mInfo2->Hmax = mInfo2->Stotal/((double)(nDivs));
+
+            Status = Lgm_TraceLine3( &v1, mInfo2->Stotal, nDivs, 1.0, 1e-7, FALSE, mInfo2 );
+
+            if ( !InitSpline( mInfo2 ) ) {
+                printf("Failed to init spline\n");
+                exit(1);
+            } else {
+
+                V = Lgm_FluxTubeVolume( mInfo2 );
+                Sb = SbIntegral_interped2( mInfo2, 0.0, mInfo2->Stotal );
+                printf( "south, north = %g %g , mlat, MLT, FluxTubeVolume, Sb = %g %g %g %g\n", mInfo2->Ssouth, mInfo2->Snorth, mlat, MLT, V, Sb );
+                fprintf( fpout, "%g %g\n", mlat, V );
+                fflush( fpout );
+                FreeSpline( mInfo2 );
+
+            }
+        } else {
+            printf("\n");
+        }
+
+        Lgm_FreeMagInfo( mInfo2 );
+
+
+    }
+
+    fclose( fpout );
+    fclose( fpout2 );
+    Lgm_FreeMagInfo( mInfo );
+
+
+
+
+}

Examples/FluxTubeVolume/Makefile

@@ -0,0 +1,9 @@
+# Very simple makefile illustrating how to use pkg-config to compile
+
+all: FluxTubeVolume
+
+FluxTubeVolume: FluxTubeVolume.c
+	gcc FluxTubeVolume.c `pkg-config --cflags --libs lgm` -o FluxTubeVolume
+
+clean:
+	rm FluxTubeVolume

Examples/KdTree/DstNNs.c

@@ -12,7 +12,7 @@ int main( ) {
     double              Time, JD, JDs, JDe, UTC, x, y, z, r, dist, *Dist, delta;
     long int            Date, KeepList[10000], nKeepList, Id;
     unsigned long int   n, *Idx, nSearches, Id_Old;
-    char                Filename[256];
+    char                Filename[256], IsoDateTime[40];
     Lgm_KdTree          *KdTree;
     Lgm_KdTreeData      *kNN;
     int                 K, Kgot, i, j, jj, d, D, k, ny, nm, nd, Keep;
@@ -30,26 +30,25 @@ int main( ) {
     /*
      * Read in Dst data from 1976 -> 2011
      */
-    Date = 19580101; UTC = 12.0; JDs = Lgm_Date_to_JD( Date, UTC, c );
-    Date = 19760101; UTC = 12.0; JDs = Lgm_Date_to_JD( Date, UTC, c );
-    Date = 20110101; UTC = 12.0; JDe = Lgm_Date_to_JD( Date, UTC, c );
-    Date = 19760102; UTC = 12.0; JDe = Lgm_Date_to_JD( Date, UTC, c );
+    Date = 19570101; UTC = 12.0; JDs = Lgm_Date_to_JD( Date, UTC, c );
+    Date = 20190630; UTC = 12.0; JDe = Lgm_Date_to_JD( Date, UTC, c );
     n = 0;
     for ( JD=JDs; JD<=JDe; JD += 1.0 ) {
         Date = Lgm_JD_to_Date( JD, &ny, &nm, &nd, &UTC );
         sprintf( Filename, "/home/mghenderson/Data/Dst/%4d/Dst_%ld.dat", ny, Date );
         printf("reading: %s\n", Filename );
-        //if ( (fp = fopen( Filename, "r" ) ) != NULL ) {
-        //    while ( fscanf( fp, "%lf %lf %lf %lf %lf\n", &Time, &dum1, &dum2, &dum3, &DstIndex ) != EOF ){
-        //        Dst[n] = DstIndex;
-        //        if ((Time > 0.0)&&(Time < 1.0)&&(Date == 20011124)) printf("n = %ld\n", n);
-        //        if ( n == 120228 ) printf("Date = %ld\n", Date);
-        //        if ( n == 250434 ) printf("Date = %ld\n", Date);
-        //        ++n;
-        //    }
-        //    fclose( fp );
-
-        //}
+        if ( (fp = fopen( Filename, "r" ) ) != NULL ) {
+            while ( fscanf( fp, "%19s %ld %lf %lf", IsoDateTime, &Date, &Time, &DstIndex ) != EOF ){
+//printf("IsoDateTime = %s\n", IsoDateTime);
+                Dst[n] = DstIndex;
+                //if ((Time > 0.0)&&(Time < 1.0)&&(Date == 20011124)) printf("n = %ld\n", n);
+                //if ( n == 120228 ) printf("Date = %ld\n", Date);
+                //if ( n == 250434 ) printf("Date = %ld\n", Date);
+                ++n;
+            }
+            fclose( fp );
+
+        }
     }
     printf("n = %ld\n", n);
 
@@ -59,14 +58,19 @@ int main( ) {
      */
     D = 96*2;
     LGM_ARRAY_2D( u, D, n, double );
+long int nj =0;
     for ( j=0; j<n-2*D; j++ ) {
+//printf("j=%d\n", j);
         for ( d=0; d<D; d++ ) {
             u[d][j] = Dst[j+d];
         }
+        ++nj;
     }
+printf("nj = %ld\n", nj);
     printf("Creating %d dimensional KdTree with %ld points\n", D, n);
     Lgm_ElapsedTimeInit( &t, 255, 150, 0 );
-    KdTree = Lgm_KdTree_Init( u, (void **)B, n, D );
+    //KdTree = Lgm_KdTree_Init( u, (void **)B, n, D );
+    KdTree = Lgm_KdTree_Init( u, (void **)B, nj, D );
     Lgm_PrintElapsedTime( &t );
 
 
@@ -75,7 +79,8 @@ int main( ) {
      */
     LGM_ARRAY_1D( q, D, double );
     for ( d=0; d<D; d++ ) {
-        q[d] = u[d][227016];
+        q[d] = u[d][227017];
+        //q[d] = u[d][226990];
     }
 
 
@@ -136,7 +141,7 @@ int main( ) {
     printf(" )\n");
 
     fp_out = fopen("DstNNs.dat", "w");
-if (nKeepList > 4) nKeepList = 4;
+if (nKeepList > 10) nKeepList = 10;
     for ( jj=0; jj<nKeepList; jj++ ) {
         j = KeepList[jj];
         fprintf( fp_out, "%d %g 3\n", 0, kNN[j].Position[0] );