CDPP / AMDA_Kernel

Commit 2cae24f09bfaa7d525399f7a587ca75e59975fe1

Authored by Hacene SI HADJ MOHAND
1 parent ebfd8474
Exists in master and in 75 other branches 10038, 10098, 5390, 5770, 6728, 7616, 7811, 7875, 8459, 8528, 8928, 9016, 9016Bis, 9017, 9017-new, 9049, 9084, 9259, 9326, 9347, 9383, 9719, 9719.new, AllLegends, CDF_ISTP, CDF_ISTP_HOPEFULLY_ALL_WORKING, FER_10400, FER_10412, FER_10663, FER_10765, FER_11000, FER_11100, FER_11242, FER_11277, FER_11500, FER_11537, FER_11588, FER_11651, FER_19836, FER_5344, FER_6031, FER_6598, FER_6966, FER_7318, FER_7618, FER_7857, FER_8455, FER_9629, FER_CR41, FER_CR42, FER_CR43, FER_TFCat, Histo1D2D, JUNO_JEDI, MAZ_10936, MAZ_6018, MAZ_9016_New, MAZ_FIX_9016, SpeasyGet, amdadev, amdadev_test, cu2022, das2, density_maps, develop, evol_9042, feature/9259, pickle, rm_6912, rm_7926, rm_8050, rm_8932, shue97, test-merge, ubuntu

shue 97 ok

Showing 3 changed files with 87 additions and 26 deletions   Show diff stats
Inline Side-by-side
src/ExternLib/Shue/ProcessShue.cc
  Diff comments   View file @2cae24f
... ... @@ -85,8 +85,8 @@ TimeStamp ProcessShue::init() {
85 85 _model98 = true;
86 86 if (_attributList.size() > 2 && _attributList[2] == "shue97")
87 87 _model98 = false;
88   - _inGSE = true;
89 88  
  89 + _inGSE = true;
90 90 if (_attributList.size() > 3 && _attributList[3] == "GSM")
91 91 _inGSE = false;
92 92  
... ...
src/ExternLib/Shue/Shue.hh
  Diff comments   View file @2cae24f
... ... @@ -209,9 +209,19 @@ namespace AMDA {
209 209 float dst;
210 210 int pos_id;
211 211 bool computed = false;
  212 + /**
  213 + * for testing to reproduce JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 103, NO. A8, PAGES 17,691-17,700, AUGUST 1, 1998
  214 + * fig 4
  215 + p_sw=14.0;
  216 + b_z_gsm = 9.3;
  217 + * */
212 218 if(_model98){
213 219 computed = ShueCompute::shue98(p_sw, b_z_gsm, sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM,
214 220 X_GSW, Y_GSW, Z_GSW, dst, pos_id);
  221 + }else{
  222 + computed = ShueCompute::shue97(p_sw, b_z_gsm, sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM,
  223 + X_GSW, Y_GSW, Z_GSW, dst, pos_id);
  224 + }
215 225 if (computed) {
216 226 if(_inGSE){
217 227 // transform results from GSM to GSE
... ... @@ -226,9 +236,6 @@ namespace AMDA {
226 236 ouputElt[2] = Z_GSW;
227 237 }
228 238 }
229   - }else{
230   - // TBD
231   - }
232 239 _paramOutput->pushTime(crtTime);
233 240 pushData(ouputElt, dst, pos_id);
234 241 }
... ...
src/ExternLib/Shue/shue_compute.hh
  Diff comments   View file @2cae24f
... ... @@ -35,7 +35,21 @@ namespace AMDA {
35 35  
36 36 static bool shue98(float Pdyn_i, float bzIMF_i, float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM,
37 37 float &x_mgnp, float &y_mgnp, float &z_mgnp, float &DIST_MGNP, int &ID_MGNP) {
38   - // We use dynamic pressure Pdyn
  38 + /**
  39 + *
  40 + * @param Pdyn_i SOLAR WIND RAM PRESSURE IN NANOPASCALS
  41 + * @param bzIMF_i IMF BZ IN NANOTESLAS
  42 + * @param sat_pos_X_GSM
  43 + * @param sat_pos_Y_GSM
  44 + * @param sat_pos_Z_GSM
  45 + * @param x_mgnp x POSITION OF THE BOUNDARY POINT
  46 + * @param y_mgnp y POSITION OF THE BOUNDARY POINT
  47 + * @param z_mgnp z POSITION OF THE BOUNDARY POINT
  48 + * @param DIST_MGNP DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSW,YGSW,ZGSW) AND THE MODEL NAGNETOPAUSE
  49 + * @param ID_MGNP POSITION FLAG: ID=+1 (-1) MEANS THAT THE OBSERVATION POINT LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
  50 + * @return true of converged or false
  51 + */
  52 +
39 53 float vel_mgnp = -1.0;
40 54  
41 55 shuetal_mgnp_08_(&Pdyn_i, &vel_mgnp, &bzIMF_i, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
... ... @@ -46,62 +60,102 @@ namespace AMDA {
46 60  
47 61 static bool shue97(float Pdyn_i, float bzIMF_i, float sat_pos_X_GSM, float sat_pos_Y_GSM, float sat_pos_Z_GSM,
48 62 float &x_mgnp, float &y_mgnp, float &z_mgnp, float &DIST_MGNP, int &ID_MGNP) {
  63 + /**
  64 + *
  65 + * @param Pdyn_i SOLAR WIND RAM PRESSURE IN NANOPASCALS
  66 + * @param bzIMF_i IMF BZ IN NANOTESLAS
  67 + * @param sat_pos_X_GSM
  68 + * @param sat_pos_Y_GSM
  69 + * @param sat_pos_Z_GSM
  70 + * @param x_mgnp x POSITION OF THE BOUNDARY POINT
  71 + * @param y_mgnp y POSITION OF THE BOUNDARY POINT
  72 + * @param z_mgnp z POSITION OF THE BOUNDARY POINT
  73 + * @param DIST_MGNP DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSW,YGSW,ZGSW) AND THE MODEL NAGNETOPAUSE
  74 + * @param ID_MGNP POSITION FLAG: ID=+1 (-1) MEANS THAT THE OBSERVATION POINT LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
  75 + * @return true of converged or false
  76 + */
49 77  
50   - float r0, r, alpha, r, rm;
51   - float x_mgnp, y_mgnp, z_mgn;
52   - float rho2, st, ct, t;
53   - int ID_MGNP;ID96, nit;
  78 + float r0, r, rm, alpha;
  79 + float x_mt96, y_mt96, z_mt96;
  80 + float rho2, rho, st, ct, t, ds, dr, dt, f, gradf, gradf_t, gradf_r;
  81 + int ID96, nit;
54 82 int iter_limit = 1000;
  83 + float eps = 1E-4;
  84 + float vel = -1.0;
55 85 /** DEFINE THE ANGLE PHI, MEASURED DUSKWARD FROM THE NOON-MIDNIGHT MERIDIAN PLANE;
56 86 IF THE OBSERVATION POINT LIES ON THE X AXIS, THE ANGLE PHI CANNOT BE UNIQUELY
57 87 C DEFINED, AND WE SET IT AT ZERO: **/
58 88 ID_MGNP = -1;
59 89 float phi=0.0;
60   - if(sat_pos_Y_GSM != 0 | sat_pos_Z_GSM !=0 )
  90 + if(sat_pos_Y_GSM != 0 || sat_pos_Z_GSM !=0 ){
61 91 phi = std::atan2(sat_pos_Y_GSM, sat_pos_Z_GSM);
  92 + }
62 93 /** FIRST, FIND OUT IF THE OBSERVATION POINT LIES INSIDE THE SHUE ET AL BDRY
63 94 AND SET THE VALUE OF THE ID FLAG: **/
64   - 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);
65   - if(bzIMF_i >0){
66   - r0 = (11.4 + 0.013*bzIMF_i)*std::pow(Pdyn_i, -1/66);
  95 + 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);
  96 + if(bzIMF_i >=0){
  97 + r0 = (11.4 + 0.013*bzIMF_i)*std::pow(Pdyn_i, -1/6.6);
67 98 }else{
68   - r0 = (11.4 + 0.14*bzIMF_i)*std::pow(Pdyn_i, -1/66);
  99 + r0 = (11.4 + 0.14*bzIMF_i)*std::pow(Pdyn_i, -1/6.6);
69 100 }
70 101 alpha= (0.58 -0.010*bzIMF_i)*(1 + 0.010*Pdyn_i);
71 102 rm = r0*std::pow((2./(1.0 + sat_pos_X_GSM/r)),alpha);
72 103 if(r<rm)
73   - id=1;
  104 + ID_MGNP =1;
74 105 /**
75 106 NOW, FIND THE CORRESPONDING T96 MAGNETOPAUSE POSITION, TO BE USED AS
76 107 A STARTING APPROXIMATION IN THE SEARCH OF A CORRESPONDING SHUE ET AL.
77 108 BOUNDARY POINT:
78 109 */
79   - // t96_mgnp_08_(&Pdyn_i, &vel_mgnp, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
80   - // &x_mgnp, &y_mgnp, &z_mgnp, &DIST_MGNP, &ID_MGNP);
81   - t96_mgnp_08_(&Pdyn_i, -1.0, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
82   - &x_mgnp, &y_mgnp, &z_mgnp, &DIST_MGNP, &ID96);
  110 + t96_mgnp_08_(&Pdyn_i, &vel, &sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
  111 + &x_mt96, &y_mt96, &z_mt96, &DIST_MGNP, &ID96);
83 112  
84   - rho2 = y_mgnp*y_mgnp + z_mgnp*z_mgnp;
85   - r = std::sqrt(rho2 + x_mgnp*x_mgnp);
  113 + rho2 = y_mt96*y_mt96 + z_mt96*z_mt96;
  114 + r = std::sqrt(rho2 + x_mt96*x_mt96);
86 115 st = std::sqrt(rho2)/r;
87   - ct = x_mgnp/r;
  116 + ct = x_mt96/r;
88 117  
89 118 /*NOW, USE NEWTON'S ITERATIVE METHOD TO FIND THE NEAREST POINT AT THE
90 119 SHUE ET AL.'S BOUNDARY:
91 120 **/
92 121 nit = 0;
93   - while(nit < iter_limit){
  122 + ds = 1000.0;
  123 + bool converged = true;
  124 + while(ds > eps){
  125 + if(nit > iter_limit){
  126 + converged = false;
  127 + break;
  128 + }
94 129 t = std::atan2(st,ct);
95 130 rm = r0*std::pow((2.0/(1.0 + ct)),alpha);
96 131  
  132 + f = r -rm;
  133 + gradf_r = 1.0;
  134 + gradf_t = -alpha/r*rm*st/(1.0+ct);
  135 + gradf=std::sqrt(gradf_r*gradf_r + gradf_t*gradf_t);
97 136  
  137 + dr = -f/(gradf*gradf);
  138 + dt = dr/r*gradf_t;
98 139  
  140 + r = r + dr;
  141 + t = t +dt;
  142 + st = sin(t);
  143 + ct = cos(t);
99 144  
  145 + ds= std::sqrt(dr*dr + (r*dt)*(r*dt));
  146 + nit += 1;
  147 + }
  148 + if (converged){
  149 + x_mgnp = r*std::cos(t);
  150 + rho = r*std::sin(t);
  151 + y_mgnp = rho*sin(phi);
  152 + z_mgnp = rho*cos(phi);
  153 + DIST_MGNP = std::sqrt((sat_pos_X_GSM -x_mgnp)*(sat_pos_X_GSM -x_mgnp) +
  154 + (sat_pos_Y_GSM -y_mgnp)*(sat_pos_Y_GSM -y_mgnp) +
  155 + (sat_pos_Z_GSM -z_mgnp)*(sat_pos_Z_GSM -z_mgnp) );
100 156 }
101 157  
102   -
103   -
104   - return ID_MGNP != 0;
  158 + return converged;
105 159 }
106 160  
107 161 };
... ...