Changeset 60828 in spip-zone

Timestamp:

Apr 30, 2012, 6:45:46 PM (9 years ago)

Author:

dorch@…

Message:

Développement des calculs d'ouvrages. Toutes les paramètres peuvent être en calcul. Lois 1,3,4,5 programmées.

File:

• _plugins_/hydraulic/branches/v0.3/hyd_inc/ouvrage.class.php (modified) (3 diffs)

_plugins_/hydraulic/branches/v0.3/hyd_inc/ouvrage.class.php

@@ -53 +53 @@
      * - CT : Coefficient de débit partie triangulaire pour les trapézoïdales
      * - CS : Coefficient de débit de la surverse
+     * - P : Précision du calcul
      */
     private $tP = array();
 
-    const G = 9.81;
-    const R2G = 4.42944; //sqrt(2*self::gP);
-    const R32 = 2.59807; // 3*sqrt(3)/2;
+    const G = 9.81; /// Constante de gravité terrestre
+    const R2G = 4.42944; /// sqrt(2*self::gP);
+    const R32 = 2.59807; /// 3*sqrt(3)/2;
+    const IDEFINT = 100; /// Pas de parcours de l'intervalle pour initialisation dichotomie
+    const IDICMAX = 100; /// Itérations maximum de la dichotomie
 
     /**
@@ -104 +107 @@
      * @return array(0=> donnée calculée, 1=> Flag d'écoulement)
      */
-    public function Calc($sCalc,$rInit=0) {
-        $rVarC = &$this->tP[$sCalc];
-
+    public function Calc($sCalc,$rInit=0.) {
+        //print_r($this->tP);
         // Calcul du débit (facile !)
-        if($sCalc=='Q') return $this->OuvrageQ();
-
-        // Sinon calcul d'une autre donnée par dichotomie
-
+        if($sCalc=='Q') {
+            return $this->OuvrageQ();
+        }
+        else {
+            // Sinon calcul d'une autre donnée par dichotomie
+            $rVarC = &$this->tP[$sCalc];
+            $QT = $this->tP['Q']; // Débit recherché (Target)
+            $XMinInit = 0;
+            $rVarC = $XMinInit;
+            list($Q1,$nFlag) = $this->OuvrageQ();
+            $XMaxInit = $rInit*10; /// @todo Boucler la valeur max sur 10,100,1000,10000
+            $rVarC = $XMaxInit;
+            list($Q2,$nFlag) = $this->OuvrageQ();
+            $DX = ($XMaxInit - $XMinInit) / floatval(self::IDEFINT);
+            $nIterMax = floor(max($XMaxInit - $rInit,$rInit - $XMinInit) / $DX + 1);
+            $Xmin = $rInit;
+            $Xmax = $rInit;
+            $X1 = $rInit;
+            $X2 = $rInit;
+            $rVarC = $rInit;
+            list($Q,$nFlag) = $this->OuvrageQ();
+            $Q1 = $Q;
+            $Q2 = $Q;
+            //echo "\n".'nIterMax='.$nIterMax.'  XMinInit='.$XMinInit.'  XMaxInit='.$XMaxInit.'  DX='.$DX;
+
+
+            for($nIter=1;$nIter<=$nIterMax;$nIter++) {
+                //Ouverture de l'intervalle des deux côtés puis à droite et à gauche
+                $Xmax = $Xmax + $DX;
+                if($Xmax > $XMaxInit xor $DX <= 0) $Xmax = $XMaxInit;
+                $rVarC = $Xmax;
+                list($Q,$nFlag) = $this->OuvrageQ();
+                if($Q1 < $Q2 xor $Q <= $Q2) {
+                    $Q2 = $Q;
+                    $X2 = $Xmax;
+                }
+                if($Q1 < $Q2 xor $Q >= $Q1) {
+                    $Q1 = $Q;
+                    $X1 = $Xmax;
+                }
+                $Xmin = $Xmin - $DX;
+                if($Xmin < $XMinInit xor $DX <= 0) {
+                    $Xmin = $XMinInit;
+                }
+                $rVarC = $Xmin;
+                list($Q,$nFlag) = $this->OuvrageQ();
+                if($Q1 < $Q2 xor $Q <= $Q2) {
+                    $Q2 = $Q;
+                    $X2 = $Xmin;
+                }
+                if($Q1 < $Q2 xor $Q >= $Q1) {
+                    $Q1 = $Q;
+                    $X1 = $Xmin;
+                }
+/*
+                echo "\n".'nIter='.$nIter.' Xmin='.$Xmin.' Xmax='.$Xmax;
+                echo "\n".'X1='.$X1.' Q1='.$Q1.' X2='.$X2.' Q2='.$Q2;
+                echo "\n".'$QT > $Q1 xor $QT >= $Q2 = '.($QT > $Q1 xor $QT >= $Q2);
+*/
+                if($QT > $Q1 xor $QT >= $Q2) {break;}
+            }
+
+            if($nIter >= self::IDEFINT) {
+                // Pas d'intervalle trouvé avec au moins une solution
+                if($Q2 < $QT and $Q1 < $QT) {
+                    // Cote de l'eau trop basse pour passer le débit il faut ouvrir un autre ouvrage
+                    $rVarC = $XmaxInit;
+                }
+                else {
+                    // Cote de l'eau trop grande il faut fermer l'ouvrage
+                    $rVarC = $XminInit;
+
+                }
+                list($Q,$nFlag) = $this->OuvrageQ();
+                $nFlag = -1;
+            }
+            else {
+                // Dichotomie
+                $X = $rInit;
+                for($nIter = 1; $nIter<=self::IDICMAX;$nIter++) {
+                    $rVarC=$X;
+                    list($Q,$nFlag) = $this->OuvrageQ();
+                    if(abs($Q/$QT-1.) <= $this->tP['P']) {break;}
+                    if($QT < $Q xor $Q1 <= $Q2) {
+                        // QT < IQ et Q(X1) > Q(X2) ou pareil en inversant les inégalités
+                        $X1=$rVarC;
+                    }
+                    else {
+                        // QT < IQ et Q(X1) < Q(X2) ou pareil en inversant les inégalités
+                        $X2=$rVarC;
+                    }
+                    $X=($X2+$X1)*0.5;
+                }
+                if($nIter == self::IDICMAX) {
+                    //IF1 <-- -10 anomalie: la dichotomie n'a pas abouti en ITER iterations
+                    $nFlag = -1;
+                }
+            }
+        }
         return array($rVarC,$nFlag);
     }
@@ -180 +277 @@
             $this->tP['ZM'] = $ZM;
         }
-
+        //echo "\n".'OuvrageQ='.$rQ.' / '.$nFlag;
         return array($rQ,$nFlag);
     }