Skip to content
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 11 additions & 11 deletions libLanlGeoMag/ComputeI_FromMltMlat.c
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ double ComputeI_FromMltMlat1( double Bm, double MLT, double mlat, double *r, dou

int reset=1, reset2;

double I, Phi, cl, sl, rat, SS1, SS2, SS, Sn, Ss, Htry, Hdid, Hnext, Bs, Be, s, sgn;
double I1, Phi, cl, sl, rat, SS1, SS2, SS, Sn, Ss, Htry, Hdid, Hnext, Bs, Be, s, sgn;
Lgm_Vector w, u, Pmirror1, Pmirror2, v1, v2, v3, Bvec, P, Ps, u_scale, Bvectmp, Ptmp;
double stmp, Btmp;

Expand All @@ -52,8 +52,8 @@ double ComputeI_FromMltMlat1( double Bm, double MLT, double mlat, double *r, dou
* Couldnt get a valid Bm. (The bracket is pretty huge,so
* we probably ought to believe there really isnt a valid one.)
*/
if (LstarInfo->VerbosityLevel > 1) printf("\t%sNo Bm found: setting I to 9e99%s\n", LstarInfo->PreStr, LstarInfo->PostStr);
I = 9e99;
if (LstarInfo->VerbosityLevel > 1) printf("\t%sNo Bm found: setting I1 to 9e99%s\n", LstarInfo->PreStr, LstarInfo->PostStr);
I1 = 9e99;

} else {

Expand Down Expand Up @@ -195,7 +195,7 @@ double ComputeI_FromMltMlat1( double Bm, double MLT, double mlat, double *r, dou
*
* Trace from Pm_North to Pm_South
*/
I = 9e99;
I1 = 9e99;
//LstarInfo->mInfo->Hmax = 10.0;
//LstarInfo->mInfo->Hmax = 0.1;
LstarInfo->mInfo->Hmax = 0.1;
Expand Down Expand Up @@ -307,17 +307,17 @@ if (0==1){
* Do I integral with interped integrand.
*/
//printf("I = %g\n", I);
I = Iinv_interped( LstarInfo->mInfo );
I1 = Iinv_interped( LstarInfo->mInfo );
Comment thread
drsteve marked this conversation as resolved.
//printf("I = %g Sm_South, Sm_North = %g %g\n", I, LstarInfo->mInfo->Sm_South, LstarInfo->mInfo->Sm_North);
// if (LstarInfo->VerbosityLevel > 1) printf("\t\t%s Integral Invariant, I (interped): %15.8g I-I0: %15.8g [a,b]: %.15g %.15g mlat: %12.8lf (nCalls = %d)%s\n", LstarInfo->PreStr, I, I-I0, LstarInfo->mInfo->Sm_South, LstarInfo->mInfo->Sm_North, mlat, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, LstarInfo->PostStr );
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: %13.6g I0: %13.6g I-I0: %13.6g [Sa,Sb]: %.8g %.8g (nCalls = %d)%s\n", LstarInfo->PreStr, mlat, I, I0, I-I0, LstarInfo->mInfo->Sm_South, LstarInfo->mInfo->Sm_North, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: %13.6g I0: %13.6g I1-I0: %13.6g [Sa,Sb]: %.8g %.8g (nCalls = %d)%s\n", LstarInfo->PreStr, mlat, I1, I0, I1-I0, LstarInfo->mInfo->Sm_South, LstarInfo->mInfo->Sm_North, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, LstarInfo->PostStr );
}
FreeSpline( LstarInfo->mInfo );

} else {

I = 9e99;
I1 = 9e99;

}

Expand All @@ -332,19 +332,19 @@ if (0==1){
/*
* Do full blown I integral.
*/
I = Iinv( LstarInfo->mInfo );
if (LstarInfo->VerbosityLevel > 1) printf("\t\t%s Integral Invariant, I (full integral): %15.8g I-I0: %15.8g mlat: %12.8lf (nCalls = %d)%s\n", LstarInfo->PreStr, I, I-I0, mlat, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, LstarInfo->PostStr );
I1 = Iinv( LstarInfo->mInfo );
if (LstarInfo->VerbosityLevel > 1) printf("\t\t%s Integral Invariant, I1 (full integral): %15.8g I1-I0: %15.8g mlat: %12.8lf (nCalls = %d)%s\n", LstarInfo->PreStr, I1, I1-I0, mlat, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, LstarInfo->PostStr );
}

} else {
I = 9e99;
I1 = 9e99;
}

}




return( I );
return( I1 );

}
40 changes: 20 additions & 20 deletions libLanlGeoMag/ComputeI_FromMltMlat2.c
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ double ComputeI_FromMltMlat2( double Bm, double MLT, double mlat, double *r, dou

int reset=1, reset2, TraceFlag;

double Bmin, I, Phi, cl, sl, rat, SS1, SS2, SS, Sn, Ss, Htry, Hdid, Hnext, Bs, Be, s, sgn;
double Bmin, I1, Phi, cl, sl, rat, SS1, SS2, SS, Sn, Ss, Htry, Hdid, Hnext, Bs, Be, s, sgn;
Lgm_Vector w, u, Pmirror1, Pmirror2, v1, v2, v3, Bvec, P, Ps, u_scale, Bvectmp, Ptmp;
double stmp, Btmp;

Expand Down Expand Up @@ -38,13 +38,13 @@ double ComputeI_FromMltMlat2( double Bm, double MLT, double mlat, double *r, dou

if ( TraceFlag != LGM_CLOSED ) {

I = 9e99;
I1 = 9e99;
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: %13.6g I0: %13.6g > Field Line not closed%s\n", LstarInfo->PreStr, mlat, I, I0, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: %13.6g I0: %13.6g > Field Line not closed%s\n", LstarInfo->PreStr, mlat, I1, I0, LstarInfo->PostStr );
}

*ErrorStatus = -1; // FL open. I undefined.
return( I );
return( I1 );

} else if ( Bmin <= Bm ) {

Expand All @@ -67,40 +67,40 @@ double ComputeI_FromMltMlat2( double Bm, double MLT, double mlat, double *r, dou
//printf("Pmirror2 = %g %g %g\n", Pmirror2.x, Pmirror2.y, Pmirror2.z);

} else {
I = 9e99;
I1 = 9e99;
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: %13.6g I0: %13.6g > Unable to find southern mirror point%s\n", LstarInfo->PreStr, mlat, I, I0, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: %13.6g I0: %13.6g > Unable to find southern mirror point%s\n", LstarInfo->PreStr, mlat, I1, I0, LstarInfo->PostStr );
}
*ErrorStatus = -3; // No valid mirror point in the south
return( I );
return( I1 );
}

// total distance between mirror point.
SS = SS1 + SS2;

// If its really small, just return 0.0 for I
if ( fabs(SS) < 1e-7) {
I = 0.0;
I1 = 0.0;
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: %13.6g I0: %13.6g > Distance between mirror points is < 1e-7Re, Assuming I=0%s\n", LstarInfo->PreStr, mlat, I, I0, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: %13.6g I0: %13.6g > Distance between mirror points is < 1e-7Re, Assuming I1=0%s\n", LstarInfo->PreStr, mlat, I1, I0, LstarInfo->PostStr );
}
*ErrorStatus = 2; // Flag that we did this
return( I );
return( I1 );
}

} else {
I = 9e99;
I1 = 9e99;
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: %13.6g I0: %13.6g > Unable to find northern mirror point%s\n", LstarInfo->PreStr, mlat, I, I0, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: %13.6g I0: %13.6g > Unable to find northern mirror point%s\n", LstarInfo->PreStr, mlat, I1, I0, LstarInfo->PostStr );
}
*ErrorStatus = -2; // No valid mirror point in the north
return( I );
return( I1 );
}

/*
* OK, we have both mirror points. Lets compute I
*/
I = 9e99;
I1 = 9e99;
LstarInfo->mInfo->Hmax = 0.1;
LstarInfo->mInfo->Hmax = SS/(double)LstarInfo->mInfo->nDivs;

Expand Down Expand Up @@ -137,16 +137,16 @@ double ComputeI_FromMltMlat2( double Bm, double MLT, double mlat, double *r, dou
/*
* Do I integral with interped integrand.
*/
I = Iinv_interped( LstarInfo->mInfo );
I1 = Iinv_interped( LstarInfo->mInfo );
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: %13.6g I0: %13.6g I-I0: %13.6g [Sa,Sb]: %.8g %.8g (nCalls = %d)%s\n", LstarInfo->PreStr, mlat, I, I0, I-I0, LstarInfo->mInfo->Sm_South, LstarInfo->mInfo->Sm_North, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: %13.6g I0: %13.6g I1-I0: %13.6g [Sa,Sb]: %.8g %.8g (nCalls = %d)%s\n", LstarInfo->PreStr, mlat, I1, I0, I1-I0, LstarInfo->mInfo->Sm_South, LstarInfo->mInfo->Sm_North, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, LstarInfo->PostStr );
}
FreeSpline( LstarInfo->mInfo );

} else {
I = 9e99;
I1 = 9e99;
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: %13.6g I0: %13.6g Couldnt initialize spline%s\n", LstarInfo->PreStr, mlat, I, I0, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: %13.6g I0: %13.6g Couldnt initialize spline%s\n", LstarInfo->PreStr, mlat, I1, I0, LstarInfo->PostStr );
}
}

Expand All @@ -160,14 +160,14 @@ double ComputeI_FromMltMlat2( double Bm, double MLT, double mlat, double *r, dou

*ErrorStatus = -10; // We didnt even try to compute I because we Field line min-B is greater than Bm.
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s mlat: %13.6g I: undefined I0: %13.6g I undefined. min-B is greater than Bm.%s\n", LstarInfo->PreStr, mlat, I0, LstarInfo->PostStr );
printf("\t\t%s mlat: %13.6g I1: undefined I0: %13.6g I1 undefined. min-B is greater than Bm.%s\n", LstarInfo->PreStr, mlat, I0, LstarInfo->PostStr );
}
return( 9e99 );

}


return( I );
return( I1 );


}
42 changes: 21 additions & 21 deletions libLanlGeoMag/ComputeLstar.c
Original file line number Diff line number Diff line change
Expand Up @@ -571,7 +571,7 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
int done2, Count, FoundShellLine, nIts, Type, retEarthTrace;
int nShabI = 0, nShabII = 0;
double rat, B, dSa, dSb, smax, SS, L, Hmax, epsabs, epsrel;
double I=-999.9, Ifound, M, MLT0, MLT, DeltaMLT, mlat, r, sa, sa2;
double I1=-999.9, Ifound, M, MLT0, MLT, DeltaMLT, mlat, r, sa, sa2;
double Phi, Phi1, Phi2, sl, cl, MirrorMLT[3*LGM_LSTARINFO_MAX_FL], MirrorMlat[3*LGM_LSTARINFO_MAX_FL], pred_mlat, mlat_try, pred_delta_mlat=0.0, mlat0, mlat1, delta;
double MirrorMLT_Old[3*LGM_LSTARINFO_MAX_FL], MirrorMlat_Old[3*LGM_LSTARINFO_MAX_FL], res;
char *PreStr, *PostStr;
Expand Down Expand Up @@ -608,7 +608,7 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
if ((LstarInfo->PitchAngle < 0.0)||(LstarInfo->PitchAngle>90.0)) return(-1);

for (k=0; k<LGM_LSTARINFO_MAX_FL; k++){
LstarInfo->I[k] = LGM_FILL_VALUE;
LstarInfo->I_data[k] = LGM_FILL_VALUE;
}


Expand Down Expand Up @@ -716,10 +716,10 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
// if FL length is small, use an approx expression for I
rat = LstarInfo->mInfo->Bmin/LstarInfo->mInfo->Bm;
if ((1.0-rat) < 0.0) {
I = 0.0;
I1 = 0.0;
} else {
// Eqn 2.66b in Roederer
I = SS*sqrt(1.0 - rat);
I1 = SS*sqrt(1.0 - rat);
}

} else if ( LstarInfo->mInfo->UseInterpRoutines ) {
Expand All @@ -737,9 +737,9 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
/*
* Do interped I integral.
*/
I = Iinv_interped( LstarInfo->mInfo );
I1 = Iinv_interped( LstarInfo->mInfo );
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t %sIntegral Invariant, I (interped): %g%s\n", PreStr, I, PostStr );
printf("\t\t %sIntegral Invariant, I1 (interped): %g%s\n", PreStr, I1, PostStr );
}

/*
Expand All @@ -759,7 +759,7 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
}
FreeSpline( LstarInfo->mInfo );
} else {
I = LGM_FILL_VALUE;
I1 = LGM_FILL_VALUE;
}

} else {
Expand All @@ -768,21 +768,21 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
* Do full blown I integral. (Integrand is evaluated by tracing to required s-values.)
* (This strategy isnt used very much anymore - 20120524, MGH)
*/
I = Iinv( LstarInfo->mInfo );
if (LstarInfo->VerbosityLevel > 1) printf("\t\t %sIntegral Invariant, I (full integral): %g%s\n", PreStr, I, PostStr );
I1 = Iinv( LstarInfo->mInfo );
if (LstarInfo->VerbosityLevel > 1) printf("\t\t %sIntegral Invariant, I1 (full integral): %g%s\n", PreStr, I1, PostStr );

}
LstarInfo->I0 = I; // save initial I in LstarInfo structure.
Ifound = I;
LstarInfo->I0 = I1; // save initial I in LstarInfo structure.
Ifound = I1;

if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t %sLgm_n_I_integrand_Calls: %d%s\n\n", PreStr, LstarInfo->mInfo->Lgm_n_I_integrand_Calls, PostStr );
printf("\t\t%sCurrent Dipole Moment, M_cd: %g%s\n", PreStr, LstarInfo->mInfo->c->M_cd, PostStr);
printf("\t\t%sReference Dipole Moment, M_cd_McIlwain: %g%s\n", PreStr, LstarInfo->mInfo->c->M_cd_McIlwain, PostStr);
printf("\t\t%sReference Dipole Moment, M_cd_2010: %g%s\n", PreStr, LstarInfo->mInfo->c->M_cd_2010, PostStr);
printf("\t\t%sDipole Moment Used, Mused: %g%s\n", PreStr, LstarInfo->Mused, PostStr);
printf("\t\t%sMcIlwain L (Hilton): %.15g%s\n", PreStr, L = LFromIBmM_Hilton( I, LstarInfo->mInfo->Bm, LstarInfo->Mused ), PostStr );
printf("\t\t%sMcIlwain L (McIlwain): %.15g%s\n", PreStr, L = LFromIBmM_McIlwain( I, LstarInfo->mInfo->Bm, LstarInfo->Mused ), PostStr );
printf("\t\t%sMcIlwain L (Hilton): %.15g%s\n", PreStr, L = LFromIBmM_Hilton( I1, LstarInfo->mInfo->Bm, LstarInfo->Mused ), PostStr );
printf("\t\t%sMcIlwain L (McIlwain): %.15g%s\n", PreStr, L = LFromIBmM_McIlwain( I1, LstarInfo->mInfo->Bm, LstarInfo->Mused ), PostStr );
}

} else {
Expand Down Expand Up @@ -1005,7 +1005,7 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
printf("\t\t%s________________________________________________________________________________________________________________________________%s\n", PreStr, PostStr );
}

FoundShellLine = FindShellLine( I, &Ifound, LstarInfo->mInfo->Bm, MLT, &mlat, &r, mlat0, mlat_try, mlat1, &nIts, 1, LstarInfo );
FoundShellLine = FindShellLine( I1, &Ifound, LstarInfo->mInfo->Bm, MLT, &mlat, &r, mlat0, mlat_try, mlat1, &nIts, 1, LstarInfo );
if (FoundShellLine > 0) {
// Found valid FL for drift shell
done2 = TRUE;
Expand Down Expand Up @@ -1040,15 +1040,15 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){

if ( (Type > 1) && (LstarInfo->ShabanskyHandling==LGM_SHABANSKY_HALVE_I) ){
if (LstarInfo->VerbosityLevel > 0) {
printf("\t\t\t%sShabansky orbit. Re-doing FL. Target I adjusted to: %g . (Original is: %g) %s\n", PreStr, I/2.0, I, PostStr );
printf("\t\t\t%sShabansky orbit. Re-doing FL. Target I1 adjusted to: %g . (Original is: %g) %s\n", PreStr, I1/2.0, I1, PostStr );
}

FoundShellLine = FindShellLine( I/2.0, &Ifound, LstarInfo->mInfo->Bm, MLT, &mlat, &r, mlat0, mlat_try, mlat1, &nIts, 1, LstarInfo );
FoundShellLine = FindShellLine( I1/2.0, &Ifound, LstarInfo->mInfo->Bm, MLT, &mlat, &r, mlat0, mlat_try, mlat1, &nIts, 1, LstarInfo );
//TODO: Do we need to test to make sure that the adjusted I is being found on a field line with multiple minima??
PredMinusActualMlat = pred_mlat - mlat;
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s________________________________________________________________________________________________________________________________%s\n\n", PreStr, PostStr );
printf("\t\t%s >> Pred/Actual/Diff mlat: %g/%g/%g MLT/MLAT: %g %g I0: %g I: %g I-I0/2: %g (SHABANSKY)%s\n", PreStr, pred_mlat, mlat, PredMinusActualMlat, MLT, mlat, I, Ifound, Ifound-I/2.0, PostStr );
printf("\t\t%s >> Pred/Actual/Diff mlat: %g/%g/%g MLT/MLAT: %g %g I0: %g I1: %g I1-I0/2: %g (SHABANSKY)%s\n", PreStr, pred_mlat, mlat, PredMinusActualMlat, MLT, mlat, I1, Ifound, Ifound-I1/2.0, PostStr );
printf("\t\t%s________________________________________________________________________________________________________________________________ %s\n\n\n", PreStr, PostStr );
}

Expand All @@ -1060,15 +1060,15 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
PredMinusActualMlat = pred_mlat - mlat;
if (LstarInfo->VerbosityLevel > 1) {
printf("\t\t%s________________________________________________________________________________________________________________________________%s\n\n", PreStr, PostStr );
printf("\t\t%s >> Pred/Actual/Diff mlat: %g/%g/%g MLT/MLAT: %g %g I0: %g I: %g I-I0: %g %s\n", PreStr, pred_mlat, mlat, PredMinusActualMlat, MLT, mlat, I, Ifound, Ifound-I, PostStr );
printf("\t\t%s >> Pred/Actual/Diff mlat: %g/%g/%g MLT/MLAT: %g %g I0: %g I1: %g I1-I0: %g %s\n", PreStr, pred_mlat, mlat, PredMinusActualMlat, MLT, mlat, I1, Ifound, Ifound-I1, PostStr );
printf("\t\t%s________________________________________________________________________________________________________________________________ %s\n\n\n", PreStr, PostStr );
}
}

} else if ( Count > 2 ) {
// Tried to find valid FL more than three times
done2 = TRUE;
if (LstarInfo->VerbosityLevel >0) { printf(" \t%sNo valid I - Drift Shell not closed: L* = undefined (FoundShellLine = %d)%s\n", PreStr, FoundShellLine, PostStr); fflush(stdout); }
if (LstarInfo->VerbosityLevel >0) { printf(" \t%sNo valid I1 - Drift Shell not closed: L* = undefined (FoundShellLine = %d)%s\n", PreStr, FoundShellLine, PostStr); fflush(stdout); }
FoundShellLine = 0;
break;//return(-3);
} else {
Expand Down Expand Up @@ -1147,7 +1147,7 @@ int Lstar( Lgm_Vector *vin, Lgm_LstarInfo *LstarInfo ){
/*
* Save individual I values
*/
LstarInfo->I[k] = Ifound;
LstarInfo->I_data[k] = Ifound;
Comment thread
drsteve marked this conversation as resolved.


MirrorMLT[k] = MLT;
Expand Down Expand Up @@ -2089,7 +2089,7 @@ int Lgm_LCDS( long int Date, double UTC, double brac1, double brac2, double Kin,
Pinner = Ptest;
LCDSv3 = v3;
LCDS = LstarInfo_test->LS;
*K = (LstarInfo_test->I[0])*sqrt(LstarInfo_test->mInfo->Bm*nTtoG);
*K = (LstarInfo_test->I_data[0])*sqrt(LstarInfo_test->mInfo->Bm*nTtoG);

//Determine the type of the orbit
LstarInfo_test->DriftOrbitType = LGM_DRIFT_ORBIT_CLOSED;
Expand Down
Loading