/************************************************************************************************
    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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// 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;    
}