GeopackWrapper.hh
4.69 KB
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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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
#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 */