/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * 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 "stopcalc.h"
#include "called.h"
#include "dense.h"
#include "co.h"
#include "phycon.h"
#include "trace.h"
#include "pressure.h"
#include "tfidle.h"
#include "converge.h"
#include "map.h"
#ifdef EPS
#	undef EPS
#endif
#define	EPS	0.005

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

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

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

	/* make sure pressure has been evaluated */
	/* this sets values of pressure.PresTotlCurr, also calls tfidle */
	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)/(float)(map.nMapStep));
		phycon.te = (float)tmin;
		tfidle(FALSE);
	}

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

		/* we want about nMapStep (=20) steps between high and low temperature */
		factor = pow(10.,log10(tmax/tmin)/(float)(map.nMapStep));

		if( thermal.lgTeHigh )
		{
			/* high te */
			factor = 1./factor;
			/* TeHighest is highest possible temperature, 1E10 */
			phycon.te = (float)(MIN2(map.RangeMap[1],StopCalc.TeHighest)/factor);
		}

		else
		{
			/* low te */
			phycon.te = (float)(tlowst/factor);
		}
		tfidle(FALSE);
	}

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

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

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

	}

	thermal.lgCNegChk = FALSE;
	map.nmap = 0;
	oldch = 0.;
	phycon.te = phycon.te *(float)factor;
	tfidle(TRUE);
	if( trace.lgTrConvg )
		fprintf(ioQQQ, "    MAP called temp range %.4e %.4e in %li stops ===============================================\n",
		tmin,
		tmax,
		map.nmap);

	while( (double)phycon.te < tmax*0.999 && (double)phycon.te > tmin*1.001 )
	{
		/* this sets values of pressure.PresTotlCurr, also calls tfidle */
		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.lgTrConvg )
			fprintf(ioQQQ, "    MAP new temp %.4e ===============================================\n",
			phycon.te );

		/* update density - temperature variables which may have changed
		tfidle(FALSE); */

		/* 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;

		/* this will force reinit of co network */
		co.iteration_co = -2;

		/* now get ionization solution for this temperature */
		ConvEdenIoniz();

		/* save results for later prints */
		map.temap[map.nmap] = phycon.te;
		map.cmap[map.nmap] = (float)thermal.ctot;
		map.hmap[map.nmap] = (float)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,"%s", PrintEfmt("%11.4e", phycon.te ) );
			fprintf(io,"%s", PrintEfmt("%11.4e", thermal.htot ) );
			fprintf(io," [%2ld][%2ld]%6.3f",
			  ksav, jsav,  
			  hfrac);
			fprintf(io,"%s", PrintEfmt("%11.4e", thermal.ctot ) );
			fprintf(io," %s %.1f%c%6.3f",
			  chLabel , 
			  wl, 
			  units, 
			  cfrac );
			fprintf(io,"%s", PrintEfmt(" %10.2e", thermal.dHeatdT ) );
			fprintf(io,"%s", PrintEfmt("%11.2e", thermal.dCooldT ) );
			fprintf(io,"%s", PrintEfmt("%11.4e", dense.eden ) );
			fprintf(io,"%s", 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);
		}

		phycon.te = phycon.te*(float)factor;
		tfidle(TRUE);
		/* increment nmap but do not exceed nMapAlloc */
		map.nmap = MIN2(map.nMapAlloc,map.nmap+1);

		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				static int kount = 0;
				factor = 1.;
				phycon.te = 8674900.f;
				tfidle(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 */
	map.lgMapOK = TRUE;
	/* >>chng 02 mar 04, lower bound had been 1, so [i-2] below was negative */
	 for( i=2; i< map.nmap-2; ++i )
	 {
		float s1,s2,s3;/* the three slopes we will use */
		s1 = map.cmap[i-2] - map.cmap[i-1];
		s2 = map.cmap[i-1] - map.cmap[i];
		s3 = map.cmap[i] - map.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.  ",
				map.temap[i]);
			fprintf( io,	" Slopes were %.2e %.2e %.2e",	s1, s2, s3);
			if( fabs(s2)/map.cmap[i] > 0.05 )
			{
				fprintf( io,
					" error large, (rel slope of %.2e).\n",
					s2 / map.cmap[i]);
				map.lgMapOK = FALSE;
			}
			else
			{
				fprintf( io,
					" error is small, (rel slope of %.2e).\n",
					s2 / map.cmap[i]);
			}
		}

		s1 = map.hmap[i-2] - map.hmap[i-1];
		s2 = map.hmap[i-1] - map.hmap[i];
		s3 = map.hmap[i]   - map.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",
				map.temap[i] );
			map.lgMapOK = FALSE;
		}
	 }

	thermal.lgCNegChk = TRUE;
	phycon.te = (float)t1;
	tfidle(FALSE);

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

	return;
}

