// 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 "IsoHyper3DOrgeas2.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 IsoHyper3DOrgeas2::indic_Verif_PoGrenobleorgeas2_et_var = 0;                            
// indicateur utilisé par Verif_Potentiel_et_var_VV  
int IsoHyper3DOrgeas2::indic_Verif_PoGrenobleorgeas2_et_var_VV = 0; 
// indicateur utilisé par CompaFormel                             
int IsoHyper3DOrgeas2::indic_Verif_CompaFormel = 0;                          
//================== fin d'initialisation des variables static ================


IsoHyper3DOrgeas2::IsoHyper3DOrgeas2 ()  : // Constructeur par defaut
 Hyper3D(ISOHYPER3DORGEAS2,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.)
 ,mu01_phi(NULL),mu02_phi(NULL),mu03_phi(NULL),Q0s_phi(NULL),Q0e_phi(NULL)
 ,mu01_nD(NULL),mu02_nD(NULL),mu03_nD(NULL),Q0s_nD(NULL),Q0e_nD(NULL)
 ,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
IsoHyper3DOrgeas2::IsoHyper3DOrgeas2 (const IsoHyper3DOrgeas2& 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)
 ,mu01_phi(loi.mu01_phi),mu02_phi(loi.mu02_phi),mu03_phi(loi.mu03_phi)
 ,Q0s_phi(loi.Q0s_phi),Q0e_phi(loi.Q0e_phi)
 ,mu01_nD(loi.mu01_nD),mu02_nD(loi.mu02_nD),mu03_nD(loi.mu03_nD),Q0s_nD(loi.Q0s_nD),Q0e_nD(loi.Q0e_nD)
 ,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);
        };
	  // pour la dépendance à la phase idem
   if (mu01_phi != NULL)
    if (mu01_phi->NomCourbe() == "_")
        {// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
         string non_courbe("_");
         mu01_phi = Courbe1D::New_Courbe1D(*loi.mu01_phi);
        };
	  if (mu02_phi != NULL)
	   if (mu02_phi->NomCourbe() == "_")
        {// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
         string non_courbe("_");   
         mu02_phi = Courbe1D::New_Courbe1D(*loi.mu02_phi);
        };
	  if (mu03_phi != NULL)
	   if (mu03_phi->NomCourbe() == "_")
        {// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
         string non_courbe("_");   
         mu03_phi = Courbe1D::New_Courbe1D(*loi.mu03_phi);
        };
	  if (Q0s_phi != NULL)
	   if (Q0s_phi->NomCourbe() == "_")
        {// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
         string non_courbe("_");   
         Q0s_phi = Courbe1D::New_Courbe1D(*loi.Q0s_phi);
        };
	  if (Q0e_phi != NULL)
	   if (Q0e_phi->NomCourbe() == "_")
        {// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
         string non_courbe("_");   
         Q0e_phi = Courbe1D::New_Courbe1D(*loi.Q0e_phi);
        };
   // idem pour les fonctions nD
   if (mu01_nD != NULL)
     if (mu01_nD->NomFonction() == "_")
         {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
          string non_courbe("_");
          mu01_nD = Fonction_nD::New_Fonction_nD(*loi.mu01_nD);
         };
    if (mu02_nD != NULL)
      if (mu02_nD->NomFonction() == "_")
          {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
           string non_courbe("_");
           mu02_nD = Fonction_nD::New_Fonction_nD(*loi.mu02_nD);
          };
    if (mu03_nD != NULL)
      if (mu03_nD->NomFonction() == "_")
          {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
           string non_courbe("_");
           mu03_nD = Fonction_nD::New_Fonction_nD(*loi.mu03_nD);
          };
    if (Q0s_nD != NULL)
      if (Q0s_nD->NomFonction() == "_")
          {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
           string non_courbe("_");
           Q0s_nD = Fonction_nD::New_Fonction_nD(*loi.Q0s_nD);
          };
    if (Q0e_nD != NULL)
      if (Q0e_nD->NomFonction() == "_")
          {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
           string non_courbe("_");
           Q0e_nD = Fonction_nD::New_Fonction_nD(*loi.Q0e_nD);
          };
	  };

IsoHyper3DOrgeas2::~IsoHyper3DOrgeas2 ()
// Destructeur 
{ if (Q0s_temperature != NULL)
      if (Q0s_temperature->NomCourbe() == "_") delete Q0s_temperature;
  if (Q0e_temperature != NULL)
      if (Q0e_temperature->NomCourbe() == "_") delete Q0e_temperature;
  if (mu01_phi != NULL) if (mu01_phi->NomCourbe() == "_") delete mu01_phi;
  if (mu02_phi != NULL) if (mu02_phi->NomCourbe() == "_") delete mu02_phi;
  if (mu03_phi != NULL) if (mu03_phi->NomCourbe() == "_") delete mu03_phi;
  if (Q0s_phi != NULL) if (Q0s_phi->NomCourbe() == "_") delete Q0s_phi;
  if (Q0e_phi != NULL) if (Q0e_phi->NomCourbe() == "_") delete Q0e_phi;
  
  if (mu01_nD != NULL) if (mu01_nD->NomFonction() == "_") delete mu01_nD;
  if (mu02_nD != NULL) if (mu02_nD->NomFonction() == "_") delete mu02_nD;
  if (mu03_nD != NULL) if (mu03_nD->NomFonction() == "_") delete mu03_nD;
  if (Q0s_nD != NULL) if (Q0s_nD->NomFonction() == "_") delete Q0s_nD;
  if (Q0e_nD != NULL) if (Q0e_nD->NomFonction() == "_") delete Q0e_nD;

  };  

// Lecture des donnees de la classe sur fichier
void IsoHyper3DOrgeas2::LectureDonneesParticulieres
    (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D,LesFonctions_nD& lesFonctionsnD)
  { string nom_class_methode("IsoHyper3DOrgeas2::LectureDonneesParticulieres");
    // 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--IsoHyper3DOrgeas2::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--IsoHyper3DOrgeas2::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--IsoHyper3DOrgeas2::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--IsoHyper3DOrgeas2::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
       
       // on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_avec_phase_"
       // ou un nouveau mot clé global auquel cas il y a pb !!
       MotCle motCle; // ref aux mots cle
       while (strstr(entreePrinc->tablcar,"fin_parametres_avec_phase_")==0) 
       {
        // si on a  un mot clé global dans la ligne courante c-a-d dans tablcar --> erreur
        if ( motCle.SimotCle(entreePrinc->tablcar))
         { cout << "\n erreur de lecture des parametre de dependance a la phase: on n'a pas trouve le mot cle "
                << " fin_parametres_avec_phase_ et par contre la ligne courante contient un mot cle global  ";
           entreePrinc->MessageBuffer
                ("** erreur5 des parametres de reglage de la loi de comportement d'IsoHyper3DOrgeas2**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
         
         // lecture d'un mot clé
         *(entreePrinc->entree) >> nom1;
         
         if ((entreePrinc->entree)->rdstate() == 0) 
           {} // lecture normale
         #ifdef ENLINUX 
         else  if ((entreePrinc->entree)->fail())
            // on a atteind la fin de la ligne et on appelle un nouvel enregistrement 
           {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
               *(entreePrinc->entree) >>nom1;
            }
         #else
         else  if ((entreePrinc->entree)->eof())
           // la lecture est bonne mais on a atteind la fin de la ligne
           { if(nom1 != "fin_parametres_avec_phase_")
               {entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom1;};
           } 
         #endif  
         else // cas d'une erreur de lecture
          { cout << "\n erreur de lecture inconnue  ";
            entreePrinc->MessageBuffer
              ("** erreur6 des parametres de reglage de la loi de comportement d'IsoHyper3DOrgeas2**");
            throw (UtilLecture::ErrNouvelleDonnee(-1));
            Sortie(1);     
          };
        
         if(nom1 == "mu01_phi=")
          {// lecture de la loi d'évolution
           *(entreePrinc->entree) >>  nom1;
           // on regarde si la courbe existe, si oui on récupère la référence
           if (lesCourbes1D.Existe(nom2))  { mu01_phi = lesCourbes1D.Trouve(nom2);}
           else { // sinon il faut la lire maintenant
              string non_courbe("_");
              mu01_phi = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom2.c_str()));
              // lecture de la courbe
              mu01_phi->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
               };
           // prepa du flot de lecture
           entreePrinc->NouvelleDonnee();
           }
         else if(nom1 == "mu02_phi=")
          {// lecture de la loi d'évolution 
           *(entreePrinc->entree) >>  nom2;    
           // on regarde si la courbe existe, si oui on récupère la référence 
           if (lesCourbes1D.Existe(nom2))  { mu02_phi = lesCourbes1D.Trouve(nom2);}
           else { // sinon il faut la lire maintenant
              string non_courbe("_");   
              mu02_phi = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom2.c_str()));
              // lecture de la courbe
              mu02_phi->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
               }; 
           // prepa du flot de lecture
           entreePrinc->NouvelleDonnee();
           }
         else if(nom1 == "mu03_phi=")
          {// lecture de la loi d'évolution 
           *(entreePrinc->entree) >>  nom2;    
           // on regarde si la courbe existe, si oui on récupère la référence 
           if (lesCourbes1D.Existe(nom2))  { mu03_phi = lesCourbes1D.Trouve(nom2);}
           else { // sinon il faut la lire maintenant
              string non_courbe("_");   
              mu03_phi = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom2.c_str()));
              // lecture de la courbe
              mu03_phi->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
               }; 
           // prepa du flot de lecture
           entreePrinc->NouvelleDonnee();
           }
         else if(nom1 == "Q0s_phi=")
          {// lecture de la loi d'évolution 
           *(entreePrinc->entree) >>  nom2;    
           // on regarde si la courbe existe, si oui on récupère la référence 
           if (lesCourbes1D.Existe(nom2))  { Q0s_phi = lesCourbes1D.Trouve(nom2);}
           else { // sinon il faut la lire maintenant
              string non_courbe("_");   
              Q0s_phi = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom2.c_str()));
              // lecture de la courbe
              Q0s_phi->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
               }; 
           // prepa du flot de lecture
           entreePrinc->NouvelleDonnee();
           }
         else if(nom1 == "Q0e_phi=")
          {// lecture de la loi d'évolution 
           *(entreePrinc->entree) >>  nom2;    
           // on regarde si la courbe existe, si oui on récupère la référence 
           if (lesCourbes1D.Existe(nom2))  { Q0e_phi = lesCourbes1D.Trouve(nom2);}
           else { // sinon il faut la lire maintenant
              string non_courbe("_");   
              Q0e_phi = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom2.c_str()));
              // lecture de la courbe
              Q0e_phi->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
               }; 
           // prepa du flot de lecture
           entreePrinc->NouvelleDonnee();
           }
         else   // cas d'une erreur
          { cout << "\n erreur en lecture de la dependance a la phase des coefficients, on aurait du lire "
                 << " le mot cle  mu01_phi= ou mu02_phi= ou mu03_phi ou Q0s_phi= ou Q0e_phi= et on a lue: "
                 << nom1;
            entreePrinc->MessageBuffer("**erreur7 IsoHyper3DOrgeas2::LectureDonneesParticulieres (...**");
            throw (UtilLecture::ErrNouvelleDonnee(-1));
            Sortie(1);
           };
       }; //-- fin du while
       avec_phase=true; // variable définit dans une classe mère, ici = vrai par défaut
       *(entreePrinc->entree) >> nom1 ; // on lit le mot clé "fin_parametres_avec_phase_"
//       entreePrinc->NouvelleDonnee();  // passage d'une ligne
     }; //--  fin lecture de la dépendance à la phase ----                       

    // ---- 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("**8--IsoHyper3DOrgeas2::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("**9--IsoHyper3DOrgeas2::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("**10--IsoHyper3DOrgeas2::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("**11--IsoHyper3DOrgeas2::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("**12--IsoHyper3DOrgeas2::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("**13--IsoHyper3DOrgeas2::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("**17--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;
          }
 	      };
		     };
     };
     
    // -------- on regarde s'il y a une dépendance à une fonction nD
    if (strstr(entreePrinc->tablcar,"avec_nD_")!=NULL)
    { string nom1,nom2;
      // lecture des  fonctions nD
      entreePrinc->NouvelleDonnee();  // passage d'une ligne
      
      // on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_avec_nD_"
      // ou un nouveau mot clé global auquel cas il y a pb !!
      MotCle motCle; // ref aux mots cle
      while (strstr(entreePrinc->tablcar,"fin_parametres_avec_nD_")==0)
      {
        // si on a  un mot clé global dans la ligne courante c-a-d dans tablcar --> erreur
        if ( motCle.SimotCle(entreePrinc->tablcar))
         { cout << "\n erreur de lecture des parametre de dependance a une fonction nD: on n'a pas trouve le mot cle "
                << " fin_parametres_avec_nD_ et par contre la ligne courante contient un mot cle global  ";
           entreePrinc->MessageBuffer
                ("** erreur51 des parametres de reglage de la loi de comportement d'IsoHyper3DOrgeas2**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
        
        // lecture d'un mot clé
        *(entreePrinc->entree) >> nom1;
        
        if ((entreePrinc->entree)->rdstate() == 0)
          {} // lecture normale
        #ifdef ENLINUX
        else  if ((entreePrinc->entree)->fail())
           // on a atteind la fin de la ligne et on appelle un nouvel enregistrement
          {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
              *(entreePrinc->entree) >>nom1;
           }
        #else
        else  if ((entreePrinc->entree)->eof())
          // la lecture est bonne mais on a atteind la fin de la ligne
          { if(nom1 != "fin_parametres_avec_nD_")
              {entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom1;};
          }
        #endif
        else // cas d'une erreur de lecture
         { cout << "\n erreur de lecture inconnue  ";
           entreePrinc->MessageBuffer
             ("** erreur61 des parametres de reglage de la loi de comportement d'IsoHyper3DOrgeas2**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
       
        if(nom1 == "mu01_nD=")
         {string nom_fonct;
          // lecture du nom de la fonction
          *(entreePrinc->entree) >>  nom_fonct;
          
          // maintenant on définit la fonction
          if (lesFonctionsnD.Existe(nom_fonct))
            {mu01_nD = lesFonctionsnD.Trouve(nom_fonct);
            }
          else
            {// sinon il faut la lire maintenant
             string non("_");
             mu01_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
             // lecture de la fonction
             mu01_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
             // maintenant on vérifie que la fonction est utilisable
             if (mu01_nD->NbComposante() != 1 )
              { cout << "\n erreur en lecture, la fonction " << nom_fonct
                     << " est une fonction vectorielle a   " << mu01_nD->NbComposante()
                     << " composante alors qu'elle devrait etre scalaire ! "
                     << " elle n'est donc pas utilisable !! ";
                string message("\n**erreur031** \n"+nom_class_methode+"(...");
                entreePrinc->MessageBuffer(message);
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);
              };
            };
          // on regarde si la fonction nD intègre la température
          const Tableau <Ddl_enum_etendu>& tab_enu = mu01_nD->Tab_enu_etendu();
          if (tab_enu.Contient(TEMP))
            thermo_dependant=true;
          // prepa du flot de lecture
          entreePrinc->NouvelleDonnee();
         }
        else if(nom1 == "mu02_nD=")
         {string nom_fonct;
          // lecture du nom de la fonction
          *(entreePrinc->entree) >>  nom_fonct;
          
          // maintenant on définit la fonction
          if (lesFonctionsnD.Existe(nom_fonct))
            {mu02_nD = lesFonctionsnD.Trouve(nom_fonct);
            }
          else
            {// sinon il faut la lire maintenant
             string non("_");
             mu02_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
             // lecture de la fonction
             mu02_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
             // maintenant on vérifie que la fonction est utilisable
             if (mu02_nD->NbComposante() != 1 )
              { cout << "\n erreur en lecture, la fonction " << nom_fonct
                     << " est une fonction vectorielle a   " << mu02_nD->NbComposante()
                     << " composante alors qu'elle devrait etre scalaire ! "
                     << " elle n'est donc pas utilisable !! ";
                string message("\n**erreur031** \n"+nom_class_methode+"(...");
                entreePrinc->MessageBuffer(message);
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);
              };
            };
          // on regarde si la fonction nD intègre la température
          const Tableau <Ddl_enum_etendu>& tab_enu = mu02_nD->Tab_enu_etendu();
          if (tab_enu.Contient(TEMP))
            thermo_dependant=true;
          // prepa du flot de lecture
          entreePrinc->NouvelleDonnee();
         }
        else if(nom1 == "mu03_nD=")
         {string nom_fonct;
          // lecture du nom de la fonction
          *(entreePrinc->entree) >>  nom_fonct;
          
          // maintenant on définit la fonction
          if (lesFonctionsnD.Existe(nom_fonct))
            {mu03_nD = lesFonctionsnD.Trouve(nom_fonct);
            }
          else
            {// sinon il faut la lire maintenant
             string non("_");
             mu03_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
             // lecture de la fonction
             mu03_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
             // maintenant on vérifie que la fonction est utilisable
             if (mu03_nD->NbComposante() != 1 )
              { cout << "\n erreur en lecture, la fonction " << nom_fonct
                     << " est une fonction vectorielle a   " << mu03_nD->NbComposante()
                     << " composante alors qu'elle devrait etre scalaire ! "
                     << " elle n'est donc pas utilisable !! ";
                string message("\n**erreur031** \n"+nom_class_methode+"(...");
                entreePrinc->MessageBuffer(message);
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);
              };
            };
          // on regarde si la fonction nD intègre la température
          const Tableau <Ddl_enum_etendu>& tab_enu = mu03_nD->Tab_enu_etendu();
          if (tab_enu.Contient(TEMP))
            thermo_dependant=true;
          // prepa du flot de lecture
          entreePrinc->NouvelleDonnee();
         }
        else if(nom1 == "Q0s_nD=")
         {string nom_fonct;
          // lecture du nom de la fonction
          *(entreePrinc->entree) >>  nom_fonct;
          
          // maintenant on définit la fonction
          if (lesFonctionsnD.Existe(nom_fonct))
            {Q0s_nD = lesFonctionsnD.Trouve(nom_fonct);
            }
          else
            {// sinon il faut la lire maintenant
             string non("_");
             Q0s_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
             // lecture de la fonction
             Q0s_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
             // maintenant on vérifie que la fonction est utilisable
             if (Q0s_nD->NbComposante() != 1 )
              { cout << "\n erreur en lecture, la fonction " << nom_fonct
                     << " est une fonction vectorielle a   " << Q0s_nD->NbComposante()
                     << " composante alors qu'elle devrait etre scalaire ! "
                     << " elle n'est donc pas utilisable !! ";
                string message("\n**erreur031** \n"+nom_class_methode+"(...");
                entreePrinc->MessageBuffer(message);
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);
              };
            };
          // on regarde si la fonction nD intègre la température
          const Tableau <Ddl_enum_etendu>& tab_enu = Q0s_nD->Tab_enu_etendu();
          if (tab_enu.Contient(TEMP))
            thermo_dependant=true;
          // prepa du flot de lecture
          entreePrinc->NouvelleDonnee();
         }
        else if(nom1 == "Q0e_nD=")
         {string nom_fonct;
          // lecture du nom de la fonction
          *(entreePrinc->entree) >>  nom_fonct;
          
          // maintenant on définit la fonction
          if (lesFonctionsnD.Existe(nom_fonct))
            {Q0e_nD = lesFonctionsnD.Trouve(nom_fonct);
            }
          else
            {// sinon il faut la lire maintenant
             string non("_");
             Q0e_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
             // lecture de la fonction
             Q0e_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
             // maintenant on vérifie que la fonction est utilisable
             if (Q0e_nD->NbComposante() != 1 )
              { cout << "\n erreur en lecture, la fonction " << nom_fonct
                     << " est une fonction vectorielle a   " << Q0e_nD->NbComposante()
                     << " composante alors qu'elle devrait etre scalaire ! "
                     << " elle n'est donc pas utilisable !! ";
                string message("\n**erreur031** \n"+nom_class_methode+"(...");
                entreePrinc->MessageBuffer(message);
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);
              };
            };
          // on regarde si la fonction nD intègre la température
          const Tableau <Ddl_enum_etendu>& tab_enu = Q0e_nD->Tab_enu_etendu();
          if (tab_enu.Contient(TEMP))
            thermo_dependant=true;
          // prepa du flot de lecture
          entreePrinc->NouvelleDonnee();
         }
        else   // cas d'une erreur
         { cout << "\n erreur en lecture de la dependance a une fonction nD des coefficients, on aurait du lire "
                << " le mot cle  mu01_nD= ou mu02_nD= ou mu03_nD ou Q0s_nD= ou Q0e_nD= et on a lue: "
                << nom1;
           entreePrinc->MessageBuffer("**erreur7 IsoHyper3DOrgeas2::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
          };
      }; //-- fin du while
      *(entreePrinc->entree) >> nom1 ; // on lit le mot clé "fin_parametres_avec_nD_"
//      entreePrinc->NouvelleDonnee();  // passage d'une ligne
    };
     
    // au cas où on passe une ligne
    if (  (strstr(entreePrinc->tablcar,"fin_parametres_avec_phase_")!=NULL)
        || (strstr(entreePrinc->tablcar,"fin_parametres_avec_nD_")!=NULL))
      entreePrinc->NouvelleDonnee();

    // 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,true);
    // 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);
    };
// affichage de la loi
void IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::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 (mu01_phi)
     {cout << " courbe mu01_phi: " << mu01_phi->NomCourbe() <<" ";
      if ( mu01_phi->NomCourbe() == "_") mu01_phi->Affiche();
     };
    if (mu02_phi)
     {cout << " courbe mu02_phi: " << mu02_phi->NomCourbe() <<" ";
      if ( mu02_phi->NomCourbe() == "_") mu02_phi->Affiche();
     };
    if (mu03_phi)
     {cout << " courbe mu03_phi: " << mu03_phi->NomCourbe() <<" ";
      if ( mu03_phi->NomCourbe() == "_") mu03_phi->Affiche();
     };
    if (Q0s_phi)
     {cout << " courbe Q0s_phi: " << Q0s_phi->NomCourbe() <<" ";
      if ( Q0s_phi->NomCourbe() == "_") Q0s_phi->Affiche();
     };
    if (Q0e_phi)
     {cout << " courbe Q0e_phi: " << Q0e_phi->NomCourbe() <<" ";
      if ( Q0e_phi->NomCourbe() == "_") Q0e_phi->Affiche();
     };
    // les dépendances à une fonction nD
    if (mu01_nD)
     {cout << " fonction nD mu01_nD: ";
      if (mu01_nD->NomFonction() != "_")
           cout << mu01_nD->NomFonction();
      else
           mu01_nD->Affiche();
     };
    if (mu02_nD)
     {cout << " fonction nD mu02_nD: ";
      if (mu02_nD->NomFonction() != "_")
           cout << mu02_nD->NomFonction();
      else
           mu02_nD->Affiche();
     };
    if (mu03_nD)
     {cout << " fonction nD mu03_nD: ";
      if (mu03_nD->NomFonction() != "_")
           cout << mu03_nD->NomFonction();
      else
           mu03_nD->Affiche();
     };
    if (Q0s_nD)
     {cout << " fonction nD Q0s_nD: ";
      if (Q0s_nD->NomFonction() != "_")
           cout << Q0s_nD->NomFonction();
      else
           Q0s_nD->Affiche();
     };
    if (Q0e_nD)
     {cout << " fonction nD Q0e_nD: ";
      if (Q0e_nD->NomFonction() != "_")
           cout << Q0e_nD->NomFonction();
      else
           Q0e_nD->Affiche();
     };
    
    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 IsoHyper3DOrgeas2::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#--------------------- dependance explicite a la phase -----------------------------------"
         << "\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 les lignes qui suivent on met des fonctions multiplicatives des parametres "
         << "\n#  initiaux, puis un mot cle signalant la fin des fonctions: fin_parametres_avec_phase_ "
         << "\n# suit un exemple de declaration (il est possible d'ommettre un groupe de 2 "
         << "\n# (nom de parametre et fonction), de meme il est possible soit d'utiliser une fonction "
         << "\n # deja defini soit de donner directement la fonction (comme pour les autres lois avec "
         << "\n # dependance a la temperature ou autre ). Enfin, apres chaque definition de fonction "
         << "\n # il faut passer a une nouvelle ligne "
         << "\n # "
         << "\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 # mu01_phi= courbe1 "
         << "\n # mu02_phi= courbe2 "
         << "\n # mu03_phi= courbe3 "
         << "\n # Q0s_phi= courbe4 "
         << "\n # Q0e_phi= courbe5 "
         << "\n # fin_parametres_avec_phase_ "
         << "\n # "
         << "\n # noter qu'il ne faut pas que mu3 devienne nulle"
         << "\n # "
         << "\n#--------------------- dependance explicite 1D a la temperature  -------------------"
         << "\n # il est egalement possible de faire dependre certains coefficients a la temperature"
         << "\n # via une fonction 1D dediee : pre-programmee ou quelconque"
         << "\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_   "
         << "\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# NB: il egalement est possible (voir si-dessous) d'utiliser la temperature via une fonction nD "
         << "\n# ceci alors pour les 5 parametres: mui et Qi"
         << "\n#"
         << "\n#--------------------- dependance explicite a une fonction nD -----------------------------------"
         << "\n# il est egalement possible de definir une dependance  des parametres a une fonction nD "
         << "\n#  dans ce cas, a la suite de la dependance explicite a la temperature via une courbr 1D "
         << "\n#  on met le mot cle: avec_nD_ "
         << "\n#  puis sur les lignes qui suivent on met des fonctions nD multiplicatives des parametres "
         << "\n#  initiaux, puis un mot cle signalant la fin des fonctions: fin_parametres_avec_nD_ "
         << "\n# suit un exemple de declaration (il est possible d'ommettre certaine fonction ! "
         << "\n# la presence de fct nD est facultative), de meme il est possible soit d'utiliser une fonction "
         << "\n # deja defini soit de donner directement la fonction nD (comme pour les autres lois) "
         << "\n # Enfin, apres chaque definition de fonction  il faut passer a une nouvelle ligne "
         << "\n # "
         << "\n#-----------------------------------------------------------------------"
         << "\n#     K    |  Q0s   |   mu01  | mu02  |  mu03  | alpha1| alpha2 | Q0e   |"
         << "\n#-----------------------------------------------------------------------"
         << "\n#  160000     150      17000    220    17000    0.002   0.004    0.01  avec_nD_ "
         << "\n # mu01_nD= fct_nD_1 "
         << "\n # mu02_nD= fct_nD_2 "
         << "\n # mu03_nD= fct_nD_3 "
         << "\n # Q0s_nD= fct_nD_4 "
         << "\n # Q0e_nD= fct_nD_5 "
         << "\n # fin_parametres_avec_nD_ "
         << "\n # "
         << "\n # Remarques: "
         << "\n # - on peut utiliser toutes les grandeurs disponibles pour alimenter la fonction nD "
         << "\n # cependant, par rapport a une dependance plus explicite (ex: phase ou temperature) "
         << "\n # l'utilisation des fonctions nD n'inclue pas la prise en compte des derives de la "
         << "\n # fonction par rapport a ses variables. Du coup, l'operateur tangent (dsig/d_ddl par ex)"
         << "\n # est moins precis que dans les cas specifiques : cf. avec_phase "
         << "\n # - les fonction nD sont des fonctions multiplicatives  "
         << "\n # des parametres initiaux (ou de ceux obtenus apres dependance a la temperature et/ou "
         << "\n # la phase) "
         << "\n # "
         << "\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 IsoHyper3DOrgeas2::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 ne sont 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 IsoHyper3DOrgeas2::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 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 IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::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;
                 }
              };
            };
           }
       else
           {ent >> toto >> Q0e ;};
       
       // puis les dépendances éventuelles à la phase
       bool test;string nom;
       ent >> nom >> test; // cas de mu01_phi
       if (test) {mu01_phi = lesCourbes1D.Lecture_pour_base_info(ent,cas,mu01_phi); };
       ent >> nom >> test; // cas de mu02_phi
       if (test) {mu02_phi = lesCourbes1D.Lecture_pour_base_info(ent,cas,mu02_phi); };
       ent >> nom >> test; // cas de mu03_phi
       if (test) {mu03_phi = lesCourbes1D.Lecture_pour_base_info(ent,cas,mu03_phi); };
       ent >> nom >> test; // cas de Q0s_phi
       if (test) {Q0s_phi = lesCourbes1D.Lecture_pour_base_info(ent,cas,Q0s_phi); };
       ent >> nom >> test; // cas de Q0e_phi
       if (test) {Q0e_phi = lesCourbes1D.Lecture_pour_base_info(ent,cas,Q0e_phi); };

       // puis les dépendances éventuelles à la phase
       ent >> nom >> test; // cas de mu01_nD
       if (test) {mu01_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mu01_nD); };
       {// on regarde si la fonction nD intègre la température
        const Tableau <Ddl_enum_etendu>& tab_enu = mu01_nD->Tab_enu_etendu();
        if (tab_enu.Contient(TEMP))
         thermo_dependant=true;
       };
       ent >> nom >> test; // cas de mu02_nD
       if (test) {mu02_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mu02_nD); };
       {// on regarde si la fonction nD intègre la température
        const Tableau <Ddl_enum_etendu>& tab_enu = mu02_nD->Tab_enu_etendu();
        if (tab_enu.Contient(TEMP))
         thermo_dependant=true;
       };
       ent >> nom >> test; // cas de mu03_nD
       if (test) {mu03_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mu03_nD); };
       {// on regarde si la fonction nD intègre la température
        const Tableau <Ddl_enum_etendu>& tab_enu = mu03_nD->Tab_enu_etendu();
        if (tab_enu.Contient(TEMP))
         thermo_dependant=true;
       };
       ent >> nom >> test; // cas de Q0s_nD
       if (test) {Q0s_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,Q0s_nD); };
       {// on regarde si la fonction nD intègre la température
        const Tableau <Ddl_enum_etendu>& tab_enu = Q0s_nD->Tab_enu_etendu();
        if (tab_enu.Contient(TEMP))
         thermo_dependant=true;
       };
       ent >> nom >> test; // cas de Q0e_nD
       if (test) {Q0e_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,Q0e_nD); };
       {// on regarde si la fonction nD intègre la température
        const Tableau <Ddl_enum_etendu>& tab_enu = Q0e_nD->Tab_enu_etendu();
        if (tab_enu.Contient(TEMP))
         thermo_dependant=true;
       };

       // calcul des constantes particulières
       alpha1_2 = Sqr(alpha1);
       alpha3_2 = Sqr(alpha3);
       avec_phase=true; // variable définit dans une classe mère, ici = vrai par défaut
       // 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 IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::Ecriture_base_info_loi(..";
                     Sortie(1);            
        		
        	 };
        }
       else
        {sort << " Q0e " << Q0e ;};
		  
		     // puis les dépendances éventuelles à la phase
       if (mu01_phi == NULL) { sort << " mu01_phi= " << false << " " ;}
       else { sort << " mu01_phi= " << true << " ";
              LesCourbes1D::Ecriture_pour_base_info(sort,cas,mu01_phi);};
       if (mu02_phi == NULL) { sort << " mu02_phi= " << false << " " ;}
       else { sort << " mu02_phi= " << true << " "; 
              LesCourbes1D::Ecriture_pour_base_info(sort,cas,mu02_phi);};
       if (mu03_phi == NULL) { sort << " mu03_phi= " << false << " " ;}
       else { sort << " mu03_phi= " << true << " "; 
              LesCourbes1D::Ecriture_pour_base_info(sort,cas,mu03_phi);};
       if (Q0s_phi == NULL) { sort << " Q0s_phi= " << false << " " ;}
       else { sort << " Q0s_phi= " << true << " "; 
              LesCourbes1D::Ecriture_pour_base_info(sort,cas,Q0s_phi);};
       if (Q0e_phi == NULL) { sort << " Q0e_phi= " << false << " " ;}
       else { sort << " Q0e_phi= " << true << " "; 
              LesCourbes1D::Ecriture_pour_base_info(sort,cas,Q0e_phi);};
              
       // puis les dépendances éventuelles à une fonction nD
       if (mu01_nD == NULL) { sort << " mu01_nD= " << false << " " ;}
       else { sort << " mu01_nD= " << true << " ";
              LesFonctions_nD::Ecriture_pour_base_info(sort,cas,mu01_nD);};
       if (mu02_nD == NULL) { sort << " mu02_nD= " << false << " " ;}
       else { sort << " mu02_nD= " << true << " ";
              LesFonctions_nD::Ecriture_pour_base_info(sort,cas,mu02_nD);};
       if (mu03_nD == NULL) { sort << " mu03_nD= " << false << " " ;}
       else { sort << " mu03_nD= " << true << " ";
              LesFonctions_nD::Ecriture_pour_base_info(sort,cas,mu03_nD);};
       if (Q0s_nD == NULL) { sort << " Q0s_nD= " << false << " " ;}
       else { sort << " Q0s_nD= " << true << " ";
              LesFonctions_nD::Ecriture_pour_base_info(sort,cas,Q0s_nD);};
       if (Q0e_nD == NULL) { sort << " Q0e_nD= " << false << " " ;}
       else { sort << " Q0e_nD= " << true << " ";
              LesFonctions_nD::Ecriture_pour_base_info(sort,cas,Q0e_nD);};
       // 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 IsoHyper3DOrgeas2::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 IsoHyper3DOrgeas2::Module_compressibilite_equivalent(Enum_dure ,const Deformation & def,SaveResul * )
   { 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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
 	  };
   // on utilise la valeur 0 pour le cos3phi, pour avoir un comportement continu
   // au cas ou la valeur correspondant à 0 ne soit pas = à 1
     {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(0.) ;
      if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(0.) ;
      if (mu01_phi) mu1 = mu01 * mu01_phi->Valeur(0.) ;
      if (mu02_phi) mu2 = mu02 * mu02_phi->Valeur(0.) ;
      if (mu03_phi) mu3 = mu03 * mu03_phi->Valeur(0.) ;
      // cas des fonctions nD éventuelles
      if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                            (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
      if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                            (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
      if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                            (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
      if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                            (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
      if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                            (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
      
     };
   // 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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
 	 };
   // on modifie les variables qui dépendent de la phase
   {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(inv_spec.cos3phi) ;
    if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(inv_spec.cos3phi) ;
    if (mu01_phi)mu2 = mu01 * mu01_phi->Valeur(inv_spec.cos3phi) ;
    if (mu02_phi)mu2 = mu02 * mu02_phi->Valeur(inv_spec.cos3phi) ;
    if (mu03_phi)mu3 = mu03 * mu03_phi->Valeur(inv_spec.cos3phi) ;
    };
   {// cas des fonctions nD éventuelles
    if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   };
   // 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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
 	 };
   // on utilise la valeur 0 pour le cos3phi, pour avoir un comportement continu
   // au cas ou la valeur correspondant à 0 ne soit pas = à 1
   {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(0.) ;
    if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(0.) ;
    if (mu01_phi)mu2 = mu01 * mu01_phi->Valeur(0.) ;
    if (mu02_phi) mu2 = mu02 * mu02_phi->Valeur(0.) ;
    if (mu03_phi) mu3 = mu03 * mu03_phi->Valeur(0.) ;
    };
   {// cas des fonctions nD éventuelles
    if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   };
   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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
 	 };
   // on modifie les variables qui dépendent de la phase
   {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(inv_spec.cos3phi) ;
    if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(inv_spec.cos3phi) ;
    if (mu01_phi)mu2 = mu01 * mu01_phi->Valeur(inv_spec.cos3phi) ;
    if (mu02_phi)mu2 = mu02 * mu02_phi->Valeur(inv_spec.cos3phi) ;
    if (mu03_phi)mu3 = mu03 * mu03_phi->Valeur(inv_spec.cos3phi) ;
    };
  {// cas des fonctions nD éventuelles
   if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
  };

   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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
    };
   // on utilise la valeur 0 pour le cos3phi, pour avoir un comportement continu
   // au cas ou la valeur correspondant à 0 ne soit pas = à 1
   {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(0.) ;
    if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(0.) ;
    if (mu01_phi)mu2 = mu01 * mu01_phi->Valeur(0.) ;
    if (mu02_phi) mu2 = mu02 * mu02_phi->Valeur(0.) ;
    if (mu03_phi) mu3 = mu03 * mu03_phi->Valeur(0.) ;
    };
   {// cas des fonctions nD éventuelles
    if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   };
   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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
 	 };
   // on modifie les variables qui dépendent de la phase
   {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(inv_spec.cos3phi) ;
    if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(inv_spec.cos3phi) ;
    if (mu01_phi)mu2 = mu01 * mu01_phi->Valeur(inv_spec.cos3phi) ;
    if (mu02_phi)mu2 = mu02 * mu02_phi->Valeur(inv_spec.cos3phi) ;
    if (mu03_phi)mu3 = mu03 * mu03_phi->Valeur(inv_spec.cos3phi) ;
    };
  {// cas des fonctions nD éventuelles
   if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
  };

   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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
 	 };
   // on utilise la valeur 0 pour le cos3phi, pour avoir un comportement continu
   // au cas ou la valeur correspondant à 0 ne soit pas = à 1
   {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(0.) ;
    if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(0.) ;
    if (mu01_phi)mu2 = mu01 * mu01_phi->Valeur(0.) ;
    if (mu02_phi) mu2 = mu02 * mu02_phi->Valeur(0.) ;
    if (mu03_phi) mu3 = mu03 * mu03_phi->Valeur(0.) ;
    };
   {// cas des fonctions nD éventuelles
    if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   };
   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);
   // 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);
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };
   
   // calcul du potentiel et de ses dérivées avec la phase
Hyper3D::PoGrenobleAvecPhaseSansVar IsoHyper3DOrgeas2::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  à 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;
 	   }
 	 };
   // on modifie les variables qui dépendent de la phase
   // et on en calcul la dérivée première
   Courbe1D::ValDer f_para;
   double dQs_dcos3phi = 0.; double dQe_dcos3phi = 0.;
   double dMu1_dcos3phi = 0.;double dMu2_dcos3phi = 0.;double dMu3_dcos3phi = 0.;
   {if (Q0s_phi )
     {f_para = Q0s_phi->Valeur_Et_derivee(i_s.cos3phi) ;
      Qs = Q0sT * f_para.valeur ;dQs_dcos3phi = Q0sT * f_para.derivee;
     };
    if (Q0e_phi )
     {f_para = Q0e_phi->Valeur_Et_derivee(i_s.cos3phi) ;
      Qe = Q0e * f_para.valeur ;dQe_dcos3phi = Q0e * f_para.derivee;
     };
    if (mu01_phi )
     {f_para = mu01_phi->Valeur_Et_derivee(i_s.cos3phi) ;
      mu1 = mu01 * f_para.valeur ;dMu1_dcos3phi = mu01 * f_para.derivee;
     };
    if (mu02_phi )
     {f_para = mu02_phi->Valeur_Et_derivee(i_s.cos3phi) ;
      mu2 = mu02 * f_para.valeur ;dMu2_dcos3phi = mu02 * f_para.derivee;
     };
    if (mu03_phi )
     {f_para = mu03_phi->Valeur_Et_derivee(i_s.cos3phi) ;
      mu3 = mu03 * f_para.valeur ;dMu3_dcos3phi = mu03 * f_para.derivee;
     };
    };
   {// cas des fonctions nD éventuelles
    if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   };
   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 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))))
           + 0.5*mu3*(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);          
//            + mu1 *(QplusQec - S2) + mu3 * (i_s.Qeps - S4 + S3);
   // cas de la variation faisant intervenir  la phase
   // 2-- variations des grandeurs intermédiaires par rapport aux grandeurs de bases
   double Qs2 = Qs * Qs;
   double dQec_Qs = - (Qs2 + Sqr(mu1) * alpha1_2) / (2. * mu1 * Qs2);
   double dS2_Qs = QplusQec * dQec_Qs / S2;
   double dS1_Qs = Qec * dQec_Qs / S1;
   double dS4_Qe = Qe / S4;
   double dS3_Qe = -(QmoinsQe) / S3;
   // 3-- variations du potentiel par rapport aux grandeurs de bases
   double dE_Qs = i_s.Qeps 
                 + 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) )
                           );
   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);
   double dE_mu2 = i_s.Qeps * i_s.Qeps;
   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);                        
   // 4-- variations du potentiel par rapport à cos3phi
   ret.Ecos3phi = dE_Qs * dQs_dcos3phi + dE_Qe * dQe_dcos3phi
                 + dE_mu2 * dMu2_dcos3phi + dE_mu3 * dMu3_dcos3phi;                        
   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 IsoHyper3DOrgeas2::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 Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance  à la phase et  température
   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 = 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;
 	   }
 	 };
   // on utilise la valeur 0 pour le cos3phi, pour avoir un comportement continu
   // au cas ou la valeur correspondant à 0 ne soit pas = à 1
     {if (Q0s_phi) Qs = Q0sT * Q0s_phi->Valeur(0.) ;
      if (Q0e_phi) Qe = Q0eT * Q0e_phi->Valeur(0.) ;
      if (mu01_phi) mu1 = mu01 * mu01_phi->Valeur(0.) ;
      if (mu02_phi) mu2 = mu02 * mu02_phi->Valeur(0.) ;
      if (mu03_phi) mu3 = mu03 * mu03_phi->Valeur(0.) ;
      };
   {// cas des fonctions nD éventuelles
    if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                          (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   };
   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 IsoHyper3DOrgeas2::PoGrenoblePhase_et_var
             (const Invariant2QepsCosphi& i_s,const Invariant & inv)
 {
//    bool affichage = (   ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 7))
//                     || (permet_affichage > 7)
//                    );

   // 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  à 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;
 	   }
 	 };
   // on modifie les variables qui dépendent de la phase
   // et on en calcul la dérivée première
   Courbe1D::ValDer2 f_para;
   double dQs_dcos3phi = 0.; double dQe_dcos3phi = 0.;
   double dMu1_dcos3phi = 0.;double dMu2_dcos3phi = 0.;double dMu3_dcos3phi = 0.;
   // 1-- calcul des variations seconde des grandeurs de bases   
   double dQs_dcos3phi2 = 0.; double dQe_dcos3phi2 = 0.;
   double dMu1_dcos3phi2 = 0.;double dMu2_dcos3phi2 = 0.;double dMu3_dcos3phi2 = 0.; 
   {if (Q0s_phi )
     {f_para = Q0s_phi->Valeur_Et_der12(i_s.cos3phi) ;
      Qs = Q0sT * f_para.valeur ;dQs_dcos3phi = Q0sT * f_para.derivee;
      dQs_dcos3phi2 = Q0sT * f_para.der_sec;
     };
    if (Q0e_phi )
     {f_para = Q0e_phi->Valeur_Et_der12(i_s.cos3phi) ;
      Qe = Q0e * f_para.valeur ;dQe_dcos3phi = Q0e * f_para.derivee;
      dQe_dcos3phi2 = Q0e * f_para.der_sec;
     };
    if (mu01_phi )
     {f_para = mu01_phi->Valeur_Et_der12(i_s.cos3phi) ;
      mu1 = mu01 * f_para.valeur ;dMu1_dcos3phi = mu01 * f_para.derivee;
      dMu1_dcos3phi2 = mu01 * f_para.der_sec;
     };
    if (mu02_phi )
     {f_para = mu02_phi->Valeur_Et_der12(i_s.cos3phi) ;
      mu2 = mu02 * f_para.valeur ;dMu2_dcos3phi = mu02 * f_para.derivee;
      dMu2_dcos3phi2 = mu02 * f_para.der_sec;
     };
    if (mu03_phi )
     {f_para = mu03_phi->Valeur_Et_der12(i_s.cos3phi) ;
      mu3 = mu03 * f_para.valeur ;dMu3_dcos3phi = mu03 * f_para.derivee;
      dMu3_dcos3phi2 = mu03 * f_para.der_sec;
     };
    };
  {// cas des fonctions nD éventuelles
   if (mu01_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu01_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (mu02_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu02_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (mu03_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (mu03_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (Q0s_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (Q0s_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
   if (Q0e_nD) mu1 *= (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                         (Q0e_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
  };

//cout << "\n cos(3phi)= " <<  i_s.cos3phi << " Qs= " << Qs << " mu2= " << mu2;     

   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 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  +++++++++++++++++++
   // a) la partie sans phase 
     
   ret.EV = K/3.*logV/inv.V;
   
// --- 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
   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);  // OK        
//   ret.EQ = Qs + 2.*mu2*i_s.Qeps
//            + mu1 *(QplusQec - S2) + mu3 * (i_s.Qeps - S4 + S3);

   // b) cas de la variation faisant intervenir  la phase
   
   // 2-- variations des grandeurs intermédiaires par rapport aux grandeurs de bases
   double Qs2 = Qs * Qs;
   double dQec_Qs = - (Qs2 + Sqr(mu1) * 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
   // 3-- variations du potentiel par rapport aux grandeurs de bases
   double dE_Qs = i_s.Qeps 
                 + 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 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_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);
   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
   // 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
   // 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
   
   // 4-- variations seconde du potentiel par rapport à cos3phi
   ret.Ecos3phi2 = dE_QsQs * dQs_dcos3phi * dQs_dcos3phi + dE_QeQe * dQe_dcos3phi * dQe_dcos3phi
                   + 2. * dE_Mu3Qe * dQe_dcos3phi * dMu3_dcos3phi
                   + dE_Qs * dQs_dcos3phi2 + dE_Qe * dQe_dcos3phi2 + dE_mu2 * dMu2_dcos3phi2
                   + dE_mu3 * dMu3_dcos3phi2; 
   // 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_QMu2*dMu2_dcos3phi + dE_QMu3*dMu3_dcos3phi;
   
    // 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), 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
//------ débug -------
/*   if (indic_Verif_CompaFormel==1)
	{ cout << "\n Qs= " << Qs << " dQs_dcos3phi= " << dQs_dcos3phi   << " dQs_dcos3phi2= " << dQs_dcos3phi2 
	       << "\n Qe= " << Qe << " dQe_dcos3phi= " << dQe_dcos3phi   << " dQe_dcos3phi2= " << dQe_dcos3phi2
	       << "\n mu2= " << mu2 << " dMu2_dcos3phi= " << dMu2_dcos3phi << " dMu2_dcos3phi2= " << dMu2_dcos3phi2
	       << "\n mu3= " << mu3 << " dMu3_dcos3phi= " << dMu3_dcos3phi << " dMu3_dcos3phi2= " << dMu3_dcos3phi2;

      cout << "\n\n dQec_Qs= " << dQec_Qs << " dQec_QsQs= " << dQec_QsQs 
           << "\n  dS1_Qs= " << dS1_Qs << " dS1_QsQs= " << dS1_QsQs
           << "\n  dS2_Qs= " << dS2_Qs << " dS2_QsQs= " << dS2_QsQs 
           << "\n  dS2_Q= " << dS2_Q << " dS2_QQ= " << dS2_QQ << " dS2_QQs= " << dS2_QQs
           << "\n  dS3_Qe= " << dS3_Qe << " dS3_QeQe= " << dS3_QeQe 
           << "\n  dS3_Q= " << dS3_Q << " dS3_QQ= " << dS3_QQ << " dS3_QQe= " << dS3_QQe
           << "\n  dS4_Qe= " << dS4_Qe << " dS4_QeQe= " << dS4_QeQe; 

      cout << "\n\n  dE_Qs= " << dE_Qs << " dE_Qe= " << dE_Qe 
           << "\n dE_mu2= " << dE_mu2 << " dE_mu3= " << dE_mu3  
           << "\n dE_QsQs= " << dE_QsQs << " dE_QeQe= " << dE_QeQe << " dE_Mu3Qe= " << dE_Mu3Qe;   
      cout << "\n dE_QQs= " << dE_QQs << " dE_QQe= " << dE_QQe << " dE_QMu2= " << dE_QMu2 << " dE_QMu3= " << dE_QMu3;
	};
    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 IsoHyper3DOrgeas2::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 
    // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_PoGrenobleorgeas2_et_var == 0)
     { cout << "\n ****verification des derivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DOrgeas2::Verif_PoGrenoble_et_var \n";
      } 
    indic_Verif_PoGrenobleorgeas2_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 IsoHyper3DOrgeas2::Verif_PoGrenoble_et_var(.."
            << " , numero d'increment = " << indic_Verif_PoGrenobleorgeas2_et_var;
       Sortie(1);
      } 
 };

    
    // vérif des dérivées du potentiels par rapport à l'invariant V, ceci par différences finies
void IsoHyper3DOrgeas2::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 
    // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_PoGrenobleorgeas2_et_var_VV == 0)
     { cout << "\n ****verification des derivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DOrgeas2::Verif_PoGrenoble_et_var \n";
      } 
    indic_Verif_PoGrenobleorgeas2_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 IsoHyper3DOrgeas2::Verif_PoGrenoble_et_var_VV(.."
            << " , numero d'increment = " << indic_Verif_PoGrenobleorgeas2_et_var_VV;
       Sortie(1);
      } 
 };

    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DOrgeas2::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_PoGrenobleorgeas2_et_var == 0)
     { cout << "\n ****vérification des dérivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DOrgeas2::Verif_PoGrenoble_et_var \n";
      } 
    indic_Verif_PoGrenobleorgeas2_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 IsoHyper3DOrgeas2::Verif_PoGrenoble_et_var(.."
            << " , numero d'increment = " << indic_Verif_PoGrenobleorgeas2_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 IsoHyper3DOrgeas2::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                        (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;};
   	  // appel du calcul des différentes valeurs
   	  PoGrenobleAvecPhaseAvecVar pot = PoGrenoblePhase_et_var(i_s,inv);
   	  // affichage
   	  cout << "\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;
   	};
		
	};