From a8c22dfe007d2dfacec85d78fcef556555437b64 Mon Sep 17 00:00:00 2001
From: ftomei <ftomei@arpae.it>
Date: Fri, 8 Mar 2024 18:26:10 +0100
Subject: [PATCH] factor of safety

---
 agrolib/soil/soil.cpp               |  2 +-
 bin/CRITERIA3D/shared/project3D.cpp | 47 ++++++++++++++++++++++++++---
 bin/CRITERIA3D/shared/project3D.h   |  2 +-
 3 files changed, 44 insertions(+), 7 deletions(-)

diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp
index b79c5a357..55018fb1c 100644
--- a/agrolib/soil/soil.cpp
+++ b/agrolib/soil/soil.cpp
@@ -759,7 +759,7 @@ namespace soil
 
         double frictionEffect =  tanFrictionAngle / tanAngle;
 
-        double unitWeight = horizonPtr->bulkDensity * GRAVITY;                   // [kN m-3]
+        double unitWeight = horizonPtr->bulkDensity * GRAVITY;                // [kN m-3]
         double cohesionEffect = 2 * (horizonPtr->effectiveCohesion + rootCohesion) / (unitWeight * depth * sin(2*slopeAngle));
 
         double suctionEffect = (suctionStress * (tanAngle + 1/tanAngle) * tanFrictionAngle) / (unitWeight * depth);
diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp
index 8a483394d..ae1dd17aa 100644
--- a/bin/CRITERIA3D/shared/project3D.cpp
+++ b/bin/CRITERIA3D/shared/project3D.cpp
@@ -1102,7 +1102,7 @@ bool Project3D::setCriteria3DMap(criteria3DVariable var, int layerIndex)
                 double value;
                 if (var == factorOfSafety)
                 {
-                    value = computeFactorOfSafety(nodeIndex);
+                    value = computeFactorOfSafety(row, col, layerIndex, nodeIndex);
                 }
                 else
                 {
@@ -1553,11 +1553,48 @@ double Project3D::assignTranspiration(int row, int col, double currentLai, doubl
 }
 
 
-float Project3D::computeFactorOfSafety(int nodeIndex)
+float Project3D::computeFactorOfSafety(int row, int col, int layerIndex, int nodeIndex)
 {
-    //TODO
-
-    return NODATA;
+    // suction stress [kPa]
+    double saturationDegree = soilFluxes3D::getDegreeOfSaturation(nodeIndex) / 100;     // [-]
+
+    // sign is just in matric potential
+    double waterPotential = soilFluxes3D::getMatricPotential(nodeIndex) * GRAVITY;      // [kPa]
+    waterPotential = std::min(0.0, waterPotential);
+    double suctionStress = waterPotential * saturationDegree;
+
+    // slope angle [rad]
+    float slopeDegree = radiationMaps->slopeMap->getValueFromRowCol(row, col);
+    float slopeAngle = slopeDegree * DEG_TO_RAD;
+
+    // friction angle [rad]
+    int soilIndex = getSoilIndex(row, col);
+    int horizonIndex = soil::getHorizonIndex(soilList[unsigned(soilIndex)], layerDepth[layerIndex]);
+    double frictionAngle = soilList[unsigned(soilIndex)].horizon[horizonIndex].frictionAngle * DEG_TO_RAD;
+
+    // effective cohesion [kPa]
+    double effectiveCohesion = soilList[unsigned(soilIndex)].horizon[horizonIndex].effectiveCohesion;
+
+    // friction effect [-]
+    double tanAngle = tan(slopeAngle);
+    double tanFrictionAngle = tan(frictionAngle);
+    double frictionEffect =  tanFrictionAngle / tanAngle;
+
+    // unit weight [kN m-3]
+    double bulkDensity = soilList[unsigned(soilIndex)].horizon[horizonIndex].bulkDensity;  // [g cm-3] --> [Mg m-3]
+    double unitWeight = bulkDensity * GRAVITY;
+
+    // cohesion effect [-]
+    double depth = layerDepth[layerIndex];   // [m]
+    // TODO root cohesion [kPa] leggere da db e assegnare in base alla ratio di root density
+    double rootCohesion = 0.;
+    double cohesionEffect = 2 * (effectiveCohesion + rootCohesion) / (unitWeight * depth * sin(2*slopeAngle));
+
+    // suction effect [-]
+    double suctionEffect = (suctionStress * (tanAngle + 1/tanAngle) * tanFrictionAngle) / (unitWeight * depth);
+
+    // factor of safety [-]
+    return frictionEffect + cohesionEffect - suctionEffect;
 }
 
 
diff --git a/bin/CRITERIA3D/shared/project3D.h b/bin/CRITERIA3D/shared/project3D.h
index 3fb14c7c9..1c8abb592 100644
--- a/bin/CRITERIA3D/shared/project3D.h
+++ b/bin/CRITERIA3D/shared/project3D.h
@@ -169,7 +169,7 @@
 
         bool setCriteria3DMap(criteria3DVariable var, int layerIndex);
 
-        float computeFactorOfSafety(int nodeIndex);
+        float computeFactorOfSafety(int row, int col, int layerIndex, int nodeIndex);
     };
 
     bool isCrit3dError(int result, QString &error);