shue_compute.hh
4.48 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
/*
* 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 */