New sigma_y and sigma_z for plume gaussian model

src/tests/plot/dispersion.hpp

@@ -27,10 +27,11 @@
 #ifndef TESTS_PLOT_DISPERSION_HPP
 #define TESTS_PLOT_DISPERSION_HPP 1
 
+#include <cmath>
+
 namespace paradevs { namespace tests { namespace plot {
 
 #define SCALE 4
-#define M_PI 3.14159265358979323846 
 
 class KleinDispersionFunction
 {
@@ -86,83 +87,48 @@ public:
         double wind_speed,
         double wind_direction) const
     {
-		paradevs::tests::boost_graph::Point destination1,destination2;
-		
-		//Modification à apporter pour prendre en compte la vitesse du vent
-		//Valeurs arbitraire à modifier en fonction du parcellaire étudié
-		double d = 400; //distance entre la source et le point d'infection maximum
-		double sigma_x = 500; //dispertion sur l'axe x
-		double moy_x   = d/std::cos(wind_direction);
-		
-		destination1  = destination;
-		
-		//Changement de repère centré en (source._x, source._y)
-        destination1._x -= source._x;
-        destination1._y -= source._y;
-        
-        //Rotation des axes d'angle wind_direction
-        destination2._x = destination1._x*std::cos(wind_direction) + 
-						  destination1._y*std::sin(wind_direction);
-		destination2._y = -destination1._x*std::sin(wind_direction) + 
-						  destination1._y*std::cos(wind_direction);
-        
-        double distance = std::sqrt(
-            (destination2._x) * (destination2._x) +
-            (destination2._y) * (destination2._y));
-        
-        std::pair<double,double> param_y;
-        int caseSwitch = 1;
-		switch (caseSwitch)
-		{
-		    case 1:
-		        param_y.first  = 0.495;
-		        param_y.second = 0.873;
-		        break;
-		    case 2:
-		        param_y.first  = 0.310;
-		        param_y.second = 0.897;
-		        break;
-		    case 3:
-		        param_y.first  = 0.197;
-		        param_y.second = 0.908;
-		        break;
-		    case 4:
-		        param_y.first  = 0.122;
-		        param_y.second = 0.916;
-		        break;
-		    case 5:
-		        param_y.first  = 0.0934;
-		        param_y.second = 0.912;
-		        break;
-		    case 6:
-		        param_y.first  = 0.0625;
-		        param_y.second = 0.911;
-		        break;
-		    default:
-		        std::cout << "Problem ! Choose a good case !" << std::endl;
-		        break;
-		}
-         
-        double sigma_y = param_y.first * std::pow(distance, param_y.second);
-        //Problème si calcul de dispertion intra parcelle -> voir modification
-        
-        double disp = 1/(M_PI*sigma_x*sigma_y) * 
-					  std::exp(-0.5*std::pow((destination2._x-moy_x)/sigma_x,2)) *
-					  std::exp(-0.5*std::pow(destination2._y/sigma_y,2));
-					  
-		/* Formule extraite du papier sur plume / par la vitesse du vent 
-		 * (mais je ne vois pas ce que ça apporte)
-		 * double disp = std::pow(10,3)/(M_PI*sigma_x*sigma_y*wind_speed) * 
-					  std::exp(-0.5*std::pow((destination2._x-moy_x)/sigma_x,2)) *
-					  std::exp(-0.5*std::pow(destination2._y/sigma_y,2));
-		 */
-		
-		//Problème avec la normalisation de disp, manque d'information sur les voisins ...
-		
-		/** ATTENTION **/ 
-		//Algorithme non validé et non testé, juste compilé.
-       
-        return disp;
+        paradevs::tests::boost_graph::Point destination1, destination2;
+
+        double delta_x = destination._x - source._x;
+
+        if (delta_x > 0) {
+            double delta_y = destination._y - source._y;
+            double distance = std::sqrt(delta_x * delta_x + delta_y * delta_y) *
+                SCALE;
+
+            destination1  = destination;
+
+            //Changement de repère centré en (source._x, source._y)
+            destination1._x -= source._x;
+            destination1._y -= source._y;
+
+            //Rotation des axes d'angle wind_direction
+            destination2._x = destination1._x * std::cos(wind_direction) +
+                destination1._y * std::sin(wind_direction);
+            destination2._y = -destination1._x * std::sin(wind_direction) +
+                destination1._y * std::cos(wind_direction);
+
+            // double sigma_y = 0.0787*x/(std::pow(1+x/707,0.135));
+            // double sigma_z = 0.0475*x/(std::pow(1+x/707, 0.465));
+
+            double sigma_y = 0.08 * destination2._x * SCALE *
+                std::pow(1 + 0.0001 * destination2._x * SCALE, -0.5);
+            double sigma_z = 0.06 * destination2._x * SCALE *
+                std::pow(1 + 0.0015 * destination2._x * SCALE , -0.5);
+
+            double h = 20;
+            double disp = 1. / (2 * M_PI * wind_speed * sigma_y * sigma_z) *
+                std::exp(-0.5 * std::pow(h / sigma_z, 2)) *
+                std::exp(-0.5 * std::pow(destination2._y * SCALE / sigma_y, 2));
+
+            // std::cout << (destination2._x * SCALE) << " "
+            //           << (destination2._y * SCALE) << " "
+            //           << (disp * 3600 * 24) << std::endl;
+
+            return disp * 3600 * 24 * 1e4;
+        } else {
+            return 0.0;
+        }
     }
 };
 

src/tests/plot/main.cpp

@@ -36,7 +36,7 @@ using namespace paradevs::common;
 using namespace std::chrono;
 using namespace paradevs::tests::plot;
 
-double plot_monothreading()
+double plot_monothreading(const std::string& shapefile)
 {
     paradevs::common::RootCoordinator <
         DoubleTime, paradevs::pdevs::Coordinator <
@@ -44,9 +44,8 @@ double plot_monothreading()
             HierarchicalGraphManager,
             paradevs::common::NoParameters,
             GraphManagerParameters >
-        > rc(0, 60, "root", paradevs::common::NoParameters(),
-             GraphManagerParameters("/home/herbez/Documents/Thèse/parcelle/oye_plage.shp",
-                                    4));
+        > rc(0, 200, "root", paradevs::common::NoParameters(),
+             GraphManagerParameters(shapefile, 4));
 
     steady_clock::time_point t1 = steady_clock::now();
 
@@ -60,7 +59,7 @@ double plot_monothreading()
     return time_span.count();
 }
 
-double plot_multithreading(int cluster_number)
+double plot_multithreading(int cluster_number, const std::string& shapefile)
 {
     paradevs::common::RootCoordinator <
         DoubleTime, paradevs::pdevs::multithreading::Coordinator <
@@ -69,8 +68,7 @@ double plot_multithreading(int cluster_number)
             paradevs::common::NoParameters,
             GraphManagerParameters >
         > rc(0, 10, "root", paradevs::common::NoParameters(),
-             GraphManagerParameters("/home/herbez/Documents/Thèse/parcelle/oye_plage.shp",
-                                    cluster_number));
+             GraphManagerParameters(shapefile, cluster_number));
 
     steady_clock::time_point t1 = steady_clock::now();
 
@@ -84,19 +82,19 @@ double plot_multithreading(int cluster_number)
     return time_span.count();
 }
 
-void plot(int n)
+void plot(int n, const std::string& shapefile)
 {
     if (n == 1) {
-        std::cout << plot_monothreading() << std::endl;
+        std::cout << plot_monothreading(shapefile) << std::endl;
     } else {
-        std::cout << plot_multithreading(n) << std::endl;
+        std::cout << plot_multithreading(n, shapefile) << std::endl;
     }
 }
 
 int main(int argc, char** argv)
 {
-    if (argc > 1) {
-        plot(atoi(argv[1]));
+    if (argc > 2) {
+        plot(atoi(argv[1]), argv[2]);
     }
     return 0;
 }

src/tests/plot/models.hpp

@@ -36,6 +36,8 @@
 
 #include <shapefil.h>
 
+#include <cmath>
+
 namespace paradevs { namespace tests { namespace plot {
 
 struct PlotParameters
@@ -137,31 +139,33 @@ public:
                     sum += _dispersion_function(
                         t,
                         paradevs::tests::boost_graph::Point(it->x, it->y),
-                        _centroid, 0., 0.) * it->ready_spore_number;
+                        _centroid, 6., 0.) * it->ready_spore_number;
 
                     // std::cout << (_index + 1) << " " << (it->index + 1) << " " << it->concentration << std::endl;
 
                 }
             }
-            sum += _dispersion_function(t, _centroid, _centroid, 0., 0.) *
+            sum += _dispersion_function(
+                t, _centroid,
+                paradevs::tests::boost_graph::Point(_centroid._x + 10,
+                                                    _centroid._y), 6., 0.) *
                 _data.ready_spore_number;
             _milsol.add_zoospore_number(sum);
 
             double p = 0;
 
             if (_milsol.get_zoospore_number() > 0) {
-                p =  std::log10(_milsol.get_zoospore_number());
+                p = std::log10(_milsol.get_zoospore_number());
             }
 
-            // std::cout << t << "\t" << get_name() << "\t"
+            // std::cout << t << "\t" << _index << "\t"
             //           << _data.ready_spore_number << "\t"
             //           << _milsol.get_zoospore_number() << "\t"
+            //           << sum << "\t"
             //           << p << std::endl;
 
             DBFWriteDoubleAttribute(_handle, _index, 0, p);
 
-            // std::cout << (_index + 1) << " " << _data.concentration << std::endl;
-
             _phase  = SEND;
             _sigma = 1;
             _received = 0;
@@ -211,7 +215,7 @@ public:
         typename common::DoubleTime::type /* t */)
     {
         if (_index == 210) {
-            _milsol.set_inoculum_primaire_number(10);
+            _milsol.set_inoculum_primaire_number(1e6);
         } else {
             _milsol.set_inoculum_primaire_number(0);
         }
@@ -263,7 +267,8 @@ private:
     std::vector < PlotData > _neighbour_data;
 
     // submodels
-    KleinDispersionFunction _dispersion_function;
+    Plume2dDispersionFunction _dispersion_function;
+    // KleinDispersionFunction _dispersion_function;
     Milsol                  _milsol;
     Climate                 _climate;