/*InitAssertResults, this must be called first, done at startup of ParseCommands*/
/*ParseAssertResults - parse input stream */
/*lgCheckAsserts, checks asserts, last thing cloudy calls, returns TRUE if all are ok, FALSE if problems */
#include <string.h>
#include <time.h>
#include "cddefines.h"
#include "input.h"
#include "getquote.h"
#include "ffmtread.h"
#include "nonumb.h"
#include "showme.h"
#include "converge.h"
#include "called.h"
#include "mappar.h"
#include "secondaries.h"
#include "radius.h"
#include "hydrogenic.h"
#include "zonecnt.h"
#include "heat.h"
#include "phycon.h"
#include "elmton.h"
#include "grainvar.h"
#include "readelement.h"
#include "cddrive.h"
#include "itercnt.h"
#include "elementnames.h"
#include "lgmatch.h"
#include "assertresults.h"

/* flag to remember that InitAssertResults was called */
static int lgInitDone=FALSE , 
	/* will be set true when space for asserts is allocated */
	lgSpaceAllocated=FALSE;

/* number of asserts we can handle, used in dim of space */
#define NASSERTS 100

/* default relative error for asserted quantities */
#define DEF_ERROR 0.05

/* dim of 5 also appears in malloc below */
#define NCHAR 5
static char **chAssertLineLabel;

/* this will be = for equal, < or > for limit */
static char *chAssertLimit;

/* this will be a two character label identifying which type of assert */
static char **chAssertType;

/* the values and error in the asserted quantity */
static double *AssertQuantity ,*AssertQuantity2 ,*AssertError;

/* this flag says where we print linear or log quantity */
static int *lgQuantityLog;
static long nAsserts=0 , *nWLAsserts ;

/*======================================================================*/
/*InitAssertResults, this must be called first, done at startup of ParseCommands*/
void InitAssertResults(void)
{
	/* set flag that init was called, and set number of asserts to zero.
	 * this is done by ParseComments for every model, even when no asserts will
	 * be done, so do not allocate space at this time */
	lgInitDone = TRUE;
	nAsserts = 0;
}

/*======================================================================*/
/* parse the assert command */
void ParseAssertResults(void)
{
	long i ,
		nelem;
	int lgEOL;
	char chLabel[80];

	if( !lgInitDone )
	{
		fprintf( ioQQQ, " ParseAssertResults called before InitAsserResults\n" );
		puts( "[Stop in ParseAssertResults]" );
		exit(1);
	}

	/* has space been allocated yet? */
	if( !lgSpaceAllocated )
	{
		/* remember that space has been allocated */
		lgSpaceAllocated = TRUE;

		/* create space for the array of labels*/
		if( (chAssertLineLabel = ((char **)malloc(NASSERTS*sizeof(char *)))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc chAssertLineLabel\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* the 2-character string saying what type of assert */
		if( (chAssertType = ((char **)malloc(NASSERTS*sizeof(char *)))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc chAssertType\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* now convert the two previous into strings */
		for( i=0; i<NASSERTS; ++i )
		{
			if( (chAssertLineLabel[i] = ((char *)malloc(NCHAR*sizeof(char )))) == NULL )
			{
				fprintf(ioQQQ," ParseAssertResults could not malloc2 chAssertLineLabel\n");
				puts( "[Stop in ParseAssertResults]" );
				exit(1);
			}
			strcpy(chAssertLineLabel[i],"unkn" );

			/* two char plus eol */
			if( (chAssertType[i] = ((char *)malloc(3*sizeof(char )))) == NULL )
			{
				fprintf(ioQQQ," ParseAssertResults could not malloc2 chAssertType\n");
				puts( "[Stop in ParseAssertResults]" );
				exit(1);
			}
			strcpy(chAssertType[i],"un" );
		}

		/* now make space for the asserted quantities  */
		if( (AssertQuantity = ((double *)malloc(NASSERTS*sizeof(double )))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc AssertQuantity\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}
		if( (AssertQuantity2 = ((double *)malloc(NASSERTS*sizeof(double )))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc AssertQuantity2\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* now the errors */
		if( (AssertError = ((double *)malloc(NASSERTS*sizeof(double )))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc AssertError\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* now the line wavelengths */
		if( (nWLAsserts = ((long *)malloc(NASSERTS*sizeof(long )))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc nWLAsserts\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* now the flag saying whether should be log */
		if( (lgQuantityLog = ((int *)malloc(NASSERTS*sizeof(int )))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc lgQuantityLog\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* the flag for upper, lower limit, or equal */
		if( (chAssertLimit = ((char *)malloc(NASSERTS*sizeof(char )))) == NULL )
		{
			fprintf(ioQQQ," ParseAssertResults could not malloc chAssertLimit\n");
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}
	}
	/* end space allocation - we are ready to roll */

	/* first say we do not know what job to do */
	strcpy( chAssertType[nAsserts] , "un" );

	/* false means print linear quantity - will be set true if entered
	 * quantity comes in as a log */
	lgQuantityLog[nAsserts] = FALSE;

	/* all asserts have option for quantity to be a limit, or the quantity itself */
	if( lgMatch("<",input.chCARDCAPS ) )
	{
		chAssertLimit[nAsserts] = '<';
	}
	else if( lgMatch(">",input.chCARDCAPS ) )
	{
		chAssertLimit[nAsserts] = '>';
	}
	else
	{
		chAssertLimit[nAsserts] = '=';
	}

	/* which quantity will we check?, first is */

	/* assert mean ionization */
	if( lgMatch("IONI",input.chCARDCAPS ) )
	{

		/* say that this will be mean ionization fraction */

		/* f will indicate average over radius, F over volumn -
		 * check whether keyword radius or volume occurs,
		 * default will be radius */
		if( lgMatch("VOLU",input.chCARDCAPS ) )
		{
			strcpy( chAssertType[nAsserts] , "F " );
		}
		else
		{
			/* this is default case, Fraction over radius */
			strcpy( chAssertType[nAsserts] , "f " );
		}

		/* first get element label and make null terminated string*/
		if( (nelem = ReadElement()) < 0 )
		{
			fprintf( ioQQQ, 
				" I could not identify an element name on this line.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem])
		 * into array that will be used to get ionization after calculation */
		strcpy( chAssertLineLabel[nAsserts], elementnames.chElementNameShort[nelem] );

		/* now get ionization stage, which will be saved into nWLAsserts */
		i = 5;
		nWLAsserts[nAsserts] = 
			(long)FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */
		if( nWLAsserts[nAsserts] < 1 || nWLAsserts[nAsserts] > nelem+2 )
		{
			fprintf( ioQQQ, 
				"  The ionization stage is inappropriate for this element.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* now get ionization fraction, log if number is negative or == 0, 
		 * linear if positive but less than or equal to 1.*/
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}

		/* now make sure we end up with positive linear ionization fraction that
		 * is greater than 0 but less than or equal to 1. */
		if( AssertQuantity[nAsserts] <= 0. )
		{
			/* log since negative or 0 */
			AssertQuantity[nAsserts] = 
				pow(10.,AssertQuantity[nAsserts] ) ;
			/* entered as a log, so print as a log too */
			lgQuantityLog[nAsserts] = TRUE;
		}
		else if( AssertQuantity[nAsserts] > 1. )
		{
			fprintf( ioQQQ, 
				"  The ionization fraction must be less than one.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
	}

	/* assert line  --  key is ine_ since linear option appears on some commands */
	else if( lgMatch("INE ",input.chCARDCAPS ) )
	{
		if(  lgMatch("LUMI",input.chCARDCAPS ) )
		{
			/* say that this is a line luminosity */
			strcpy( chAssertType[nAsserts] , "Ll" );

			/* entered as a log, so print as a log too */
			lgQuantityLog[nAsserts] = TRUE;
		}
		else
		{
			/* say that this is an intensity - this is the default */
			strcpy( chAssertType[nAsserts] , "Li" );

			/* entered as a log, so print as a log too */
			lgQuantityLog[nAsserts] = FALSE;
		}

		/* this will check a line intensity, first get labels within quotes,
		 * will be null terminated */
		GetQuote( chLabel , input.chCARDCAPS );
		/* check that there were exactly 4 characters in the label*/
		if( chLabel[4] != '\0' )
		{
			fprintf( ioQQQ, 
				" The label must be exactly 4 char long, between double quotes.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}
		/* copy string into array */
		strcpy( chAssertLineLabel[nAsserts], chLabel );

		/* now get line wavelength */
		i = 5;
		nWLAsserts[nAsserts] = 
			(long)FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}

		/* now get intensity or luminosity - intensity is linear and luminosity is log */
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* luminosity was entered as a log */
		if( lgQuantityLog[nAsserts] )
		{
			AssertQuantity[nAsserts] = 
				pow(10.,AssertQuantity[nAsserts] ) ;
		}

		/* optional error, default available */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
	}

	/* assert departure coefficients */
	else if( lgMatch("DEPA",input.chCARDCAPS ) )
	{
		/* for now only do case where atom is in STE */
		AssertQuantity[nAsserts] = 1.;
		/* label for all cases is just dept */
		strcpy( chAssertLineLabel[nAsserts] , "dept" );
		/* this is relative error, max departure from unity of any level */
		i = 5;
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}

		/* H-like key means do one of the hydrogenic ions */
		if( lgMatch("H-LI",input.chCARDCAPS ) )
		{
			/* remember this is departure coef for some element */
			strcpy( chAssertType[nAsserts] , "D " );
			/* now get element number for h ion */
			if( (	nWLAsserts[nAsserts] = ReadElement()) < 0 )
			{
				fprintf( ioQQQ, 
					" I could not identify an element name on this line.\n");
				fprintf( ioQQQ, " Sorry.\n" );
				puts( "[Stop in ParseAssertResults]" );
				exit(1);
			}
		}

		/* this is the large FeII ion */
		else if( lgMatch("FEII",input.chCARDCAPS ) )
		{
			/* remember this is departure coef for feii */
			strcpy( chAssertType[nAsserts] , "d " );
		}
		else
		{
			fprintf( ioQQQ, 
				" There must be a second key, FEII or h-like followed by element.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}
	}

	/* assert some results from map */
	else if( lgMatch(" MAP",input.chCARDCAPS ) )
	{

		/* must have heating or cooling, since will check one or the other */
		/* check heating cooling results from map at some temperature */
		if( lgMatch("HEAT",input.chCARDCAPS ) )
		{
			strcpy( chAssertType[nAsserts] , "mh" );
		}
		else if( lgMatch("COOL",input.chCARDCAPS ) )
		{
			strcpy( chAssertType[nAsserts] , "mc" );
		}
		else
		{
			fprintf( ioQQQ, 
				" There must be a second key, HEATing or COOLing.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* no ionization stage*/
		nWLAsserts[nAsserts] = 0;

		/* i was set above for start of scan */
		/* now get temperature for AssertQuantity2 array*/
		i = 5;
		AssertQuantity2[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}

		if( AssertQuantity2[nAsserts] <= 10. )
		{
			/* entered as log, but we will compare with linear */
			AssertQuantity2[nAsserts] = 
				pow(10.,AssertQuantity2[nAsserts] ) ;
		}

		/* heating or cooling, both log, put into error */
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}

		/* AssertQuantity array will have heating or cooling */
		AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts]);

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}

		/* entered as a log, so print as a log too */
		lgQuantityLog[nAsserts] = TRUE;
	}

	/* assert column density of something */
	else if( lgMatch("COLU",input.chCARDCAPS ) )
	{

		/* this is default case, Fraction over radius */
		strcpy( chAssertType[nAsserts] , "cd" );

		/* we want to remember where the match occurred within the string
		 * since do not want to count the 2 as the first number */
		i=lgMatch(" H2 ",input.chCARDCAPS );
		if( i )
		{
			strcpy( chAssertLineLabel[nAsserts], "H2  " );
			/* increment to get past the label */
			i += 3;
		}
		else if( lgMatch(" CO ",input.chCARDCAPS ) )
		{
			strcpy( chAssertLineLabel[nAsserts], "CO  " );
			i = 5;
		}
		else
		{
			fprintf( ioQQQ, 
				" I could not identify CO or H2 on this line.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* this says to look for molecular column density, also could be ion stage */
		nWLAsserts[nAsserts] = 0;

		/* no ionization stage*/
		nWLAsserts[nAsserts] = 0;

		/* i was set above for start of scan */
		/* now get log of column density */
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* entered as log, but we will compare with linear */
		AssertQuantity[nAsserts] = 
			pow(10.,AssertQuantity[nAsserts] ) ;

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
		/* entered as a log, so print as a log too */
		lgQuantityLog[nAsserts] = TRUE;
	}

	/* assert mean temperature */
	else if( lgMatch("TEMP",input.chCARDCAPS ) )
	{
		/* say that this will be mean temperature, electron or grain */

		/* t will indicate temperature average over radius, T over volume -
		 * check whether keyword radius or volume occurs,
		 * default will be radius */
		if( lgMatch("VOLU",input.chCARDCAPS ) )
		{
			strcpy( chAssertType[nAsserts] , "T ");
		}
		else
		{
			/* this is default case, Fraction over radius */
			strcpy( chAssertType[nAsserts] , "t ");
		}

		/* first look for keyword Grains, since label silicate may be on it,
		 * and this would trigger the element search */
		if( lgMatch("GRAI",input.chCARDCAPS ) )
		{
			/* grains, copy 4-char string "grai"
			 * into array that will be used to get temperature after calculation */
			strcpy( chAssertLineLabel[nAsserts], "GRAI" );
			/* this is to make sure that pointer to grain type is valid, we check
			 * that is is less than this below */
			nelem = NDUST-2;
			/* the first numerical ary on the temper grain command is the grain
			 * type as defined in hazy, counting from 1.  When it is used to
			 * to get mean temps below we will subtract 1 from this to get onto
			 * the c scale.  but must leave on physical scale here to pass sanity
			 * checks that occur below */
		}
		/* look for element label within quotes,
		 * will be null terminated */
		else if( (nelem = ReadElement()) >= 0 )
		{
			/* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem])
			 * into array that will be used to get ionization after calculation */
			strcpy( chAssertLineLabel[nAsserts], elementnames.chElementNameShort[nelem] );
		}
		else
		{
			fprintf( ioQQQ, 
				" I could not identify an element name on this line.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* now get ionization stage (or grain type), which will be saved into nWLAsserts */
		i = 5;
		nWLAsserts[nAsserts] = 
			(long)FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */
		if( nWLAsserts[nAsserts] < 1 || nWLAsserts[nAsserts] > nelem+2 )
		{
			fprintf( ioQQQ, 
				"  The ionization stage is inappropriate for this element.\n");
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		/* now get temperature, log if number is <= 10, else linear */
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}

		/* now make sure we end up with positive linear temperature
		 * number is log if <=10 unless linear keyword appears */
		if( AssertQuantity[nAsserts] <= 10. && !lgMatch( "LINE" ,input.chCARDCAPS ) )
		{
			/* log since negative or 0 */
			AssertQuantity[nAsserts] = 
				pow(10.,AssertQuantity[nAsserts] ) ;
			/* entered as a log, so print as a log too */
			lgQuantityLog[nAsserts] = TRUE;
		}

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
	}

	/* assert (probably upper limit to) number of zones */
	else if( lgMatch("NZON",input.chCARDCAPS ) )
	{
		/* this flag will mean number of zones */
		strcpy( chAssertType[nAsserts] , "z " );
		/* say that this is zone numbering */
		strcpy( chAssertLineLabel[nAsserts], "zone" );

		/* now get number of zones, which will be saved into nWLAsserts */
		i = 5;
		nWLAsserts[nAsserts] = 
			(long)FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* copy it here although we will not use it */
		AssertQuantity[nAsserts] = (double)nWLAsserts[nAsserts];

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
	}

	/* assert secondary ionization rate, csupra */
	else if( lgMatch("CSUP",input.chCARDCAPS ) )
	{
		/* this flag will mean secondary ionization, entered as log */
		strcpy( chAssertType[nAsserts] , "sc" );
		/* say that this is sec ioniz */
		strcpy( chAssertLineLabel[nAsserts], "sion" );

		/* now get rate, saved into assert quantity */
		i = 5;
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* entered as log, make linear */
		AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts] );

		/* no wavelength */
		nWLAsserts[nAsserts] = 0;

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}

		/* we want to print the log of eden, not linear value */
		lgQuantityLog[nAsserts] = TRUE;
	}

	/* assert heating rate, erg/cm3/s, htot */
	else if( lgMatch("HTOT",input.chCARDCAPS ) )
	{
		/* this flag will mean heating, entered as log */
		strcpy( chAssertType[nAsserts] , "ht" );

		/* say that this is heating rate */
		strcpy( chAssertLineLabel[nAsserts], "heat" );

		/* now get rate, saved into assert quantity */
		i = 5;
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* entered as log, make linear */
		AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts] );

		/* no wavelength */
		nWLAsserts[nAsserts] = 0;

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}

		/* we want to print the log of the heating, not linear value */
		lgQuantityLog[nAsserts] = TRUE;
	}

	/* assert number of iterations per zone, a test of convergence */
	else if( lgMatch("ITRZ",input.chCARDCAPS ) )
	{
		/* this flag will mean number of iterations per zone */
		strcpy( chAssertType[nAsserts] , "iz" );
		/* say that this is iterations per zone  */
		strcpy( chAssertLineLabel[nAsserts], "itrz" );

		/* now get quantity */
		i = 5;
		AssertQuantity[nAsserts] = 
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* wavelength is meaningless */
		nWLAsserts[nAsserts] = 0;

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
	}

	/* assert electron density */
	else if( lgMatch("EDEN",input.chCARDCAPS ) )
	{
		/* this flag will mean electron density */
		strcpy( chAssertType[nAsserts] , "e " );
		/* say that this is electron density */
		strcpy( chAssertLineLabel[nAsserts], "eden" );

		/* now get electron density, which is a log */
		i = 5;
		AssertQuantity[nAsserts] = 
			pow(10., FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL) );
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
		nWLAsserts[nAsserts] = 0;

		/* we want to print the log of eden, not linear value */
		lgQuantityLog[nAsserts] = TRUE;
	}

	/* assert thickness of model */
	else if( lgMatch("THIC",input.chCARDCAPS ) )
	{
		/* this flag will mean electron density */
		strcpy( chAssertType[nAsserts] , "th" );
		/* say that this is electron density */
		strcpy( chAssertLineLabel[nAsserts], "thic" );

		/* now get thickness, which is a log */
		i = 5;
		AssertQuantity[nAsserts] = 
			pow(10., FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL) );
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}

		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
		nWLAsserts[nAsserts] = 0;

		/* we want to print the log of eden, not linear value */
		lgQuantityLog[nAsserts] = TRUE;
	}

	/* assert (probably upper limit to) number of iterations */
	else if( lgMatch("NITE",input.chCARDCAPS ) )
	{
		/* this flag will mean number of iterations */
		strcpy( chAssertType[nAsserts] , "Z " );
		/* say that this is iteration numbering */
		strcpy( chAssertLineLabel[nAsserts], "iter" );

		/* now get number of iterations, which will be saved into nWLAsserts */
		i = 5;
		nWLAsserts[nAsserts] = 
			(long)FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			NoNumb(input.chOrgCard);
		}
		/* copy it here although we will not use it */
		AssertQuantity[nAsserts] = (double)nWLAsserts[nAsserts];


		/* optional error, default available (cannot do befor loop since we
		 * do not know how many numbers are on line */
		AssertError[nAsserts] =
			FFmtRead(input.chCARDCAPS ,&i, 76,&lgEOL);
		if( lgEOL )
		{
			/* default error was set in define above */
			AssertError[nAsserts] = DEF_ERROR;
		}
	}
	else
	{
		/* did not recognize a command */
		fprintf( ioQQQ, 
			"  Unrecognized command.  The line image was==%s==\n", 
		  input.chOrgCard );
		fprintf( ioQQQ, 
			"  The options I know about are: lines, ionization, departure coef, column, "
			"temperature, nzone, csupre, htot, itrz, eden, niter, \n");
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in ParseAssertResults]" );
		exit(1);
	}

	/* increment number of asserts and confirm that limit not exceeded */
	++nAsserts;
	if( nAsserts >= NASSERTS )
	{
		fprintf(ioQQQ,
			" ParseAssertResults: too many asserts, limit is NASSERT=%d\n",
			NASSERTS );
		puts( "[Stop in ParseAssertResults]" );
		exit(1);
	}
}

/*============================================================================*/
/*lgCheckAsserts, checks asserts, last thing cloudy calls, returns TRUE if all are ok, FALSE if problems */
int lgCheckAsserts(
	/* this is the file we will write this to, usually standard output, 
	 * but can be punch */
	FILE * ioASSERT )
{
	double PredQuan[NASSERTS] , RelError[NASSERTS];
	double relint , absint ;
	int lg1OK[NASSERTS];
	long i,j;

	/* this is a global variable in assertresults.h, and can be checked by
	 * other routines to see if asserts are ok - (most runs will not use asserts) */
	lgAssertsOK = TRUE;

	/* first check all asserts, there probably are none */
	for(i=0; i<nAsserts; ++i )
	{
		lg1OK[i] = TRUE;
		PredQuan[i] = 0.;
		RelError[i] = 0.;

		/* which type of assert? */
		/* is it intensity? */
		if( strcmp( chAssertType[i] , "Li")==0 )
		{
			/* this is an intensity, get the line, returns false if could not find it */
			if( !cdLine( chAssertLineLabel[i] , nWLAsserts[i] , &relint,&absint ) )
			{
				fprintf( ioASSERT, 
					" assert error: lgCheckAsserts could not find a line with label %s %ld \n",
					chAssertLineLabel[i] , nWLAsserts[i] );
				fprintf( ioASSERT, 
					" assert error: The punch line labels command is a good way to get a list of line labels.\n\n");
				/* go to next line */
				lg1OK[i] = FALSE;
				RelError[i] = 0;
				PredQuan[i] = 0;
				continue;
			}
			PredQuan[i] = relint;
			RelError[i] = 1. -  PredQuan[i]/AssertQuantity[i]  ;

		}

		/*this is line luminosity */
		else if( strcmp( chAssertType[i] , "Ll")==0 )
		{
			/* this is a luminosity, get the line, returns false if could not find it */
			if( !cdLine( chAssertLineLabel[i] , nWLAsserts[i] , &relint,&absint ) )
			{
				fprintf( ioASSERT, 
					" assert error: lgCheckAsserts could not find a line with label %s %ld \n",
					chAssertLineLabel[i] , nWLAsserts[i] );
				fprintf( ioASSERT, 
					" assert error: The punch line labels command is a good way to get a list of line labels.\n\n");
				/* go to next line */
				lg1OK[i] = FALSE;
				RelError[i] = 0;
				PredQuan[i] = 0;
				continue;
			}
			/* absint was returned as the log */
			PredQuan[i] = pow(10.,absint);
			RelError[i] = 1. -  PredQuan[i]/AssertQuantity[i]  ;

		}
		else if( strcmp( chAssertType[i] , "z " ) == 0)
		{
			/* this is (probably an upper limit) to the number of zones */
			PredQuan[i] = ZoneCnt.nzone;
			/* two integers in difference */
			RelError[i] = nWLAsserts[i] - ZoneCnt.nzone ;
		}

		else if( strcmp( chAssertType[i] , "iz") == 0 )
		{
			/* this is number of iterations per zone, a test of convergence properties */
			PredQuan[i] = (double)(conv.nTotalIoniz)/(double)(ZoneCnt.nzone);
			/* this is relative error */
			RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
		}

		else if( strcmp( chAssertType[i] , "e ") == 0 )
		{
			/* this is electron density */
			PredQuan[i] = phycon.eden;
			/* this is relative error */
			RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
		}

		else if( strcmp( chAssertType[i] , "th") == 0 )
		{
			/* this is thickness */
			PredQuan[i] = radius.depth;
			/* this is relative error */
			RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
		}

		else if( strcmp( chAssertType[i] , "sc") == 0 )
		{
			/* this is secondary ionization rate */
			PredQuan[i] = Secondaries.csupra;
			/* this is relative error */
			RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
		}

		else if( strcmp( chAssertType[i] , "ht") == 0 )
		{
			/* this is heating rate */
			PredQuan[i] = heat.htot;
			/* this is relative error */
			RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
		}

		else if( strcmp( chAssertType[i] , "Z ") == 0 )
		{
			/* this is (probably an upper limit) to the number of iterations */
			PredQuan[i] = IterCnt.iter;
			/* two integers in difference */
			RelError[i] = nWLAsserts[i] - IterCnt.iter ;
		}

		else if( strcmp( chAssertType[i] , "D ") == 0 )
		{
			long ipZ , n;
			/* this is departure coef for H-like ion given by nWLAsserts */
			/* stored number was element number on C scale */
			ipZ = nWLAsserts[i];
			if( !elmton.lgElmtOn[ipZ] )
			{
				fprintf(ioQQQ," asserted element %ld is not turned on!\n",ipZ);
				PredQuan[i] = ZoneCnt.nzone;
				RelError[i] = nWLAsserts[i] - ZoneCnt.nzone ;
			}
			else
			{
				RelError[i] = 0.;
				PredQuan[i] = 0.;
				for( n=IP1S; n<=nhlevel; ++n )
				{
					PredQuan[i] += hydro.hbn[ipZ][n];
					/* relerror is the largest deviation from unity */
					RelError[i] = MAX2( RelError[i] , hydro.hbn[ipZ][n] -1.);
					RelError[i] = MAX2( RelError[i] , 1. - hydro.hbn[ipZ][n] );
				}
				PredQuan[i] /= (double)(nhlevel);
			}
		}
		else if( strcmp( chAssertType[i] , "d ") == 0 )
		{
			/* this is depature coef for FeII */
			fprintf(ioQQQ," FeII not yet in \n");
			PredQuan[i] = ZoneCnt.nzone;
			RelError[i] = nWLAsserts[i] - ZoneCnt.nzone ;
		}

		/* this would be ionization fraction */
		else if( (strcmp( chAssertType[i] , "f ") == 0) || 
			(strcmp( chAssertType[i] , "F ") == 0) )
		{
			char chWeight[7];
			if( strcmp( chAssertType[i] , "F ") == 0 )
			{
				strcpy( chWeight , "VOLUME" );
			}
			else
			{
				/* this is default case, Fraction over radius */
				strcpy( chWeight , "RADIUS" );
			}
			/* get ionization fraction, returns false if could not find it */
			if( !cdIonFrac(
				/* four char string, null terminzed, giving the element name */
				chAssertLineLabel[i], 
				/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
				nWLAsserts[i], 
				/* will be fractional ionization */
				&relint, 
				/* how to weight the average, must be "VOLUME" or "RADIUS" */
				chWeight )  )
			{
				fprintf( ioASSERT, 
					" assert error: lgCheckAsserts could not find a line with label %s %ld \n",
					chAssertLineLabel[i] , nWLAsserts[i] );
				fprintf( ioASSERT, 
					" assert error: The punch line labels command is a good way to get a list of line labels.\n\n");
				/* go to next line */
				lg1OK[i] = FALSE;
				RelError[i] = 0;
				PredQuan[i] = 0;
				continue;
			}
			/* this is ionization fraction */
			PredQuan[i] = relint;
			RelError[i] = 1. -  PredQuan[i]/AssertQuantity[i]  ;
		}

		/* this would be column density of CO or H2 */
		else if( strcmp( chAssertType[i] , "cd") == 0) 
		{
			/* get ionization fraction, returns false if could not find it */
			if( !cdColm(
				/* four char string, null terminzed, giving the element name */
				chAssertLineLabel[i], 
				/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
				 * zero for molecule*/
				nWLAsserts[i], 
				/* will be fractional ionization */
				&relint) )
			{
				fprintf( ioASSERT, 
					" assert error: lgCheckAsserts could not find a molecule with label %s %ld \n",
					chAssertLineLabel[i] , nWLAsserts[i] );
				fprintf( ioASSERT, 
					" assert error: The punch line labels command is a good way to get a list of line labels.\n\n");
				/* go to next line */
				lg1OK[i] = FALSE;
				RelError[i] = 0;
				PredQuan[i] = 0;
				continue;
			}
			/* this is ionization fraction */
			PredQuan[i] = relint;
			RelError[i] = 1. -  PredQuan[i]/AssertQuantity[i]  ;
		}

		/* check heating/cooling at some temperature in a thermal map */
		else if( strcmp( chAssertType[i] , "mh") == 0 || strcmp( chAssertType[i] , "mc") == 0) 
		{
			/* check heating or cooling (stored in error array) at temperature in assert results */
			/* check that map was done, and arrays have nmap elements */
			if( MapPar.nMapAlloc == 0 )
			{
				/* this happens if map not done and space for h/c not allocated */
				fprintf( ioASSERT, 
					" assert error: lgCheckAsserts cannot check map since map not done.\n");
				/* go to next line */
				lg1OK[i] = FALSE;
				RelError[i] = 0;
				PredQuan[i] = 0;
				continue;
			}
			/* now check that requested temperature is within the range of the map we computed */
			if( AssertQuantity2[i]< MapPar.temap[0] || AssertQuantity2[i]> MapPar.temap[MapPar.nmap-1] )
			{
				fprintf( ioASSERT, 
					" assert error: lgCheckAsserts cannot check map since temperature not within range.\n");
				/* go to next line */
				lg1OK[i] = FALSE;
				RelError[i] = 0;
				PredQuan[i] = 0;
				continue;
			}

			/* we should have valid data - find closest tempeature >- requested temperature */
			j = 0;
			while( AssertQuantity2[i]>MapPar.temap[j]*1.001 && j < MapPar.nmap )
			{
				++j; 
			}

			/* j points to correct cell in heating cooling array */
			/* we will not interpolate, just use this value, and clobber te to prove it*/
			if( strcmp( chAssertType[i] , "mh") == 0 )
			{
				/* heating */
				PredQuan[i] = MapPar.hmap[j];
				strcpy(chAssertLineLabel[i],"MapH" );
			}
			else if( strcmp( chAssertType[i] , "mc") == 0)
			{
				/* cooling */
				PredQuan[i] = MapPar.cmap[j];
				strcpy(chAssertLineLabel[i],"MapC" );
			}
			RelError[i] = 1. -  PredQuan[i]/AssertQuantity[i]  ;
		}

		/* this will be an average temperatre */
		else if( (strcmp( chAssertType[i] , "t ") == 0) || 
			(strcmp( chAssertType[i] , "T ") == 0) )
		{
			char chWeight[7];
			if( strcmp( chAssertType[i] , "T ") == 0 )
			{
				strcpy( chWeight , "VOLUME" );
			}
			else
			{
				/* this is default case, Fraction over radius */
				strcpy( chWeight , "RADIUS" );
			}

			/* two options are average Te for ion or temp for grain */
			if( strcmp( chAssertLineLabel[i], "GRAI" ) == 0 )
			{
				/* the minus one is because the grain types are counted from one,
				 * but stuffed into the c array, that counts from zero */
				relint = GrainVar.avdust[nWLAsserts[i]-1]/GrainVar.reftot;
			}
			else
			{
				/* get temperature, returns false if could not find it */
				if( !cdTemp(
					/* four char string, null terminzed, giving the element name */
					chAssertLineLabel[i], 
					/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
					nWLAsserts[i], 
					/* will be mean temperatue */
					&relint, 
					/* how to weight the average, must be "VOLUME" or "RADIUS" */
					chWeight )  )
				{
					fprintf( ioASSERT, 
						" assert error: lgCheckAsserts could not find an ion with label %s %ld \n",
						chAssertLineLabel[i] , nWLAsserts[i] );
					fprintf( ioASSERT, 
						" assert error: The punch line labels command is a good way to get a list of line labels.\n\n");
					/* go to next line */
					lg1OK[i] = FALSE;
					RelError[i] = 0;
					PredQuan[i] = 0;
					continue;
				}
			}
			/* this is the temperature */
			PredQuan[i] = relint;
			RelError[i] = 1. -  PredQuan[i]/AssertQuantity[i]  ;
		}

		else
		{
			fprintf( ioASSERT, 
				" assert error: lgCheckAsserts received an insane chAssertType=%s, impossible\n",
				chAssertType[i] );
			ShowMe();
			puts( "[Stop in ParseAssertResults]" );
			exit(1);
		}

		if( chAssertLimit[i] == '=' )
		{
			/* predicted quantity should be within error of expected */
			if( fabs(RelError[i]) > AssertError[i] )
			{
				lg1OK[i] = FALSE;
				lgAssertsOK = FALSE;
			}
		}
		else if( chAssertLimit[i] == '<' )
		{
			/* expected is an upper limit, so PredQuan/AssertQuantity should
			 * be less than one, and so RelError should be postive */
			if( RelError[i] <= 0.-AssertError[i])
			{
				lg1OK[i] = FALSE;
				lgAssertsOK = FALSE;
			}
		}
		else if( chAssertLimit[i] == '>' )
		{
			/* expected is a lower limit, so PredQuan/AssertQuantity should
			 * be greater than one, and so RelError should be negative */
			if( RelError[i] >= 0.+AssertError[i])
			{
				lg1OK[i] = FALSE;
				lgAssertsOK = FALSE;
			}
		}
	}

	/* only print summary if we are talking */
	if( called.lgTalk && nAsserts>0 )
	{ 
		char chVer[10];
		time_t now;
		cdVersion(chVer ) ;
		/* write start of title and version number of code */
		fprintf( ioASSERT, "=============Results of asserts: C%s ",chVer );
		/* now add date of this run */
		now = time(NULL);
		/* now print this time at the end of the string.  the system put cr at the end of the string */
		fprintf(ioASSERT,"%s", ctime(&now) );

		if( lgAssertsOK )
		{
			fprintf( ioASSERT, " No errors were found.  Summary follows.\n");
		}
		else
		{
			fprintf( ioASSERT, " Errors were found.  Summary follows.\n");
		}

		fprintf( ioASSERT, 
			"                  Label line  computed    asserted Rel Err Set err\n");
		/* now print a summary */
		/*lint -e771 vars conceivably not initialized */
		for( i=0; i<nAsserts; ++i )
		{
			double prtPredQuan, prtAssertQuantity;
			/* this is option to print log of quantity rather than linear.  default is
			 * linear, and log only for a few such as ionization to see small numbers */
			if( lgQuantityLog[i] )
			{
				prtPredQuan = log10( MAX2( SMALLFLOAT , PredQuan[i] ) );
				prtAssertQuantity = log10( MAX2( SMALLFLOAT , AssertQuantity[i] ) );
			}
			else
			{
				prtPredQuan = PredQuan[i] ;
				prtAssertQuantity = AssertQuantity[i];
			}
			/* start of line, flag problems */
			if( lg1OK[i] )
			{
				/* the ChkAssert is a unique label so that we can grep on all outpu
				 * and see what happened, without picking up input stream */
				fprintf( ioASSERT, " ChkAssert        ");
			}
			else
			{
				fprintf( ioASSERT, " ChkAssert botch>>");
			}
			fprintf( ioASSERT , "%4s %5ld %9.3f %c %9.3f %7.3f %7.3f \n",
				chAssertLineLabel[i] , nWLAsserts[i] ,
				prtPredQuan ,  chAssertLimit[i]  , prtAssertQuantity , 
				/*PredQuan[i] ,  chAssertLimit[i]  , AssertQuantity[i] , */
				RelError[i] , AssertError[i]);
		}
		/*lint +e771 vars conceivably not initialized */
		fprintf( ioASSERT , " \n");
		if( !lgAssertsOK )
		{
			fprintf( ioASSERT, "  BOTCHED ASSERTS!!!   Botched asserts!!! \n\n");
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->PrtFinal()\n", debug_fp );
#	endif
	/* return the value */
	return(lgAssertsOK);
}