#include "geopack.hh"


namespace AMDA {
namespace Parameters {
namespace Tsyganenko96 {
namespace geopack {
   
#define LMAX_VAL 2000
 
class GeopackWrapper {
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 isInMagnetopause(float Pdyn_i, float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM) {
        // We use dynamic pressure Pdyn
        float vel_mgnp = -1.0;
        
        float x_mgnp = 0.;
        float y_mgnp = 0.;
        float z_mgnp = 0.;
        float DIST_MGNP = 0.;
        int ID_MGNP = 0;
        
        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);
            
        return (ID_MGNP == 1);
    }

    static bool computeGeomagneticFieldInGSM(float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM,
        float Pdyn, float Dst, float B_Y_GSM, float B_Z_GSM,
	float &B_X_GSM_RES, float& B_Y_GSM_RES, float& B_Z_GSM_RES) {

	// COMPUTES T96 GEOMAGNETIC FIELD MODEL
	int IOPT = 0;
	float PARMOD[10];
	memset(PARMOD, 0, 10 * sizeof(float));
	PARMOD[0] = Pdyn;
	PARMOD[1] = Dst;
	PARMOD[2] = B_Y_GSM;
	PARMOD[3] = B_Z_GSM;

	float BX = 0;
	float BY = 0;
	float BZ = 0;
	t96_01_modified_(&IOPT, PARMOD, &sat_pos_X_GSM, &sat_pos_Y_GSM,
            &sat_pos_Z_GSM, &BX, &BY, &BZ);

	// COMPUTES GEODIPOLE MAGNETIC field
	float HX = 0.;
	float HY = 0.;
	float HZ = 0.;
	igrf_gsw_08_(&sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM, &HX, &HY, &HZ);

	B_X_GSM_RES = BX + HX;
	B_Y_GSM_RES = BY + HY;
	B_Z_GSM_RES = BZ + HZ;

	return true;
    }
    
    static bool computeGeomagneticFieldInGSE(float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM,
        float Pdyn, float Dst, float B_Y_GSM, float B_Z_GSM,
        float &B_X_GSE_RES, float& B_Y_GSE_RES, float& B_Z_GSE_RES) {
	float B_X_GSM_RES, B_Y_GSM_RES, B_Z_GSM_RES;
	// COMPUTE GEOMAGNETIC FIELD IN GSM
	if (!GeopackWrapper::computeGeomagneticFieldInGSM(sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM, Pdyn, Dst, B_Y_GSM, B_Z_GSM, B_X_GSM_RES, B_Y_GSM_RES, B_Z_GSM_RES)) {
		return false;
	}
	// APPLY TRANSFORMATION FROM GSM TO GSE
	int transform_flag = 1;
        gswgse_08_(&B_X_GSM_RES, &B_Y_GSM_RES, &B_Z_GSM_RES , &B_X_GSE_RES, &B_Y_GSE_RES, &B_Z_GSE_RES, &transform_flag);
	return true;
    }


    static bool computeFootprint(bool south, float altitude, float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM,
        float Pdyn, float Dst, float B_Y_GSM, float B_Z_GSM, float& FP_X_GSM, float& FP_Y_GSM, float& FP_Z_GSM) {
            
        float DIR = south ? 1.0 : -1.0;
        float DSMAX = 0.5;     // Upper limit of the stepsize
        float ERR   = 0.0001;  // Permissible step error
        float RLIM  = 50.0;    // Outter boundary radii
        float R0    = 1. + altitude;  // Inner boundary radii
        int   IOPT  = 0;       // Model index
        
        float PARMOD[10];
        memset(PARMOD, 0, 10 * sizeof(float));
        PARMOD[0] = Pdyn;
        PARMOD[1] = Dst;
        PARMOD[2] = B_Y_GSM;
        PARMOD[3] = B_Z_GSM;
        
        int LMAX = LMAX_VAL;
        
        float XX[LMAX_VAL];
        float YY[LMAX_VAL];
        float ZZ[LMAX_VAL];
        
        memset(XX, 0, LMAX * sizeof(float));
        memset(YY, 0, LMAX * sizeof(float));
        memset(ZZ, 0, LMAX * sizeof(float));
        
        int L = 0;
        
        trace_08_modified_(&sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM, &DIR,
            &DSMAX, &ERR, &RLIM, &R0, &IOPT, PARMOD, 
            &FP_X_GSM, &FP_Y_GSM, &FP_Z_GSM, XX, YY, ZZ, &L, &LMAX);
            
        return true;
    }
    
    static void toRLatLong(float FP_X_GSM, float FP_Y_GSM, float FP_Z_GSM, float& r, float& lat, float& lon) {
        int transform_flag = -1;
        
        float FP_GEO_X_i = 0.;
        float FP_GEO_Y_i = 0.;
        float FP_GEO_Z_i = 0.;
    
        geogsw_08_(&FP_GEO_X_i, &FP_GEO_Y_i, &FP_GEO_Z_i, 
                  &FP_X_GSM, &FP_Y_GSM, &FP_Z_GSM, 
                  &transform_flag);
                  
        float SQ = 0.;
        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);
        lon=atan2(FP_GEO_Y_i, FP_GEO_X_i);
        lat=atan2(FP_GEO_Z_i, SQ);
      
        if (lon < 0.) 
        {
            lon += 2.*M_PI;
        }
        
        lon *= (180./M_PI);
        lat *= (180./M_PI);
    }
};

} /* geopack */
} /* Tsyganenko96 */
} /* Parameters */
} /* AMDA */