/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

/* 
 * File:   shue_compute.hh
 * Author: hacene
 *
 * Created on December 31, 2020, 2:50 PM
 */

#ifndef SHUE_COMPUTE_HH
#define SHUE_COMPUTE_HH

#include "geopack.hh"

namespace AMDA {
    namespace Parameters {
        namespace Shue {
#define LMAX_VAL 2000

            class ShueCompute {
            public:

                static void initInGSM(int iyr, int iday, int ihour, int min, int isec) {
                    // We use GSM coordinate frames
                    float v_default_x = -400.0;
                    float v_default_y = 0.0;
                    float v_default_z = 0.0;

                    recalc_08_(&iyr, &iday, &ihour, &min, &isec, &v_default_x, &v_default_y, &v_default_z);
                }

                static bool shue98(float Pdyn_i, float bzIMF_i, float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM,
                        float &x_mgnp, float &y_mgnp, float &z_mgnp, float &DIST_MGNP, int &ID_MGNP) {
                    // We use dynamic pressure Pdyn
                    float vel_mgnp = -1.0;

                    shuetal_mgnp_08_(&Pdyn_i, &vel_mgnp, &bzIMF_i, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
                            &x_mgnp, &y_mgnp, &z_mgnp, &DIST_MGNP, &ID_MGNP);

                    return ID_MGNP != 0;
                }
                
              static bool shue97(float Pdyn_i, float bzIMF_i, float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM,
                    float &x_mgnp, float &y_mgnp, float &z_mgnp, float &DIST_MGNP, int &ID_MGNP) {
                  
                  float r0, r, alpha, r, rm;
                  float x_mgnp, y_mgnp, z_mgn;
                  float rho2, st, ct, t;
                  int ID_MGNP;ID96, nit;
                  int iter_limit = 1000;
                  /** DEFINE THE ANGLE PHI, MEASURED DUSKWARD FROM THE NOON-MIDNIGHT MERIDIAN PLANE;
                    IF THE OBSERVATION POINT LIES ON THE X AXIS, THE ANGLE PHI CANNOT BE UNIQUELY
                    C  DEFINED, AND WE SET IT AT ZERO: **/
                  ID_MGNP = -1;
                  float phi=0.0;
                  if(sat_pos_Y_GSM != 0 | sat_pos_Z_GSM !=0 )
                      phi = std::atan2(sat_pos_Y_GSM, sat_pos_Z_GSM);
                  /** FIRST, FIND OUT IF THE OBSERVATION POINT LIES INSIDE THE SHUE ET AL BDRY
                   AND SET THE VALUE OF THE ID FLAG: **/
                  r=std::sqrt(sat_pos_X_GSM*sat_pos_X_GSM, sat_pos_Y_GSM*sat_pos_Y_GSM,  sat_pos_Z_GSM*sat_pos_Z_GSM);
                  if(bzIMF_i >0){
                      r0 = (11.4 + 0.013*bzIMF_i)*std::pow(Pdyn_i, -1/66);
                  }else{
                      r0  = (11.4 + 0.14*bzIMF_i)*std::pow(Pdyn_i, -1/66);
                  }
                  alpha= (0.58 -0.010*bzIMF_i)*(1 + 0.010*Pdyn_i);
                  rm = r0*std::pow((2./(1.0 + sat_pos_X_GSM/r)),alpha);
                  if(r<rm)
                      id=1;
                  /**
                    NOW, FIND THE CORRESPONDING T96 MAGNETOPAUSE POSITION, TO BE USED AS
                    A STARTING APPROXIMATION IN THE SEARCH OF A CORRESPONDING SHUE ET AL.
                    BOUNDARY POINT: 
                   */
                  //  t96_mgnp_08_(&Pdyn_i, &vel_mgnp, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
                 //               &x_mgnp, &y_mgnp, &z_mgnp, &DIST_MGNP, &ID_MGNP);
                  t96_mgnp_08_(&Pdyn_i, -1.0, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
                                &x_mgnp, &y_mgnp, &z_mgnp, &DIST_MGNP, &ID96);
                  
                  rho2 = y_mgnp*y_mgnp + z_mgnp*z_mgnp;
                  r = std::sqrt(rho2 + x_mgnp*x_mgnp);
                  st = std::sqrt(rho2)/r;
                  ct = x_mgnp/r;
                  
                  /*NOW, USE NEWTON'S ITERATIVE METHOD TO FIND THE NEAREST POINT AT THE
                      SHUE ET AL.'S BOUNDARY:
                   **/
                  nit = 0;
                  while(nit < iter_limit){
                  t = std::atan2(st,ct);
                  rm = r0*std::pow((2.0/(1.0 + ct)),alpha);
                  
                  
                  
                  
                  }
                  
                  

                    return ID_MGNP != 0;
                }
                
            };
            
            
        }

    }
}


#endif /* SHUE_COMPUTE_HH */