/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*map_do produce map of heating-cooling space for specified zone, called as result of
 * map command  */
#include "cddefines.h"
#include "thermal.h"
#include "cooling.h"
#include "called.h"
#include "dense.h"
#include "mole.h"
#include "phycon.h"
#include "trace.h"
#include "pressure.h"
#include "conv.h"
#include "hcmap.h"
#ifdef EPS
#	undef EPS
#endif
#define	EPS	0.005

t_hcmap hcmap;

void map_do(
			FILE *io, 
			const char *chType)
{
	char chLabel[NCOLNT_LAB_LEN+1];
	char units;
	long int i, 
	  ksav, 
	  j, 
	  jsav, 
	  k,
	  nelem;
	realnum wl;
	double cfrac, 
	  ch, 
	  fac, 
	  factor, 
	  hfrac, 
	  oldch, 
	  ratio, 
	  strhet, 
	  strong, 
	  t1, 
	  tlowst, 
	  tmax, 
	  tmin, 
	  torg;

	DEBUG_ENTRY( "map_do()" );

	t1 = phycon.te;
	torg = phycon.te;
	hcmap.lgMapOK = true;
	/* flag indicating that we have computed a map */
	hcmap.lgMapDone = true;

	/* make sure pressure has been evaluated */
	/* this sets values of pressure.PresTotlCurr */
	PresTotCurrent();

	/* print out all coolants if all else fails */
	if( called.lgTalk )
	{
		fprintf( io, " Cloudy punts, Te=%10.3e HTOT=%10.3e CTOT=%10.3e nzone=%4ld\n", 
		  phycon.te, thermal.htot, thermal.ctot, nzone );
		fprintf( io, " COOLNG array is\n" );

		if( thermal.ctot > 0. )
		{
			coolpr(io, "ZERO",1,0.,"ZERO");
			for( i=0; i < thermal.ncltot; i++ )
			{
				ratio = thermal.cooling[i]/thermal.ctot;
				if( ratio>EPS )
				{
					coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
					  ratio,"DOIT");
				}
			}

			coolpr(io, "DONE",1,0.,"DONE");
			fprintf( io, " Line heating array follows\n" );
			coolpr(io, "ZERO",1,0.,"ZERO");

			for( i=0; i < thermal.ncltot; i++ )
			{
				ratio = thermal.heatnt[i]/thermal.ctot;
				if( ratio>EPS )
				{
					coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
					  ratio,"DOIT");
				}
			}

			coolpr(io,"DONE",1,0.,"DONE");
		}
	}

	/* map out te-ionization-cooling space before punching out. */
	if( called.lgTalk )
	{
		fprintf( io, " map of heating, cooling, vs temp, follows.\n");
		fprintf( io, 
			"     Te     Heat--------------------> Cool--------------------->   dH/dT     dC/DT       Ne         NH      HII       Helium \n" );
	}

	if( strcmp(chType,"punt") == 0 )
	{
		/* this is the original use of punt, we are punting
		 * only map within factor of two of final temperature
		 * fac will the range to either side of punted temperature */
		fac = 1.5;
		tmin = torg/fac;
		tmax = torg*fac;

		/* we want about 20 steps between high and low temperature
		 * default of nMapStep is 20, set with set nmaps command */
		factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
		TempChange(tmin , false);
	}

	else if( strcmp(chType," map") == 0 )
	{
		/* create some sort of map of heating-cooling */
		tlowst = MAX2(hcmap.RangeMap[0],phycon.TEMP_LIMIT_LOW);
		tmin = tlowst*0.998;
		tmax = MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)*1.002;

		/* we want about nMapStep (=20) steps between high and low temperature */
		factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
		double TeNew;
		if( thermal.lgTeHigh )
		{
			/* high te */
			factor = 1./factor;
			/* TeHighest is highest possible temperature, 1E10 */
			TeNew = (MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)/factor);
		}

		else
		{
			/* low te */
			TeNew = (tlowst/factor);
		}
		TempChange(TeNew , false);
	}

	else
	{
		/* don't know what to do */
		fprintf( ioQQQ, " PUNT called with insane argument,=%4.4s\n", 
		  chType );
		cdEXIT(EXIT_FAILURE);
	}

	/* now allocate space for te, c, h vectors in map, if not already done */
	if( hcmap.nMapAlloc==0 )
	{
		/* space not allocated, do so now */
		hcmap.nMapAlloc = hcmap.nMapStep+4;

		/* now make the space */
		hcmap.temap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
		if( hcmap.temap == NULL )
		{ 
			printf( " not enough memory to allocate hcmap.temap in map_do\n" );
			cdEXIT(EXIT_FAILURE);
		}
		hcmap.cmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
		if( hcmap.cmap == NULL )
		{ 
			printf( " not enough memory to allocate hcmap.cmap in map_do\n" );
			cdEXIT(EXIT_FAILURE);
		}
		hcmap.hmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
		if( hcmap.hmap == NULL )
		{ 
			printf( " not enough memory to allocate hcmap.hmap in map_do\n" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	thermal.lgCNegChk = false;
	hcmap.nmap = 0;
	oldch = 0.;
	TempChange(phycon.te *factor , true);
	if( trace.nTrConvg )
		fprintf(ioQQQ, "    MAP called temp range %.4e %.4e in %li stops ===============================================\n",
		tmin,
		tmax,
		hcmap.nmap);

	while( (double)phycon.te < tmax*0.999 && (double)phycon.te > tmin*1.001 )
	{
		/* this sets values of pressure.PresTotlCurr */
		PresTotCurrent();

		/* must reset upper and lower bounds for ionization distributions */
		/* fix number of stages of ionization */
		for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				dense.IonLow[nelem] = 0;
				/* 
				* IonHigh[n] is the highest stage of ionization present
				* the IonHigh array index is on the C scake, so [0] is hydrogen
				* the value is also on the C scale, so element [nelem] can range
				* from 0 to nelem+1 
				*/
				dense.IonHigh[nelem] = nelem + 1;
			}
			else
			{
				/* this element is turned off, make stages impossible */
				dense.IonLow[nelem] = -1;
				dense.IonHigh[nelem] = -1;
			}
		}

		/* this turns on constant reevaluation of everything */
		conv.lgSearch = true;

		if( trace.nTrConvg )
			fprintf(ioQQQ, "    MAP new temp %.4e ===============================================\n",
			phycon.te );

		/* this counts how many times ionize is called in this model after startr,
		 * and is flag used by ionize to understand it is being called the first time*/
		conv.nTotalIoniz = 0;

		/* now get ionization solution for this temperature */
		if( ConvEdenIoniz() )
			lgAbort = true;

		/* save results for later prints */
		hcmap.temap[hcmap.nmap] = (realnum)phycon.te;
		hcmap.cmap[hcmap.nmap] = (realnum)thermal.ctot;
		hcmap.hmap[hcmap.nmap] = (realnum)thermal.htot;

		wl = 0.f;
		strong = 0.;

		for( j=0; j < thermal.ncltot; j++ )
		{
			if( thermal.cooling[j] > strong )
			{
				strcpy( chLabel, thermal.chClntLab[j] );
				strong = thermal.cooling[j];
				wl = thermal.collam[j];
			}
		}

		cfrac = strong/thermal.ctot;
		strhet = 0.;
		/* these will be reset in following loop*/
		ksav = -INT_MAX;
		jsav = -INT_MAX;

		for( k=0; k < LIMELM; k++ )
		{
			for( j=0; j < LIMELM; j++ )
			{
				if( thermal.heating[k][j] > strhet )
				{
					strhet = thermal.heating[k][j];
					jsav = j;
					ksav = k;
				}
			}
		}

		ch = thermal.ctot - thermal.htot;
		/* use ratio to check for change of sign since product
		 * can underflow at low densities */
		if( oldch/ch < 0. && called.lgTalk )
		{
			fprintf( io, " ----------------------------------------------- Probable thermal solution here. --------------------------------------------\n" );
		}

		oldch = ch;
		hfrac = strhet/thermal.htot;
		if( called.lgTalk )
		{
			/* convert to micros if IR transition */
			if( wl < 100000.f )
			{
				units = 'A';
			}

			else
			{
				wl /= 10000.f;
				units = 'm';
			}

			if( trace.lgTrace )
			{
				fprintf( io, "TRACE: te, htot, ctot%11.3e%11.3e%11.3e\n", 
				  phycon.te, thermal.htot, thermal.ctot );
			}

			/*fprintf( io, 
				"%10.4e%11.4e%4ld%4ld%6.3f%11.4e %4.4s %4ld%c%6.3f%10.2e%11.4e%11.4e%6.2f%6.2f%6.2f%6.2f\n",;*/
			fprintf(io, PrintEfmt("%11.4e", phycon.te ) );
			fprintf(io, PrintEfmt("%11.4e", thermal.htot ) );
			fprintf(io," [%2ld][%2ld]%6.3f",
			  ksav, jsav,  
			  hfrac);
			fprintf(io, PrintEfmt("%11.4e", thermal.ctot ) );
			fprintf(io," %s %.1f%c%6.3f",
			  chLabel , 
			  wl, 
			  units, 
			  cfrac );
			fprintf(io, PrintEfmt(" %10.2e", thermal.dHeatdT ) );
			fprintf(io, PrintEfmt("%11.2e", thermal.dCooldT ) );
			fprintf(io, PrintEfmt("%11.4e", dense.eden ) );
			fprintf(io, PrintEfmt("%11.4e", dense.gas_phase[ipHYDROGEN] ) );
			if( dense.lgElmtOn[ipHELIUM] )
			{
				fprintf(io,"%6.2f%6.2f%6.2f%6.2f\n",
				log10(MAX2(1e-9,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])), 
				log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][0]/dense.gas_phase[ipHELIUM])), 
				log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM])), 
				log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][2]/dense.gas_phase[ipHELIUM])) );
			}
			fflush(io);
		}

		TempChange(phycon.te*factor , true);
		/* increment nmap but do not exceed nMapAlloc */
		hcmap.nmap = MIN2(hcmap.nMapAlloc,hcmap.nmap+1);

		{
			enum {DEBUG_LOC=false};
			if( DEBUG_LOC )
			{
				static int kount = 0;
				factor = 1.;
				TempChange(8674900. , true);
				++kount;
				if( kount >=100 )
				{
					fprintf(ioQQQ," exiting in map_do\n");
					break;
				}
			}
		}
	}

	/* now check whether sharp inflections occurred, and also find the biggest jump
	 * in the heating and cooling */
	hcmap.lgMapOK = true;
	/* >>chng 02 mar 04, lower bound had been 1, so [i-2] below was negative */
	 for( i=2; i< hcmap.nmap-2; ++i )
	 {
		realnum s1,s2,s3;/* the three slopes we will use */
		s1 = hcmap.cmap[i-2] - hcmap.cmap[i-1];
		s2 = hcmap.cmap[i-1] - hcmap.cmap[i];
		s3 = hcmap.cmap[i] - hcmap.cmap[i+1];
		if( s1*s3 > 0. && s2*s3 < 0. )
		{
			 /* of the three points, the outer had the same slope 
			  * (their product was positive) but there was an inflection
			  * between them (the negative product).  The data chain looked like
			  *     2 4
			  *    1 3  or vice versa, either case is wrong,
			  * with the logic in the above test, the problem point will aways be s2 */
			fprintf( io,
				" cooling curve had double inflection at T=%.2e.  ",
				hcmap.temap[i]);
			fprintf( io,	" Slopes were %.2e %.2e %.2e",	s1, s2, s3);
			if( fabs(s2)/hcmap.cmap[i] > 0.05 )
			{
				fprintf( io,
					" error large, (rel slope of %.2e).\n",
					s2 / hcmap.cmap[i]);
				hcmap.lgMapOK = false;
			}
			else
			{
				fprintf( io,
					" error is small, (rel slope of %.2e).\n",
					s2 / hcmap.cmap[i]);
			}
		}

		s1 = hcmap.hmap[i-2] - hcmap.hmap[i-1];
		s2 = hcmap.hmap[i-1] - hcmap.hmap[i];
		s3 = hcmap.hmap[i]   - hcmap.hmap[i+1];
		if( s1*s3 > 0. && s2*s3 < 0. )
		{
			 /* of the three points, the outer had the same slope 
			  * (their product was positive) but there was an inflection
			  * between them (the negative product).  The data chain looked like
			  *     2 4
			  *    1 3  or vice versa, either case is wrong */
			fprintf( io,
				" heating curve had double inflection at T=%.2e.\n",
				hcmap.temap[i] );
			hcmap.lgMapOK = false;
		}
	 }

	thermal.lgCNegChk = true;
	TempChange(t1 , false);
	return;
}