// kzk.cpp
// This version: 16 February 2001
//
// This file can be freely distributed, however, all header
// information must be retained.
//
// Robin Cleveland
// Department of Aerospace and Mechanical Engineering
// Boston University
// Boston, MA 02215
// robinc@bu.edu
//
// Mark Hamilton
// Department of Mechanical Engineering
// The University of Texas at Austin
// Austin, TX
// hamilton@mail.utexas.edu
//
// Sebastien Manneville
// Department of Aerospace and Mechanical Engineering
// Boston University
// Boston, MA 02215
// sebm@bu.edu 
// sebastien.manneville@espci.fr

// This program is a C++ version of the fortran KZKtexas2.f code.
// The fortran source code for a time-domain computer code was developed at
// The University of Texas at Austin to model axisymmetric focused
// sound beams in fluids. The code is based on an augmented KZK equation
// that accounts for nonlinearity, diffraction, thermoviscous absorption,
// and absorption and dispersion due to an arbitrary number of relaxation
// phenomena. The principal references for this code are the following:
//
// UNFOCUSED CODE:
//   REFER to: "Time-domain modeling of pulsed finite-amplitude
//   sound beams", Yang-Sub Lee and Mark F. Hamilton, J. Acoust.
//   Soc. Am. 97, 906-917 (1995).
//
//   Modified by Robin Cleveland, Applied Physics Laboratory
//   University of Washington at Seattle, January 1996
//
//   REFER to: "Time-domain modeling of finite-amplitude waves in relaxing fluids",
//   Robin O. Cleveland, Mark F. Hamilton, and David T. Blackstock,
//   J. Acoust. Soc. Am. 99, 3312-3318 (1996).
//
// FOCUSING CODE. Modified as per Lee and Hamilton paper to solve
//   the KZK equation on a "flat" grid which allows the code to be
//   used for moderately focused fields.
//
//   REFER to: "Nonlinear distortion of short pulses radiated by plane
//   and focused circular pistons", Michalakis A. Averkiou and Mark F. Hamilton, 
//   J. Acoust. Soc. Am. 102, 2539-2548 (1997).
//
//

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <malloc.h>
#include <math.h>
#include <process.h>
#include "declar.h"

#define	PI		3.14159265358979

extern int InputConfig(LPCONFIG lpConfig, char *sConfigFileName);
extern int Init(LPCONFIG lpConfig, double **P);
extern void IBFDinitMat(LPCONFIG lpConfig, LPMAT_1 lpMatriceAbsorb, LPMAT_1 lpMatriceDiffract, LPMAT_2 lpMatriceRelax);
extern void CNFDinitMat(LPCONFIG lpConfig, LPMAT_1 lpMatriceAbsorb, LPMAT_1 lpMatriceDiffract, LPMAT_2 lpMatriceRelax);
extern int NonLinear(LPCONFIG lpConfig, double **P, double *Pold, double *TauDisto);
extern void Tridiag(LPCONFIG lpConfig, LPMAT_1 lpMatrice, double *Rhs, double *Sol, double *Beta, double *Gamma);


//
// MAIN PROGRAM : KZK CODE FOR AXISYMMETRIC FOCUSED BEAMS
//


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

{
LPCONFIG lpConfig;
LPMAT_1 lpMatriceAbsorb, lpMatriceDiffract, lpMatriceRelaxTmp;
LPMAT_2 lpMatriceRelax;
char sConfigFileName[100];
int  i, j, jTmp, n1, n2, iTargetID, iResInit, iResNL;
double dTarget, dDeltaSigmaOld, dRhoTmp, dDummy;
double *RhsRho, *SumRho, *BetaRho, *GammaRho, *RhsTau, *SolRho, *SolTau, *BetaTau, *GammaTau, *Pold, *TauDisto;
double **P;
FILE *fp;

printf("KZKpulse program\n");
printf("________________\n\n");

if (argc != 2) {
	printf("\nPE.exe requires 1 input parameter !\n");
	printf("\nSimulation aborted !\n");
	return(-1);
}
strcpy(sConfigFileName,argv[1]);


//
// Initializations and memory allocations
//


lpConfig			= (LPCONFIG)  malloc(sizeof(CONFIG));
lpMatriceAbsorb		= (LPMAT_1)malloc(sizeof(MAT_1));
lpMatriceDiffract	= (LPMAT_1)malloc(sizeof(MAT_1));
lpMatriceRelax		= (LPMAT_2)malloc(sizeof(MAT_2));
lpMatriceRelaxTmp	= (LPMAT_1)malloc(sizeof(MAT_1));

iResInit=InputConfig(lpConfig, sConfigFileName);
if ( iResInit == -1) {
	printf("\nSimulation aborted !\n");
	return(-1);
}

P = (double **) malloc((lpConfig->NbTau+1)*sizeof(double *));
for (i=0; i<=lpConfig->NbTau; i++)
	P[i]  = (double *) malloc((lpConfig->NbRho+1)*sizeof(double));

RhsRho = (double *) malloc((lpConfig->NbRho+1)*sizeof(double));
SumRho = (double *) malloc((lpConfig->NbRho+1)*sizeof(double));
SolRho = (double *) malloc((lpConfig->NbRho+1)*sizeof(double));
BetaRho = (double *) malloc((lpConfig->NbRho+1)*sizeof(double));
GammaRho = (double *) malloc((lpConfig->NbRho+1)*sizeof(double));
RhsTau = (double *) malloc((lpConfig->NbTau+1)*sizeof(double));
SolTau = (double *) malloc((lpConfig->NbTau+1)*sizeof(double));
BetaTau = (double *) malloc((lpConfig->NbTau+1)*sizeof(double));
GammaTau = (double *) malloc((lpConfig->NbTau+1)*sizeof(double));
Pold = (double *) malloc((lpConfig->NbTau+1)*sizeof(double));
TauDisto = (double *) malloc((lpConfig->NbTau+1)*sizeof(double));

iResInit=Init(lpConfig, P);
if ( iResInit == -1) {
	printf("\nSimulation aborted !\n");
	return(-1);
}

if ( lpConfig->ThermoViscous ) {
	lpMatriceAbsorb->Size=lpConfig->NbTau+1;
	lpMatriceAbsorb->a=(double *)malloc((lpMatriceAbsorb->Size)*sizeof(double));
	lpMatriceAbsorb->b=(double *)malloc((lpMatriceAbsorb->Size)*sizeof(double));
	lpMatriceAbsorb->c=(double *)malloc((lpMatriceAbsorb->Size)*sizeof(double));
}

if ( lpConfig->Diffraction ) {
	lpMatriceDiffract->Size=lpConfig->NbRho+1;
	lpMatriceDiffract->a=(double *)malloc((lpMatriceDiffract->Size)*sizeof(double));
	lpMatriceDiffract->b=(double *)malloc((lpMatriceDiffract->Size)*sizeof(double));
	lpMatriceDiffract->c=(double *)malloc((lpMatriceDiffract->Size)*sizeof(double));
}

if ( lpConfig->Relaxation ) {
	lpMatriceRelax->NbMat=lpConfig->NbRelax;
	lpMatriceRelax->Size=lpConfig->NbTau+1;
	lpMatriceRelax->Lambda = (double *) malloc((lpMatriceRelax->NbMat)*sizeof(double));
	lpMatriceRelax->Mu = (double *) malloc((lpMatriceRelax->NbMat)*sizeof(double));
	lpMatriceRelax->a = (double **) malloc((lpMatriceRelax->NbMat)*sizeof(double *));
	lpMatriceRelax->b = (double **) malloc((lpMatriceRelax->NbMat)*sizeof(double *));
	lpMatriceRelax->c = (double **) malloc((lpMatriceRelax->NbMat)*sizeof(double *));
	for (j=0; j<lpMatriceRelax->NbMat; j++) {
		lpMatriceRelax->a[j]=(double *)malloc(lpMatriceRelax->Size*sizeof(double));
		lpMatriceRelax->b[j]=(double *)malloc(lpMatriceRelax->Size*sizeof(double));
		lpMatriceRelax->c[j]=(double *)malloc(lpMatriceRelax->Size*sizeof(double));
	}

}

fp=fopen(lpConfig->OutputFileName,"wb");
if( fp == NULL ) {
    printf("\nCould not open output file %s\n", lpConfig->OutputFileName);
	printf("\nSimulation aborted !\n");
	return(-1);
}
fwrite(&(lpConfig->NbTau), sizeof(int), 1, fp);
fwrite(lpConfig->Tau, sizeof(double), lpConfig->NbTau+1, fp);
fwrite(&(lpConfig->NbWatch), sizeof(int), 1, fp);

//
// IBFD iterations
//


printf("\n   Beginning IBFD iterations...\n");
iTargetID=0;
dTarget=lpConfig->SigmaWatch[0];
lpConfig->Sigma=0;
lpConfig->DeltaSigma=__min(lpConfig->DeltaSigmaIB,lpConfig->DeltaSigmaShock/2.0);
IBFDinitMat(lpConfig, lpMatriceAbsorb, lpMatriceDiffract, lpMatriceRelax);

for (n1=0; n1<lpConfig->NbStepsIB; n1++) {

	dDeltaSigmaOld=lpConfig->DeltaSigma;
	lpConfig->DeltaSigma=__min(lpConfig->DeltaSigmaIB,lpConfig->DeltaSigmaShock/2.0);
	if ( dDeltaSigmaOld != lpConfig->DeltaSigma) {
		IBFDinitMat(lpConfig, lpMatriceAbsorb, lpMatriceDiffract, lpMatriceRelax);
		printf("         DeltaSigma changed to %lf\n",lpConfig->DeltaSigma);
	}
	lpConfig->Sigma += lpConfig->DeltaSigma;
	dDummy = 100.0 * (lpConfig->Sigma / lpConfig->SigmaMax);
	printf("      Current Sigma = %lf\t\t Done %3.2lf %%\n",lpConfig->Sigma, dDummy);

	if ( lpConfig->Diffraction ) {		// IBFD DIFFRACTION
		for (j=lpMatriceDiffract->Start; j<=lpMatriceDiffract->Finish; j++)
			SumRho[j]=0.0;
		for (i=1; i<=lpConfig->NbTau-1; i++) {
			for (j=lpMatriceDiffract->Start; j<=lpMatriceDiffract->Finish; j++)
				SumRho[j] += P[i-1][j];
			RhsRho[0]=P[i][0] + lpMatriceDiffract->Lambda*(SumRho[1]-SumRho[0]);
			for (j=1; j<=lpConfig->NbRho-2; j++)
				RhsRho[j]=P[i][j] + (1.0 - 1.0/(2.0*j)) * lpMatriceDiffract->Lambda * SumRho[j-1]/4.0
								- lpMatriceDiffract->Lambda * SumRho[j]/2.0 +
								(1.0 + 1.0/(2.0*j)) * lpMatriceDiffract->Lambda * SumRho[j+1]/4.0;
			RhsRho[lpConfig->NbRho-1]=P[i][lpConfig->NbRho-1] +
								(1.0 - 1.0/(2.0*(lpConfig->NbRho-1))) * lpMatriceDiffract->Lambda * SumRho[lpConfig->NbRho-2]/4.0
								- lpMatriceDiffract->Lambda * SumRho[lpConfig->NbRho-1]/2.0;
			Tridiag(lpConfig, lpMatriceDiffract, RhsRho, SolRho, BetaRho, GammaRho);
			for (j=lpMatriceDiffract->Start; j<=lpMatriceDiffract->Finish; j++)
				P[i][j]=SolRho[j];
		}
	}

	if ( lpConfig->ThermoViscous ) {	// IBFD ABSORPTION
		for (j=0; j<=lpConfig->NbRho-1; j++) {
			for (i=lpMatriceAbsorb->Start; i<=lpMatriceAbsorb->Finish; i++)
				RhsTau[i]=P[i][j];
			Tridiag(lpConfig, lpMatriceAbsorb, RhsTau, SolTau, BetaTau, GammaTau);
			for (i=lpMatriceAbsorb->Start; i<=lpMatriceAbsorb->Finish; i++)
				P[i][j]=SolTau[i];
		}
	}
	
	if ( lpConfig->Relaxation ) {		// IBFD RELAXATION
		for (n2=0; n2<lpConfig->NbRelax; n2++) {
			lpMatriceRelaxTmp->a=lpMatriceRelax->a[n2];
			lpMatriceRelaxTmp->b=lpMatriceRelax->b[n2];
			lpMatriceRelaxTmp->c=lpMatriceRelax->c[n2];
			lpMatriceRelaxTmp->Size=lpMatriceRelax->Size;
			lpMatriceRelaxTmp->Start=lpMatriceRelax->Start;
			lpMatriceRelaxTmp->Finish=lpMatriceRelax->Finish;
			for (j=0; j<=lpConfig->NbRho-1; j++) {
				for (i=lpMatriceRelaxTmp->Start; i<=lpMatriceRelaxTmp->Finish; i++)
					RhsTau[i]=-lpMatriceRelax->Mu[n2]*P[i-1][j] + P[i][j] +
								lpMatriceRelax->Mu[n2]*P[i+1][j];				
				Tridiag(lpConfig, lpMatriceRelaxTmp, RhsTau, SolTau, BetaTau, GammaTau);
				for (i=lpMatriceRelaxTmp->Start; i<=lpMatriceRelaxTmp->Finish; i++)
					P[i][j]=SolTau[i];
			}
		}
	}

	if ( lpConfig->NonLinearity ) {
		iResNL=NonLinear(lpConfig, P, Pold, TauDisto);
		if ( iResNL == -1) {
			printf("\nSimulation aborted !\n");
			return(-1);
		}
	}

	if ( lpConfig->Sigma >= dTarget ) {
		printf("         Writing data to output file...\n");
		while ( lpConfig->SigmaWatch[iTargetID] == dTarget ) {
			fwrite(&(lpConfig->Sigma), sizeof(double), 1, fp);
			jTmp = (int) floor( lpConfig->RhoWatch[iTargetID] * lpConfig->NbPiston );
			dRhoTmp = jTmp * lpConfig->DeltaRho;
			fwrite(&(dRhoTmp), sizeof(double), 1, fp);
			for (i=0; i<=lpConfig->NbTau; i++)
				fwrite(&(P[i][jTmp]), sizeof(double), 1, fp);
			printf("            Sigma = %lf Rho = %lf\n",lpConfig->Sigma,dRhoTmp);
			iTargetID++;
		}
		dTarget=lpConfig->SigmaWatch[iTargetID];
	}
}

printf("   Ended IBFD iterations at Sigma = %lf\n", lpConfig->Sigma);


//
// CNFD iterations
//


printf("\n   Beginning CNFD iterations...\n");
lpConfig->DeltaSigma=__min(lpConfig->DeltaSigmaCN,lpConfig->DeltaSigmaShock/2.0);
if ( dDeltaSigmaOld != lpConfig->DeltaSigma) {
	printf("         DeltaSigma changed to %lf\n",lpConfig->DeltaSigma);
}
CNFDinitMat(lpConfig, lpMatriceAbsorb, lpMatriceDiffract, lpMatriceRelax);

while ( lpConfig->Sigma <= lpConfig->SigmaMax ) {

	dDeltaSigmaOld=lpConfig->DeltaSigma;
	lpConfig->DeltaSigma=__min(lpConfig->DeltaSigmaCN,lpConfig->DeltaSigmaShock/2.0);
	if ( dDeltaSigmaOld != lpConfig->DeltaSigma) {
		CNFDinitMat(lpConfig, lpMatriceAbsorb, lpMatriceDiffract, lpMatriceRelax);
		printf("         DeltaSigma changed to %lf\n",lpConfig->DeltaSigma);
	}
	lpConfig->Sigma += lpConfig->DeltaSigma;
	dDummy = 100.0 * (lpConfig->Sigma / lpConfig->SigmaMax) ;
	printf("      Current Sigma = %lf\t\t Done %3.2lf %%\n",lpConfig->Sigma, dDummy);

	if ( lpConfig->Diffraction ) {		// CNFD DIFFRACTION
		for (j=lpMatriceDiffract->Start; j<=lpMatriceDiffract->Finish; j++) {
			SumRho[j]=0.0;
			SolRho[j]=0.0;
		}
		for (i=1; i<=lpConfig->NbTau-1; i++) {
			for (j=lpMatriceDiffract->Start; j<=lpMatriceDiffract->Finish; j++)
				SumRho[j] += SolRho[j];
			RhsRho[0]=2.0*P[i][0] + lpMatriceDiffract->Lambda*(SumRho[1]-SumRho[0])/2.0;
			for (j=1; j<=lpConfig->NbRho-2; j++)
				RhsRho[j]=2.0*P[i][j] + (1.0 - 1.0/(2.0*j)) * lpMatriceDiffract->Lambda * SumRho[j-1]/8.0
								- lpMatriceDiffract->Lambda * SumRho[j]/4.0 +
								(1.0 + 1.0/(2.0*j)) * lpMatriceDiffract->Lambda * SumRho[j+1]/8.0;
			RhsRho[lpConfig->NbRho-1]=2.0*P[i][lpConfig->NbRho-1] +
								(1.0 - 1.0/(2.0*(lpConfig->NbRho-1))) * lpMatriceDiffract->Lambda * SumRho[lpConfig->NbRho-2]/8.0
								- lpMatriceDiffract->Lambda * SumRho[lpConfig->NbRho-1]/4.0;
			Tridiag(lpConfig, lpMatriceDiffract, RhsRho, SolRho, BetaRho, GammaRho);
			for (j=lpMatriceDiffract->Start; j<=lpMatriceDiffract->Finish; j++)
				P[i][j]=SolRho[j]-P[i][j];
		}
	}

	if ( lpConfig->ThermoViscous ) {	// CNFD ABSORPTION
		for (j=0; j<=lpConfig->NbRho-1; j++) {
			for (i=lpMatriceAbsorb->Start; i<=lpMatriceAbsorb->Finish; i++)
				RhsTau[i] = lpMatriceAbsorb->Lambda * P[i-1][j]/2.0
								+ (1-lpMatriceAbsorb->Lambda) * P[i][j]
								+ lpMatriceAbsorb->Lambda * P[i+1][j]/2.0;
			Tridiag(lpConfig, lpMatriceAbsorb, RhsTau, SolTau, BetaTau, GammaTau);
			for (i=lpMatriceAbsorb->Start; i<=lpMatriceAbsorb->Finish; i++)
				P[i][j]=SolTau[i];
		}
	}
	
	if ( lpConfig->Relaxation ) {		// CNFD RELAXATION
		for (n2=0; n2<lpConfig->NbRelax; n2++) {
			lpMatriceRelaxTmp->a=lpMatriceRelax->a[n2];
			lpMatriceRelaxTmp->b=lpMatriceRelax->b[n2];
			lpMatriceRelaxTmp->c=lpMatriceRelax->c[n2];
			lpMatriceRelaxTmp->Size=lpMatriceRelax->Size;
			lpMatriceRelaxTmp->Start=lpMatriceRelax->Start;
			lpMatriceRelaxTmp->Finish=lpMatriceRelax->Finish;
			for (j=0; j<=lpConfig->NbRho-1; j++) {
				for (i=lpMatriceRelaxTmp->Start; i<=lpMatriceRelaxTmp->Finish; i++)
					RhsTau[i] = -lpMatriceRelaxTmp->c[1] * P[i-1][j]
								+ (1.0-2.0*lpMatriceRelax->Lambda[n2]) * P[i][j]
								- lpMatriceRelaxTmp->a[1] * P[i+1][j];				
				Tridiag(lpConfig, lpMatriceRelaxTmp, RhsTau, SolTau, BetaTau, GammaTau);
				for (i=lpMatriceRelaxTmp->Start; i<=lpMatriceRelaxTmp->Finish; i++)
					P[i][j]=SolTau[i];
			}
		}
	}

	if ( lpConfig->NonLinearity ) {
		iResNL=NonLinear(lpConfig, P, Pold, TauDisto);
		if ( iResNL == -1) {
			printf("\nSimulation aborted !\n");
			return(-1);
		}
	}

	if ( lpConfig->Sigma >= dTarget ) {
		printf("         Writing data to output file...\n");
		while ( lpConfig->SigmaWatch[iTargetID] == dTarget ) {
			fwrite(&(lpConfig->Sigma), sizeof(double), 1, fp);
			jTmp = (int) floor( lpConfig->RhoWatch[iTargetID] * lpConfig->NbPiston );
			dRhoTmp = jTmp * lpConfig->DeltaRho;
			fwrite(&(dRhoTmp), sizeof(double), 1, fp);
			for (i=0; i<=lpConfig->NbTau; i++)
				fwrite(&(P[i][jTmp]), sizeof(double), 1, fp);
			printf("            Sigma = %lf Rho = %lf\n",lpConfig->Sigma,dRhoTmp);
			iTargetID++;
		}
		dTarget=lpConfig->SigmaWatch[iTargetID];
	}
}

printf("   Ended CNFD iterations at Sigma = %lf\n\n", lpConfig->Sigma);

fclose(fp);
return 0;
}


//////////////////////////////////////////////////////////////////////////

//
// Load input file
//


int InputConfig(LPCONFIG lpConfig, char *sConfigFileName)
{
FILE *fp;
int i, iDummy;
double dDummy1, dDummy2;

fp=fopen(sConfigFileName,"rt");
if( fp == NULL ) {
    printf("\nParameter file %s not found !\n", sConfigFileName);
	return(-1);
}
else {
	printf("   Input file is %s\n", sConfigFileName );
	printf("   Loading simulation parameters...\n");
	fscanf( fp, "Nonlinear effects ?\t%d\n", &((bool)lpConfig->NonLinearity) );
	fscanf( fp, "Diffraction effects ?\t%d\n", &((bool)lpConfig->Diffraction) );
	fscanf( fp, "Thermoviscous effects ?\t%d\n", &((bool)lpConfig->ThermoViscous) );
	fscanf( fp, "Relaxation effects ?\t%d\n", &((bool)lpConfig->Relaxation) );
	fscanf( fp, "Source file ?\t\t%d\n\n", &((bool)lpConfig->SourceFile) );
	fscanf( fp, "Output filename:\t\t%s\n",  &(lpConfig->OutputFileName) );

	fscanf( fp, "\nSource filename:\t\t%s\n",  &(lpConfig->SourceFileName) );
	fscanf(fp, "Number of cycles:\t\t%lf\n", &(lpConfig->Cycles) );
	fscanf(fp, "Envelope exponent:\t%d\n\n", &(lpConfig->Exponent) );

	fscanf( fp, "\nGain:\t\t\t\t%lf\n", &(lpConfig->Gain) );
	fscanf( fp, "\nAbsorption coeff:\t\t%lf\n", &(lpConfig->AbsorptionCoeff) );
	fscanf( fp, "\nNonlinearity coeff:\t%lf\n\n", &(lpConfig->NonLinearityCoeff) );
	fscanf( fp, "\nNumber of relaxation processes:\t%d\n", &(lpConfig->NbRelax) );
	if ( lpConfig->NbRelax > 0 ) {
		lpConfig->CNu=(double *)malloc(lpConfig->NbRelax*sizeof(double));
		lpConfig->ThetaNu=(double *)malloc(lpConfig->NbRelax*sizeof(double));
		for (i=0; i<lpConfig->NbRelax; i++) {
			fscanf( fp, "%lf\t%lf\n", &dDummy1, &dDummy2 );
			lpConfig->CNu[i]=dDummy1;
			lpConfig->ThetaNu[i]=dDummy2;
		}
	}
	fscanf( fp, "\nNumber of watchpoints:\t%d\n", &(lpConfig->NbWatch) );
	lpConfig->SigmaWatch=(double *)malloc(lpConfig->NbWatch*sizeof(double));
	lpConfig->RhoWatch=(double *)malloc(lpConfig->NbWatch*sizeof(double));
	for (i=0; i<lpConfig->NbWatch; i++) {
		fscanf( fp, "%lf\t%lf\n", &dDummy1, &dDummy2 );
		lpConfig->SigmaWatch[i]=dDummy1;
		lpConfig->RhoWatch[i]=dDummy2;
		if ( (i>0) && ((lpConfig->SigmaWatch[i]-lpConfig->SigmaWatch[i-1]) < 0) ) {
			printf("\nWarning: SigmaWatch values must be in ascending order !\n");
			return(-1);
		}
		if ( lpConfig->RhoWatch[i] < 0)  {
			printf("\nWarning: RhoWatch values must be positive !\n");
			return(-1);
		}
	} 
	lpConfig->SigmaMax=lpConfig->SigmaWatch[lpConfig->NbWatch-1];
	fscanf( fp, "\nTauMin:\t%lf\n", &(lpConfig->TauMin) );
	fscanf( fp, "TauMax:\t%lf\n", &(lpConfig->TauMax) );
	if ( lpConfig->TauMin >= lpConfig->TauMax)  {
		printf("\nWarning: TauMin must be smaller than TauMax !\n");
		return(-1);
	}
	fscanf( fp, "Number of steps per cycle:\t%d\n", &iDummy );
	lpConfig->NbTau=(int) (ceil((double) (iDummy * (lpConfig->TauMax-lpConfig->TauMin))));
	lpConfig->TauMin *= (2*PI);
	lpConfig->TauMax *= (2*PI);
	lpConfig->DeltaTau=(lpConfig->TauMax-lpConfig->TauMin) / (double) lpConfig->NbTau;
	lpConfig->Tau=(double *)malloc((lpConfig->NbTau+1)*sizeof(double));
	for (i=0; i<=lpConfig->NbTau; i++)
		lpConfig->Tau[i]=lpConfig->TauMin + i * lpConfig->DeltaTau;

	fscanf( fp, "RhoMax:\t%lf\n", &(lpConfig->RhoMax) );	
	fscanf( fp, "Number of points across piston:\t%d\n", &(lpConfig->NbPiston) );
	lpConfig->NbRho=(int) (lpConfig->RhoMax * lpConfig->NbPiston);
	lpConfig->DeltaRho=1.0 / (double) lpConfig->NbPiston;
	fscanf( fp, "DeltaSigmaIB:\t%lf\n", &(lpConfig->DeltaSigmaIB) );
	fscanf( fp, "DeltaSigmaCN:\t%lf\n", &(lpConfig->DeltaSigmaCN) );
	fscanf( fp, "Number of steps for IB:\t%d\n", &(lpConfig->NbStepsIB) );
	lpConfig->DeltaSigmaShock=11*lpConfig->DeltaSigmaIB;
	if ( lpConfig->SigmaMax <= (lpConfig->NbStepsIB * lpConfig->DeltaSigmaIB) )  {
		printf("\nWarning: Last SigmaWatch must fall into CNFD range !\n");
		return(-1);
	}
	lpConfig->SigmaMax += __max(lpConfig->DeltaSigmaIB, lpConfig->DeltaSigmaCN);
}

fclose(fp);
return(0);
}


//////////////////////////////////////////////////////////////////////////

//
// Initialize P array		/////////////////////////////
//


int Init(LPCONFIG lpConfig, double **P)

{
double dDummy, dRho, dTauIn, dTauNext, dTauLast, dDeltaTauIn, dPNext, dPLast, dDeltaP;
int i, j, k, iRhoIndex, iNbRhoIn, iNbTauIn, jRho, jRhoLast, jDeltaRho;
FILE *fp;

if ( lpConfig->SourceFile ) {
	fp=fopen(lpConfig->SourceFileName,"rb");
	if( fp == NULL ) {
		printf("\nSource file %s not found !\n", lpConfig->SourceFileName);
		return(-1);
	}
	else {	
		printf("\n   Initial data file is %s\n", lpConfig->SourceFileName);
		printf("   Loading initial data...\n");
	}
	fread( &iNbRhoIn, sizeof(int), 1, fp );
	jRhoLast=-1;
	for (iRhoIndex=0; iRhoIndex<iNbRhoIn; iRhoIndex++) {
		fread( &dRho, sizeof(double), 1, fp );
		fread( &dTauIn, sizeof(double), 1, fp );
		fread( &dDeltaTauIn, sizeof(double), 1, fp );
		fread( &iNbTauIn, sizeof(int), 1, fp );
		dDummy = dRho * lpConfig->NbPiston;
		jRho = (int) floor( dDummy );
		if ( (double) (dDummy - jRho) > 0.5 )
			jRho++;
		if ( ( iRhoIndex == 0 ) && ( jRho != 0 ) ) {
			printf("\nWarning: First Rho in source file has to be zero !\n");
			return(-1);
		}
		i=0;
		while ( lpConfig->Tau[i] < dTauIn ) {
			P[i][jRho] = 0.0;
			i++;
		}
		fread( &dPNext, sizeof(double), 1, fp );
		dTauNext=dTauIn;
		k=1;
		while ( k <= iNbTauIn ) {
			dPLast=dPNext;
			dTauLast=dTauNext;
			dTauNext=dTauIn + k*dDeltaTauIn;
			fread( &dPNext, sizeof(double), 1, fp );
			dDeltaP=(dPNext-dPLast) / dDeltaTauIn;
			k++;
			while ( ( lpConfig->Tau[i] > dTauNext ) && ( k <= iNbTauIn ) ){
				dPLast=dPNext;
				dTauLast=dTauNext;
				dTauNext=dTauIn + k*dDeltaTauIn;
				fread( &dPNext, sizeof(double), 1, fp );
				dDeltaP=(dPNext-dPLast) / dDeltaTauIn;
				k++;
			}
			if ( k <= iNbTauIn ) {
				while ( lpConfig->Tau[i] <= dTauNext ) {
					P[i][jRho] = dPLast + ( lpConfig->Tau[i] - dTauLast ) * dDeltaP;
					i++;
				}
			}
		}
		while ( i <= lpConfig->NbTau ) {
			P[i][jRho] = 0.0;
			i++;
		}
		jDeltaRho=jRho-jRhoLast;
		if ( jDeltaRho > 1 ) {
			for (i=0; i<=lpConfig->NbTau; i++) {
				dDeltaP = ( P[i][jRho] - P[i][jRhoLast] ) / (double) jDeltaRho;
				for (j=1; j<= jDeltaRho-1; j++) 
					P[i][jRhoLast+j] = P[i][jRhoLast] + dDeltaP * (double) j;
			}
		}
		if ( jDeltaRho < 0 ) {	
			printf("\nWarning: Rho values in source file are not in ascending order !\n");
			return(-1);
		}
		jRhoLast=jRho;	
	}
	fclose(fp);
	if ( jRhoLast < lpConfig->NbPiston ) {
		printf("\nWarning: Initial data is not specified at edge of piston !\n");
		printf("         Assuming uniform source from last rho point to piston edge.\n");
		for (j=jRhoLast+1; j<=lpConfig->NbPiston; j++) {
			dDummy = (double) lpConfig->Gain * ( j*j - jRhoLast*jRhoLast ) * lpConfig->DeltaRho*lpConfig->DeltaRho;
			k = (int) floor( lpConfig->NbTau*dDummy/(2*PI) - 0.5 );
			dDummy = (dDummy - 2*PI*k/lpConfig->NbTau) / lpConfig->DeltaTau;
			for (i=k; i<=(lpConfig->NbTau-1); i++ )
				P[i-k][j] = P[i][jRhoLast] + (P[i+1][jRhoLast]-P[i][jRhoLast])*dDummy;
			for (i=(lpConfig->NbTau-k); i<=lpConfig->NbTau; i++)
				P[i][j] = 0.0;
		}
		jRhoLast=lpConfig->NbPiston;
	}
	if ( jRhoLast > lpConfig->NbPiston ) {
		printf("\nWarning: Initial data is specified beyond edge of piston !\n");
		printf("         Simulation parameters may not make sense with this source file.\n");
	}
	for (j=jRhoLast+1; j<=lpConfig->NbRho; j++) {
		for (i=0; i<=lpConfig->NbTau; i++)
			P[i][j] = 0.0;
	}
}
else {
	for (i=0; i<=lpConfig->NbTau; i++) {
		for (j=0; j<=lpConfig->NbPiston; j++) {
			dDummy = lpConfig->Tau[i] + lpConfig->Gain*j*j*lpConfig->DeltaRho*lpConfig->DeltaRho;
			P[i][j]=(double) exp( - pow(dDummy/(PI*lpConfig->Cycles),2*lpConfig->Exponent) )
									* sin( dDummy );
		}
		for (j=lpConfig->NbPiston+1; j<=lpConfig->NbRho; j++) {
			P[i][j]=0.0;
		}
	}
}

return(0);
}


//////////////////////////////////////////////////////////////////////////

//
// Initialize IB matrices		/////////////////////////////
//


void IBFDinitMat(LPCONFIG lpConfig, LPMAT_1 lpMatriceAbsorb, LPMAT_1 lpMatriceDiffract, LPMAT_2 lpMatriceRelax)

{
int i, j;

if ( lpConfig->ThermoViscous ) {
	lpMatriceAbsorb->Start=1;
	lpMatriceAbsorb->Finish=lpConfig->NbTau-1;
	lpMatriceAbsorb->Lambda=lpConfig->AbsorptionCoeff * lpConfig->DeltaSigma / (lpConfig->DeltaTau*lpConfig->DeltaTau);
	for (i=lpMatriceAbsorb->Start; i<=lpMatriceAbsorb->Finish; i++) {
		lpMatriceAbsorb->a[i]=-lpMatriceAbsorb->Lambda;
		lpMatriceAbsorb->b[i]=1.0+2.0*lpMatriceAbsorb->Lambda;
		lpMatriceAbsorb->c[i]=-lpMatriceAbsorb->Lambda;
	}
}

if ( lpConfig->Diffraction ) {
	lpMatriceDiffract->Start=0;
	lpMatriceDiffract->Finish=lpConfig->NbRho-1;
	lpMatriceDiffract->Lambda=lpConfig->DeltaTau * lpConfig->DeltaSigma / (lpConfig->Gain*lpConfig->DeltaRho*lpConfig->DeltaRho);
	lpMatriceDiffract->a[0]=0.0;
	lpMatriceDiffract->b[0]=1.0+lpMatriceDiffract->Lambda/2.0;
	lpMatriceDiffract->c[0]=-lpMatriceDiffract->Lambda/2.0;
	for (j=lpMatriceDiffract->Start+1; j<=lpMatriceDiffract->Finish-1; j++) {
		lpMatriceDiffract->a[j]=(-1.0+1.0/(2.0*j))*lpMatriceDiffract->Lambda/8.0;
		lpMatriceDiffract->b[j]=1.0+lpMatriceDiffract->Lambda/4.0;
		lpMatriceDiffract->c[j]=-(1.0+1.0/(2.0*j))*lpMatriceDiffract->Lambda/8.0;
	}
	lpMatriceDiffract->a[lpMatriceDiffract->Finish]=(-1.0+1.0/(2.0*lpMatriceDiffract->Finish))*lpMatriceDiffract->Lambda/8.0;
	lpMatriceDiffract->b[lpMatriceDiffract->Finish]=1.0+lpMatriceDiffract->Lambda/4.0;
	lpMatriceDiffract->c[lpMatriceDiffract->Finish]=0.0;
}

if ( lpConfig->Relaxation ) {
	lpMatriceRelax->Start=1;
	lpMatriceRelax->Finish=lpConfig->NbTau-1;
	for (j=0; j<lpMatriceRelax->NbMat; j++) {
		lpMatriceRelax->Mu[j] = 0.5 * lpConfig->ThetaNu[j] / lpConfig->DeltaTau;
		lpMatriceRelax->Lambda[j] = lpConfig->CNu[j] * lpConfig->DeltaSigma / (lpConfig->DeltaTau*lpConfig->DeltaTau);
		for (i=lpMatriceRelax->Start; i<=lpMatriceRelax->Finish; i++) {
			lpMatriceRelax->a[j][i]=-lpMatriceRelax->Mu[j]-lpMatriceRelax->Lambda[j];
			lpMatriceRelax->b[j][i]=1.0+2.0*lpMatriceRelax->Lambda[j];
			lpMatriceRelax->c[j][i]=lpMatriceRelax->Mu[j]-lpMatriceRelax->Lambda[j];
		}
	}
}

}


//////////////////////////////////////////////////////////////////////////

//
// Initialize CN matrices		/////////////////////////////
//


void CNFDinitMat(LPCONFIG lpConfig, LPMAT_1 lpMatriceAbsorb, LPMAT_1 lpMatriceDiffract, LPMAT_2 lpMatriceRelax)

{
int i, j;

if ( lpConfig->ThermoViscous ) {
	lpMatriceAbsorb->Start=1;
	lpMatriceAbsorb->Finish=lpConfig->NbTau-1;
	lpMatriceAbsorb->Lambda=lpConfig->AbsorptionCoeff * lpConfig->DeltaSigma / (lpConfig->DeltaTau*lpConfig->DeltaTau);
	for (i=lpMatriceAbsorb->Start; i<=lpMatriceAbsorb->Finish; i++) {
		lpMatriceAbsorb->a[i]=-lpMatriceAbsorb->Lambda/2.0;
		lpMatriceAbsorb->b[i]=1.0+lpMatriceAbsorb->Lambda;
		lpMatriceAbsorb->c[i]=-lpMatriceAbsorb->Lambda/2.0;
	}
}

if ( lpConfig->Diffraction ) {
	lpMatriceDiffract->Start=0;
	lpMatriceDiffract->Finish=lpConfig->NbRho-1;
	lpMatriceDiffract->Lambda=lpConfig->DeltaTau * lpConfig->DeltaSigma / (lpConfig->Gain*lpConfig->DeltaRho*lpConfig->DeltaRho);
	lpMatriceDiffract->a[0]=0.0;
	lpMatriceDiffract->b[0]=1.0+lpMatriceDiffract->Lambda/4.0;
	lpMatriceDiffract->c[0]=-lpMatriceDiffract->Lambda/4.0;
	for (j=lpMatriceDiffract->Start+1; j<=lpMatriceDiffract->Finish-1; j++) {
		lpMatriceDiffract->a[j]=-(1.0-1.0/(2.0*j))*lpMatriceDiffract->Lambda/16.0;
		lpMatriceDiffract->b[j]=1.0+lpMatriceDiffract->Lambda/8.0;
		lpMatriceDiffract->c[j]=-(1.0+1.0/(2.0*j))*lpMatriceDiffract->Lambda/16.0;
	}
	lpMatriceDiffract->a[lpMatriceDiffract->Finish]=-(1.0-1.0/(2.0*lpMatriceDiffract->Finish))*lpMatriceDiffract->Lambda/16.0;
	lpMatriceDiffract->b[lpMatriceDiffract->Finish]=1.0+lpMatriceDiffract->Lambda/8.0;
	lpMatriceDiffract->c[lpMatriceDiffract->Finish]=0.0;
}

if ( lpConfig->Relaxation ) {
	lpMatriceRelax->Start=1;
	lpMatriceRelax->Finish=lpConfig->NbTau-1;
	for (j=0; j<lpMatriceRelax->NbMat; j++) {
		lpMatriceRelax->Mu[j] = 0.5 * lpConfig->ThetaNu[j] / lpConfig->DeltaTau;
		lpMatriceRelax->Lambda[j] = 0.5 * lpConfig->CNu[j] * lpConfig->DeltaSigma / (lpConfig->DeltaTau*lpConfig->DeltaTau);
		for (i=lpMatriceRelax->Start; i<=lpMatriceRelax->Finish; i++) {
			lpMatriceRelax->a[j][i]=-lpMatriceRelax->Mu[j]-lpMatriceRelax->Lambda[j];
			lpMatriceRelax->b[j][i]=1.0+2.0*lpMatriceRelax->Lambda[j];
			lpMatriceRelax->c[j][i]=lpMatriceRelax->Mu[j]-lpMatriceRelax->Lambda[j];
		}
	}
}

}


//////////////////////////////////////////////////////////////////////////

//
// Tridiag Subroutine		/////////////////////////////
//


void Tridiag(LPCONFIG lpConfig, LPMAT_1 lpMatrice, double *Rhs, double *Sol, double *Beta, double *Gamma)

{
int i;

Beta[lpMatrice->Start] = lpMatrice->b[lpMatrice->Start];
Gamma[lpMatrice->Start] = Rhs[lpMatrice->Start] / Beta[lpMatrice->Start];

for (i=lpMatrice->Start+1; i<=lpMatrice->Finish; i++) {
	Beta[i]=lpMatrice->b[i] - lpMatrice->a[i] * lpMatrice->c[i-1] / Beta[i-1];
	Gamma[i]=(Rhs[i] - lpMatrice->a[i] * Gamma[i-1]) / Beta[i];
}

Sol[lpMatrice->Finish]=Gamma[lpMatrice->Finish];

for (i=lpMatrice->Finish-1; i>=lpMatrice->Start; i--)
	Sol[i]=Gamma[i] - lpMatrice->c[i] * Sol[i+1] / Beta[i];
}


//////////////////////////////////////////////////////////////////////////

//
// Nonlinear Subroutine		/////////////////////////////
//


int NonLinear(LPCONFIG lpConfig, double **P, double *Pold, double *TauDisto)

{
int i, j, k;
double dDisto, dDeltaPmax, dDeltaPdt;

dDisto = lpConfig->NonLinearityCoeff * lpConfig->DeltaSigma;
TauDisto[0]=lpConfig->Tau[0];
TauDisto[lpConfig->NbTau]=lpConfig->Tau[lpConfig->NbTau];
dDeltaPmax=0.0;

for (j=0; j<=lpConfig->NbRho-1; j++) {
	Pold[0]=P[0][j];
	for (i=1; i<=lpConfig->NbTau; i++) {
		Pold[i]=P[i][j];
		TauDisto[i]=lpConfig->Tau[i]-Pold[i]*dDisto;
		if ( TauDisto[i] <= TauDisto[i-1] ) {
			printf("\nWarning: Discontinuity at Sigma = %lf and time step %d\n",lpConfig->Sigma,i);
			printf("         Smaller DeltaSigma required !\n");
			return(-1);
		}
	}
	k=1;
	for (i=1; i<=lpConfig->NbTau-1; i++) {
		while (lpConfig->Tau[i] > TauDisto[k] ) {
			k++;
			if ( k > lpConfig->NbTau ) {
				for (k=i; k<=lpConfig->NbTau; k++)
					P[k][j]=P[lpConfig->NbTau][j];
			}
		}
		dDeltaPdt = (Pold[k]-Pold[k-1]) / (TauDisto[k]-TauDisto[k-1]);
		if ( dDeltaPdt > dDeltaPmax )
			dDeltaPmax = dDeltaPdt;
		P[i][j] = Pold[k-1] + dDeltaPdt * (lpConfig->Tau[i] - TauDisto[k-1]); 
	}
}

lpConfig->DeltaSigmaShock = 1.0 / (lpConfig->NonLinearityCoeff * dDeltaPmax);

return(0);
}