shue_compute.hh
7.25 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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
/*
* 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 */