/************************************************************************************************ Name : makeFootprints Author : Arnaud BIEGUN Version : 2.0 (10/07/2015) Brief : Computes geomagnetic field lines feet on a specified spherical surface *************************************************************************************************/ /**************** HEADERS ******************/ // C LIB #include #include #include #include // CONSTANTS // #define PLASMA_FILENAME "inputs.txt" // FORTRAN BINDING extern void recalc_08_(); extern void trace_08_(); extern void t96_01_(); extern void geogsw_08_(); extern void gswgse_08_(); /**********END OF HEADERS ******************/ int main(int argc,char *argv[]) { /****************** LOCAL VARIABLES ******************/ // TIME int iyr, iday, ihour, min__, isec; // SATELLITE float sat_pos_GSE_X_i, sat_pos_GSE_Y_i, sat_pos_GSE_Z_i; float sat_pos_GSW_X_i, sat_pos_GSW_Y_i, sat_pos_GSW_Z_i; // PLASMA float Pdyn_i, sym_h_i; float V_GSE_X_i, V_GSE_Y_i, V_GSE_Z_i; float B_GSE_X_i, B_GSE_Y_i, B_GSE_Z_i; float B_GSW_X_i, B_GSW_Y_i, B_GSW_Z_i; // TSYGANENKO MODEL float DIR, R0, RLIM, DSMAX; float PSI; float BT96_GSW_X, BT96_GSW_Y, BT96_GSW_Z; float BT96_GSE_X, BT96_GSE_Y, BT96_GSE_Z; float x_mgnp, y_mgnp, z_mgnp, vel_mgnp; float DIST_MGNP; float XX[2000], YY[2000], ZZ[2000], PARMOD[10]; float FP_GSW_X_i, FP_GSW_Y_i, FP_GSW_Z_i; float FP_GEO_X_i, FP_GEO_Y_i, FP_GEO_Z_i; float ERR; int L, LMAX, IOPT, IMAP, ISWq = 0; int ID_MGNP; // +1 : inside; -1 : outside // RADIANS TO DEGREES float SQ, R, PHI, THETA; // OTHERS FILE *plasma_file = NULL; // FILE *fp_file = NULL; char line[1000]; char errMsg[256]; int SWFLAG_i; int transform_flag; int i; /****************** END OF LOCAL VARIABLES ******************/ /*********************CHECK I/O FILES **********************/ char* plasmaFilename = argv[1]; plasma_file = fopen(plasmaFilename, "r"); if(plasma_file == NULL) { exit(EXIT_FAILURE); } /****************** END OF CHECK I/O FILES ********************/ /****************** FOOTPRINTS COMPUTING ******************/ // T96 MODEL CONFIGURATION R0 = 1.02; // ~ 120 km // Inner boundary radii RLIM = 50.0; // Outter boundary radii DSMAX = 0.5; // Upper limit of the stepsize LMAX = 2000; // Max number of field line points ERR = 0.0001; // Permissible step error IOPT = 0; // Model index PSI = -0.47; // MAIN LOOP while ( fgets(line, 1000, plasma_file) != NULL ) { // GET VALUES FROM INPUT FILE sscanf(line, // "%4d %3d %2d %2d %f %f %f %f %f %f %f %f %f %f %f %d", "%d %d %d %d %f %f %f %f %f %f %f %f %f %f %f %d", &iyr, &iday, &ihour, &min__, &sat_pos_GSE_X_i, &sat_pos_GSE_Y_i, &sat_pos_GSE_Z_i, &B_GSE_X_i, &B_GSE_Y_i, &B_GSE_Z_i, &V_GSE_X_i, &V_GSE_Y_i, &V_GSE_Z_i, &Pdyn_i, &sym_h_i, &SWFLAG_i ); // CONFIGURATION OF GEOMAGNETIC FIELD MODELS recalc_08_(&iyr, &iday, &ihour, &min__, &isec, &V_GSE_X_i, &V_GSE_Y_i, &V_GSE_Z_i); printf("%d %d %d %d ", iyr, iday, ihour, min__); // fprintf(fp_file, "%d %d %d %d ", iyr, iday, ihour, min__); // GSE --> GSW ( GSM ) transform_flag = -1; gswgse_08_(&sat_pos_GSW_X_i, &sat_pos_GSW_Y_i, &sat_pos_GSW_Z_i, &sat_pos_GSE_X_i, &sat_pos_GSE_Y_i, &sat_pos_GSE_Z_i, &transform_flag); gswgse_08_(&B_GSW_X_i, &B_GSW_Y_i, &B_GSW_Z_i, &B_GSE_X_i, &B_GSE_Y_i, &B_GSE_Z_i, &transform_flag); // COMPUTES T96 GEOMAGNETIC FIELD MODEL PARMOD[0] = Pdyn_i; PARMOD[1] = sym_h_i; PARMOD[2] = B_GSW_Y_i; PARMOD[3] = B_GSW_Z_i; t96_01_(&IOPT, PARMOD, &PSI, &sat_pos_GSW_X_i, &sat_pos_GSW_Y_i, &sat_pos_GSW_Z_i, &BT96_GSW_X, &BT96_GSW_Y, &BT96_GSW_Z); transform_flag = 1; gswgse_08_(&BT96_GSW_X, &BT96_GSW_Y, &BT96_GSW_Z, &BT96_GSE_X, &BT96_GSE_Y, &BT96_GSE_Z, &transform_flag); // IS IN MAGNETOPAUSE ? vel_mgnp = -1.0; // We use dynamic pressure Pdyn // shuetal_mgnp_08__(&Pdyn_i, &vel_mgnp, // &B_GSW_Z_i, // &sat_pos_GSW_X_i, // &sat_pos_GSW_Y_i, // &sat_pos_GSW_Z_i, // &x_mgnp, &y_mgnp, &z_mgnp, // &DIST_MGNP, &ID_MGNP); t96_mgnp_08_(&Pdyn_i, &vel_mgnp, &sat_pos_GSW_X_i, &sat_pos_GSW_Y_i, &sat_pos_GSW_Z_i, &x_mgnp, &y_mgnp, &z_mgnp, &DIST_MGNP, &ID_MGNP); // FOOTPRINTS -> STEP 1 : NORTH DIR = -1.0; // Parallel to the total field vector trace_08_(&sat_pos_GSW_X_i, &sat_pos_GSW_Y_i, &sat_pos_GSW_Z_i, &DIR, &DSMAX, &ERR, &RLIM, &R0, &IOPT, PARMOD, &FP_GSW_X_i, &FP_GSW_Y_i, &FP_GSW_Z_i, XX, YY, ZZ, &L, &LMAX); transform_flag = -1; geogsw_08_(&FP_GEO_X_i, &FP_GEO_Y_i, &FP_GEO_Z_i, &FP_GSW_X_i, &FP_GSW_Y_i, &FP_GSW_Z_i, &transform_flag); SQ=FP_GEO_X_i*FP_GEO_X_i + FP_GEO_Y_i*FP_GEO_Y_i; R=sqrt(SQ+FP_GEO_Z_i*FP_GEO_Z_i); SQ=sqrt(SQ); PHI=atan2(FP_GEO_Y_i, FP_GEO_X_i); THETA=atan2(FP_GEO_Z_i, SQ); if (PHI < 0.) { PHI += 2.*M_PI; } if ( (ID_MGNP > 0) && (R <= R0) ) { printf("%.2f %.2f %.2f ", R, THETA*(180./M_PI), PHI*(180./M_PI)); // fprintf(fp_file, "%.2f %.2f %.2f ", R, THETA*(180./M_PI), PHI*(180./M_PI)); } else { printf("%.2f %.2f %.2f ", nanf(" "), nanf(" "), nanf(" ")); } // FOOTPRINTS -> STEP 2 : SOUTH DIR = 1.0; // Antiparallel to the total field vector trace_08_(&sat_pos_GSW_X_i, &sat_pos_GSW_Y_i, &sat_pos_GSW_Z_i, &DIR, &DSMAX, &ERR, &RLIM, &R0, &IOPT, PARMOD, &FP_GSW_X_i, &FP_GSW_Y_i, &FP_GSW_Z_i, XX, YY, ZZ, &L, &LMAX); transform_flag = -1; geogsw_08_(&FP_GEO_X_i, &FP_GEO_Y_i, &FP_GEO_Z_i, &FP_GSW_X_i, &FP_GSW_Y_i, &FP_GSW_Z_i, &transform_flag); SQ=FP_GEO_X_i*FP_GEO_X_i + FP_GEO_Y_i*FP_GEO_Y_i; R=sqrt(SQ+FP_GEO_Z_i*FP_GEO_Z_i); SQ=sqrt(SQ); PHI=atan2(FP_GEO_Y_i, FP_GEO_X_i); THETA=atan2(FP_GEO_Z_i, SQ); if (PHI < 0.) { PHI += 2.*M_PI; } if ( (ID_MGNP > 0) && (R <= R0) ) { printf("%.2f %.2f %.2f ", R, THETA*(180./M_PI), PHI*(180./M_PI)); // fprintf(fp_file, "%.2f %.2f %.2f ", R, THETA*(180./M_PI), PHI*(180./M_PI)); } else { printf("%.2f %.2f %.2f ", nanf(" "), nanf(" "), nanf(" ")); } // Print SW flag index printf("%d ", SWFLAG_i); // fprintf(fp_file, "%d ", SWFLAG_i); // Print distance to magnetopause printf("%.2f ", DIST_MGNP*ID_MGNP); // fprintf(fp_file, "%f ", DIST_MGNP*ID_MGNP); // Print positions (GSE) and T96 output (GSE) printf("%.2f %.2f %.2f ", sat_pos_GSE_X_i, sat_pos_GSE_Y_i, sat_pos_GSE_Z_i); printf("%.2f %.2f %.2f\n", BT96_GSE_X, BT96_GSE_Y, BT96_GSE_Z); // fprintf(fp_file, "%.2f %.2f %.2f ", sat_pos_GSE_X_i, sat_pos_GSE_Y_i, sat_pos_GSE_Z_i); // fprintf(fp_file, "%.2f %.2f %.2f\n", BT96_GSE_X, BT96_GSE_Y, BT96_GSE_Z); } /****************** END OF FOOTPRINTS COMPUTING ******************/ /************************** FREE I/O FILES ***********************/ // fclose(fp_file); fclose(plasma_file); /********************* END OF FREE I/O FILES *********************/ return 0; }