9719 ok

Hacene SI HADJ MOHAND
1 parent 88605249
Showing 2 changed files with 122 additions and 28 deletions Show diff stats
src/Parameters/DataTypeMath.hh
test/data/functions.xml
@@ -317,7 +317,7 @@ Type1 alfvenVelocity(Type1 density, Type2 mag) {
     if (density <= 0) {
         val << NotANumber();
     } else {
-        val = 21.83 *mag / std::sqrt(density);
+        val = 21.83 * mag / std::sqrt(density);
     }
     return val;
 }
@@ -342,8 +342,8 @@ Type1 fairfield70(Type1  y_sm, Type2 z_sm,Type3 x_ss) {
  * @return delta Z 
  */
 template <typename Type1, typename Type2>
-Type1 fairfield70(std::vector<Type1>  xyz_sm, Type2 x_ss) {
-    Type1 val = xyz_sm[2]-std::sqrt(121.0 - 0.537778 * pow(xyz_sm[1], 2.0)) * std::sin(x_ss*DEG2RAD);
+Type1 fairfield70(std::vector<Type1> xyz_sm, Type2 x_ss) {
+    Type1 val = xyz_sm[2] - std::sqrt(121.0 - 0.537778 * pow(xyz_sm[1], 2.0)) * std::sin(x_ss * DEG2RAD);
     return val;
 }
  
@@ -354,8 +354,8 @@ Type1 fairfield70(std::vector&lt;Type1&gt;  xyz_sm, Type2 x_ss) {
  * @return delta Z 
  */
 template <typename Type1, typename Type2>
-Type1 fairfield80(std::vector<Type1>  xyz_sm, Type2 x_ss) {
-    Type1 val = xyz_sm[2]-10.5*std::sqrt(1.0-pow(xyz_sm[1], 2.0)/289.0)* std::sin(x_ss * DEG2RAD);
+Type1 fairfield80(std::vector<Type1> xyz_sm, Type2 x_ss) {
+    Type1 val = xyz_sm[2] - 10.5 * std::sqrt(1.0 - pow(xyz_sm[1], 2.0) / 289.0) * std::sin(x_ss * DEG2RAD);
     return val;
 }
  
@@ -370,8 +370,8 @@ Type1 fairfield80(std::vector&lt;Type1&gt;  xyz_sm, Type2 x_ss) {
  *            denotes the latitudinal angle of the IMF vector measured from the solar equatorial plane.
  */
 template <typename Type1, typename Type2, typename Type3>
-Type1 perreault78(Type1  E_mag, Type2 B_mag,Type3 theta_p) {
-    Type1 val = E_mag*B_mag/(4*std::atan(1)*4)*std::pow((7*RE*std::pow(sin(theta_p/2),2)),2);
+Type1 perreault78(Type1 E_mag, Type2 B_mag, Type3 theta_p) {
+    Type1 val = E_mag * B_mag / (4 * std::atan(1)*4) * std::pow((7 * RE * std::pow(sin(theta_p / 2), 2)), 2);
     return val;
 }
  
@@ -383,11 +383,11 @@ Type1 perreault78(Type1  E_mag, Type2 B_mag,Type3 theta_p) {
  *            theta IMF clock angle arctan(By/Bz)
  */
 template <typename Type1, typename Type2>
-Type1 newell2007(Type1  mu, std::vector<Type2> Bxyz) {
-    Type2 B_mag = std::sqrt(Bxyz[0]*Bxyz[0] + Bxyz[1]*Bxyz[1] + Bxyz[2]*Bxyz[2]);
-    Type2 theta = std::atan2(Bxyz[1]/Bxyz[2]);
-    Type1 val = std::pow(mu , 4.0/3.0)*std::pow(B_mag , 2.0/3.0)*pow(std::sin(theta/2),8.0/3.0);
-    return val; 
+Type1 newell2007(Type1 mu, std::vector<Type2> Bxyz) {
+    Type2 B_mag = std::sqrt(Bxyz[0] * Bxyz[0] + Bxyz[1] * Bxyz[1] + Bxyz[2] * Bxyz[2]);
+    Type2 theta = std::atan2(Bxyz[1] / Bxyz[2]);
+    Type1 val = std::pow(mu, 4.0 / 3.0) * std::pow(B_mag, 2.0 / 3.0) * pow(std::sin(theta / 2), 8.0 / 3.0);
+    return val;
 }
  
 /**
@@ -398,12 +398,81 @@ Type1 newell2007(Type1  mu, std::vector&lt;Type2&gt; Bxyz) {
  * R is the radial distance in RE and DTA is the dipole tilt angle 
  */
 template <typename Type1, typename Type2, typename Type3, typename Type4>
-Type1 lopez90(Type1  Kp, Type2 phi,Type3 R, Type4  DTA) {
-    Type1 val = -(0.14*Kp+0.69)*std::pow(phi, 0.3333)*(0.065*R-0.16)*DTA*DEG2RAD;
+Type1 lopez90(Type1 Kp, Type2 phi, Type3 R, Type4 DTA) {
+    Type1 val = -(0.14 * Kp + 0.69) * std::pow(phi, 0.3333)*(0.065 * R - 0.16) * DTA*DEG2RAD;
     return val;
 }
  
+template <typename Type1, typename Type2>
+Type1 pv_ionos_te(Type1 alt, Type2 sza) {
+
+    Type1 te;
+    if (alt > 3000 || alt < 150) {
+        te << NotANumber();
+        return te;
+    } else {
+
+        std::vector<double> pb;
+        pb.resize(4);
+        // 0 was added to fdsmodt in order to keep indexes as in fortran original source 
+        const std::vector< double> fdsmodt{0,
+            3.567296, 1.3836078e-02, 1.5086544e-02, 6.0729804e-03, 4.0306677e-03,
+            -1.3217842e-02, 7.8089219e-03, -2775.023, -123.7310, 8.828382,
+            -96.94563, 15.84634, 138.9812, -139.7511, -84.03277,
+            0.2649297, -1.444215, 0.6347786, -2.192531, -0.5396207,
+            0.6054179, 1.0569388e-04, -8.0908003e-06, 1.3877957e-05, -1.2327257e-05,
+            -1.1256760e-05, 1.3830228e-05, -1.4350858e-05};
+        double s = std::sin(sza);
+        double c = std::cos(sza);
+        double s2 = std::sin(2 * sza);
+        double c2 = std::cos(2 * sza);
+        double s3 = std::sin(3 * sza);
+        double c3 = std::cos(3 * sza);
+        for (int j = 1; j < 5; j++) {
+            int k = 7 * j;
+            pb[j - 1] = fdsmodt[k - 6] + fdsmodt[k - 5] * s + fdsmodt[k - 4] * c + fdsmodt[k - 3] * s2 + fdsmodt[k - 2] * c2 +
+                    fdsmodt[k - 1] * s3 + fdsmodt[k] * c3;
+        }
+        te = (Type1) pb[0] + pb[1] / ((alt + pb[2])*(alt + pb[2])) + pb[3] * alt;
+        return te;
+    }
+}
+
+template <typename Type1, typename Type2>
+Type1 pv_ionos_dn(Type1 alt, Type2 sza) {
  
+    Type1 dn;
+    if (alt > 3000 || alt < 150) {
+        dn << NotANumber();
+        return dn;
+    } else {
+        std::vector<double> bp;
+        bp.resize(4);
+        // 0 was added to fdsmodde in order to keep indexes as in fortran original source 
+        const std::vector<double> fdsmodde{0,
+            5.334022, 252.5587, 120.4444, 9.1163822E-02, -6.195264,
+            37.38898, 3.331818, 245.5529, 0.7587564, 0.1453792,
+            -0.9461438, 127.1485, -0.3391595, 2.4936652E-02, 0.1879423,
+            -5.3126566E-02, 2.9439656E-02, 1.654282, 415.7859, 2.6281483E-02,
+            145.3951, 2.045347, 0.4421963, 3.036170, 9.1009791E-04,
+            141.5824};
+        double s2 = std::sin(sza) * std::sin(sza);
+        double  szap = sza * (1. - fdsmodde[23] * s2 * std::exp(-fdsmodde[24] * s2 - fdsmodde[25]*(alt - fdsmodde[26])*(alt - fdsmodde[26])));
+        double sinz = std::sin(fdsmodde[18] + szap);
+        double s = fdsmodde[1] + fdsmodde[2] / (alt - fdsmodde[3])
+                + sinz * fdsmodde[13]
+                + fdsmodde[19] * std::exp(-fdsmodde[20]*((alt - fdsmodde[21])*(alt - fdsmodde[21])) - fdsmodde[22] * szap * szap);
+        double o = fdsmodde[4] + fdsmodde[5] / (alt - fdsmodde[6]) + sinz * fdsmodde[14];
+        double c = fdsmodde[7] + fdsmodde[8] / (alt - fdsmodde[9]) + sinz * fdsmodde[15];
+        double m = fdsmodde[10] + fdsmodde[11] / (alt - fdsmodde[12]) + sinz * fdsmodde[16];
+        double ss = s * cos(fdsmodde[17] + szap);
+        double calc = (1.77245 * ss * erfc(-ss) + std::exp(-ss * ss) + o);
+        if (calc <= 0)
+            calc = 0.00001;
+        dn = (Type1) c + m * log(calc);
+        return dn;
+    }
+}
  
 /**
  * @brief Operator cross product between  two  vector
@@ -1391,18 +1460,20 @@ std::vector&lt;Type&gt; total(const AMDA::Parameters::Tab2DData&lt;Type&gt;&amp; a, int mode) {
 }
  
 // Reshape a Tab2D to a vector
+
 template <typename Type>
 std::vector<Type> reshapeToVector(const AMDA::Parameters::Tab2DData<Type>& a) {
     std::vector<Type> v;
-    v.resize(a.getDim1Size()*a.getDim2Size());
+    v.resize(a.getDim1Size() * a.getDim2Size());
     for (int r = 0; r < a.getDim1Size(); ++r) {
         for (int c = 0; c < a.getDim2Size(); ++c) {
-            v[r*a.getDim2Size() + c] = a[r][c];
+            v[r * a.getDim2Size() + c] = a[r][c];
         }
     }
  
     return v;
 }
+
 /*
  * Evaluation of SuperMAG auroral electrojet indices as indicators
  * of substorms and auroral power
@@ -1410,7 +1481,7 @@ std::vector&lt;Type&gt; reshapeToVector(const AMDA::Parameters::Tab2DData&lt;Type&gt;&amp; a) {
  */
 template <typename Type>
 Type newell2011(Type a) {
-    return 0.048*a + 0.241*std::sqrt(a) ;
+    return 0.048 * a + 0.241 * std::sqrt(a);
 }
  
 template <typename Type1, typename Type2>
@@ -1490,15 +1561,14 @@ std::vector&lt;Type1&gt; operator/(AMDA::Parameters::TabRow&lt;Type1&gt; a, std::vector&lt;Type
 }
  
 template<typename Type>
-bool isValid(const std::string& num) { 
-   bool flag = true; 
-   try { 
-      Type tmp = boost::lexical_cast<Type>(num); 
-   } 
-   catch (boost::bad_lexical_cast &e) { 
-      flag = false; 
-   } 
-   return flag; 
+bool isValid(const std::string& num) {
+    bool flag = true;
+    try {
+        Type tmp = boost::lexical_cast<Type>(num);
+    } catch (boost::bad_lexical_cast &e) {
+        flag = false;
+    }
+    return flag;
 }
  
 template <typename Type>
@@ -1520,8 +1590,8 @@ std::pair&lt;Type, Type&gt; operator+(std::pair&lt;Type, Type&gt; a, std::pair&lt;Type, Type&gt; b
 template <typename Type>
 std::pair<Type, Type> operator*(std::pair<Type, Type> a, std::pair<Type, Type> b) {
     std::pair<Type, Type> c;
-    c.first = a.first* b.firest;
-    c.second = a.second* b.second;
+    c.first = a.first * b.firest;
+    c.second = a.second * b.second;
     return c;
 }
  
@@ -189,6 +189,30 @@
         </info_brief>
         <new_kernel>lopez90</new_kernel>
     </function>
+    <function name="pv_ionos_dn(,)" args="2" kind="model" group="space">
+        <prompt_param></prompt_param>
+        <default_args>alt(Km), sza(Radians)</default_args>
+        <info_brief>
+            Electron Density around Venus &lt;br/&gt;
+            dn = pv_ionos_dn(alt, sza) &lt;br/&gt;
+	    alt is the altitude in Km 
+	    sza is the solat zenith angle in radians 
+            dn is the base 10 log of the model electron density in no/cc
+        </info_brief>
+        <new_kernel>pv_ionos_dn</new_kernel>
+    </function>
+    <function name="pv_ionos_te(,)" args="2" kind="model" group="space">
+        <prompt_param></prompt_param>
+        <default_args>alt(Km), sza(Radians)</default_args>
+        <info_brief>
+            Electron temperature around Venus &lt;br/&gt;
+            te = pv_ionos_dn(alt, sza) &lt;br/&gt;
+            alt is the altitude in Km
+            sza is the solat zenith angle in radians
+            te is in log10 deg K
+        </info_brief>
+        <new_kernel>pv_ionos_te</new_kernel>
+    </function>
     <function name="framesTransformation(,,,)" args="1" kind="frames" group="space">
 	<prompts>
         	<prompt type="list" subtype="frames" default="GSM">Input frame:</prompt>