/*
 * 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) {

                    /**
                   * 
                   * @param Pdyn_i SOLAR WIND RAM PRESSURE IN NANOPASCALS
                   * @param bzIMF_i IMF BZ IN NANOTESLAS
                   * @param sat_pos_X_GSM
                   * @param sat_pos_Y_GSM
                   * @param sat_pos_Z_GSM
                   * @param x_mgnp x POSITION OF THE BOUNDARY POINT
                   * @param y_mgnp y POSITION OF THE BOUNDARY POINT
                   * @param z_mgnp z POSITION OF THE BOUNDARY POINT
                   * @param DIST_MGNP DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSW,YGSW,ZGSW) AND THE MODEL NAGNETOPAUSE
                   * @param ID_MGNP POSITION FLAG:  ID=+1 (-1) MEANS THAT THE OBSERVATION POINT LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
                   * @return true of converged or false 
                   */

                    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) {

                  /**
                   * 
                   * @param Pdyn_i SOLAR WIND RAM PRESSURE IN NANOPASCALS
                   * @param bzIMF_i IMF BZ IN NANOTESLAS
                   * @param sat_pos_X_GSM
                   * @param sat_pos_Y_GSM
                   * @param sat_pos_Z_GSM
                   * @param x_mgnp x POSITION OF THE BOUNDARY POINT
                   * @param y_mgnp y POSITION OF THE BOUNDARY POINT
                   * @param z_mgnp z POSITION OF THE BOUNDARY POINT
                   * @param DIST_MGNP DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSW,YGSW,ZGSW) AND THE MODEL NAGNETOPAUSE
                   * @param ID_MGNP POSITION FLAG:  ID=+1 (-1) MEANS THAT THE OBSERVATION POINT LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
                   * @return true of converged or false 
                   */

                  float r0, r, rm, alpha;
                  float x_mt96, y_mt96, z_mt96;
                  float rho2, rho, st, ct, t, ds, dr, dt, f, gradf, gradf_t, gradf_r;
                  int ID96, nit;

                  int iter_limit = 1000;

                  float eps = 1E-4;
                  float vel = -1.0;

                  /** 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/6.6);

                  }else{

                      r0  = (11.4 + 0.14*bzIMF_i)*std::pow(Pdyn_i, -1/6.6);

                  }
                  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_MGNP =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, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
                                &x_mt96, &y_mt96, &z_mt96, &DIST_MGNP, &ID96);

                  rho2 = y_mt96*y_mt96 + z_mt96*z_mt96;
                  r = std::sqrt(rho2 + x_mt96*x_mt96);

                  st = std::sqrt(rho2)/r;

                  ct = x_mt96/r;

                  
                  /*NOW, USE NEWTON'S ITERATIVE METHOD TO FIND THE NEAREST POINT AT THE
                      SHUE ET AL.'S BOUNDARY:
                   **/
                  nit = 0;

                  ds = 1000.0;
                  bool converged = true;
                  while(ds > eps){
                      if(nit >  iter_limit){
                          converged = false;
                          break; 
                      }

                  t = std::atan2(st,ct);
                  rm = r0*std::pow((2.0/(1.0 + ct)),alpha);
                  

                  f = r -rm;                  
                  gradf_r = 1.0;
                  gradf_t = -alpha/r*rm*st/(1.0+ct);
                  gradf=std::sqrt(gradf_r*gradf_r + gradf_t*gradf_t);

                  dr = -f/(gradf*gradf);
                  dt = dr/r*gradf_t;

                  r = r + dr;                  
                  t = t +dt;
                  st = sin(t);
                  ct = cos(t);

                  ds= std::sqrt(dr*dr + (r*dt)*(r*dt));  
                  nit += 1;
                  }
                  if (converged){
                      x_mgnp = r*std::cos(t);
                      rho = r*std::sin(t);
                      y_mgnp = rho*sin(phi);
                      z_mgnp = rho*cos(phi);
                      DIST_MGNP = std::sqrt((sat_pos_X_GSM -x_mgnp)*(sat_pos_X_GSM -x_mgnp) + 
                              (sat_pos_Y_GSM -y_mgnp)*(sat_pos_Y_GSM -y_mgnp)  +
                              (sat_pos_Z_GSM -z_mgnp)*(sat_pos_Z_GSM -z_mgnp) );

                  }
                  

                    return converged;

                }
                

            };

        }

    }
}


#endif /* SHUE_COMPUTE_HH */