// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

#include "IsoHyper3DOrgeas1.h"
//#include "ComLoi_comp_abstraite.h"

# include <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "CharUtil.h"

//================== initialisation des variables static ======================
// indicateur utilisé par Verif_Potentiel_et_var  
int IsoHyper3DOrgeas1::indic_Verif_PoGrenobleorgeas1_et_var = 0;                            
// indicateur utilisé par Verif_Potentiel_et_var_VV  
int IsoHyper3DOrgeas1::indic_Verif_PoGrenobleorgeas1_et_var_VV = 0; 
// indicateur utilisé par CompaFormel                             
int IsoHyper3DOrgeas1::indic_Verif_CompaFormel = 0;                          
//================== fin d'initialisation des variables static ================


IsoHyper3DOrgeas1::IsoHyper3DOrgeas1 ()  : // Constructeur par defaut
 Hyper3D(ISOHYPER3DORGEAS1,CAT_MECANIQUE,false),K(ConstMath::tresgrand),Q0s(ConstMath::tresgrand)
 ,mu01(ConstMath::tresgrand),mu02(ConstMath::tresgrand),mu03(ConstMath::tresgrand)
 ,alpha1(ConstMath::tresgrand),alpha3(ConstMath::tresgrand),Q0e(ConstMath::tresgrand)
 ,alpha1_2(0.),alpha3_2(0.) //,mu1_2(0.)
 ,nQs(1.),gammaQs(0.),nQe(1.),gammaQe(0.),nMu1(1.),gammaMu1(0.),nMu2(1.),gammaMu2(0.),nMu3(1.),gammaMu3(0.)
 ,cas_Q0s_thermo_dependant(0),cas_Q0e_thermo_dependant(0)
 ,T0rQe(ConstMath::tresgrand),Qe_T0rQe(ConstMath::tresgrand),h1(ConstMath::tresgrand)
 ,T0r(ConstMath::tresgrand),Gr(ConstMath::tresgrand)
 ,Qrmax(ConstMath::tresgrand),ar(ConstMath::tresgrand),Tc(0.),Ts(0.),Qs_T0rQe(0.)
 ,Q0s_temperature(NULL),Q0e_temperature(NULL)
       {};
// Constructeur de copie
IsoHyper3DOrgeas1::IsoHyper3DOrgeas1 (const IsoHyper3DOrgeas1& loi) :
  Hyper3D (loi),K(loi.K),Q0s(loi.Q0s),mu01(loi.mu01),mu02(loi.mu02)
  ,mu03(loi.mu03),alpha1(loi.alpha1),alpha3(loi.alpha3),Q0e(loi.Q0e)
 ,alpha1_2(loi.alpha1_2)
 ,alpha3_2(loi.alpha3_2) //,mu1_2(loi.mu1_2)
 ,nQs(loi.nQs),gammaQs(loi.gammaQs),nQe(loi.nQe),gammaQe(loi.gammaQe)
 ,nMu1(loi.nMu2),gammaMu1(loi.gammaMu2),nMu2(loi.nMu2),gammaMu2(loi.gammaMu2),nMu3(loi.nMu3),gammaMu3(loi.gammaMu3)
 ,cas_Q0s_thermo_dependant(loi.cas_Q0s_thermo_dependant),cas_Q0e_thermo_dependant(loi.cas_Q0e_thermo_dependant)
 ,T0rQe(loi.T0rQe),Qe_T0rQe(loi.Qe_T0rQe),h1(loi.h1)
 ,T0r(loi.T0r),Gr(loi.Gr),Qrmax(loi.Qrmax),ar(loi.ar)
 ,Tc(loi.Tc),Ts(loi.Ts),Qs_T0rQe(loi.Qs_T0rQe)
 ,Q0s_temperature(loi.Q0s_temperature),Q0e_temperature(loi.Q0e_temperature)
	 {// on regarde s'il s'agit d'une courbe locale ou d'une courbe globale pour Qs(T)
	  if (Q0s_temperature != NULL)
	   if (Q0s_temperature->NomCourbe() == "_")
        {// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
         string non_courbe("_");   
         Q0s_temperature = Courbe1D::New_Courbe1D(*loi.Q0s_temperature);
        };
	  // on regarde s'il s'agit d'une courbe locale ou d'une courbe globale pour Qe(T)
	  if (Q0e_temperature != NULL)
	   if (Q0e_temperature->NomCourbe() == "_")
        {// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
         string non_courbe("_");   
         Q0e_temperature = Courbe1D::New_Courbe1D(*loi.Q0e_temperature);
        };
  };

IsoHyper3DOrgeas1::~IsoHyper3DOrgeas1 ()
// Destructeur 
{ if (Q0s_temperature != NULL)
      if (Q0s_temperature->NomCourbe() == "_") delete Q0s_temperature;
  if (Q0e_temperature != NULL)
      if (Q0e_temperature->NomCourbe() == "_") delete Q0e_temperature;
};
  

// Lecture des donnees des paramètres matériau sur fichier
void IsoHyper3DOrgeas1::LectureDonneesParticulieres
      (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D,LesFonctions_nD& lesFonctionsnD)
  { // lecture des quatres coefficients de la loi
    *(entreePrinc->entree) >> K ; 
    // on regarde s'il y a une thermo-dépendance pour Q0s
    if (strstr(entreePrinc->tablcar,"Q0s_thermo_dependant_")!=NULL)
     { string nom; 
       *(entreePrinc->entree) >> nom;
       if (nom != "Q0s_thermo_dependant_") // vérification
        {  cout << "\n erreur en lecture du parametre Q0s, on attendait le mot cle "
    	           << " Q0s_thermo_dependant_ et on a lue: " << nom;
           entreePrinc->MessageBuffer("**1--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
        };
        Q0s=0.; // init néanmoins de la valeur
        cas_Q0s_thermo_dependant =1; // on signale qu'il y a thermo dépendance (= 1 n'est pas définitif)
     }
    else // sinon cas sans thermo-dépendance
     { *(entreePrinc->entree) >> Q0s; };
    
    // puis lecture des autres paramètres
     
    *(entreePrinc->entree) >> mu01 >> mu02 >> mu03 
                           >> alpha1 >> alpha3; 
    if ((alpha1 <= 0.) || (alpha3 <= 0.) )
       { cout << "\n erreur en lecture des parametres alpha1= " << alpha1 << " et alpha3 = " << alpha3
    	      << " , ces deux parametres doivent etre strictement positifs pour eviter des / par 0 ";
         entreePrinc->MessageBuffer("**2--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
         throw (UtilLecture::ErrNouvelleDonnee(-1));
         Sortie(1);     
       };
    if (mu01 < ConstMath::petit)    
       { cout << "\n erreur en lecture du parametre mu1 " << mu01 << " ce potentiel n'est pas prevu pour "
    	      << " fonctionner avec un mu1 quasi-nul ";
         entreePrinc->MessageBuffer("**3--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
         throw (UtilLecture::ErrNouvelleDonnee(-1));
         Sortie(1);     
       };
    if (mu03 < ConstMath::petit)    
       { cout << "\n erreur en lecture du parametre mu3 " << mu03 << " ce potentiel n'est pas prevu pour "
    	      << " fonctionner avec un mu3 quasi-nul ";
         entreePrinc->MessageBuffer("**4--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
         throw (UtilLecture::ErrNouvelleDonnee(-1));
         Sortie(1);     
       };
		  
    // on regarde s'il y a une thermo-dépendance pour Q0e
    if (strstr(entreePrinc->tablcar,"Q0e_thermo_dependant_")!=NULL)
     { string nom; 
       *(entreePrinc->entree) >> nom;
       if (nom != "Q0e_thermo_dependant_") // vérification
        {  cout << "\n erreur en lecture du parametre Q0e, on attendait le mot cle "
    	           << " Q0e_thermo_dependant_ et on a lue: " << nom;
           entreePrinc->MessageBuffer("**1_1--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
        };
       Q0e=0.; // init néanmoins de la valeur
       cas_Q0e_thermo_dependant =1; // on signale qu'il y a thermo dépendance (= 1 n'est pas définitif)
     }
    else // sinon cas sans thermo-dépendance
     { *(entreePrinc->entree) >> Q0e; };
                               
    // on regarde s'il y a dépendance à la phase                       
    if (strstr(entreePrinc->tablcar,"avec_phase")!=NULL)
     { string nom1,nom2;
       // lecture des  paramètres de phase
       entreePrinc->NouvelleDonnee();  // passage d'une ligne
       // cas de Qsigrev
       if (strstr(entreePrinc->tablcar,"nQs=")!=NULL)
        { *(entreePrinc->entree) >> nom1 >> nQs >> nom2 >> gammaQs;
          if ((nom1 != "nQs=") || (nom2 != "gammaQs=")) 
           	{ cout << "\n erreur en lecture des parametres de phase pour Q0sig_rev, on attendait les mots cles "
    	           << " nQs= et gammaQs= et on a lue: " << nom1 << " " << nom2;
              entreePrinc->MessageBuffer("**5--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
    	    };
        };
       // cas de Q0epsrev
       if (strstr(entreePrinc->tablcar,"nQe=")!=NULL)
        { *(entreePrinc->entree) >> nom1 >> nQe >> nom2 >> gammaQe;
          if ((nom1 != "nQe=") || (nom2 != "gammaQe=")) 
           	{ cout << "\n erreur en lecture des parametres de phase pour Qeps_rev, on attendait les mots cles "
    	           << " nQe= et gammaQe= et on a lue: " << nom1 << " " << nom2;
              entreePrinc->MessageBuffer("**6--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
    	    };
        };
       // cas de mu1rev
       if (strstr(entreePrinc->tablcar,"nMu1=")!=NULL)
        { *(entreePrinc->entree) >> nom1 >> nMu1 >> nom2 >> gammaMu1;
          if ((nom1 != "nMu1=") || (nom2 != "gammaMu1=")) 
           	{ cout << "\n erreur en lecture des parametres de phase pour mu1_rev, on attendait les mots cles "
    	           << " nMu1= et gammaMu1= et on a lue: " << nom1 << " " << nom2;
              entreePrinc->MessageBuffer("**71--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
    	    };
        };
       // cas de mu2rev
       if (strstr(entreePrinc->tablcar,"nMu2=")!=NULL)
        { *(entreePrinc->entree) >> nom1 >> nMu2 >> nom2 >> gammaMu2;
          if ((nom1 != "nMu2=") || (nom2 != "gammaMu2=")) 
           	{ cout << "\n erreur en lecture des parametres de phase pour mu2_rev, on attendait les mots cles "
    	           << " nMu2= et gammaMu2= et on a lue: " << nom1 << " " << nom2;
              entreePrinc->MessageBuffer("**7--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
    	    };
        };
       // cas de mu3rev
       if (strstr(entreePrinc->tablcar,"nMu3=")!=NULL)
        { *(entreePrinc->entree) >> nom1 >> nMu3 >> nom2 >> gammaMu3;
          if ((nom1 != "nMu3=") || (nom2 != "gammaMu3=")) 
           	{ cout << "\n erreur en lecture des parametres de phase pour mu3_rev, on attendait les mots cles "
    	           << " nMu3= et gammaMu3= et on a lue: " << nom1 << " " << nom2;
              entreePrinc->MessageBuffer("**8--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
    	    };
        };
       avec_phase=true;
     };
    // ---- maintenant on traite le cas des thermo-dépendances
	 // .... cas de Q0s .......
    if (cas_Q0s_thermo_dependant != 0)
     { thermo_dependant=true; // on indique que la loi est thermo-dépendant
       string nom; 
       entreePrinc->NouvelleDonnee();  // passage d'une ligne
       *(entreePrinc->entree) >> nom ;
       if (nom != "cas_Q0s_thermo_dependant_") // vérification
        {  cout << "\n erreur en lecture du parametre cas_Q0s_thermo_dependant, on attendait le mot cle "
    	           << " cas_Q0s_thermo_dependant_ et on a lue: " << nom;
           entreePrinc->MessageBuffer("**9--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
        };
       *(entreePrinc->entree) >> cas_Q0s_thermo_dependant;
       if (cas_Q0s_thermo_dependant == 1)
         {// cas de T0r
          *(entreePrinc->entree) >> nom ;
          if (nom != "T0r=") // vérification
           {  cout << "\n erreur en lecture du parametre T0r, on attendait le mot cle "
       	           << " T0r= et on a lue: " << nom;
              entreePrinc->MessageBuffer("**10--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
           };
          *(entreePrinc->entree) >> T0r;
          // cas de Gr
          *(entreePrinc->entree) >> nom ;
          if (nom != "Gr=") // vérification
           {  cout << "\n erreur en lecture du parametre Gr, on attendait le mot cle "
       	           << " Gr= et on a lue: " << nom;
              entreePrinc->MessageBuffer("**11--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
           };
          *(entreePrinc->entree) >> Gr;
          // cas de Qrmax
          *(entreePrinc->entree) >> nom ;
          if (nom != "Qrmax=") // vérification
           {  cout << "\n erreur en lecture du parametre Qrmax, on attendait le mot cle "
       	           << " Qrmax= et on a lue: " << nom;
              entreePrinc->MessageBuffer("**12--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
           };
          *(entreePrinc->entree) >> Qrmax;
          // cas de ar
          *(entreePrinc->entree) >> nom ;
          if (nom != "ar=") // vérification
           {  cout << "\n erreur en lecture du parametre ar, on attendait le mot cle "
       	           << " ar= et on a lue: " << nom;
              entreePrinc->MessageBuffer("**13--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
           };
          *(entreePrinc->entree) >> ar;
          // on calcul maintenant les variables intermédiaires fixes
          Tc=T0r+(Qrmax/Gr/2.);
          Ts=(ar*ar- Sqr(Qrmax/Gr) )/2./(Qrmax/Gr);
         }
       else if (cas_Q0s_thermo_dependant == 2)
       	{// lecture de la loi d'évolution de Q0s en fonction de la température
         *(entreePrinc->entree) >>  nom;    
         // on regarde si la courbe existe, si oui on récupère la référence 
         if (lesCourbes1D.Existe(nom)) 
           { Q0s_temperature = lesCourbes1D.Trouve(nom);
            }
         else
          { // sinon il faut la lire maintenant
            string non_courbe("_");   
            Q0s_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
            // lecture de la courbe
            Q0s_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
           } 
       	}
       else
       	{ cout << "\n erreur en lecture du type de dependance de Q0s de la temperature "
       	       << " actuellement seules les cas 1 et 2 sont possible et on a lue: " << cas_Q0s_thermo_dependant;
              entreePrinc->MessageBuffer("**14--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
       	};	
     };

	 // .... cas de Q0e .......
    if (cas_Q0e_thermo_dependant != 0)
     { thermo_dependant=true; // on indique que la loi est thermo-dépendant
       string nom; 
       entreePrinc->NouvelleDonnee();  // passage d'une ligne
       *(entreePrinc->entree) >> nom ;
       if (nom != "cas_Q0e_thermo_dependant_") // vérification
        {  cout << "\n erreur en lecture du parametre cas_Q0e_thermo_dependant, on attendait le mot cle "
    	           << " cas_Q0e_thermo_dependant_ et on a lue: " << nom;
           entreePrinc->MessageBuffer("**15--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
        };
       *(entreePrinc->entree) >> cas_Q0e_thermo_dependant;
       if (cas_Q0e_thermo_dependant == 1)
         {// cas de T0rQe
          *(entreePrinc->entree) >> nom ;
          if (nom != "T0rQe=") // vérification
           {  cout << "\n erreur en lecture du parametre T0rQe, on attendait le mot cle "
       	           << " T0rQe= et on a lue: " << nom;
              entreePrinc->MessageBuffer("**16--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
           };
          *(entreePrinc->entree) >> T0rQe;
          // cas de Qe_T0rQe
          *(entreePrinc->entree) >> nom ;
          if (nom != "Qe_T0rQe=") // vérification
           {  cout << "\n erreur en lecture du parametre Qe_T0rQe, on attendait le mot cle "
       	           << " Qe_T0rQe= et on a lue: " << nom;
              entreePrinc->MessageBuffer("**17--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
           };
          *(entreePrinc->entree) >> Qe_T0rQe;
          // cas de h1
          *(entreePrinc->entree) >> nom ;
          if (nom != "h1=") // vérification
           {  cout << "\n erreur en lecture du parametre h1, on attendait le mot cle "
       	           << " h1= et on a lue: " << nom;
              entreePrinc->MessageBuffer("**17--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
           };
          *(entreePrinc->entree) >> h1;
         }
       else if (cas_Q0e_thermo_dependant == 2)
       	{// lecture de la loi d'évolution de Q0e en fonction de la température
         *(entreePrinc->entree) >>  nom;    
         // on regarde si la courbe existe, si oui on récupère la référence 
         if (lesCourbes1D.Existe(nom)) 
           { Q0e_temperature = lesCourbes1D.Trouve(nom);
            }
         else
          { // sinon il faut la lire maintenant
            string non_courbe("_");   
            Q0e_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
            // lecture de la courbe
            Q0e_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
           } 
       	}
       else
       	{ cout << "\n erreur en lecture du type de dependance de Q0e de la temperature "
       	       << " actuellement seules les cas 1 et 2 sont possible et on a lue: " << cas_Q0e_thermo_dependant;
              entreePrinc->MessageBuffer("**18--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
       	};	
       // maintenant on calcule la valeur de Qs_T0rQr
       {Qs_T0rQe = Q0s; // au cas où il n'y a pas de dépendance de Qs à la température
          switch (cas_Q0s_thermo_dependant)
            {case 1:// cas d'une dépendance à la température selon fct laurent
              { if (T0rQe >= Tc)
                  {Qs_T0rQe=0.5*Gr*( (T0rQe-Tc+Ts)-sqrt(Sqr(T0rQe-Tc+Ts) + ar*ar)) + Qrmax;}
                else
                  {Qs_T0rQe=0.5*Gr*( (T0rQe-Tc-Ts)+sqrt(Sqr(T0rQe-Tc-Ts) + ar*ar));};
                break;
              }
             case 2:// cas d'une dépendance à la température selon une fonction d'évolution
              { Qs_T0rQe = Q0s_temperature->Valeur(T0rQe);
                break;
              }
            };
       };	
     };
     
    // Lecture des paramètre specifique sur fichier (exe: regularisation sortie_post...
    Hyper3D::LectParaSpecifiques(entreePrinc,lesCourbes1D,lesFonctionsnD);
   
    // appel au niveau de la classe mère
    Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
                 (*entreePrinc,lesFonctionsnD);
    // mise à jour et même calcul des variables intermédiaires
    // les grandeurs de base ici sans phases pour une première vérification
    alpha1_2 = Sqr(alpha1);
    alpha3_2 = Sqr(alpha3);
//    mu1_2 = Sqr(mu1);
    };
// affichage de la loi
void IsoHyper3DOrgeas1::Affiche() const 
  { cout << " \n loi de comportement 3D hyperelastique isotrope Orgeas 1 : " << Nom_comp(id_comp)
         << " parametres : \n"; 
    cout << " K= " << K ;
    if (cas_Q0s_thermo_dependant != 0)
      { cout << " Q0s thermo dependant, cas= " << cas_Q0s_thermo_dependant;
        switch (cas_Q0s_thermo_dependant)
          { case 1: cout << " T0r= " << T0r <<  " Gr= " << Gr << " Qrmax= " << Qrmax << " ar= " << ar << " ";
                    break;
            case 2: { cout << " courbe Q0s=f(T): " << Q0s_temperature->NomCourbe() <<" ";
                      if ( Q0s_temperature->NomCourbe() == "_") Q0s_temperature->Affiche();
                     }
                    break;
            default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
                          << " cas non pris en compte !!!! dans Affiche() " 
                          << "\n IsoHyper3DOrgeas1::Affiche()";
                     Sortie(1);            
        		
        	};
    	}
    else
     {cout << " ; Q0s = " << Q0s ;};
    // puis les autres para
    cout << " ; mu01 = " << mu01 <<  " ; mu02 = " << mu02 <<  " ; mu03 = " << mu03 
         << " ; alpha1 = " << alpha1 << " ; alpha3 = " << alpha3 ;
    // Qe thermo dépendant ou pas
    if (cas_Q0e_thermo_dependant != 0)
      { cout << " Q0e thermo dependant, cas= " << cas_Q0e_thermo_dependant;
        switch (cas_Q0e_thermo_dependant)
          { case 1: cout << " T0rQe= " << T0rQe << " Qe_T0rQe= " << Qe_T0rQe << " h1= " << h1 << " ";
                    break;
            case 2: { cout << " courbe Q0e=f(T): " << Q0e_temperature->NomCourbe() <<" ";
                      if ( Q0e_temperature->NomCourbe() == "_") Q0e_temperature->Affiche();
                     }
                    break;
            default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
                          << " cas non pris en compte !!!! dans Affiche() " 
                          << "\n IsoHyper3DOrgeas1::Affiche()";
                     Sortie(1);            
        		
        	};
    	}
    else
     {cout << " ; Q0e = " << Q0e ;};
    //  les dépendances à la phase
    if (avec_phase)
        { cout << "\n nQs= "  << nQs  << " gammaQs= "  << gammaQs 
               << "\n nQe= "  << nQe  << " gammaQe= "  << gammaQe 
               << "\n nMu1= " << nMu1 << " gammaMu1= " << gammaMu1 
               << "\n nMu2= " << nMu2 << " gammaMu2= " << gammaMu2 
               << "\n nMu3= " << nMu3 << " gammaMu3= " << gammaMu3; 
         };
    cout << endl;
    // appel de la classe mère
    Loi_comp_abstraite::Affiche_don_classe_abstraite();
  };
            
// affichage et definition interactive des commandes particulières à chaques lois
void IsoHyper3DOrgeas1::Info_commande_LoisDeComp(UtilLecture& entreePrinc)
 {  ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
	   cout << "\n definition standart (rep o) ou exemples exhaustifs (rep n'importe quoi) ? ";
	   string rep = "_";
    // procédure de lecture avec prise en charge d'un retour chariot
    rep = lect_return_defaut(true,"o");

    sort << "\n# .......  loi de comportement 3D hyperelastique isotrope Orgeas 1  ........"
         << "\n#-----------------------------------------------------------------------------------------"
         << "\n#     K      |     Q0s     |   mu01     |   mu02  |  mu03    | alpha1 |  alpha2 |    Q0e   |"
         << "\n#-----------------------------------------------------------------------------------------"
         << "\n     160000        150       17000       220     17000     0.002      0.004     0.01   ";
    if ((rep != "o") && (rep != "O" ) && (rep != "0") )
      { sort << "\n# il est egalement possible de definir une dependance a la phase des parametres"
         << "\n#  dans ce cas, a la suite du parametre Q0e on met le mot cle: avec_phase "
         << "\n#  puis sur la ligne qui suit on met les 10 coefficients de controle de la dependance "
         << "\n# suit un exemple de declaration (il est possible d'ommettre un groupe de 2 "
         << "\n#  mais il faut respecter l'ordre ) "
         << "\n#-----------------------------------------------------------------------"
         << "\n#     K    |  Q0s   |   mu01  | mu02  |  mu03  | alpha1| alpha2 | Q0e   |"
         << "\n#-----------------------------------------------------------------------"
         << "\n#  160000     150      17000    220    17000    0.002   0.004    0.01  avec_phase "
         << "\n # nQs= 0.2 gammaQs= 0.5  nQe= 0.2 gammaQe= 0.5  nMu1= 0.2 gammaMu1= 0.5 nMu2= 0.2 gammaMu2= 0.5  nMu3= 0.2 gammaMu3= 0.5  "
         << "\n # "
         << "\n # noter qu'il ne faut pas que mu3 devienne nulle"
         << "\n # "
         << "\n # il est egalement possible de faire dependre certains coefficients a la temperature"
         << "\n # pour l'instant seule les parametres Q0s et Q0e peuvent-etre thermo-dependant, un exemple:"
         << "\n#-----------------------------------------------------------------------"
         << "\n#     K    |  Q0s   |   mu01  | mu02  |  mu03  | alpha1| alpha2 | Q0e   |"
         << "\n#-----------------------------------------------------------------------"
         << "\n#  160000     Q0s_thermo_dependant_   17000    220    17000    0.002   0.004    Q0e_thermo_dependant_  avec_phase "
         << "\n # nQs= 1.1 gammaQs= 1.  nQe= 1.1 gammaQe= 1.  nMu2= 1.1 gammaMu2= 1.  nMu3= 1.1 gammaMu3= 1.  "
         << "\n #  cas_Q0s_thermo_dependant_  1  T0r= 308  Gr= 7.5   Qrmax=  570    ar= 9 "
         << "\n #  cas_Q0e_thermo_dependant_  1  T0rQe= 308 Qe_T0rQe= 0.01 h1= 1.2"
         << "\n # "
         << "\n # la valeur de Q0s est remplassee par le mot cle: Q0s_thermo_dependant_,"
         << "\n # idem (eventuellement) pour Q0e qui est remplace par Q0e_thermo_dependant_,  ensuite apres la definition"
         << "\n # de tous les autres parametres, sur la ligne qui suit, on definit la thermodependance:"
         << "\n # 1) tout d'abord pour Q0s"
         << "\n # cas_Q0s_thermo_dependant_ : indique le cas de la thermodependance, "
         << "\n #  = 1: indique que l'on a une fonction pre-programmee (voir doc) donnee par Laurent orgeas, "
         << "\n # dependant des parametre qui suivent "
         << "\n #  = 2: indique Q0s depend d'une fonction Q0s(T) que l'on indique ensuite ex: "
         << "\n #  cas_Q0s_thermo_dependant_  2  courbe_evol_Q0s "
         << "\n #    ici courbe_evol_Q0s est le nom d'une courbe qui doit avoir ete definit auparavant "
         << "\n # "
         << "\n # 2) puis eventuellement pour Q0e, dans ce cas on passe une ligne puis: "
         << "\n # cas_Q0e_thermo_dependant_ : indique le cas de la thermodependance, "
         << "\n #  = 1: indique que l'on a la fonction pre-programmee suivante: " 
			      << "\n # la thermo-dépendance est fonction de la pente mu3+mu2 -> Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2)) "
         << "\n # dependant des parametres  T0rQe  Qe_T0rQe et h1 "
         << "\n # ex:  cas_Q0e_thermo_dependant_  1  T0rQe= 273  Qe_T0rQe= 0.07  h1= 1.2 "
         << "\n # second cas: = 2: indique Q0e depend d'une fonction Q0e(T) que l'on indique ensuite ex: "
         << "\n #  cas_Q0e_thermo_dependant_  2  courbe_evol_Q0e "
         << "\n #    ici courbe_evol_Q0e est le nom d'une courbe qui doit avoir ete definit auparavant "
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  il est possible d'indiquer 2 parametres specifiques de limite inf "
         << "\n#  limite_inf_Qeps_ : si Qeps < limite_inf_Qeps on considere que l'angle de phase et ses derivees sont nulles "
         << "\n#  limite_inf_bIIb_ : si Dabs(inv.bIIb) < limite_inf_bIIb   "
         << "\n#    les derivees du potentiel sont calculees pas difference fini, sauf ceux relatif a la phase qui sont mises a 0 "
         << "\n#  les  mots cle limite_inf_Qeps_ suivi du facteur voulu et limite_inf_bIIb_ suivi du facteur "
         << "\n#  doivent etre dans cet ordre "
         << "\n#  ex: "
         << "\n#    limite_inf_Qeps_  8.5e-5  limite_inf_bIIb_ 36.e-10 "
         << "\n#  ces mots cle doivent  se situer avant le mot cle avec_regularisation_ "
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  il est possible d'indiquer un facteur de regularisation qui permet d'eviter "
         << "\n#  de potentiels  problemes de NaN, de type division par 0 par exemple "
         << "\n#  1/a est remplace par 1/(a+fact_regularisation), par defaut fact_regularisation = 1.e-12 "
         << "\n#  pour indiquer un facteur de regulation non nul on indique en dernier parametre "
         << "\n#  le mot cle avec_regularisation_ suivi du facteur voulu "
         << "\n#  ex: "
         << "\n#    avec_regularisation_  1.e-12 "
         << "\n#  ce mot cle doit se situer avant le mot cle sortie_post_ "
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  il est possible de recuperer differentes grandeurs de travail par exemple "
         << "\n#  l'intensite du potentiel, comme ces grandeurs sont calculees au moment de la resolution "
         << "\n#  et ne sont pas stockees, il faut indiquer le mot cle sortie_post_ suivi de 1 (par defaut = 0) "
         << "\n#  ensuite au moment de la constitution du .CVisu on aura acces aux grandeurs de travail "
         << "\n#  ex: "
         << "\n#    sortie_post_  1 "
         << "\n#  ce mot cle est le dernier des parametres specifiques de la loi il doit se situe "
         << "\n#  a la fin de la derniere ligne de donnees "
         << "\n#"
         << "\n#------------------------------------------------------------------------------------"
         << "\n # ";
      };
    sort << endl;
    // appel de la classe Hyper3D
    Hyper3D::Info_commande_LoisDeComp_hyper3D(entreePrinc);
    // appel de la classe mère
    Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);     
  };  		  	  

// test si la loi est complete
int IsoHyper3DOrgeas1::TestComplet()
    { int ret = LoiAbstraiteGeneral::TestComplet();
      if ((K == ConstMath::tresgrand) || (Q0s == ConstMath::tresgrand) || 
          (mu01 == ConstMath::tresgrand) || (mu02 == ConstMath::tresgrand) 
          || (mu03 == ConstMath::tresgrand) || (alpha1 == ConstMath::tresgrand)
          || (alpha3 == ConstMath::tresgrand) || (Q0e == ConstMath::tresgrand))
       { cout << " \n Les parametres ne sont pas defini pour  la loi " << Nom_comp(id_comp)
              << '\n';
         Affiche();     
         ret = 0;
       }
      if (cas_Q0s_thermo_dependant != 0)
        { switch (cas_Q0s_thermo_dependant)
            { case 1: if ((T0r == ConstMath::tresgrand) || (Gr == ConstMath::tresgrand) || 
                          (Qrmax == ConstMath::tresgrand) || (ar == ConstMath::tresgrand)) 
                      { cout << " \n un des parametres T0r, Gr, Qrmax ou ar n'est pas defini pour  la loi " 
                             << Nom_comp(id_comp) << '\n';
                         Affiche();ret = 0;
                      }
                      break;
              case 2: if (Q0s_temperature == NULL)
                      { cout << " \n la loi Q0s(T) n'est pas defini pour  la loi " 
                             << Nom_comp(id_comp) << '\n';
                         Affiche();ret = 0;
                      }
                      break;
              default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
                            << " cas non pris en compte !!!! dans Affiche() " 
                            << "\n IsoHyper3DOrgeas1::TestComplet()";
                       ret = 0.; Sortie(1);            
          	};
        };
      if (cas_Q0e_thermo_dependant != 0)
        { switch (cas_Q0e_thermo_dependant)
            { case 1: if ((T0rQe == ConstMath::tresgrand) || (Qe_T0rQe == ConstMath::tresgrand) 
				             || (h1 == ConstMath::tresgrand) ) 
                      { cout << " \n un des parametres T0rQe, Qe_T0rQe ou h1 n'est pas defini pour  la loi " 
                             << Nom_comp(id_comp) << '\n';
                         Affiche();ret = 0;
                      }
                      break;
              case 2: if (Q0e_temperature == NULL)
                      { cout << " \n la loi Q0e(T) n'est pas defini pour  la loi " 
                             << Nom_comp(id_comp) << '\n';
                         Affiche();ret = 0;
                      }
                      break;
              default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
                            << " cas non pris en compte !!!! dans Affiche() " 
                            << "\n IsoHyper3DOrgeas1::TestComplet()";
                       ret = 0.; Sortie(1);            
          	};
        };
       return ret; 
    }; 
  
	   //----- lecture écriture de restart -----
	   // cas donne le niveau de la récupération
       // = 1 : on récupère tout
       // = 2 : on récupère uniquement les données variables (supposées comme telles)
void IsoHyper3DOrgeas1::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
                                             ,LesFonctions_nD& lesFonctionsnD)
{ string toto;
  if (cas == 1)
   { ent >> toto >> K >> toto >> cas_Q0s_thermo_dependant;
     if (cas_Q0s_thermo_dependant != 0)
         { thermo_dependant=true; // on indique que la loi est thermo-dépendant
           switch (cas_Q0s_thermo_dependant)
           { case 1:{ent >> toto >> T0r >> toto >> Gr >> toto >> Qrmax >> toto >> ar ;
                     // on calcul maintenant les variables intermédiaires fixes
                     Tc=T0r+(Qrmax/Gr/2.);
                     Ts=(ar*ar- Sqr(Qrmax/Gr) )/2./(Qrmax/Gr);
                     break;
                    }
             case 2:{ ent >> toto; 
                      if (toto != " Q0s(T) ")
                       { cout << "\n erreur en lecture de la fonction temperature, on attendait Q0s(T) et on a lue " 
                              << toto
                              << "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(...";
                         Sortie(1);
                       };
                       Q0s_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,Q0s_temperature);
                     break;
                    };
             default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
                           << " cas non pris en compte !!!! dans Lecture_base_info_loi(.. " 
                           << "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(..";
                      Sortie(1);            
           
           };
         }
     else
         {ent >> toto >> Q0s ;};
     // puis les autres para sauf Qe
     ent >> toto >> Q0s >> toto >> mu01 >> toto >> mu02 >> toto >> mu03
         >> toto >> alpha1 >> toto >> alpha3 ;
         // Qe dépendant éventuellement de la température  			 
     ent  >> toto >> cas_Q0e_thermo_dependant;
     if (cas_Q0e_thermo_dependant != 0)
       { thermo_dependant=true; // on indique que la loi est thermo-dépendant
         switch (cas_Q0e_thermo_dependant)
          { case 1:{ent >> toto >> T0rQe >> toto >> Qe_T0rQe >> toto >> h1  ;
                      break;
                   }
            case 2:{ ent >> toto; 
                     if (toto != " Q0e(T) ")
                      { cout << "\n erreur en lecture de la fonction temperature, on attendait Q0e(T) et on a lue " 
                               << toto
                               << "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(...";
                          Sortie(1);
                      };
                     Q0e_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,Q0e_temperature);
                     break;
                   };
            default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
                          << " cas non pris en compte !!!! dans Lecture_base_info_loi(.. " 
                          << "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(..";
                     Sortie(1);
          };

         // maintenant on calcule la valeur de Qs_T0rQr
         {Qs_T0rQe = Q0s ; // au cas où il n'y a pas de dépendance de Qs à la température
          switch (cas_Q0s_thermo_dependant)
              {case 1:// cas d'une dépendance à la température selon fct laurent
                { if (T0rQe >= Tc)
                    {Qs_T0rQe=0.5*Gr*( (T0rQe-Tc+Ts)-sqrt(Sqr(T0rQe-Tc+Ts) + ar*ar)) + Qrmax;}
                  else
                    {Qs_T0rQe=0.5*Gr*( (T0rQe-Tc-Ts)+sqrt(Sqr(T0rQe-Tc-Ts) + ar*ar));};
                  break;
                }
               case 2:// cas d'une dépendance à la température selon une fonction d'évolution
                { Qs_T0rQe = Q0s_temperature->Valeur(T0rQe);
                  break;
                }
               default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
                           << " cas non pris en compte !!!! dans Lecture_base_info_loi(.. " 
                           << "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(..";
                      Sortie(1);
              };
         };
       }
     else
       {ent >> toto >> Q0e ;};
     
     // puis les dépendances éventuelles à la phase
     
     ent >> toto >> avec_phase;
     if (avec_phase)
         { ent >> toto >> nQs  >> toto >> gammaQs 
               >> toto >> nQe  >> toto >> gammaQe
               >> toto >> nMu1 >> toto >> gammaMu1
               >> toto >> nMu2 >> toto >> gammaMu2
               >> toto >> nMu3 >> toto >> gammaMu3;
         };
     // calcul des constantes particulières
     alpha1_2 = Sqr(alpha1);
     alpha3_2 = Sqr(alpha3);
 //       mu1_2 = Sqr(mu1);
     // prise en compte éventuelle de paramètres particulier: ex: régularisation
     // géré par HyperD et Hyper3D
     Hyper3D::Lecture_para_specifiques(ent,cas);
     // gestion du post-traitement
     ent >> toto >> sortie_post ;
   };
  // appel class mère
  Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD); 
};
       // cas donne le niveau de sauvegarde
       // = 1 : on sauvegarde tout
       // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void IsoHyper3DOrgeas1::Ecriture_base_info_loi(ofstream& sort,const int cas)
  { if (cas == 1) 
     { sort << "  module_dilatation " << K << " cas_Q0s_thermo_dependant " << cas_Q0s_thermo_dependant;
       if (cas_Q0s_thermo_dependant != 0)
        { sort << " cas_Q0s_thermo_dependant " << cas_Q0s_thermo_dependant;
          switch (cas_Q0s_thermo_dependant)
          { case 1: sort << " T0r= " << T0r <<  " Gr= " << Gr << " Qrmax= " << Qrmax << " ar= " << ar << " ";
                    break;
            case 2: {sort << " Q0s(T) ";
                     LesCourbes1D::Ecriture_pour_base_info(sort,cas,Q0s_temperature);
                    };
                    break;
            default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
                          << " cas non pris en compte !!!! dans Ecriture_base_info_loi(.. " 
                          << "\n IsoHyper3DOrgeas1::Ecriture_base_info_loi(..";
                     Sortie(1);            
        		
        	};
        }
       else
        {sort << " Q0s " << Q0s ;};
       // puis les autres para sauf Qe
       sort << "  mu01 " << mu01 << "  mu02 " << mu02 << "  mu03 " << mu03
            << " alpha1 " << alpha1 << " alpha3 " << alpha3 << " ";
		     // Qe	 eventuellement thermo dépendant
       sort <<  " cas_Q0e_thermo_dependant " << cas_Q0e_thermo_dependant;
       if (cas_Q0e_thermo_dependant != 0)
        { sort << " cas_Q0e_thermo_dependant " << cas_Q0e_thermo_dependant;
          switch (cas_Q0e_thermo_dependant)
          { case 1: sort << " T0rQe= " << T0rQe <<  " Qe_T0rQe= " << Qe_T0rQe  << " h1= " << h1 << " ";
                    break;
            case 2: {sort << " Q0e(T) ";
                     LesCourbes1D::Ecriture_pour_base_info(sort,cas,Q0e_temperature);
                    };
                    break;
            default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
                          << " cas non pris en compte !!!! dans Ecriture_base_info_loi(.. " 
                          << "\n IsoHyper3DOrgeas1::Ecriture_base_info_loi(..";
                     Sortie(1);            
        		
        	 };
        }
       else
        {sort << " Q0e " << Q0e ;};
		  
		     // puis les dépendances éventuelles à la phase
		     sort  << " avec_phase " << avec_phase;
       if (avec_phase)
        { sort << "\n nQs= "  << nQs  << " gammaQs= "  << gammaQs 
               << "\n nQe= "  << nQe  << " gammaQe= "  << gammaQe 
               << "\n nMu1= " << nMu1 << " gammaMu1= " << gammaMu1 
               << "\n nMu2= " << nMu2 << " gammaMu2= " << gammaMu2 
               << "\n nMu3= " << nMu3 << " gammaMu3= " << gammaMu3; 
         };
       // prise en compte éventuelle de paramètres particulier: ex: régularisation
       // géré par HyperD et Hyper3D
       Hyper3D::Ecriture_para_specifiques(sort,cas);
         
       // gestion du post-traitement
       sort << " sortie_post= "<< sortie_post << " ";
     };
    // appel de la classe mère
    Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);  
   }; 
        // calcul d'un module d'young équivalent à la loi
double IsoHyper3DOrgeas1::Module_young_equivalent(Enum_dure ,const Deformation &,SaveResul * )
   { return (9.*K*(mu01+mu02)/((mu01+mu02)+3.*K));
   };
// récupération d'un module de compressibilité équivalent à la loi, ceci pour un chargement nul
// il s'agit ici de la relation -pression = sigma_trace/3. = module de compressibilité * I_eps
double IsoHyper3DOrgeas1::Module_compressibilite_equivalent(Enum_dure ,const Deformation & def)
   { return K/3.;};

// ===========  METHODES Protégées dérivant de virtuelles : ==============

//  METHODES internes spécifiques à l'hyperélasticité isotrope découlant de 
//  méthodes virtuelles de Hyper3D

    // calcul du potentiel tout seul sans la phase car Qeps est  nul
    // ou très proche de 0
double IsoHyper3DOrgeas1::PoGrenoble
             (const double & Qeps,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01;double mu2 = mu02; double mu3 = mu03;
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe = Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = Qeps * Qeps;

   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);

   double QplusQec = (Qeps + Qec);

   double S1 = sqrt(alpha1_2 + Qec2);

   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   double E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(Qeps*(Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4);
   // retour
   return E;
   };

    // calcul du potentiel tout seul avec la phase donc dans le cas où Qeps est non nul
double IsoHyper3DOrgeas1::PoGrenoble
             (const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et  température
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	};
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   // dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
   if (avec_phase)
     {Qs = Q0sT / pow((1.+gammaQs * inv_spec.cos3phi),nQs);
      Qe = Q0eT / pow((1.+gammaQe * inv_spec.cos3phi),nQe);
      mu1 = mu01 / pow((1.+gammaMu1 * inv_spec.cos3phi),nMu1);
      mu2 = mu02 / pow((1.+gammaMu2 * inv_spec.cos3phi),nMu2);
      mu3 = mu03 / pow((1.+gammaMu3 * inv_spec.cos3phi),nMu3);
      };
 
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = inv_spec.Qeps * inv_spec.Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (inv_spec.Qeps + Qec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = inv_spec.Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);

   // le potentiel 
   double E = K/6. * (logV)*(logV) + Qs * inv_spec.Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* inv_spec.Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(inv_spec.Qeps*(inv_spec.Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4);
   // retour
   return E;
   };

    // calcul du potentiel tout seul sans la phase car Qeps est  nul
    // ou très proche de 0, et de sa variation suivant V uniquement
Hyper3D::PoGrenoble_V IsoHyper3DOrgeas1::PoGrenoble_et_V
             (const double & Qeps,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	};
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe =  Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe =  Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   PoGrenoble_V ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = Qeps * Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (Qeps + Qec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(Qeps*(Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4);
   // dérivées premières
   ret.EV = K/3.*logV/inv.V;
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

    // calcul du potentiel et de sa variation suivant V uniquement
Hyper3D::PoGrenoble_V IsoHyper3DOrgeas1::PoGrenoble_et_V
             (const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et  température
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	};
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   // dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
   if (avec_phase)
     {Qs = Q0sT / pow((1.+gammaQs * inv_spec.cos3phi),nQs);
      Qe = Q0eT / pow((1.+gammaQe * inv_spec.cos3phi),nQe);
      mu1 = mu01 / pow((1.+gammaMu1 * inv_spec.cos3phi),nMu1);
      mu2 = mu02 / pow((1.+gammaMu2 * inv_spec.cos3phi),nMu2);
      mu3 = mu03 / pow((1.+gammaMu3 * inv_spec.cos3phi),nMu3);
      };
   PoGrenoble_V ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps = inv_spec.Qeps; // pour simplifier la notation
   double Qeps2 = Qeps * Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (Qeps + Qec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(Qeps*(Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4);
   // dérivées premières
   ret.EV = K/3.*logV/inv.V;
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

    // calcul du potentiel tout seul sans la phase car Qeps est  nul
    // ou très proche de 0, et de ses variations première et seconde suivant V uniquement
Hyper3D::PoGrenoble_VV IsoHyper3DOrgeas1::PoGrenoble_et_VV
             (const double & Qeps,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	};
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe =  Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe =  Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   PoGrenoble_VV ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = Qeps * Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (Qeps + Qec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(Qeps*(Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4);
/* double toto = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2;
 toto =  + mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(QplusQec/alpha1_2 * S2) - log(Qec/alpha1_2 * S1)));
 cout << "toto = " << toto;                      
 toto =  + mu3/2.*(Qeps*(Qeps-2. * S4));
 cout << "toto = " << toto;                      
 toto=  + alpha3_2*(log(QmoinsQe/alpha3_2 * S3) - log(-Qe/alpha3_2 * S4));
 cout << "toto = " << toto;                      
 toto =  + QmoinsQe * S3 + Qe * S4;                           
 cout << "toto = " << toto;  */                    
   // dérivées premières
   ret.EV = K/3.*logV/inv.V;
   // dérivées secondes
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
//++++verif*** 
//Verif_PoGrenoble_et_var_VV(Qeps,inv,ret );
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

// calcul du potentiel et de sa variation première et seconde suivant V uniquement
// on prend en compte une dépendance éventuelle à la phase    
Hyper3D::PoGrenoble_VV IsoHyper3DOrgeas1::PoGrenoble_et_VV
             (const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et  température
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	};
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   // dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
   if (avec_phase)
     {Qs = Q0sT / pow((1.+gammaQs * inv_spec.cos3phi),nQs);
      Qe = Q0eT / pow((1.+gammaQe * inv_spec.cos3phi),nQe);
      mu1 = mu01 / pow((1.+gammaMu1 * inv_spec.cos3phi),nMu1);
      mu2 = mu02 / pow((1.+gammaMu2 * inv_spec.cos3phi),nMu2);
      mu3 = mu03 / pow((1.+gammaMu3 * inv_spec.cos3phi),nMu3);
      };
   PoGrenoble_VV ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps = inv_spec.Qeps; // pour simplifier la notation
   double Qeps2 = Qeps * Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (Qeps + Qec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
	// le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(Qeps*(Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4);
   // dérivées premières
   ret.EV = K/3.*logV/inv.V;
   // dérivées secondes
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

 // calcul du potentiel et de ses dérivées premières non compris la phase
Hyper3D::PoGrenobleSansPhaseSansVar IsoHyper3DOrgeas1::PoGrenoble
            (const InvariantQeps & i_s,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe = Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   PoGrenobleSansPhaseSansVar ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = i_s.Qeps * i_s.Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (i_s.Qeps + Qec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = i_s.Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4); // ok
   // dérivées premières
   ret.EV = K/3.*logV/inv.V;
   double dS2_Q = QplusQec / S2;
   double dS3_Q = QmoinsQe / S3;
   ret.EQ = Qs + 2.*mu2*i_s.Qeps
            + 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2 
                      -alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
            + 0.5*mu3*(2.*(i_s.Qeps - S4)
                      +alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
                      +S3+QmoinsQe*dS3_Q);          

   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };
   
   // calcul du potentiel et de ses dérivées avec la phase
Hyper3D::PoGrenobleAvecPhaseSansVar IsoHyper3DOrgeas1::PoGrenoblePhase
             (const InvariantQepsCosphi& i_s,const Invariant & inv) 
 { // les grandeurs de base,
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et  température
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	};
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   // dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
   if (avec_phase)
     {Qs = Q0sT / pow((1.+gammaQs * i_s.cos3phi),nQs);
      Qe = Q0eT / pow((1.+gammaQe * i_s.cos3phi),nQe);
      mu1 = mu01 / pow((1.+gammaMu1 * i_s.cos3phi),nMu1);
      mu2 = mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2);
      mu3 = mu03 / pow((1.+gammaMu3 * i_s.cos3phi),nMu3);
      };
	  double mu1_2 = Sqr(mu1);
   PoGrenobleAvecPhaseSansVar ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = i_s.Qeps * i_s.Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (i_s.Qeps + Qec);
	  double QplusQec_2 = Sqr(QplusQec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = i_s.Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4); // ok
   // dérivées premières
   ret.EV = K/3.*logV/inv.V;
   double dS2_Q = QplusQec / S2;
   double dS3_Q = QmoinsQe / S3;
   ret.EQ = Qs + 2.*mu2*i_s.Qeps
            + 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2 
                      -alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
            + 0.5*mu3*(2.*(i_s.Qeps - S4)
                      +alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
                      +S3+QmoinsQe*dS3_Q);          
//            + mu1 *(QplusQec - S2) + mu3 * (i_s.Qeps - S4 + S3);
   // cas de la variation faisant intervenir  la phase
   // 1-- calcul des variations des grandeurs de bases   
   double dQs_dcos3phi = - nQs * gammaQs * Qs / (1.+gammaQs * i_s.cos3phi);
   double dQe_dcos3phi = - nQe * gammaQe * Qe / (1.+gammaQe * i_s.cos3phi);
   double dMu1_dcos3phi = - nMu1 * gammaMu1 * mu1 / (1.+gammaMu1 * i_s.cos3phi);
   double dMu2_dcos3phi = - nMu2 * gammaMu2 * mu2 / (1.+gammaMu2 * i_s.cos3phi);
   double dMu3_dcos3phi = - nMu3 * gammaMu3 * mu3 / (1.+gammaMu3 * i_s.cos3phi);
   // 2-- variations des grandeurs intermédiaires par rapport aux grandeurs de bases
   double Qs2 = Qs * Qs;
   double dQec_Qs = - (Qs2 + mu1_2 * alpha1_2) / (2. * mu1 * Qs2); // OK
   double dS2_Qs = QplusQec * dQec_Qs / S2; // OK
   double dS1_Qs = Qec * dQec_Qs / S1; // OK
   double dS4_Qe = Qe / S4;  // ok
   double dS3_Qe = -(QmoinsQe) / S3;  // OK
	
   double dQec_Mu1 = (Qs2 + mu1_2 * alpha1_2) / (2. * mu1_2 * Qs);
   double dQec2_Mu1 = 2. * Qec * dQec_Mu1;
   double dQplusQec_Mu1 = dQec_Mu1;  // pourra être supprimé par la suite mais pour la mise au point ... c'est plus cohérent
   double dS1_Mu1 = dQec_Mu1 * Qec/ S1; 
   double dS2_Mu1 = dQec_Mu1 * QplusQec/ S2;
  
   // 3-- variations du potentiel par rapport aux grandeurs de bases
   // A_x,B_x,C_x sont des variables intermédiaires
	  double A_x = (dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
                           -alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
                           );
   double dE_Qs = i_s.Qeps + 0.5*mu1*A_x;
//                 + 0.5*mu1*(dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
//                           -alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
//                           ); // OK
   double dE_Qe = 0.5*mu3*(-2.*i_s.Qeps*dS4_Qe
                           +alpha3_2*((-1.+dS3_Qe)/(QmoinsQe + S3)-(-1.+dS4_Qe)/(-Qe+S4))
                           - S3 + S4 + Qe*dS4_Qe + QmoinsQe * dS3_Qe); // OK

	  double B_x = 2.* i_s.Qeps * dQec_Mu1 - (dQplusQec_Mu1 * S2 + QplusQec * dS2_Mu1) + dQec_Mu1* S1 + Qec* dS1_Mu1
						         - alpha1_2*( ((dQplusQec_Mu1+dS2_Mu1 )) / (QplusQec+S2)
									            - ( (dQec_Mu1+dS1_Mu1)) / (Qec+S1)      
													);
	  double C_x = Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
				  - alpha1_2*(log(Dabs(QplusQec + S2)) - log(Dabs(Qec +  S1)));				   											
   double dE_Mu1 = 0.5 * ( C_x + mu1 * B_x ); // ok

   double dE_Mu2 = i_s.Qeps * i_s.Qeps;  // OK
   double dE_Mu3 = 0.5 * (i_s.Qeps*(i_s.Qeps-2. * S4)
                          + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
                          + QmoinsQe * S3 + Qe * S4);    // OK                     
   // 4-- variations du potentiel par rapport à cos3phi
   ret.Ecos3phi = dE_Qs * dQs_dcos3phi + dE_Qe * dQe_dcos3phi
                 + dE_Mu1 * dMu1_dcos3phi + dE_Mu2 * dMu2_dcos3phi + dE_Mu3 * dMu3_dcos3phi;   // OK                     
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };
      

// calcul  du potentiel sans phase et dérivées avec  ses variations par rapport aux invariants
Hyper3D::PoGrenobleSansPhaseAvecVar IsoHyper3DOrgeas1::PoGrenoble_et_var
             (const Invariant2Qeps& i_s,const Invariant & inv)
 { // les grandeurs de base, qui ici ne dépendent pas de la phase 
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	  };
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe =  Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe =  Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	  };
   PoGrenobleSansPhaseAvecVar ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = i_s.Qeps * i_s.Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (i_s.Qeps + Qec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
   double QmoinsQe = i_s.Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4);
   // dérivées premières
   ret.EV = K/3.*logV/inv.V;
   double dS2_Q = QplusQec / S2;
   double dS3_Q = QmoinsQe / S3;
   ret.EQ = Qs + 2.*mu2*i_s.Qeps
            + 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2 
                      -alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
            + 0.5*mu3*(2.*(i_s.Qeps - S4)
                      +alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
                      +S3+QmoinsQe*dS3_Q);          
//   ret.EQ = Qs + 2.*mu2*i_s.Qeps
//            + mu1 *(QplusQec - S2) + mu3 * (i_s.Qeps - S4 + S3);
   // dérivées secondes
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
   ret.EQV = 0.;
//   double dS2_QQ = (S2-dS2_Q)/(S2*S2);
//   double dS3_QQ = (S3-dS3_Q)/(S3*S3);
   double dS2_QQ = (1. - dS2_Q*dS2_Q)/S2;   // OK
   double dS3_QQ = (1. - dS3_Q*dS3_Q)/S3;   // OK
   ret.EQQ = 2. * mu2 
           + 0.5*mu1*( 2.*(1-dS2_Q) - QplusQec * dS2_QQ 
                      -alpha1_2*(dS2_QQ/(QplusQec + S2)-Sqr((1.+dS2_Q)/(QplusQec + S2))))
           + 0.5*mu3*(2. 
                      +alpha3_2*(dS3_QQ/(QmoinsQe + S3)-Sqr((1.+dS3_Q)/(QmoinsQe + S3)))
                      +2.* dS3_Q + QmoinsQe * dS3_QQ);           
           
//   ret.EQQ = 2. * mu2 + mu1*(1.-QplusQec*S2) + mu3*(1. + QmoinsQe*S3);
    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
//++++verif*** 
//   Verif_PoGrenoble_et_var(i_s.Qeps,inv,ret );
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };
    
    // calcul  du potentiel avec phase et dérivées avec  ses variations par rapport aux invariants
Hyper3D::PoGrenobleAvecPhaseAvecVar IsoHyper3DOrgeas1::PoGrenoblePhase_et_var
             (const Invariant2QepsCosphi& i_s,const Invariant & inv)
 { // les grandeurs de base,
   double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
   double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et  température
   switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	};
   switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	 };
   // dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
   if (avec_phase)
     {Qs = Q0sT / pow((1.+gammaQs * i_s.cos3phi),nQs);
      Qe = Q0eT / pow((1.+gammaQe * i_s.cos3phi),nQe);
      mu1 = mu01 / pow((1.+gammaMu1 * i_s.cos3phi),nMu1);
      mu2 = mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2);
      mu3 = mu03 / pow((1.+gammaMu3 * i_s.cos3phi),nMu3);
//cout << "\n cos(3phi)= " <<  i_s.cos3phi << " Qs= " << Qs << " mu2= " << mu2;     
      };
   double mu1_2 = Sqr(mu1);
	  double Qs_2 = Sqr(Qs);
   PoGrenobleAvecPhaseAvecVar ret;
   // des variables intermédiaires
   double logV = log(inv.V);
   double Qeps2 = i_s.Qeps * i_s.Qeps;
   double Qe2 = Qe * Qe;
   double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
   double Qec2 = Sqr(Qec);
   double QplusQec = (i_s.Qeps + Qec);
	  double QplusQec_2 = Sqr(QplusQec);
   double S1 = sqrt(alpha1_2 + Qec2);
   double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
	  double S2_2 = Sqr(S2);
   double QmoinsQe = i_s.Qeps - Qe;
   double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
   double S4 = sqrt(alpha3_2+Qe2);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
           + mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
                       - alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec +  S1))))
           + mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
           + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
           + QmoinsQe * S3 + Qe * S4); // ok
           
   // +++++++++++++++++++  dérivées premières  +++++++++++++++++++
   // a) la partie sans phase 
     
   ret.EV = K/3.*logV/inv.V; // ok
   
// --- pour le débug   
//   switch(fpclassify(ret.EV))
//    { case FP_NAN: case FP_INFINITE:
//        { // cas d'un nombre infini 
//          cout << "\n erreur valeur nan pour EV! ";
//        };
//        break;
//    };
    
// --- fin débug
    
   double dS2_Q = QplusQec / S2; // OK
   double dS3_Q = QmoinsQe / S3; // OK
   // D_x est une variable intermédiaire, utilisée par la suite
	  double D_x = QplusQec*(2.-dS2_Q) -S2 - alpha1_2*((1.+dS2_Q)/(QplusQec + S2));
   ret.EQ = Qs + 2.*mu2*i_s.Qeps
            + 0.5 * mu1 * D_x 
//            + 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2 
//                      -alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
            + 0.5*mu3*(2.*(i_s.Qeps - S4)
                      +alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
                      +S3+QmoinsQe*dS3_Q);  // OK        
							 
   // b) cas de la variation faisant intervenir  la phase
   
   // 1-- calcul des variations des grandeurs de bases   
   double dQs_dcos3phi = - nQs * gammaQs * Qs / (1.+gammaQs * i_s.cos3phi); // 0K
   double dQe_dcos3phi = - nQe * gammaQe * Qe / (1.+gammaQe * i_s.cos3phi); // OK
   double dMu1_dcos3phi = - nMu1 * gammaMu1 * mu1 / (1.+gammaMu1 * i_s.cos3phi); // OK
   double dMu2_dcos3phi = - nMu2 * gammaMu2 * mu2 / (1.+gammaMu2 * i_s.cos3phi); // OK
   
// double toto  =  - nMu2 * gammaMu2 * mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2+1);  
//   if (Abs(toto -dMu2_dcos3phi) > ConstMath::unpeupetit)
//      cout << "\n pb dMu2_dcos3phi" << toto << " " << dMu2_dcos3phi << endl;
   
   
   double dMu3_dcos3phi = - nMu3 * gammaMu3 * mu3 / (1.+gammaMu3 * i_s.cos3phi); // OK
   // 2-- variations des grandeurs intermédiaires par rapport aux grandeurs de bases
   double Qs2 = Qs * Qs;
   double dQec_Qs = - (Qs2 + mu1_2 * alpha1_2) / (2. * mu1 * Qs2); // OK
   double dS2_Qs = QplusQec * dQec_Qs / S2; // OK
   double dS1_Qs = Qec * dQec_Qs / S1; // OK
   double dS4_Qe = Qe / S4;  // ok
   double dS3_Qe = -(QmoinsQe) / S3;  // OK
	
   double dQec_Mu1 = (Qs2 + mu1_2 * alpha1_2) / (2. * mu1_2 * Qs);
   double dQec2_Mu1 = 2. * Qec * dQec_Mu1;
   double dQplusQec_Mu1 = dQec_Mu1;  // pourra être supprimé par la suite mais pour la mise au point ... c'est plus cohérent
   double dS1_Mu1 = dQec_Mu1 * Qec/ S1; 
   double dS2_Mu1 = dQec_Mu1 * QplusQec/ S2;
  
   // 3-- variations du potentiel par rapport aux grandeurs de bases
   // A_x,B_x,C_x sont des variables intermédiaires
	  double A_x = (dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
                           -alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
                           );
   double dE_Qs = i_s.Qeps + 0.5*mu1*A_x;
//                 + 0.5*mu1*(dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
//                           -alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
//                           ); // OK
   double dE_Qe = 0.5*mu3*(-2.*i_s.Qeps*dS4_Qe
                           +alpha3_2*((-1.+dS3_Qe)/(QmoinsQe + S3)-(-1.+dS4_Qe)/(-Qe+S4))
                           - S3 + S4 + Qe*dS4_Qe + QmoinsQe * dS3_Qe); // OK

	  double B_x = 2.* i_s.Qeps * dQec_Mu1 - (dQplusQec_Mu1 * S2 + QplusQec * dS2_Mu1) + dQec_Mu1* S1 + Qec* dS1_Mu1
						         - alpha1_2*( ((dQplusQec_Mu1+dS2_Mu1 )) / (QplusQec+S2)
									            - ( (dQec_Mu1+dS1_Mu1)) / (Qec+S1)      
													);
	  double C_x = Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
				  - alpha1_2*(log(Dabs(QplusQec + S2)) - log(Dabs(Qec +  S1)));				   											
   double dE_Mu1 = 0.5 * ( C_x + mu1 * B_x ); // ok

   double dE_Mu2 = i_s.Qeps * i_s.Qeps;  // OK
   double dE_Mu3 = 0.5 * (i_s.Qeps*(i_s.Qeps-2. * S4)
                          + alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
                          + QmoinsQe * S3 + Qe * S4);    // OK                     
   // 4-- variations du potentiel par rapport à cos3phi
   ret.Ecos3phi = dE_Qs * dQs_dcos3phi + dE_Qe * dQe_dcos3phi
                 + dE_Mu1 * dMu1_dcos3phi + dE_Mu2 * dMu2_dcos3phi + dE_Mu3 * dMu3_dcos3phi;   // OK                     
// --- pour le débug   
//   switch(fpclassify(ret.Ecos3phi))
//    { case FP_NAN: case FP_INFINITE:
//        { // cas d'un nombre infini 
//          cout << "\n erreur valeur nan pour EV! ";
//        };
//        break;
//    };
    
// --- fin débug

   // +++++++++++++++++++   dérivées secondes +++++++++++++++++++  
   // a) la partie sans phase 
   
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V); // ok
   ret.EQV = 0.;
   double dS2_QQ = (1. - dS2_Q*dS2_Q)/S2;   // OK
   double dS3_QQ = (1. - dS3_Q*dS3_Q)/S3;   // OK
   ret.EQQ = 2. * mu2 
           + 0.5*mu1*( 2.*(1-dS2_Q) - QplusQec * dS2_QQ
                      -alpha1_2*(dS2_QQ/(QplusQec + S2)-Sqr((1.+dS2_Q)/(QplusQec + S2))))
           + 0.5*mu3*(2. 
                      +alpha3_2*(dS3_QQ/(QmoinsQe + S3)-Sqr((1.+dS3_Q)/(QmoinsQe + S3)))
                      +2.* dS3_Q + QmoinsQe * dS3_QQ);  // OK 
							 
//   ret.EQQ = 2. * mu2 + mu1*(1.-QplusQec*S2) + mu3*(1. + QmoinsQe*S3);

    // b) la partie dépendant de la phase
   // 1-- calcul des variations seconde des grandeurs de bases   
   double dQs_dcos3phi2 = - (nQs+1) * gammaQs * dQs_dcos3phi / (1.+gammaQs * i_s.cos3phi); // OK
   double dQe_dcos3phi2 = - (nQe+1) * gammaQe * dQe_dcos3phi / (1.+gammaQe * i_s.cos3phi); // OK
   double dMu1_dcos3phi2 = - (nMu1+1) * gammaMu1 * dMu1_dcos3phi / (1.+gammaMu1 * i_s.cos3phi); // OK
   double dMu2_dcos3phi2 = - (nMu2+1) * gammaMu2 * dMu2_dcos3phi / (1.+gammaMu2 * i_s.cos3phi); // OK
   
//   toto =  nMu2 *(nMu2+1)* gammaMu2*gammaMu2 * mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2+2);  
   
//   if (Abs(toto -dMu2_dcos3phi2) > ConstMath::unpeupetit)
//      cout << "\n pb dMu2_dcos3phi2" << toto << " " << dMu2_dcos3phi2 << endl;
   
   double dMu3_dcos3phi2 = - (nMu3+1) * gammaMu3 * dMu3_dcos3phi / (1.+gammaMu3 * i_s.cos3phi); // OK
   // 2-- variations seconde des grandeurs intermédiaires par rapport aux grandeurs de bases
   double dQec_QsQs = mu1 * alpha1_2 / (Qs2 * Qs); // OK
   double dS2_QsQs = ((dQec_Qs*dQec_Qs + QplusQec * dQec_QsQs) - Sqr(QplusQec*dQec_Qs/S2))/S2; // OK
   double dS1_QsQs = ((dQec_Qs*dQec_Qs + Qec * dQec_QsQs) - Sqr(Qec*dQec_Qs/S1))/S1; // OK
   double dS4_QeQe = (1.-Sqr(Qe/S4))/S4; // OK
   double dS3_QeQe = (1.-Sqr(QmoinsQe/S3))/S3; // OK
   double dS2_QQs = (dQec_Qs - (QplusQec) * dS2_Qs / S2) / S2; // OK
   double dS3_QQe = (-1. - QmoinsQe * dS3_Qe / S3) / S3 ; // OK
	
   double dQec_Mu1Mu1 = -Qs / (mu1 * mu1_2); // ok
   double dQec_QsMu1 = -alpha1_2/Qs_2 + (Qs_2+mu1_2*alpha1_2) / (2.*mu1_2*Qs_2); // OK
   double dQec_Mu1Q = 0.; 
   double dS1_Mu1Mu1 = ((dQec_Mu1*dQec_Mu1 + Qec*dQec_Mu1Mu1) -Sqr(Qec*dQec_Mu1)/Sqr(S1)) / S1; // ok
   double dS1_QsMu1 = ((dQec_Mu1*dQec_Qs+Qec*dQec_QsMu1)-dQec_Mu1*dQec_Qs*Sqr(Qec/S1)) / S1; // ok
   double dS2_Mu1Mu1 = dQec_Mu1 * dQec_Mu1 * (1.-QplusQec_2/S2_2)/S2 + QplusQec/S2 * dQec_Mu1Mu1; // ok
   double dS2_QsMu1 = dQec_Mu1 * dQec_Qs * (1.-QplusQec_2/S2_2)/S2 + QplusQec/S2 * dQec_QsMu1; // ok
   double dS2_QMu1 = dQec_Mu1 * (1.-QplusQec_2/S2_2)/S2;	// ok
	
   // 3-- variations seconde du potentiel par rapport aux grandeurs de bases
   double dE_QsQs = 0.5*mu1*(dQec_QsQs*(2.*i_s.Qeps-S2+S1) - 2.*dQec_Qs*dS2_Qs -QplusQec*dS2_QsQs
                             +2.*dQec_Qs*dS1_Qs + Qec * dS1_QsQs
                             -alpha1_2*((dQec_QsQs+dS2_QsQs)/(QplusQec+S2)-Sqr((dQec_Qs+dS2_Qs)/(QplusQec+S2)) 
                                        -(dQec_QsQs+dS1_QsQs)/(Qec+S1)+Sqr((dQec_Qs+dS1_Qs)/(Qec+S1))
                                       )
                             );  // OK         
   double dE_QeQe = 0.5*mu3*(-2.*i_s.Qeps*dS4_QeQe
                             +alpha3_2*((dS3_QeQe)/(QmoinsQe + S3)-Sqr((-1.+dS3_Qe)/(QmoinsQe + S3))
                                       -(dS4_QeQe)/(-Qe + S4)+Sqr((-1.+dS4_Qe)/(-Qe + S4))
                                      )
                             -2.*dS3_Qe + 2.*dS4_Qe + Qe*dS4_QeQe + QmoinsQe*dS3_QeQe
                             );  // OK        
   double dE_Mu3Qe = dE_Qe/mu3; // OK
   // puis pour la dérivée croisée
   double dE_QQs = 1. + 0.5*mu1*(dQec_Qs*(2.-dS2_Q)-dS2_Qs-QplusQec*dS2_QQs
                               -alpha1_2*( (dS2_QQs-(1.+dS2_Q)*(dQec_Qs+dS2_Qs)/(QplusQec+S2))/(QplusQec+S2))
                               ); // OK
   double dE_QQe = 0.5*mu3*(-2.*dS4_Qe
                            +alpha3_2*((dS3_QQe-(1.+dS3_Q)*(-1.+dS3_Qe)/(QmoinsQe+S3))/(QmoinsQe+S3))
                            +dS3_Qe - dS3_Q + QmoinsQe*dS3_QQe
                            );   // OK                        
   double dE_QMu2 = 2.* i_s.Qeps;  // OK
   double dE_QMu3 = 0.5 * (2.*i_s.Qeps - 2. * S4 
                           +alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
                           +S3+QmoinsQe*dS3_Q);   // OK 
									
   double dE_QsMu1 = 0.5 * A_x 
	                 + 0.5*mu1 * (dQec_QsMu1*(2.*i_s.Qeps-S2+S1)+dQec_Qs*(-dS2_Mu1+dS1_Mu1)-dQec_Mu1*dS2_Qs-QplusQec*dS2_QsMu1
						           +dQec_Mu1*dS1_Qs+Qec*dS1_QsMu1
								-alpha1_2*( (dQec_QsMu1+dS2_QsMu1)/(QplusQec + S2)-(dQec_Qs+dS2_Qs)*(dQec_Mu1+dS2_Mu1)/Sqr(QplusQec + S2)
								       -	(dQec_QsMu1+dS1_QsMu1)/(Qec + S1) + (dQec_Qs+dS1_Qs)*(dQec_Mu1+dS1_Mu1)/Sqr(Qec+S1))
										 ); // ok
   double dE_QsQe = 0.;
   double dE_Mu1Mu1 =  B_x  // en fait ici il y a la somme de deux termes 0.5 * B 
						  + 0.5 * mu1 * ( dQec_Mu1Mu1 * (2.* i_s.Qeps - S2 + S1) - 2.*dQec_Mu1*(dS2_Mu1 - dS1_Mu1)
												- QplusQec * dS2_Mu1Mu1 + Qec * dS1_Mu1Mu1
												- alpha1_2*( (dQec_Mu1Mu1+dS2_Mu1Mu1 ) / (QplusQec+S2) - Sqr( (dQec_Mu1+dS2_Mu1) / (QplusQec+S2))      
													 - (dQec_Mu1Mu1+dS1_Mu1Mu1)/(Qec + S1) + Sqr((dQec_Mu1+dS1_Mu1)/(Qec+S1)))
											 ); // ok
   double dE_QMu1 =  0.5 * D_x	  
	                + 0.5 * mu1	 * (dQec_Mu1 * (2.-dS2_Q) - dS2_Mu1 - QplusQec * dS2_QMu1	
						               -alpha1_2*( (dS2_QMu1)/(QplusQec+S2)
									                  -(dQec_Mu1+dS2_Mu1)*(1.+dS2_Q)/Sqr(QplusQec+S2)
														  )
											 ); // ok																	
   
   // 4-- variations seconde du potentiel par rapport à cos3phi
	// on tiens compte que dE_QsMu2=0, dE_QsMu3=0, dE_QeMu1=, dE_QeMu2=0, dE_Mu1Mu2=0, dE_Mu1Mu3=0, 
	// et également que dE_Mu2Mu2=0, dE_Mu3Mu3=0,
   ret.Ecos3phi2 = dE_QsQs * Sqr(dQs_dcos3phi) + dE_Qs * dQs_dcos3phi2
	                + dE_QeQe * Sqr(dQe_dcos3phi) + dE_Qe * dQe_dcos3phi2
                   + dE_Mu1Mu1 * Sqr(dMu1_dcos3phi) + dE_Mu1 * dMu1_dcos3phi2	
//                   + dE_Mu2Mu2 * Sqr(dMu2_dcos3phi) + dE_Mu2 * dMu2_dcos3phi2	
                                                    + dE_Mu2 * dMu2_dcos3phi2	
//                   + dE_Mu3Mu3 * Sqr(dMu3_dcos3phi) + dE_Mu3 * dMu3_dcos3phi2	
						                                  + dE_Mu3 * dMu3_dcos3phi2	
                   + 2. * dE_Mu3Qe * dQe_dcos3phi * dMu3_dcos3phi
                   + 2. * dE_QsQe * dQs_dcos3phi * dQe_dcos3phi
                   + 2. * dE_QsMu1 * dQs_dcos3phi * dMu1_dcos3phi; // ok
   // 5-- variations coisée seconde du potentiel par rapport à cos3phi et Qepsilon
   ret.EVcos3phi= 0.; 
   ret.EQcos3phi= dE_QQs*dQs_dcos3phi + dE_QQe*dQe_dcos3phi + dE_QMu1*dMu1_dcos3phi+ dE_QMu2*dMu2_dcos3phi + dE_QMu3*dMu3_dcos3phi; // ok
   
    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
//++++verif*** ++++++
//  la vérification a été effectuée selon trois étapes:
// 1) développement des formules avec maxima, et comparaison des formules
// 2) vérification par rapport aux dérivées numériques (pas satisfaisant: ne permet pas de voir où est l'erreur)
// 3) calcul de la dérivée totale du potentiel par rapport aux invariants V,Q,cos(3phi) (en calcul formel), et évaluation numérique
//    du résultat pour différentes valeurs de (V,Q,cos(3phi)) et comparaison avec le calcul donné par herezh
//    --> extrémement puissant, à permis de trouver toutes les erreurs, et on peut être a peut-près certain
//    qu'actuellement il n'y a plus d'erreur (jamais complètement certain ...)
//    la méthode CompaFormel() est utilisée pour cela

//   Verif_PoGrenoble_et_var(i_s.Qeps,inv,ret );
//     Verif_PoGrenoble_et_var(i_s.Qeps,inv,i_s.cos3phi,ret );
// ---   Invariant2Qeps toto(inv_spec.Qeps,inv_spec.dQepsdbIIb,inv_spec.dQepsdbIIb2);
// ---   Hyper3D::PoGrenobleSansPhaseAvecVar truc = PoGrenoble_et_var(toto,inv);
// utilitaire de vérif pour comparaison avec des calculs effectués avec un programme
// de calcul formel: par exemple maxima
// pour la vérification des dérivées / mu1, on utilise l'outil yacas: les informations sont dans les fichiers remarques_verif_yacas
// associées à plusieurs autres fichiers... ça marche super !!!

//------ débug -------
/*
   if (indic_Verif_CompaFormel==1)
	{ cout << "\n\n----------------------------------------------------------------";
	  int prec = ParaGlob::NbdigdoCA();
     cout << "\n V= " << setprecision(prec) << inv.V << " Qeps= " 
	                   << setprecision(prec) << i_s.Qeps << " cos3phi= " 
							 << setprecision(prec) << i_s.cos3phi;
	  cout << "\n Qs= " << setprecision(prec) << Qs 
	       << " dQs_dcos3phi= " << setprecision(prec) << dQs_dcos3phi   
			 << " dQs_dcos3phi2= " << setprecision(prec) << dQs_dcos3phi2 
	       << "\n Qe= " << setprecision(prec) << Qe 
			 << " dQe_dcos3phi= " << setprecision(prec) << dQe_dcos3phi   
			 << " dQe_dcos3phi2= " << setprecision(prec) << dQe_dcos3phi2
	       << "\n mu1= " << setprecision(prec) << mu1 
			 << " dMu1_dcos3phi= " << setprecision(prec) << dMu1_dcos3phi 
			 << " dMu1_dcos3phi2= " << setprecision(prec) << dMu1_dcos3phi2
	       << "\n mu2= " << setprecision(prec) << mu2 
			 << " dMu2_dcos3phi= " << setprecision(prec) << dMu2_dcos3phi 
			 << " dMu2_dcos3phi2= " << setprecision(prec) << dMu2_dcos3phi2
	       << "\n mu3= " << setprecision(prec) << mu3 
			 << " dMu3_dcos3phi= " << setprecision(prec) << dMu3_dcos3phi 
			 << " dMu3_dcos3phi2= " << setprecision(prec) << dMu3_dcos3phi2;

      cout << "\n\n Qec= " << setprecision(prec) << Qec << " dQec_Qs= " << setprecision(prec) << dQec_Qs 
			  << " dQec_QsQs= " << setprecision(prec) << dQec_QsQs 
		     << "\n  dQec_Mu1= " << setprecision(prec) << dQec_Mu1 << " dQec_Mu1Mu1= " << setprecision(prec) << dQec_Mu1Mu1 
			  << " dQec_QsMu1= " << setprecision(prec) << dQec_QsMu1
           << "\n  S1= " << setprecision(prec) << S1 << " dS1_Qs= " << setprecision(prec) << dS1_Qs 
			  << " dS1_QsQs= " << setprecision(prec) << dS1_QsQs
			  << "\n  dS1_Mu1= " << setprecision(prec) << dS1_Mu1 << " dS1_Mu1Mu1= " << setprecision(prec) << dS1_Mu1Mu1 
			  << " dS1_QsMu1= " << setprecision(prec) << dS1_QsMu1
           << "\n  S2= " << setprecision(prec) << S2 << "  dS2_Qs= " << setprecision(prec) << dS2_Qs 
			  << " dS2_QsQs= " << setprecision(prec) << dS2_QsQs 
           << "\n  dS2_Q= " << setprecision(prec) << dS2_Q << " dS2_QQ= " << setprecision(prec) << dS2_QQ 
			  << " dS2_QQs= " << setprecision(prec) << dS2_QQs
			  << "\n  dS2_Mu1= " << setprecision(prec) << dS2_Mu1 << " dS2_Mu1Mu1= " << setprecision(prec) << dS2_Mu1Mu1 
			  << " dS2_QsMu1= " << setprecision(prec) << dS2_QsMu1 << " dS2_QMu1= " << setprecision(prec) << dS2_QMu1
           << "\n  S3= " << setprecision(prec) << S3 << "  dS3_Qe= " << setprecision(prec) << dS3_Qe 
			  << " dS3_QeQe= " << setprecision(prec) << dS3_QeQe 
           << "\n  dS3_Q= " << setprecision(prec) << dS3_Q << " dS3_QQ= " << setprecision(prec) << dS3_QQ 
			  << " dS3_QQe= " << setprecision(prec) << dS3_QQe
           << "\n  S4= " << S4 << "  dS4_Qe= " << dS4_Qe << " dS4_QeQe= " << dS4_QeQe; 

      cout << "\n\n  dE_Qs= " << setprecision(prec) << dE_Qs << " dE_Qe= " << setprecision(prec) << dE_Qe 
		     << "\n dE_Mu1= " << setprecision(prec) << dE_Mu1 << "\n dE_Mu2= " << setprecision(prec) << dE_Mu2 
			  << " dE_mu3= " << setprecision(prec) << dE_Mu3  
           << "\n dE_QsQs= " << setprecision(prec) << dE_QsQs << " dE_QeQe= " << setprecision(prec) << dE_QeQe 
			  << " dE_Mu3Qe= " << setprecision(prec) << dE_Mu3Qe
			  << "\n dE_QsMu1= " << setprecision(prec) << dE_QsMu1 << " dE_Mu1Mu1= " << setprecision(prec) << dE_Mu1Mu1;   
      cout << "\n dE_QQs= " << setprecision(prec) << dE_QQs << " dE_QQe= " << setprecision(prec) << dE_QQe 
		     << " dE_QMu2= " << setprecision(prec) << dE_QMu2 << " dE_QMu3= " << setprecision(prec) << dE_QMu3
		     << "\n dE_QMu1= " << setprecision(prec) << dE_QMu1;
	   indic_Verif_CompaFormel++; // pour éviter une boucle infinie
	};
    if (indic_Verif_CompaFormel==0) 
       CompaFormel();
*/
// ------- fin débug --------       
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

             
///=============== fonctions pour la vérification et la mise au point ===============
    
    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var
            (const double & Qeps,const Invariant & inv
               ,const PoGrenobleSansPhaseAvecVar& potret )
 {  // les grandeurs de base, qui ici ne dépendent pas de la phase 
    double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
    switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	  };
    switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe =  Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe =  Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	  };
    // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_PoGrenobleorgeas1_et_var == 0)
     { cout << "\n ****verification des derivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var \n";
     };
    indic_Verif_PoGrenobleorgeas1_et_var++;
    // potret contiend le potentiel et ses variations première et seconde
    // on va calculer ces mêmes variations par différences finies et comparer les deux résultats
    // calcul des invariants sous la nouvelle forme
    
    double toto = potret.E; // pour que le débugger affiche potret
    
    Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl 
    double E_n = PoGrenoble(Qeps,inv_n0);  // la valeur du potentiel de référence
    // ici on est sans phase donc deux invariants indépendant V et Qeps,
    // pour calculer les variations on définit des points distants d'un incrément puis de deux incréments
    // pour les dérivées premières et secondes
    double delta = ConstMath::unpeupetit;
    double Qeps_et_dQeps = Qeps+delta;
    double E_et_dQeps = PoGrenoble(Qeps_et_dQeps,inv_n0);
    
    double Qeps_et_dQeps2 = Qeps + 2.*delta;
    double E_et_dQeps2 = PoGrenoble(Qeps_et_dQeps2,inv_n0);
    
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    double E_et_dV = PoGrenoble(Qeps,inv_et_dV);
    
    Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
    double E_et_dV2 = PoGrenoble(Qeps,inv_et_dV2);
    
    Invariant inv_et_dVdQeps = inv_n0;
    inv_et_dVdQeps.V += delta;
    double Qeps_et_dVdQeps = Qeps+delta;
    double E_et_dVdQeps = PoGrenoble(Qeps_et_dVdQeps,inv_et_dVdQeps);
    // calcul des dérivées premières
    double E_V = (E_et_dV  - E_n)/delta;
    double E_Qeps = (E_et_dQeps - E_n)/delta;
    // calcul des dérivées secondes
    double E_V2 = (E_et_dV2 - 2.*E_et_dV  + E_n ) /(delta * delta);
    double E_Qeps2 = (E_et_dQeps2 - 2.*E_et_dQeps + E_n)/(delta * delta);
    double E_V_a_dQeps = (E_et_dVdQeps  - E_et_dQeps )/delta;
    double E_VQeps = ( E_V_a_dQeps - E_V)/delta;
    // comparaison avec les valeurs de dérivées analytiques
    bool erreur = false;
    if (diffpourcent(potret.EV,E_V,MaX(Dabs(E_V),Dabs(potret.EV)),0.05))
     {if (MiN(Dabs(E_V),Dabs(potret.EV)) == 0.)
        {if ( MaX(Dabs(E_V),Dabs(potret.EV)) > 50.*delta) erreur = true;}
      else erreur = true;  
     };
    if (diffpourcent(potret.EQ,E_Qeps,MaX(Dabs(E_Qeps),Dabs(potret.EQ)),0.05))
     {if (MiN(Dabs(E_Qeps),Dabs(potret.EQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps),Dabs(potret.EQ)) > 50.*delta) erreur = true;}
      else erreur = true;  
     };
    if (diffpourcent(potret.EQQ,E_Qeps2,MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)),0.05))
     {if (MiN(Dabs(E_Qeps2),Dabs(potret.EQQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)) > 50.*delta) erreur = true;}
      else erreur = true;  
     };
    if (diffpourcent(potret.EVV,E_V2,MaX(Dabs(E_V2),Dabs(potret.EVV)),0.05))
     {if (MiN(Dabs(E_V2),Dabs(potret.EVV)) == 0.)
        {if ( MaX(Dabs(E_V2),Dabs(potret.EVV)) > 50.*delta) erreur = true;}
      else erreur = true;
     };
    if (diffpourcent(potret.EQV,E_VQeps,MaX(Dabs(E_VQeps),Dabs(potret.EQV)),0.05))
     {if (MiN(Dabs(E_VQeps),Dabs(potret.EQV)) == 0.)
       {if (Dabs(potret.EQV) == 0.)
          // si la valeur de la dérivée numérique est nulle cela signifie peut-être
          // que l'on est à un extréma, la seule chose que le l'on peut faire est de
          // vérifier que l'on tend numériquement vers cette extréma ici un minima
          {
            Invariant inv_et_ddVddQeps = inv_n0;
            inv_et_ddVddQeps.V += delta/10.;
            double Qeps_et_ddVddQeps = Qeps+delta/10.;
            double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
            double Qeps_et_ddQeps = Qeps+delta/10.;
            double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
            double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
            if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps)) erreur = true;
           }
        else if ( MaX(Dabs(E_VQeps),Dabs(potret.EQV)) > 150.*delta) erreur = true;
        }
      else erreur = true;
     };
    if (erreur)
     { cout << "\n erreur dans le calcul analytique des dérivees du potentiel";
       cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var(.."
            << " , numero d'increment = " << indic_Verif_PoGrenobleorgeas1_et_var;
       Sortie(1);
     };
 };

    
    // vérif des dérivées du potentiels par rapport à l'invariant V, ceci par différences finies
void IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var_VV
            (const double & Qeps,const Invariant & inv
               ,const PoGrenoble_VV& potret )
 {  // les grandeurs de base, qui ici ne dépendent pas de la phase 
    double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
    switch (cas_Q0s_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon fct laurent
 	   { if ((*temperature) >= Tc)
 		   {Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
 	     else	
 		   {Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
 	     break;
 	   }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qs = Q0s_temperature->Valeur(*temperature);
 	     break;
 	   }
 	  };
    switch (cas_Q0e_thermo_dependant)
    {case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
 	   { Qe =  Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
 	     break;
 	    }
     case 2:// cas d'une dépendance à la température selon une fonction d'évolution
 	   { Qe =  Q0e_temperature->Valeur(*temperature);
 	     break;
 	   }
 	  };
    // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_PoGrenobleorgeas1_et_var_VV == 0)
     { cout << "\n ****verification des derivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var \n";
     };
    indic_Verif_PoGrenobleorgeas1_et_var_VV++;
    // potret contiend le potentiel et ses variations première et seconde
    // on va calculer ces mêmes variations par différences finies et comparer les deux résultats
    // calcul des invariants sous la nouvelle forme
    
    double toto = potret.E; // pour que le débugger affiche potret
    
    Invariant inv_n0(inv); // recopie des invariants  
    double E_n = PoGrenoble(Qeps,inv_n0);  // la valeur du potentiel de référence
    // ici on est sans phase donc deux invariants indépendant V et Qeps,
    // pour calculer les variations on définit des points distants d'un incrément puis de deux incréments
    // pour les dérivées premières et secondes
    double delta = ConstMath::unpeupetit;
    
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    double E_et_dV = PoGrenoble(Qeps,inv_et_dV);
    
    Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
    double E_et_dV2 = PoGrenoble(Qeps,inv_et_dV2);
    
    // calcul des dérivées premières
    double E_V = (E_et_dV  - E_n)/delta;
    // calcul des dérivées secondes
    double E_V2 = (E_et_dV2 - 2.*E_et_dV  + E_n ) /(delta * delta);
    // comparaison avec les valeurs de dérivées analytiques
    bool erreur = false;
    if (diffpourcent(potret.EV,E_V,MaX(Dabs(E_V),Dabs(potret.EV)),0.05))
     {if (MiN(Dabs(E_V),Dabs(potret.EV)) == 0.)
        {if ( MaX(Dabs(E_V),Dabs(potret.EV)) > 50.*delta) erreur = true;}
      else erreur = true;  
     };
    if (diffpourcent(potret.EVV,E_V2,MaX(Dabs(E_V2),Dabs(potret.EVV)),0.05))
     {if (MiN(Dabs(E_V2),Dabs(potret.EVV)) == 0.)
        {if ( MaX(Dabs(E_V2),Dabs(potret.EVV)) > 50.*delta) erreur = true;}
      else erreur = true;
     };
    if (erreur)
     { cout << "\n erreur dans le calcul analytique des derivees du potentiel selon V et VV";
       cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var_VV(.."
            << " , numero d'increment = " << indic_Verif_PoGrenobleorgeas1_et_var_VV;
       Sortie(1);
     };
 };

    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var
            (const double & Qeps,const Invariant & inv,const double& cos3phi
               ,const PoGrenobleAvecPhaseAvecVar& potret )
  { // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_PoGrenobleorgeas1_et_var == 0)
     { cout << "\n ****vérification des dérivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var \n";
     };
    indic_Verif_PoGrenobleorgeas1_et_var++;
    // potret contiend le potentiel et ses variations première et seconde
    // on va calculer ces mêmes variations par différences finies et comparer les deux résultats
    // calcul des invariants sous la nouvelle forme
    
    double toto = potret.E; // pour que le débugger affiche potret
    
    Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl 
    Invariant0QepsCosphi invcos3phi(Qeps,cos3phi);
    double E_n = PoGrenoble(invcos3phi,inv_n0);  // la valeur du potentiel de référence
    // ici on est sans phase donc deux invariants indépendant V et Qeps,
    // pour calculer les variations on définit des points distants d'un incrément puis de deux incréments
    // pour les dérivées premières et secondes
    double delta = 100.*ConstMath::unpeupetit;
    double Qeps_et_dQeps = Qeps+delta;
    Invariant0QepsCosphi invcos3phi1(Qeps_et_dQeps,cos3phi);
    double E_et_dQeps = PoGrenoble(invcos3phi1,inv_n0);
    
    double Qeps_moins_dQeps = Qeps - delta;
    Invariant0QepsCosphi invcos3phi2(Qeps_moins_dQeps,cos3phi);
    double E_moins_dQeps = PoGrenoble(invcos3phi2,inv_n0);
    
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    double E_et_dV = PoGrenoble(invcos3phi,inv_et_dV);
    
    Invariant inv_moins_dV = inv_n0;inv_moins_dV.V -=  delta;
    double E_moins_dV = PoGrenoble(invcos3phi,inv_moins_dV);
    
    Invariant inv_et_dVdQeps = inv_n0;
    inv_et_dVdQeps.V += delta;
    double Qeps_et_dVdQeps = Qeps+delta;
    Invariant0QepsCosphi invcos3phi3(Qeps_et_dVdQeps,cos3phi);
    double E_et_dVdQeps = PoGrenoble(invcos3phi3,inv_et_dVdQeps);
    
    // cas des variations avec cos3phi    
    double cos3phi_et_dcos3phi = cos3phi+delta;
    Invariant0QepsCosphi invcos3phi_et_dcos3phi(Qeps,cos3phi_et_dcos3phi);
    double E_et_dcos3phi = PoGrenoble(invcos3phi_et_dcos3phi,inv_n0);
    
    double cos3phi_moins_dcos3phi = cos3phi-delta;
    Invariant0QepsCosphi invcos3phi_moins_dcos3phi(Qeps,cos3phi_moins_dcos3phi);
    double E_moins_dcos3phi = PoGrenoble(invcos3phi_moins_dcos3phi,inv_n0);
    
    Invariant0QepsCosphi invcos3phi_et_dcos3phi_a_dQeps(Qeps_et_dQeps,cos3phi_et_dcos3phi);
    double E_et_dcos3phi_a_dQeps = PoGrenoble(invcos3phi_et_dcos3phi_a_dQeps,inv_n0);
    
    double E_et_dcos3phi_a_dV = PoGrenoble(invcos3phi_et_dcos3phi,inv_et_dV);
    
    
    // calcul des dérivées premières
    double E_V = (E_et_dV  - E_n)/delta;
    double E_Qeps = (E_et_dQeps - E_n)/delta;
    double E_cos3phi = (E_et_dcos3phi -E_n)/delta;
    // calcul des dérivées secondes
    double E_V2 = (E_et_dV - 2.*E_n  + E_moins_dV ) /(delta * delta);
    double E_Qeps2 = (E_et_dQeps - 2.*E_n + E_moins_dQeps)/(delta * delta);
    double E_V_a_dQeps = (E_et_dVdQeps  - E_et_dQeps )/delta;
    double E_VQeps = ( E_V_a_dQeps - E_V)/delta;
    double E_cos3phi_2 = (E_et_dcos3phi - 2.*E_n  + E_moins_dcos3phi ) /(delta * delta);
    double E_dQeps_a_dcos3phi = (E_et_dcos3phi_a_dQeps  - E_et_dcos3phi )/delta;
    double E_dcos3phi_dQeps = (E_dQeps_a_dcos3phi  - E_Qeps )/delta;
    double E_dV_a_dcos3phi = (E_et_dcos3phi_a_dV  - E_et_dcos3phi )/delta;
    double E_dcos3phi_dV = (E_dV_a_dcos3phi  - E_V )/delta;

 // dans les dérivées secondes par rapport à cos3phi on a des erreurs, sans doute à cause des petites valeurs ??
  
    // comparaison avec les valeurs de dérivées analytiques
    int erreur = 0;
    if (diffpourcent(potret.EV,E_V,MaX(Dabs(E_V),Dabs(potret.EV)),0.05))
     {if (MiN(Dabs(E_V),Dabs(potret.EV)) == 0.)
        {if ( MaX(Dabs(E_V),Dabs(potret.EV)) > 50.*delta) erreur = 1;}
      else 
       erreur = 2; } 
    if (diffpourcent(potret.EQ,E_Qeps,MaX(Dabs(E_Qeps),Dabs(potret.EQ)),0.05)) 
     {if (MiN(Dabs(E_Qeps),Dabs(potret.EQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps),Dabs(potret.EQ)) > 50.*delta) erreur = 3;}
      else 
        erreur = 4;}        
    if (diffpourcent(potret.Ecos3phi,E_cos3phi,MaX(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)),0.05)) 
     {if (MiN(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)) == 0.)
        {if ( MaX(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)) > 50.*delta) erreur = 31;}
      else 
        erreur = 41;}
    if (diffpourcent(potret.EQQ,E_Qeps2,MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)),0.05)) 
     {if (MiN(Dabs(E_Qeps2),Dabs(potret.EQQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)) > 50.*delta) erreur = 5;}
      else 
        erreur = 6;}  
    if (diffpourcent(potret.EVV,E_V2,MaX(Dabs(E_V2),Dabs(potret.EVV)),0.05)) 
     {if (MiN(Dabs(E_V2),Dabs(potret.EVV)) == 0.)
        {if ( MaX(Dabs(E_V2),Dabs(potret.EVV)) > 50.*delta) erreur = 7;}
      else 
        erreur = 8;}        
    if (diffpourcent(potret.Ecos3phi2,E_cos3phi_2,MaX(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)),0.05)) 
     {if (MiN(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)) == 0.)
        {if ( MaX(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)) > 50.*delta) erreur = 71;}
      else 
        {erreur = 81;cout << "  " << potret.Ecos3phi2 << " " <<E_cos3phi_2; }
     }
    if (diffpourcent(potret.EQV,E_VQeps,MaX(Dabs(E_VQeps),Dabs(potret.EQV)),0.05)) 
     {if (MiN(Dabs(E_VQeps),Dabs(potret.EQV)) == 0.)
      {if (Dabs(potret.EQV) == 0.)
         // si la valeur de la dérivée numérique est nulle cela signifie peut-être
         // que l'on est à un extréma, la seule chose que le l'on peut faire est de
         // vérifier que l'on tend numériquement vers cette extréma ici un minima
         {
           Invariant inv_et_ddVddQeps = inv_n0;
           inv_et_ddVddQeps.V += delta/10.;
           double Qeps_et_ddVddQeps = Qeps+delta/10.;
           double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
           double Qeps_et_ddQeps = Qeps+delta/10.;
           double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
           double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
           if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps)) erreur = 9;
          }
       else 
         if ( MaX(Dabs(E_VQeps),Dabs(potret.EQV)) > 150.*delta) erreur = 10;
        } 
      else 
        erreur = 11;
      }
      
    if (diffpourcent(potret.EVcos3phi,E_dcos3phi_dV,MaX(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)),0.05)) 
     {if (MiN(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)) == 0.)
      {if (Dabs(potret.EVcos3phi) == 0.)
         // si la valeur de la dérivée numérique est nulle cela signifie peut-être
         // que l'on est à un extréma, la seule chose que le l'on peut faire est de
         // vérifier que l'on tend numériquement vers cette extréma ici un minima
         {
           /*Invariant inv_et_ddVddQeps = inv_n0;
           inv_et_ddVddQeps.V += delta/10.;
           double Qeps_et_ddVddQeps = Qeps+delta/10.;
           double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
           double Qeps_et_ddQeps = Qeps+delta/10.;
           double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
           double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
           if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps))*/ erreur = 91;
          }
       else 
         if ( MaX(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)) > 150.*delta) erreur = 101;
        } 
      else 
        erreur = 111;
      }
      
    if (diffpourcent(potret.EQcos3phi,E_dcos3phi_dQeps,MaX(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)),0.05)) 
     {if (MiN(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)) == 0.)
      {if (Dabs(potret.EQcos3phi) == 0.)
         // si la valeur de la dérivée numérique est nulle cela signifie peut-être
         // que l'on est à un extréma, la seule chose que le l'on peut faire est de
         // vérifier que l'on tend numériquement vers cette extréma ici un minima
         {
           /*Invariant inv_et_ddVddQeps = inv_n0;
           inv_et_ddVddQeps.V += delta/10.;
           double Qeps_et_ddVddQeps = Qeps+delta/10.;
           double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
           double Qeps_et_ddQeps = Qeps+delta/10.;
           double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
           double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
           if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps))*/ erreur = 92;
          }
       else 
         if ( MaX(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)) > 150.*delta) erreur = 102;
        } 
      else 
        erreur = 112;
      }
      
        
    if (erreur > 0)
     { cout << "\n erreur dans le calcul analytique des derivees du potentiel";
       cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var(.."
            << " , numero d'increment = " << indic_Verif_PoGrenobleorgeas1_et_var
            << " , numero d'erreur : " << erreur ;
 //      Sortie(1);
      } 
 };
 
// utilitaire de vérif pour comparaison avec des calculs effectués avec un programme
// de calcul formel: par exemple maxima
void IsoHyper3DOrgeas1::CompaFormel()
 { indic_Verif_CompaFormel++; // pour éviter une boucle infinie
   // l'objectif ici est d'obtenir des résultats numériques à partir de données entrées interactivement
   bool fin_interactif = false;
   // on définit des invariants initialisées à 0 par défaut
   Invariant2QepsCosphi i_s;
   Invariant  inv;
   while (!fin_interactif)
   	{ cout << "\n variation relative de volume (rep V   ) ? "
   	       << "\n Qeps                         (rep Qeps) ? "
   	       << "\n cos(3phi)                    (rep cos)  ? "
   	       << "\n affichage                    (rep a )   ? "
   	       << "\n arret verif                  (rep f )   ? ";
   	  string rep;
      // procédure de lecture avec prise en charge d'un retour chariot
      rep = lect_return_defaut(false,"o");
   	  if (rep == "V")
   	  	{ cout << "\n valeur de V "; inv.V=lect_double();}
   	  else if (rep == "Qeps")	
   	  	{ cout << "\n valeur de Qeps "; i_s.Qeps=lect_double();}
   	  else if (rep == "cos")	
   	  	{ cout << "\n valeur de cos3phi "; i_s.cos3phi=lect_double();}
   	  else if (rep=="f")
   	  	{fin_interactif=true;}
		  else if (rep=="a")
		   {	
   	  // appel du calcul des différentes valeurs
   	  PoGrenobleAvecPhaseAvecVar pot = PoGrenoblePhase_et_var(i_s,inv);
   	  // affichage
   	  cout << "\n\n----------------------------------------------------------------"
		       << "\n V= " << inv.V << " Qeps= " << i_s.Qeps << " cos3phi= " << i_s.cos3phi << "\n\n"; 
   	  
      cout << "\n  E=          "<< pot.E; // valeur du potentiel
      cout << "\n  EV=         "<< pot.EV; // variation du potentiel par rapport à V
      cout << "\n  EQ=         "<< pot.EQ; // variation du potentiel par rapport à Q
      cout << "\n  Ecos3phi=   "<< pot.Ecos3phi; // variation du potentiel par rapport à cos3phi
      cout << "\n  EVV=        "<< pot.EVV;   // variation seconde par rapport à V
      cout << "\n  EQQ=        "<< pot.EQQ; // variation seconde par rapport à Q
      cout << "\n  Ecos3phi2=  "<< pot.Ecos3phi2; // variation seconde par rapport à cos3phi
      cout << "\n  EQV=        "<< pot.EQV; // variation seconde par rapport à Q et V
      cout << "\n  EVcos3phi=  "<< pot.EVcos3phi; // variation seconde par rapport à V cos3phi
      cout << "\n  EQcos3phi=  "<< pot.EQcos3phi; // variation seconde par rapport à Q cos3phi
      cout << "\n\n\n" << endl;
		  }
   	};
		
	};