/*CheckSanity, check that various parts of the code still work */
#include <float.h>
#include <string.h>
#include "cddefines.h"
#include "expnf.h"
#include "ee1.h"
#include "matin1.h"
#include "elmton.h"
#include "linpack.h"
#include "veclib.h"
#include "showme.h"
#include "taulines.h"
#include "e2.h"
#include "bevington.h"
#include "typmatrx.h"
#include "opacpoint.h"
#include "nsshells.h"
#include "trace.h"
#include "checksanity.h"

void CheckSanity(void)
{
	int lgOK=TRUE;
	long ner , ipiv[3] , job , i , j , 
		ipZ ,
		nelem , 
		nion ,
		nshells;
	/* this will be charge tot he 4th power */
	double Z4;

#	define NDIM 10
	double x , ans1 , ans2  , xMatrix[NDIM][NDIM] , yVector[NDIM] ,
		rcond, 
	  work[NDIM];

#	ifdef DEBUG_FUN
	fputs( "<+>CheckSanity()\n", debug_fp );
#	endif

	/*********************************************************
	 *                                                       *
	 * confirm that various part of cloudy still work        *
	 *                                                       *
	 *********************************************************/

	/* if this is no longer true at end, we have a problem */
	lgOK = TRUE;

	/*********************************************************
	 *                                                       *
	 * confirm that hydrogen einstein As are still valid     *
	 *                                                       *
	 *********************************************************/
	for( ipZ=0; ipZ<2; ++ipZ )
	{
		/* this element may be turned off */
		if( elmton.lgElmtOn[ipZ] )
		{ 
			/*Z4 = (double)(POW2(ipZ+1)*POW2(ipZ+1));*/
			/* form charge to the 4th power */
			Z4 = (double)(ipZ+1);
			Z4 *= Z4;
			Z4 *= Z4;
			/* H Lya */
			ans1 = HydroLines[ipZ][IP2P][IP1S].Aul;
			ans2 = 6.265e8*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-3 )
			{
				fprintf(ioQQQ , "CheckSanity finds insane A for H Lya %g %g ipZ=%li\n",
					ans1 , ans2 , ipZ );
				lgOK = FALSE;
			}

			/* H two photon */
			ans1 = HydroLines[ipZ][IP2S][IP1S].Aul;
			ans2 = 8.23*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-3 )
			{
				fprintf(ioQQQ , "CheckSanity finds insane A for H 2-phot %g %g ipZ=%li\n",
					ans1 , ans2 , ipZ );
				lgOK = FALSE;
			}

			/* Balmer alpha - only way out of 3 for case b*/
			ans1 = HydroLines[ipZ][3][IP2S].Aul;
			ans2 = 1.10318e7*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-3 )
			{
				fprintf(ioQQQ , "CheckSanity finds insane A for H 2s %g %g ipZ=%li\n",
					ans1 , ans2 , ipZ );
				lgOK = FALSE;
			}
			ans1 = HydroLines[ipZ][3][IP2P].Aul;
			ans2 = 3.3095e7*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-3 )
			{
				fprintf(ioQQQ , "CheckSanity finds insane A for H 2s %g %g ipZ=%li\n",
					ans1 , ans2 , ipZ );
				lgOK = FALSE;
			}

			/* for higher line the branching ratio depends on l-mixing and so no close
			 * test is possible.  However all decays under case b should add up to
			 * hlife for that level */
			/* Balmer beta 4 - 2 plus Paschen alpha, 4-3*/
			ans1 = HydroLines[ipZ][4][IP2S].Aul + 
				HydroLines[ipZ][4][IP2P].Aul + HydroLines[ipZ][4][3].Aul;
			ans2 = (8.424e6+ 8.991e6)*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-2 )
			{
				fprintf(ioQQQ , 
					"CheckSanity finds insane A for summed H=4 decays found=%g correct=%g ipZ=%li\n",
					ans1 , ans2 , ipZ );
				lgOK = FALSE;
			}

			/* Balmer gamma 5 - 2 + 5-4 + 5-3*/
			ans1 = HydroLines[ipZ][5][IP2S].Aul + 
				HydroLines[ipZ][5][IP2P].Aul + 
				HydroLines[ipZ][5][3].Aul + HydroLines[ipZ][5][4].Aul;
			ans2 = (2.532e6+2.20e6+2.70e6)*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-2 )
			{
				fprintf(ioQQQ , 
					"CheckSanity finds insane A for H5-2 found=%g correct=%g ipZ=%li\n",
					ans1 , ans2 , ipZ );
				lgOK = FALSE;
			}

		}
	}

	/*********************************************************
	 *                                                       *
	 * confirm that exponential integral routines still work *
	 *                                                       *
	 *********************************************************/

	/* check that first and second exponential integrals are ok,
	 * step through range of values, beginning with following */
	x = 1e-3;
	do
	{
		/* check that fast e1 routine is ok */
		ans1 = ee1(x);
		ans2 = expnf( 1 , x );
		if( fabs(ans1-ans2)/(ans1+ans2) > 1e-2 )
		{
			fprintf(ioQQQ , "CheckSanity finds insane E1 %g %g %g\n",
				x , ans1 , ans2 );
			lgOK = FALSE;
		}

		/* check that e2 is ok */
		ans1 = e2(x,exp(-x));
		ans2 = expnf( 2 , x );
		if( fabs(ans1-ans2)/(ans1+ans2) > 1e-2 )
		{
			fprintf(ioQQQ , "CheckSanity finds insane E2 %g %g %g\n",
				x , ans1 , ans2 );
			lgOK = FALSE;
		}

		/* now increment x */
		x *= 2.;
		/* following limit set by sexp returning zero, used in ee1 */
	} while( x < 64. );


	/*********************************************************
	 *                                                       *
	 * confirm that matrix inversion routine still works     *
	 *                                                       *
	 *********************************************************/

	/* these are the answer, chosen to get xvec 1,2,3 */
	yVector[0] = 1.;
	yVector[1] = 3.;
	yVector[2] = 3.;

	/* zero out the main matrix */
	for(i=0;i<3;++i)
	{
		for( j=0;j<3;++j )
		{
			xMatrix[i][j] = 0.;
		}
	}

	/* remember that order is column, row, alphabetical order, rc */
	xMatrix[0][0] = 1.;
	xMatrix[0][1] = 1.;
	xMatrix[1][1] = 1.;
	xMatrix[2][2] = 1.;

	/* which matrix solver? */
	if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 )
	{
		/*matin1();*/
		ner = 1;
		if( ner != 0 )
		{
			fprintf( ioQQQ, " CheckSanity MATRIX ERROR.\n" );
			puts( "[Stop in CheckSanity]" );
			exit(1);
		}
	}

#	if 0
	/* this test is for 2-d form of xMartrix, now commmented out */
	/* this is the default matrix solver */
	else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 )
	{
		/*void DGETRF(long,long,double*,long,long[],long*);*/
		DGETRF(3,3,(double*)xMatrix,NDIM,ipiv,&ner);
		if( ner != 0 )
		{
			fprintf( ioQQQ, " CheckSanity DGETRF error\n" );
			puts( "[Stop in CheckSanity]" );
			exit(1);
		}

		/* usage DGETRS, 'N' = no transpose
		 * order of matrix,
		 * number of cols in bvec, =1
		 * array
		 * leading dim of array */
		DGETRS('N',3,1,(double*)xMatrix,NDIM,ipiv,yVector,3,&ner);

		if( ner != 0 )
		{
			fprintf( ioQQQ, " CheckSanity DGETRS error\n" );
			puts( "[Stop in CheckSanity]" );
			exit(1);
		}
	}
#endif

	/* this is the default matrix solver */
	/* this test is the 1-d matrix with 2-d macro simulation */
	else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 )
	{
		double *A;
	  /* LDA is right dimension of matrix */
#		define AA(I_,J_)	(*(A+(I_)*(NDIM)+(J_)))
		/* malloc space for the  1-d array */
		if( (A=(double*) malloc( (sizeof(double)*NDIM*NDIM ))) == NULL )
		{
			fprintf( ioQQQ, " CheckSanity malloc A error\n" );
			puts( "[Stop in CheckSanity]" );
			exit(1);
		}
		/* copy over the main matrix */
		for(i=0;i<3;++i)
		{
			for( j=0;j<3;++j )
			{
				AA(i,j) = xMatrix[i][j];
			}
		}

		/*void DGETRF(long,long,double*,long,long[],long*);*/
		DGETRF(3,3,A,NDIM,ipiv,&ner);
		if( ner != 0 )
		{
			fprintf( ioQQQ, " CheckSanity DGETRF error\n" );
			puts( "[Stop in CheckSanity]" );
			exit(1);
		}

		/* usage DGETRS, 'N' = no transpose
		 * order of matrix,
		 * number of cols in bvec, =1
		 * array
		 * leading dim of array */
		DGETRS('N',3,1,A,NDIM,ipiv,yVector,3,&ner);

		if( ner != 0 )
		{
			fprintf( ioQQQ, " CheckSanity DGETRS error\n" );
			puts( "[Stop in CheckSanity]" );
			exit(1);
		}
		/* release the vector */
		free( A );
	}

	else if( strcmp(TypMatrx.chMatrix,"veclib ") == 0 )
	{
		/* Jason found this one on the Exemplar, distributed source just stops */
		fprintf( ioQQQ, " this has not been checked since H atom conv\n" );
		/*666 error! in following one is dimension, second is number of levels */
		rcond = 0.;
		job = 0;
		dgeco((double*)xMatrix,NDIM,NDIM,ipiv,rcond,work);
		dgesl((double*)xMatrix,NDIM,NDIM,ipiv,yVector,job);
		/* now put results back into z so rest of code treates only
		 * one case - as if matin1 had been used */

		puts( "[Stop in CheckSanity]" );
		exit(1);
	}

	else if( strcmp(TypMatrx.chMatrix,"Bevingt") == 0 )
	{
		/* this came from Bevington */
		matinv(xMatrix,3,&rcond);
		for( j=0; j < 3; j++ )
		{
			work[j] = 0.;
			for( i=0; i < 3; i++ )
			{
				work[j] += yVector[i]*xMatrix[i][j];
			}
		}
		for( j=0; j < 3; j++ )
		{
			yVector[j] = work[j];
		}
	}

	else
	{
		fprintf( ioQQQ, " chMatrix type insane in CheckSanity, was%7.7s\n", 
		  TypMatrx.chMatrix );
		puts( "[Stop in CheckSanity]" );
		exit(1);
	}

	/* now check on validity of the solution, demand that this
	 * simple problem have come within a few epsilons of the
	 * correct answer */

	/* find largest deviation */
	rcond = 0.;
	for(i=0;i<3;++i)
	{
		x = fabs( yVector[i]-i-1.);
		rcond = MAX2( rcond, x );
		/*printf(" %g ", yVector[i]);*/
	}

	if( rcond>DBL_EPSILON)
	{
		fprintf(ioQQQ,
			"CheckSanity found too large a deviation in matrix solver = %g \n", 
			rcond);
		/* set flag saying that things are not ok */
		lgOK = FALSE;
	}
	/* end matrix inversion check */


	/* these pointers were set to INT_MIN in CreatePoint,
	 * then set to valid numbers in ipShells and MakeOpacity
	 * this checks that all values have been properly filled */
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		if( elmton.lgElmtOn[nelem] )
		{
			for( nion=0; nion<=nelem; ++nion )
			{
				/* loop over all shells that are defined */
				for( nshells=0; nshells<nsShellsCom.nsShells[nion][nelem]; ++nshells )
				{
					for( j=0; j<3; ++j )
					{
						if( OpacPoint.ipElement[j][nshells][nion][nelem] <0 )
						{
							/* this is not possible */
							fprintf(ioQQQ,
								"CheckSanity found insane ipElement for nelem=%li nion=%li nshells=%li j=%li \n", 
								nelem , nion , nshells, j );
							fprintf(ioQQQ,
								"value was %li  \n", OpacPoint.ipElement[j][nshells][nion][nelem] );
							/* set flag saying that things are not ok */
							lgOK = FALSE;
						}
					}
				}
			}
		}
	}


	/*********************************************************
	 *                                                       *
	 * everything is done, all checks make, did we pass them?*
	 *                                                       *
	 *********************************************************/

	if( lgOK )
	{
		/*return if ok */
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " CheckSanity returns OK\n");
		}

#	ifdef DEBUG_FUN
		fputs( " <->CheckSanity()\n", debug_fp );
#	endif
		return;
	}

	else
	{
		/* stop since problem encountered, lgEOF set FALSE */
		fprintf(ioQQQ , "CheckSanity finds insanity so exiting\n");
		ShowMe();
		exit(1);
	}
}