//12 decembre 2003
// novembre 2004
#include <stdio.h>
#include <stdlib.h>
#include "math.h"
 
/* calcul de l'expression du frottement DANS L'ESPACE DE FOURIER en fonction de la forme de la bosse DANS L'ESPACE DE FOURIER et du regime  */
/*
  TF bosse  ,TF pente bosse, TF pression, TF tau       , nbre , taille ,regime */
void hydro(
double TFf[],double TFfp[],double TFp[],double TFtau[],int nx,double Lx,int reg)
{ extern double Pi;
  static int i;
  static double a,a13,*Tzr,*Tzi;
  
Tzr =(double*)calloc(nx+2,sizeof(double)); 
Tzi =(double*)calloc(nx+2,sizeof(double)); 
  
  i=0;
  TFp[2*i+1]  = 0;
  TFp[2*i+2]  = 0;
  TFtau[2*i+1]= 0;
  TFtau[2*i+2]= 0;
  TFtau[2*(nx-i)+1]= TFtau[2*i+1];
  TFtau[2*(nx-i)+2]=-TFtau[2*i+2];
  TFp[2*(nx-i)+1]  = TFp[2*i+1];
  TFp[2*(nx-i)+2]  =-TFp[2*i+2];
  
  
  
  //valcoeffzrzi3.nb
  ///Volumes/MacintoshHD/DOKUMENTS/Documents/2003/helec/graphs 
  //
  /*
   <<"Algebra`ReIm`"

ap=N[AiryAiPrime[0]];
  ap3=3 ap;
  a0=N[AiryAi[0]];
  Bsta1=  N[(-I )^(1/3)/ap3]   a13;
  a13/:Im[a13]=0;
  Im[a]^=0;
  a23/:Im[a23]=0;
  taufDIVpf=a13^2  N[(-I )^(2/3) (a0/ap)  ]   ;
  {CForm[ Simplify[Re[1/Bsta1 ]]],
     CForm[Simplify[Im[1/Bsta1 ]]]}
     {CForm[ Simplify[Re[taufDIVpf]]],
        CForm[Simplify[Im[taufDIVpf]]]}
*/
   //cas tuyau //-0.672433 - 0.388229 \[ImaginaryI]
if (reg==0){
 for (i=1; i<=nx/2;i++)
 {  
 a=(i)*(Pi/Lx);
 a13=pow(a,1./3);
 Tzr[i]=(-0.672433)/a13;
 Tzi[i]=(-0.388229)/a13;}}
 
   //cas subsonique, transformee de Hilbert
   /*
   B = 1/a;
   pf =1/(Bsta1- B) ;
   {CForm[ Simplify[Re[pf]]],
   CForm[Simplify[Im[pf]]]}
   */
if (reg==2){
 for (i=1; i<=nx/2;i++)
 {  
 a=(i)*(Pi/Lx);
 a13=pow(a,1./3);
 Tzr[i]=(a*(-1 - 1.1153535259122473*a*a13))/
   (1 + 2.2307070518244947*a*a13 + 
     1.6586846503531767*pow(a,2)*pow(a13,2));
 Tzi[i]=(-0.6439496584270346*pow(a,2)*a13)/
   (1 + 2.2307070518244947*a*a13 + 
     1.6586846503531767*pow(a,2)*pow(a13,2));}}
     
    //cas subcritique
    /*
     B = 1;
     pf =1/(Bsta1- B) ;
    */
if (reg==-1){
 for (i=1; i<=nx/2;i++)
 {  
 a=(i)*(Pi/Lx);
 a13=pow(a,1./3);
  Tzr[i]=(1 + 1.1153535259122473*a13)/
   (-1 - 2.2307070518244947*a13 - 
     1.6586846503531767*pow(a13,2));
  Tzi[i]=(-0.6439496584270346*a13)/
   (1 + 2.2307070518244947*a13 + 
     1.6586846503531767*pow(a13,2)) ;}}

    //cas supercritique 
if (reg==1){
 for (i=1; i<=nx/2;i++)
 {  
 a=(i)*(Pi/Lx);
 a13=pow(a,1./3);
 Tzr[i]=(1 - 1.1153535259122473*a13)/
   (1 - 2.2307070518244947*a13 + 
     1.6586846503531767*pow(a13,2));
  Tzi[i]=(-0.6439496584270346*a13)/
   (1 - 2.2307070518244947*a13 + 
     1.6586846503531767*pow(a13,2));}}


// cas tuyau bis    
if (reg==3){
for (i=1; i<=nx/2;i++)
{  
a=(i)*(Pi/Lx);
a13=pow(a,1./3);
Tzr[i]=(0.8108021448049955 + 0.11052934823644561*pow(a, 2)*a13 + 
     0.005977178098051019*pow(a, 4)*pow(a13, 2) + 
     0.00014814814814814812*pow(a, 6)*pow(a13, 3) + 
     1.3769796616200583*.000001 
     *pow(a, 8)*pow(a13, 4))/(-2.411549416067905*a13 - 
     0.2689730144122959*pow(a, 2)*pow(a13, 2) - 
     0.012222222222222221*pow(a, 4)*pow(a13, 3) - 
     0.0002478563390916105*pow(a, 6)*pow(a13, 4) - 2.0477588275965144*.000001*pow(a, 8)*pow(a13, 5));
 Tzi[i]=(0.46811683656269026 + 0.05221157641990838*pow(a, 2)*a13 + 
     0.0025881940379280694*pow(a, 4)*pow(a13, 2) + 
     0.00006415002990995844*pow(a, 6)*pow(a13, 3) + 7.94999578304981*.0000001*pow(a, 8)*pow(a13, 4))/(-2.411549416067905*a13 - 
     0.2689730144122959*pow(a, 2)*pow(a13, 2) - 
     0.012222222222222221*pow(a, 4)*pow(a13, 3) - 
     0.0002478563390916105*pow(a, 6)*pow(a13, 4) - 2.0477588275965144*.000001*pow(a, 8)*pow(a13, 5));}}             

     
 //super sonique     
 /*
  B = 1/((I a ));
  pf = - 1/(B-Bsta1);  {CForm[\ Simplify[Re[pf]]],   CForm[Simplify[Im[pf]]]}
 */ 
if (reg==10){
 for (i=1; i<=nx/2;i++)
 {  
 a=(i)*(Pi/Lx);
 a13=pow(a,1./3);
 Tzr[i]=(-1.1153535259122473*pow(a, 2)*a13)/(1 + 1.287899316854069*a*a13 + 
1.6586846503531767*pow(a, 2)*pow(a13, 2));
  Tzi[i]=(a*(-1 - 0.6439496584270346*a*a13))/(1 + 1.287899316854069*
            a*a13 + 1.6586846503531767*pow(a, 2)*pow(a13, 2));}}
          
/* fin de la premiere etape dependant du regime */

 for (i=1; i<=nx/2;i++)
 {  
//pression a partir de la bosse et de Tz
  TFp[2*i+1]= ( Tzr[i] * TFf[2*i+1] - Tzi[i] * TFf[2*i+2]) ;
  TFp[2*i+2]= ( Tzr[i] * TFf[2*i+2] + Tzi[i] * TFf[2*i+1]) ;
  
// frottement
//-0.685861 + 1.18795 \[ImaginaryI]
  TFtau[2*i+1]= ((-0.685861) * TFp[2*i+1] - ( 1.18795) * TFp[2*i+2])*pow(i*(Pi/Lx) ,2./3);
  TFtau[2*i+2]= ((-0.685861) * TFp[2*i+2] + ( 1.18795) * TFp[2*i+1])*pow(i*(Pi/Lx) ,2./3);
  
//   TFtau[2*i+1]=  TFf[2*i+2]* i*(Pi/Lx);
//   TFtau[2*i+2]= -TFf[2*i+1]* i*(Pi/Lx);
	   
	   
	   
/* derivee de yw: pente de la bosse TFp */
  
    TFfp[2*i+1]= TFf[2*i+2]* i*(Pi/Lx);
    TFfp[2*i+2]=-TFf[2*i+1]* i*(Pi/Lx);
	   
	   
	   
	   
 /* autres modes -> conjugue! car on travaille sur une fction Reelle  */
   TFtau[2*(nx-i)+1]= TFtau[2*i+1];
   TFtau[2*(nx-i)+2]=-TFtau[2*i+2];
   
   TFp[2*(nx-i)+1]= TFp[2*i+1];
   TFp[2*(nx-i)+2]=-TFp[2*i+2];
	

   TFfp[2*(nx-i)+1]= TFfp[2*i+1];	   
   TFfp[2*(nx-i)+2]=-TFfp[2*i+2];
  
  } 
  
/*=========================================CAS PARTICULIERS====================*/ 
/* CAS EXNER */
 if (reg==100){
 for (i=1; i<=nx/2;i++)
 {  
 a=(i)*(Pi/Lx);
 Tzr[i]= - a;
 Tzi[i]= 0;
  TFp[2*i+1]= ( Tzr[i] * TFf[2*i+1] - Tzi[i] * TFf[2*i+2]) ;
  TFp[2*i+2]= ( Tzr[i] * TFf[2*i+2] + Tzi[i] * TFf[2*i+1]) ;
//-1   cas exner pur du CRAS tau = - u 
  TFtau[2*i+1]= -TFp[2*i+1] ;
  TFtau[2*i+2]= -TFp[2*i+2] ;
 /* derivee de yw */
    TFfp[2*i+1]= TFf[2*i+2]* i*(Pi/Lx);
    TFfp[2*i+2]=-TFf[2*i+1]* i*(Pi/Lx);
 /* autres modes -> conjugue! car on travaille sur une fction Reelle  */
   TFtau[2*(nx-i)+1]= TFtau[2*i+1];
   TFtau[2*(nx-i)+2]=-TFtau[2*i+2];
   //
   TFp[2*(nx-i)+1]= TFp[2*i+1];
   TFp[2*(nx-i)+2]=-TFp[2*i+2];
   //
   TFfp[2*(nx-i)+1]= TFfp[2*i+1];	   
   TFfp[2*(nx-i)+2]=-TFfp[2*i+2];   
   }}	
   
   
/* CAS cas Hele Chaud juin 2003   */
 
   if (reg==3){
   for (i=1; i<=nx/2;i++)
   {  
   a=(i)*(Pi/Lx);
   Tzr[i]= -6*fabs(a)/5.;
   Tzi[i]=  - fabs(a)/a;
	TFp[2*i+1]= ( Tzr[i] * TFf[2*i+1] - Tzi[i] * TFf[2*i+2]) ;
	TFp[2*i+2]= ( Tzr[i] * TFf[2*i+2] + Tzi[i] * TFf[2*i+1]) ;
  //-1 
  //
	Tzr[i]=  (100 + 45*a*a)*fabs(a) /(100 + 9*a*a);
	Tzi[i]=     -120*a*fabs(a)      /(100 + 9*a*a);
	TFtau[2*i+1]= ( Tzr[i] * TFf[2*i+1] - Tzi[i] * TFf[2*i+2]) ;
	TFtau[2*i+2]= ( Tzr[i] * TFf[2*i+2] + Tzi[i] * TFf[2*i+1]) ;
		 /* derivee de yw */
	  TFfp[2*i+1]= TFf[2*i+2]* i*(Pi/Lx);
	  TFfp[2*i+2]=-TFf[2*i+1]* i*(Pi/Lx);
	/* autres modes -> conjugue! car on travaille sur une fction Reelle  */
	 TFtau[2*(nx-i)+1]= TFtau[2*i+1];
	 TFtau[2*(nx-i)+2]=-TFtau[2*i+2];
	 //
	 TFp[2*(nx-i)+1]= TFp[2*i+1];
	 TFp[2*(nx-i)+2]=-TFp[2*i+2];
	 //
	 TFfp[2*(nx-i)+1]= TFfp[2*i+1];	   
	 TFfp[2*(nx-i)+2]=-TFfp[2*i+2];   
	 }}	
	   
   
   
   
   
/* CAS ALLER RETOUR SINUS */
 if (reg==111){
 for (i=1; i<=nx/2;i++)
 {  
 a= (i)*(Pi/Lx) ;
 Tzr[i]= 1;
 Tzi[i]= 0;
  TFp[2*i+1]= ( Tzr[i] * TFf[2*i+1] - Tzi[i] * TFf[2*i+2]) ;
  TFp[2*i+2]= ( Tzr[i] * TFf[2*i+2] + Tzi[i] * TFf[2*i+1]) ; 
// 
  TFtau[2*i+1]= TFf[2*i+2]  ;
  TFtau[2*i+2]= TFf[2*i+1] ;
	   /* derivee de yw */
    TFfp[2*i+1]= TFf[2*i+2]* i*(Pi/Lx);
    TFfp[2*i+2]=-TFf[2*i+1]* i*(Pi/Lx);
  /* autres modes -> conjugue! car on travaille sur une fction Reelle  */
   TFtau[2*(nx-i)+1]= TFtau[2*i+1];
   TFtau[2*(nx-i)+2]=-TFtau[2*i+2];
   //
   TFp[2*(nx-i)+1]= TFp[2*i+1];
   TFp[2*(nx-i)+2]=-TFp[2*i+2];
   //
   TFfp[2*(nx-i)+1]= TFfp[2*i+1];	   
   TFfp[2*(nx-i)+2]=-TFfp[2*i+2];   
   }}	

/* ----------------------------------------------------------*/ 
/* liberez les tableaux !!                                   */
        free(Tzr);
        free(Tzi);

}