CDPP / AMDA_Kernel

Blame view

src/ExternLib/Tsyganenko96/GeopackWrapper.hh 4.69 KB
Edit Raw Normal View History
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
 
#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;
    }
d5d85a78   Benjamin Renard   Add process to co... 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
 
    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;
    }
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
 
    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 */