/* tr2.f -- translated by f2c cc*/
//OK fev 2006
/*fev 98  slop val*/ 

/* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */
/* Subroutine */
/* void tridiag3fr_ */

import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
import java.io.*;

 
 
 
public class tridiago    
 {
   static double diagi[];
   static double diagp[];
   static double diags[];
   static double rhs[];
   static double u[];
   static double uo[];
   static double a[];
   static double b[];
   static double Dy[];
   static double y[];
   static double nut[];
   static double nu0;
   static double longp[];
   static int n;
   static double dy;
   static double ry;
   static double dt;
   static double u0;
   static double un;
   static double px; 
   
   static double q;
   static double Bi;
   static boolean slope=false;
   static boolean mixte=false;
   static boolean val  =false;
   static boolean sym  =false;
   static double turb=1;


 
  public static void prems()
 {  int j;
   System.out.println(" initialize tridiago "); 
   diagi = new double[2*n+1];
   diagp = new double[2*n+1];
   diags = new double[2*n+1];
   rhs   = new double[2*n+1]; 
   u     = new double[2*n+1];
   uo    = new double[2*n+1];  
   a     = new double[2*n+1];
   b     = new double[2*n+1];
   Dy    = new double[2*n+1];
   y     = new double[2*n+1];
   nut   = new double[2*n+1];
   longp = new double[2*n+1];
   Dy[0]=dy;
   y[0]=0;
   y[2*n]=2;
   Dy[2*n]=dy;
   
   for (j=1; j<=n;j++)
   {Dy[j]=Dy[j-1]*ry;
	y[j]=y[j-1]+Dy[j-1];}	
	//y[n]=1;
	//Dy[n-1]=1-y[n-1];
	
   for (j=0; j<n;j++)	
	{y[2*n-j]=2-y[j];	}
   for (j=n; j<=2*n;j++)		
	{Dy[j]=y[j]-y[j-1];}
	
	System.out.println(" ymax "+y[n] + " " + y[2*n] +  " " + Dy[n-1]);


   }
   
   
   
public static void nuturb()
 { int j;
   double umy2;
   for (j=0; j<=n;j++)
   {/* Van Driest C p 363 */
   if(turb==0){longp[j]=0;}
   if(turb==1){longp[j]= 0.41 * y[j];}
   if(turb==2){longp[j]= 0.41 * y[j] *(1-Math.exp(-(y[j])/nu0/26));}
   
	umy2=(1-y[j])*(1-y[j]);   // formule de Nikuradse Schlichting 
	if(turb==3){longp[j]= 0.14 - 0.08 * umy2 - 0.06 * umy2*umy2;}
	//if(longp[j]>0.05){longp[j]=0.05;}
	}
	for (j=1; j<n;j++)
   {longp[2*n-j]=longp[j];}
	
	 //for (j=1; j<=2*n;j++) //
	 for (j=1; j<=n;j++)
   {
	 double lpp,lpm;       
	 //lpp =(longp[j]+longp[j+1])/2 ;
	//nut[j]=(lpp*lpp*(uo[j+1]-uo[j])/Dy[j]);
	 lpm =(longp[j]+longp[j-1])/2  ; 
	  nut[j-1]=(lpm*lpm*(uo[j]-uo[j-1])/Dy[j-1]);
	}
	
	for (j=1; j<2*n;j++){ nut[2*n-j]=nut[j-1]; }
     
	 nut[2*n]=0; 
	 
 }    

 public static void inverse()
 { int j, nn;
   double bb,cc,dd;
   double qn=0;
/*     inversion de matrice tridiago                              */
/*     ai ui-1 + bi ui +ci ui+1 = rhs i                          */
/*     diagi[j] u[j-1] +  diagp[j] u[j] + diags[j] u[j+1]= rhs[j]*/
/*    	u[j] = a[j] * u[j + 1] + b[j];                            */
/*     14 nov 90 -  23 nov 92 -  aout 95 - 8 mars 96 Aout 05  */
/*     Peyret Taylor p 20 */
/*     l indice varie de 0 a N */
	   
   nn=2*n;
   for (j=1; j<=nn;j++)
   {
   /*     (nu u')'   +  b u' + c u    =  d                          */ 
/********************************************************************/
		  
          bb = 0;				                       
		  cc =  - (1./dt) ; 
		  dd =  -(uo[j])/dt + px ;  
		 				  
	 diags[j] =    (nut[j]  + nu0)/(Dy[j]  *(Dy[j]+Dy[j-1])/2)         + bb/ (Dy[j]+Dy[j-1]);
	 diagp[j] =  -((nut[j]  + nu0)/(Dy[j]  *(Dy[j]+Dy[j-1])/2) 
				  +(nut[j-1]+ nu0)/(Dy[j-1]*(Dy[j]+Dy[j-1])/2))        + cc  ;
	 diagi[j] =    (nut[j-1]+ nu0)/(Dy[j-1]*(Dy[j]+Dy[j-1])/2)         - bb /(Dy[j]+Dy[j-1]);
	 rhs[j] =           dd  ;
		}
    j = 0;
    /*a[0] = 0.f;  b[0] = *u0;  donne'e de u0*/
    /*a[0] = 1.f;  b[0] = 0; modif du/dy=0 en 0*/
    /*a[0] = 1.f;  b[0] = dy;du/dy=1 en 0*/
     
  if(slope){  
   a[0] = 1.;  b[0] = (q)*Dy[0];
   }
  if(val){
   a[0] = 0.;  b[0] = u0; }
   
  if(mixte){  double Te=0;
   a[0]= (Bi/Dy[0])/(1+Bi/Dy[0]); b[0] = Te/(1+Bi/Dy[0]);}
 
    for (j = 1; j < nn; ++j) {
	dd = 1./(diagp[j] + diagi[j] * a[j - 1]);
	a[j] = -diags[j] * dd ;
	b[j] = (rhs[j] - diagi[j] * b[j - 1]) * dd ;
    }

if(sym)
    {
	 u[nn-1] = (rhs[nn-1]- b[nn-2]*diagi[nn-1] - diags[nn-1]*qn*Dy[nn-1])/(diagp[nn-1] + diagi[nn-1] * a[nn - 2] + diags[nn-1]);
	 u[nn] = u[nn-1] + qn*Dy[nn-1];
	 for (j = nn - 2; j >= 0; --j) {   
	u[j] = a[j] * u[j + 1] + b[j];
    }
	}
else
	{
	u[nn] = un;
    for (j = nn - 1; j >= 0; --j) {   
	u[j] = a[j] * u[j + 1] + b[j];
    }
    }
 /* tridiago2_ */
 
  }
}//