/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*SanityCheck, check that various parts of the code still work, called by Cloudy after continuum
 * and optical depth arrays are set up, but before initial temperature and ionization */
#include "cddefines.h"
#include "physconst.h"
#include "math_complex.h"
#include "dense.h"
#include "lapack.h"
#include "helike.h"
#include "helike_recom.h"
#include "rfield.h"
#include "taulines.h"
#include "hypho.h"
#include "iso.h"
#include "opacity.h"
#include "hydro_bauman.h"
#include "hydrogenic.h"
#include "heavy.h"
#include "kshllenr.h"
#include "trace.h"
#include "gaunt.h"
#include "cloudy.h"

/* NB - this routine must not change any global variables - any that are changed as part of
 * a test must be reset, so that the code retains state */

static void SanityCheckBegin(void );
static void SanityCheckFinal(void );
/* chJob is either "begin" or "final" 
 * "begin is before code starts up
 * "final is after model is complete */
void SanityCheck( const char *chJob )
{
	if( strcmp(chJob,"begin") == 0 )
	{
		SanityCheckBegin();
	}
	else if( strcmp(chJob,"final") == 0 )
	{
		SanityCheckFinal();
	}
	else
	{
		fprintf(ioQQQ,"SanityCheck called with insane argument.\n");
		puts( "[Stop in SanityCheck]" );
		cdEXIT(EXIT_FAILURE);
	}
}

static void SanityCheckFinal(void )
{
	/* PrtComment also has some ending checks on sanity */
}

static void SanityCheckBegin(void )
{
	int lgOK=TRUE;
	int lgFlag;
	int32 ner, ipiv[3];
	long	i , 
		j , 
		nelem , 
		ion ,
		nshells;
	double *A;

	/* this will be charge to the 4th power */
	double Aul ,
		error,
		Z4, gaunt;

	long n, logu, loggamma2;
	long ipISO;

#	define NDIM 10
	double x , ans1 , ans2  , xMatrix[NDIM][NDIM] , yVector[NDIM] ,
		rcond;
	float *fvector;
	long int *ipvector;

#	ifdef DEBUG_FUN
	fputs( "<+>SanityCheck()\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;

	/*********************************************************
	 *                                                       *
	 * check that all the Lyas As are ok                     *
	 *                                                       *
	 *********************************************************/
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		/* this element may be turned off */
		if( dense.lgElmtOn[nelem] )
		{ 
			/* H_Einstein_A( n, l, np, lp, iz ) - all are on physics scale */
			Aul = H_Einstein_A( 2, 1, 1, 0, nelem+1 );
			/*fprintf(ioQQQ,"%li\t%.4e\n", nelem+1, EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Aul );*/
			if( fabs( Aul - EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Aul ) /Aul > 0.01 )
			{
				fprintf(ioQQQ," SanityCheck found insane H-like As.\n");
				lgOK = FALSE;
			}
		}
	}

	/*********************************************************
	 *                                                       *
	 * check that gaunt factors are good                     *
	 *                                                       *
	 *********************************************************/
	/* Uncommenting each of the four print statements here
	 * will produce a nice table comparable to Sutherland 98, Table 2.	*/
	/* fprintf(ioQQQ,"u\t-4\t-3\t-2\t-1\t0\t1\t2\t3\t4\n");*/
	for( logu=-4; logu<=4; logu++)
	{
		/*fprintf(ioQQQ,"%li\t", logu);*/
		for(loggamma2=-4; loggamma2<=4; loggamma2++)
		{ 
			double SutherlandGff[9][9]=
			{	{5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
				{4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
				{3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
				{1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
				{0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
				{0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
				{0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
				{0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
				{0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};

			gaunt = CalcThermAveGaunt( TE1RYD/pow(10.,(double)loggamma2), 1., pow(10.,(double)(logu-loggamma2)) );
			error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
			/*fprintf(ioQQQ,"%1.3f\t", gaunt);*/
			if( error>0.11 || ( loggamma2<2 && error>0.05 ) )
			{
				fprintf(ioQQQ," SanityCheck found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
					logu, loggamma2, error );
				lgOK = FALSE;
			}
		}
		/*fprintf(ioQQQ,"\n");*/
	}

	/*********************************************************
	 *                                                       *
	 * check some transition probabililties for he-like ions *
	 *                                                       *
	 *********************************************************/
	for( nelem=1; nelem<LIMELM; ++nelem )
	{
		/* the helike 9-1 transition, A(3^3P to 2^3S) */
		double as[]={
		 /* updated with Johnson values */
		 0.       ,9.47e+006 ,3.44e+008 ,1.74e+009 ,5.51e+009 ,1.34e+010 ,
		2.79e+010 ,5.32E+010 ,8.81e+010 ,1.46E+011 ,2.15e+011 ,3.15e+011 ,
		4.46e+011 ,6.39E+011 ,8.26e+011 ,1.09e+012 ,1.41e+012 ,1.86E+012 ,
		2.26e+012 ,2.80e+012 ,3.44e+012 ,4.18e+012 ,5.04e+012 ,6.02e+012 ,
		7.14e+012 ,8.40e+012 ,9.83e+012 ,1.14e+013 ,1.32e+013 ,1.52e+013
		};

		if( iso.numLevels_max[ipHE_LIKE][nelem] > 8 && dense.lgElmtOn[nelem])
		{	
			/* following used to print current values of As */
			/*fprintf(ioQQQ,"%.2e\n", EmisLines[ipHE_LIKE][nelem][9][1].Aul );*/
			if( fabs( as[nelem] - EmisLines[ipHE_LIKE][nelem][9][1].Aul ) /as[nelem] > 0.025 )
			{
				fprintf(ioQQQ,
					" SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n",
					nelem,
					as[nelem] , 
					EmisLines[ipHE_LIKE][nelem][9][1].Aul );
				lgOK = FALSE;
			}
		}
	}

	/* only do this test if case b is not in effect */
	if( !opac.lgCaseB )
	{

		for( i = 0; i <=110; i++ )
		{
			double DrakeTotalAuls[111] = {
				-1.0000E+00, -1.0000E+00, -1.0000E+00, 1.02160E+07,
				1.02160E+07, 1.02160E+07, 1.80090E+09, 2.78530E+07,
				1.82990E+07, 1.05480E+07, 7.07210E+07, 6.37210E+07,
				5.79960E+08, 1.60330E+07, 1.13640E+07, 7.21900E+06,
				3.11920E+07, 2.69830E+07, 1.38380E+07, 1.38330E+07,
				2.52270E+08, 9.20720E+06, 6.82220E+06, 4.56010E+06,
				1.64120E+07, 1.39290E+07, 7.16030E+06, 7.15560E+06,
				4.25840E+06, 4.25830E+06, 1.31150E+08, 5.62960E+06,
				4.29430E+06, 2.95570E+06, 9.66980E+06, 8.12340E+06,
				4.19010E+06, 4.18650E+06, 2.48120E+06, 2.48120E+06,
				1.64590E+06, 1.64590E+06, 7.65750E+07, 3.65330E+06,
				2.84420E+06, 1.99470E+06, 6.16640E+06, 5.14950E+06,
				2.66460E+06, 2.66200E+06, 1.57560E+06, 1.57560E+06,
				1.04170E+06, 1.04170E+06, 7.41210E+05, 7.41210E+05,
				4.84990E+07, 2.49130E+06, 1.96890E+06, 1.39900E+06,
				4.16900E+06, 3.46850E+06, 1.79980E+06, 1.79790E+06,
				1.06410E+06, 1.06410E+06, 7.02480E+05, 7.02480E+05,
				4.98460E+05, 4.98460E+05, -1.0000E+00, -1.0000E+00,
				3.26190E+07, 1.76920E+06, 1.41440E+06, 1.01460E+06,
				2.94830E+06, 2.44680E+06, 1.27280E+06, 1.27140E+06,
				7.52800E+05, 7.52790E+05, 4.96740E+05, 4.96740E+05,
				3.51970E+05, 3.51970E+05, -1.0000E+00, -1.0000E+00,
				-1.0000E+00, -1.0000E+00, 2.29740E+07, 1.29900E+06,
				1.04800E+06, 7.57160E+05, 2.16090E+06, 1.79030E+06,
				9.33210E+05, 9.32120E+05, 5.52310E+05, 5.52310E+05,
				3.64460E+05, 3.64460E+05, 2.58070E+05, 2.58070E+05,
				-1.0000E+00, -1.0000E+00, -1.0000E+00, -1.0000E+00,
				-1.0000E+00, -1.0000E+00, 1.67840E+07};

			if( DrakeTotalAuls[i] > 0. && 
				i < iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM])
			{
				if( fabs( DrakeTotalAuls[i] - (1./helike.Lifetime[ipHELIUM][i]) ) /DrakeTotalAuls[i] > 0.001 )
				{
					fprintf(ioQQQ,
						" SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n",
						i,
						DrakeTotalAuls[i], 
						(1./helike.Lifetime[ipHELIUM][i]) );
					lgOK = FALSE;
				}
			}
		}
	}

	/*********************************************************
	 *                                                       *
	 * check the threshold photoionization cs for He I       *
	 *                                                       *
	 *********************************************************/
	/* Don't check for sanity in these if running the Benjamin test...	*/
	if(  dense.lgElmtOn[ipHELIUM] && !helike.lgSetBenjamin )
	{
		/* HeI photoionization cross sections at threshold for lowest 20 levels */
#		define NHE1CS	20
		double he1cs[NHE1CS] = 
		{
			5.480E-18 , 9.253E-18 , 1.598E-17 , 1.598E-17 , 1.598E-17 , 1.348E-17 , 
			8.025E-18 , 1.449E-17 , 2.852E-17 , 1.848E-17 , 1.813E-17 , 2.699E-17 , 
			1.077E-17 , 2.038E-17 , 4.159E-17 , 3.670E-17 , 3.575E-17 , 1.900E-17 , 
			1.900E-17 , 4.175E-17 
		};

		/* loop over levels and check on photo cross section */
		j = MIN2( NHE1CS+1 , iso.numLevels_max[ipHE_LIKE][ipHELIUM] -iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] );
		for( n=1; n<j; ++n )
		{
			/* above list of levels does not include the ground */
			i = iso.ipOpac[ipHE_LIKE][ipHELIUM][n];
			ASSERT( i>0 );
			/*fprintf(ioQQQ,"%li\t%lin", n , i );*/
			/* >>chng 02 apr 10, from 0.01 to 0.02, values stored
			 * where taken from calc at low contin resolution, when continuum
			 * resolution changed this changes too */
			/*fprintf(ioQQQ,"%li %.2e\n", n,( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] );*/
			/* >>chng 02 jul 16, limt from 0.02 to 0.04, so that "set resolution 4" will work */
			/* >>chng 04 may 18, levels 10 and 11 are about 12% off - because of energy binning, chng from 0.08 to 0.15 */ 
			if( fabs( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] > 0.15 )
			{
				fprintf(ioQQQ,
					" SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n",
					n,
					he1cs[n-1] , 
					opac.OpacStack[i - 1]);
				fprintf(ioQQQ,
					" n=%li, l=%li, s=%li\n",
					iso.quant_desig[ipHE_LIKE][ipHELIUM][n].n ,
					iso.quant_desig[ipHE_LIKE][ipHELIUM][n].l ,
					iso.quant_desig[ipHE_LIKE][ipHELIUM][n].s);
				lgOK = FALSE;
			}
		}

		/* Check for agreement between on-the-fly and interpolation calculations
		 * of recombination, but only if interpolation is turned on. */
		if( !helike.lgNoRecombInterp )
		{
			/* check the recombination coefficients for he ground */
			error = fabs( HelikeCheckRecomb( ipHELIUM , 0 , 7500. ) );
			if( error > 0.01 )
			{
				fprintf(ioQQQ,
					" SanityCheck found insane1 HeI recom coef: expected, n=%i error: %.2e \n",
					0,
					error );
				lgOK = FALSE;
			}


			/* check the recombination coefficients for he 23P */
			error = fabs( HelikeCheckRecomb( ipHELIUM , 1 , 12500. ) );
			if( error > 0.01 )
			{
				fprintf(ioQQQ,
					" SanityCheck found insane2 HeI recom coef: expected, n=%i error: %.2e \n",
					1,
					error );
				lgOK = FALSE;
			}
		}

#		undef NHE1CS
	}


	/*********************************************************
	 *                                                       *
	 * check out the sorting routine                         *
	 *                                                       *
	 *********************************************************/

#	define NSORT 100 

	if( (fvector = (float *)MALLOC((NSORT)*sizeof(float) )) ==NULL )
		BadMalloc();

	if( (ipvector = (long *)MALLOC((NSORT)*sizeof(long int) )) ==NULL )
		BadMalloc();

	nelem = 1;
	/* make up some unsorted values */
	for( i=0; i<NSORT; ++i )
	{
		nelem *= -1;
		fvector[i] = (float)nelem * ((float)NSORT-i);
	}

	/*spsort netlib routine to sort array returning sorted indices */
	spsort(fvector, 
		   NSORT, 
		  ipvector, 
		  /* flag saying what to do - 1 sorts into increasing order, not changing
		   * the original routine */
		  1, 
		  &lgFlag);

	if( lgFlag ) lgOK = FALSE;

	for( i=1; i<NSORT; ++i )
	{
		/*fprintf(ioQQQ," %li %li %.0f\n", 
			i, ipvector[i],fvector[ipvector[i]] );*/
		if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] )
		{
			fprintf(ioQQQ," SanityCheck found insane sort\n");
			lgOK = FALSE;
		}
	}

	free( fvector );
	free( ipvector);
#	undef NSORT 

#	if 0
	ttemp = (float)sqrt(phycon.te);
	/* check that the temperatures make sense */
	if( fabs(ttemp - phycon.sqrte )/ttemp > 1e-5 )
	{
		fprintf(ioQQQ , "SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n",
			phycon.te , 
			sqrt(phycon.te) , 
			phycon.sqrte , 
			fabs(sqrt(phycon.te) - phycon.sqrte ) );
		lgOK = FALSE;
	}
#	endif

	/*********************************************************
	 *                                                       *
	 * confirm that widflx and anu arrays correspond         *
	 * to one another                                        *
	 *                                                       *
	 *********************************************************/

#	if 0
	/* this check on widflx can't be used since some sharpling curved continua, like laser,
	 * totally fail due to non-linear nature of widflx and anu relationship */
#	if !defined(NDEBUG)
	x = 0.;
	for( i=1; i<rfield.nupper-1; ++i )
	{
		if( fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.) > 0.02 )
		{
			ans1 = fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.) ;
			fprintf(ioQQQ," SanityCheck found insane widflx anu[i+1]=%e anu[i]=%e widflx=%e delta=%e rel err %e\n",
			rfield.anu[i+1] , rfield.anu[i] , rfield.widflx[i] , rfield.anu[i+1] -rfield.anu[i] , ans1 );
			lgOK = FALSE;
			x = MAX2( ans1 , x);
		}
		/* problems when at energy where resolution of grid changes dramatically */
		/* this is resolution at current energy */
		ans1 = rfield.widflx[i] / rfield.anu[i];
		if( (rfield.anu[i]+rfield.widflx[i]/2.)*(1.-ans1/10.) > rfield.anu[i+1] - rfield.widflx[i+1]/2.) 
		{
			fprintf(ioQQQ," SanityCheck found insane overlap1 widflx %e %e %e %e %e %e\n",
			rfield.anu[i] , rfield.widflx[i], rfield.anu[i] + rfield.widflx[i]/2. , rfield.anu[i+1], 
			rfield.widflx[i+1], rfield.anu[i+1] -rfield.widflx[i+1]/2. );
			lgOK = FALSE;
		}
		if( !lgOK )
		{
			fprintf(ioQQQ," big error was %e\n", x);
		}
	}
#	endif
#	endif


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

#			if 0
			/*must disable since, at this time, induced is included in Aul */
			/* H two photon */
			ans1 = EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].Aul;
			ans2 = 8.226*powi( (1.+nelem) , 6 );
			if( fabs(ans1-ans2)/ans2 > 1e-3 )
			{
				fprintf(ioQQQ , "SanityCheck finds insane A for H 2-phot %g %g nelem=%li\n",
					ans1 , ans2 , nelem );
				lgOK = FALSE;
			}
#			endif

#			if 0
			/* Balmer alpha - only way out of 3 for case b*/
			ans1 = EmisLines[ipH_LIKE][nelem][3][ipH2s].Aul;
			ans2 = 1.10318e7*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-3 )
			{
				fprintf(ioQQQ , "SanityCheck finds insane A for H 2s %g %g nelem=%li\n",
					ans1 , ans2 , nelem );
				lgOK = FALSE;
			}
			ans1 = EmisLines[ipH_LIKE][nelem][3][ipH2p].Aul;
			ans2 = 3.3095e7*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-3 )
			{
				fprintf(ioQQQ , "SanityCheck finds insane A for H 2s %g %g nelem=%li\n",
					ans1 , ans2 , nelem );
				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 = EmisLines[ipH_LIKE][nelem][4][ipH2s].Aul + 
				EmisLines[ipH_LIKE][nelem][4][ipH2p].Aul + EmisLines[ipH_LIKE][nelem][4][3].Aul;
			ans2 = (8.424e6+ 8.991e6)*Z4;
			if( fabs(ans1-ans2)/ans2 > 1e-2 )
			{
				fprintf(ioQQQ , 
					"SanityCheck finds insane A for summed H=4 decays found=%g correct=%g nelem=%li\n",
					ans1 , ans2 , nelem );
				lgOK = FALSE;
			}

#			endif
			if( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > 5 )
			{
				/* Balmer gamma 5 - 2 + 5-4 + 5-3*/
				ans1 = EmisLines[ipH_LIKE][nelem][5][ipH2s].Aul + 
					EmisLines[ipH_LIKE][nelem][5][ipH2p].Aul + 
					EmisLines[ipH_LIKE][nelem][5][3].Aul + EmisLines[ipH_LIKE][nelem][5][4].Aul;
				ans2 = (2.532e6+2.20e6+2.70e6)*Z4;
				if( fabs(ans1-ans2)/ans2 > 1e-2 )
				{
					fprintf(ioQQQ , 
						"SanityCheck finds insane A for H5-2 found=%g correct=%g nelem=%li\n",
						ans1 , ans2 , nelem );
					lgOK = FALSE;
				}
			}
		}
	}

	/* check that hydrogenic branching ratios add up to unity */
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		int ipHi, ipLo;
		for( ipHi=4; ipHi<=15; ++ipHi )
		{
			double sum = 0.;
			for( ipLo=2; ipLo<ipHi;++ipLo )
			{
				sum += HydroBranch(ipHi,ipLo,nelem+1);
			}
			if( fabs(sum-1.)>0.01 ) 
			{
				fprintf(ioQQQ , 
					"SanityCheck H branching ratio sum not unity for nelem=%li upper n=%i sum=%.3e\n",
					nelem, ipHi, sum );
				lgOK = FALSE;
			}
		}
	}

	/* check photo cross sections for H */
	ipISO = 0;
	nelem = 0;
	/* loop starts at 3, the first level with n = n and full l */
	for( n=3; n<MIN2(100., iso.numLevels_max[ipISO][nelem] ) ; ++n )
	{
		float anu[1]={1.} , cs[1] ;
		double energy;

		/* photon energy where cross section will be evaluated,
		 * this is in Ryd */
		energy = iso.xIsoLevNIonRyd[ipISO][nelem][n];
		anu[0] = (float)(energy*1.01);

		hypho(
			/* Z-Nelec, the residual charge, 1 for hydrogen, 2 for helium */
			(double)(nelem+1), 
			/* principle quantum number */
		  n, 
		  /* lowest angular momentum */
		  0, 
		  /* highest angular momentum */
		  n-1, 
		  /* scaled lowest photon energy, 
		   * => incorrect?? in units of zed^2/n^2, 
		   * at which the cs will be done */
		  energy, 
		  /* number of points to do */
		  1, 
		  /* array of energies (in units given above) */
		  anu , 
		  /* calculated photoionization cross section in cm^-2 */
		  cs);

		error = fabs(cs[0] - opac.OpacStack[iso.ipOpac[ipISO][nelem][n]-1] )/
			(        (cs[0] + opac.OpacStack[iso.ipOpac[ipISO][nelem][n]-1] )/2.);
		/*fprintf(ioQQQ,"z=%ld n=%ld error %g old %e new %e\n",nelem,n, error,
			opac.OpacStack[iso.ipOpac[ipISO][nelem][n]-1] ,cs[0] );*/
		if( error > 0.05 )
		{
			fprintf(ioQQQ , "SanityCheck finds insane H photo cs\n");
			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 , "SanityCheck 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 , "SanityCheck 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.;

	/* this is the default matrix solver */
	/* this test is the 1-d matrix with 2-d macro simulation */
	/* 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, " SanityCheck MALLOC A error\n" );
		puts( "[Stop in SanityCheck]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* copy over the main matrix */
	for(i=0;i<3;++i)
	{
		for( j=0;j<3;++j )
		{
			AA(i,j) = xMatrix[i][j];
		}
	}

	ner = 0;

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

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

	if( ner != 0 )
	{
		fprintf( ioQQQ, " SanityCheck DGETRS error\n" );
		puts( "[Stop in SanityCheck]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* release the vector */
	free( A );
#		undef AA

	/* 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,
			"SanityCheck 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 ContCreatePointers,
	 * then set to valid numbers in ipShells and OpacityCreate1Element
	 * this checks that all values have been properly filled */
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		/* must reset state of code after tests performed, remember state here */
		float xIonF[NISO][LIMELM];
		double hbn[NISO][LIMELM], hn[NISO][LIMELM];

		if( dense.lgElmtOn[nelem] )
		{
			/* set these abundances so that opacities can be checked below */
			hbn[ipH_LIKE][nelem] = iso.DepartCoef[ipH_LIKE][nelem][ipH1s];
			hn[ipH_LIKE][nelem] = iso.Pop2Ion[ipH_LIKE][nelem][ipH1s];
			xIonF[ipH_LIKE][nelem] = dense.xIonDense[nelem][nelem+1];

			iso.Pop2Ion[ipH_LIKE][nelem][ipH1s] = 1.;
			iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = 0.;
			dense.xIonDense[nelem][nelem+1] = 1.;

			if( nelem > ipHYDROGEN )
			{

				hbn[ipHE_LIKE][nelem] = iso.DepartCoef[ipHE_LIKE][nelem][ipHe1s1S];
				hn[ipHE_LIKE][nelem] = iso.Pop2Ion[ipHE_LIKE][nelem][ipHe1s1S];
				xIonF[ipHE_LIKE][nelem] = dense.xIonDense[nelem][nelem];

				/* this does not exist for hydrogen itself */
				iso.Pop2Ion[ipHE_LIKE][nelem][ipH1s] = 1.;
				iso.DepartCoef[ipHE_LIKE][nelem][ipH1s] = 0.;
				dense.xIonDense[nelem][nelem] = 1.;
			}

			for( ion=0; ion<=nelem; ++ion )
			{
				/* loop over all shells that are defined */
				for( nshells=0; nshells<Heavy.nsShells[nelem][ion]; ++nshells )
				{
					for( j=0; j<3; ++j )
					{
						/* >>chng 00 apr 05, array index is on fortran scale so must be
						 * >= 1.  This test had been <0, correct for C.  Caught by Peter van Hoof */
						if( opac.ipElement[nelem][ion][nshells][j] <=0 )
						{
							/* this is not possible */
							fprintf(ioQQQ,
								"SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n", 
								nelem , ion , nshells, j );
							fprintf(ioQQQ,
								"value was %li  \n", opac.ipElement[nelem][ion][nshells][j] );
							/* set flag saying that things are not ok */
							lgOK = FALSE;
						}
					}
				}

				if( nelem > 1 )
				{
					float saveion[LIMELM+3];
					/* check that photoionization cross sections are ok */
					for( j=1; j <= (nelem + 2); j++ )
					{
						saveion[j] = dense.xIonDense[nelem][j-1];
						dense.xIonDense[nelem][j-1] = 0.;
					}

					dense.xIonDense[nelem][ion] = 1.;

					OpacityZero();
					opac.lgRedoStatic = TRUE;

					/* generate opacity with standard routine - this is the one
					 * called in OpacityAddTotal to make opacities in usual calculations */
					OpacityAdd1Element(nelem);

					/* this starts one beyond energy of threshold since cs may be zero there */
					for( j=Heavy.ipHeavy[nelem][ion]; j < MIN2(rfield.nflux,KshllEnr.KshellLimit) ; j++ )
					{
						if( opac.opacity_abs[j]+opac.OpacStatic[j] < FLT_MIN )
						{
							/* this is not possible */
							fprintf(ioQQQ,
								"SanityCheck found non-positive photo cs for nelem=%li ion=%li \n", 
								nelem , ion );
							fprintf(ioQQQ,
								"value was %.2e + %.2e nelem %li ion %li at energy %.2e\n", 
								opac.opacity_abs[j] ,
								opac.OpacStatic[j] ,
								nelem , 
								ion , 
								rfield.anu[j]);
							/* set flag saying that things are not ok */
							lgOK = FALSE;
							break;/**/
						}
					}
					/* reset the ionization distribution */
					for( j=1; j <= (nelem + 2); j++ )
					{
						/*lint -e771 savion not init */
						dense.xIonDense[nelem][j-1] = saveion[j];
						/*lint +e771 savion not init */
					}

				}
			}
			iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = hbn[ipH_LIKE][nelem];
			iso.Pop2Ion[ipH_LIKE][nelem][ipH1s] = hn[ipH_LIKE][nelem];
			dense.xIonDense[nelem][nelem+1] = xIonF[ipH_LIKE][nelem];

			if( nelem > ipHYDROGEN )
			{
				iso.DepartCoef[ipHE_LIKE][nelem][ipHe1s1S] = hbn[ipHE_LIKE][nelem];
				/* >>chng 03 april 17, this Pop2Ion had H_LIKE instead of HE_LIKE */
                 iso.Pop2Ion[ipHE_LIKE][nelem][ipHe1s1S] = hn[ipHE_LIKE][nelem];
				dense.xIonDense[nelem][nelem] = xIonF[ipHE_LIKE][nelem];
			}
		}
	}


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

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

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

	else
	{
		/* stop since problem encountered, lgEOF set FALSE */
		fprintf(ioQQQ , "SanityCheck finds insanity so exiting\n");
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}
}
#undef NDIM
