/* copyright ECMWF, Nils Peter Wedi, May 2001 ECMWF training course */

/* computer program which enables the user to play with the parameters
   of the various numerical schemes examined during the seminar. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "svd.h"

#define PI 3.1415927

#define gamma 1.0 /* implicitness parameter */
#define dt 0.166  /* time step */
#define IMAX 290 /* number of iterations , e.g. t=i*dt, i<IMAX */

/* for a description of the parameters see the text */

/* parameters for stiffness */
#define Kstiff 100.0
/* #define IMAX 20 */

/* parameters for boundaries */
#define NMAX 16
#define dz 100.0
#define T00 16.0
#define T0z 0.0
#define kappa 20.0
/* #define IMAX 36*5 */

/* parameters for non-linear terms */
#define Knonl 10
#define Pnonl 3


/* ************************************ */

double exact_stiff( double );
void stiffness();

void time_step_boundaries(double [NMAX], double [NMAX]);
void boundaries();

void non_linear();
void non_linear_predictor_corrector();

/* ************************************ */

int main(int argc,char *argv[]) {


  /* stiffness */
  /* stiffness(); */

  /* boundaries */
  /* boundaries(); */

  /* non linear terms */
  /* non_linear(); */
  non_linear_predictor_corrector();

return(0);

}

/* ************************ */

double exact_stiff(double t) {

  /* u(0)==0 */

  return(-0.5*exp(-Kstiff*t)+t/10.0+0.5);

}

/* ************************ */

void stiffness() {

  int i=0;
  double t;
  double T0=0.0,T1=0.0;


  /* stiffness */
  printf("gamma %f K %f dt %f\n",gamma,Kstiff,dt);
  printf("t %f exact %f num %f\n",0.0,exact_stiff(0.0),0.0);
  for (i=1;i<IMAX;i++) {
    t=i*dt;
    T1=T0 - Kstiff*dt*(1.-gamma)*T0+Kstiff*dt*(t/10.+0.5)+dt*1./10.;
    T1*=1./(1.+Kstiff*dt*gamma);
    printf("t %f exact %f num %f\n",t,exact_stiff(t),T1);
    T0=T1;
  }

}

/* ************************ */

void time_step_boundaries(double T9[NMAX], double T1[NMAX]) {

/* (C) Copr. 1986-92 Numerical Recipes Software #1--3Y".*/
/* ************************ */
  void lubksb(double **a, int n, int *indx, double b[]);
  void ludcmp(double **a, int n, int *indx, double *d);
/* ************************ */

  /* initialize */

  double **A,*b,d,D;
  int *indx,i,k;

  indx=ivector(1,NMAX);
  b=dvector(1,NMAX);
  A=dmatrix(1,NMAX,1,NMAX);

  for( i=1;i<NMAX-1;i++ )
    b[i+1]=T9[i]+dt*(1-gamma)*kappa*(T9[i+1]-2*T9[i]+T9[i-1])/(dz*dz);
  b[1]=T0z;
  b[NMAX]=T0z;
  
  k=0;
  D=kappa*gamma*dt/(dz*dz);
  for(i=1;i<NMAX-1;i++) {
    A[i+1][k+1]=-D;
    A[i+1][k+2]=1+2*D;
    A[i+1][k+3]=-D;
    k++;
  }
  A[1][1]=1.0;
  A[16][16]=1.0;
  for(k=2;k<NMAX;k++) {
    A[1][k]=0.0;
    A[16][k]=0.0;
  }
    
  /* solve A*T = b */

  ludcmp(A,NMAX,indx,&d);
  lubksb(A,NMAX,indx,b);

  /* solution is in b , A is destroyed */

  for(i=0;i<NMAX;i++)
    T1[i]=b[i+1];

}

/* ************************ */

void boundaries() {

  int i,k;
  double T9[NMAX],T1[NMAX];
  double t;

  /* initialize */

  for( i=1;i<NMAX-1;i++ )
    T9[i]=T00;
  T9[0]=0.0;
  T9[NMAX-1]=0.0;

  printf("z %f kappa %f gamma %f dt %f cfl %f\n",
	 100.0,kappa,gamma,dt,kappa*dt/(dz*dz));

  printf("t %f T %f\n",0.0,T00);
  for (i=1;i<IMAX;i++) {
    t=i*dt;
    time_step_boundaries(T9,T1);
    for(k=1;k<NMAX-1;k++)
      T9[k]=T1[k];
    printf("t %f T %f\n",t,T1[1]);
  }

}

/* ************************ */

void non_linear() {

  int i=0;
  double t;
  double T0=0.5,T1=0.0;

  /* non linear terms */
  printf("gamma %f K %d P %d dt %f\n",gamma,Knonl,Pnonl,dt);

  printf("t %f T %f\n",0.0,T0);
  for (i=1;i<IMAX;i++) {
    t=i*dt;
    T1=T0*(1.- Knonl*dt*(1.-gamma)*pow(T0,Pnonl))+dt-dt*sin(2.*PI*t/24.);
    T1*=1./(1.+Knonl*dt*gamma*pow(T0,Pnonl));
    printf("t %f T %f\n",t,T1);
    T0=T1;
  }

}

void non_linear_predictor_corrector() {

  int i=0;
  double t;
  double T0=0.5,T1=0.0,T11=0.0;

  /* gamma should be one for the example of McDonald */

  /* non linear terms */
  printf("gamma %f K %d P %d dt %f\n",gamma,Knonl,Pnonl,dt);

  printf("t %f T %f\n",0.0,T0);
  for (i=1;i<IMAX;i++) {
    t=i*dt;
    T11=T0*(1.- Knonl*dt*(1.-gamma)*pow(T0,Pnonl))+dt-dt*sin(2.*PI*t/24.);
    T11*=1./(1.+Knonl*dt*gamma*pow(T0,Pnonl));

    T1=T0*(1.- Knonl*dt*(1.-gamma)*pow(T11,Pnonl))+dt-dt*sin(2.*PI*t/24.);
    T1*=1./(1.+Knonl*dt*gamma*pow(T11,Pnonl));

    printf("t %f T %f\n",t,T1);
    T0=T1;
  }

}