/*punt produce map of heating-cooling space for specified zone */
#include <string.h>
#include <limits.h>
#include "cddefines.h"
#ifdef EPS
#	undef EPS
#endif
#define	EPS	0.005
#include "coolants.h"
#include "mappar.h"
#include "cnegtk.h"
#include "heating.h"
#include "stopcalc.h"
#include "printe82.h"
#include "he.h"
#include "thigh.h"
#include "called.h"
#include "heat.h"
#include "phycon.h"
#include "zonecnt.h"
#include "h.h"
#include "tcool.h"
#include "trace.h"
#include "coolpr.h"
#include "pressure.h"
#include "tfidle.h"
#include "boltgn.h"
#include "nsset.h"
#include "converge.h"
#include "punt.h"

void punt(char *chType)
{
	char chLabel[5];
	char units;
	long int i, 
	  ih, 
	  j, 
	  jh, 
	  k, 
	  lambda;
	double cfrac, 
	  ch, 
	  fac, 
	  factor, 
	  hfrac, 
	  oldch, 
	  ratio, 
	  strhet, 
	  strong, 
	  t1, 
	  tlowst, 
	  tmax, 
	  tmin, 
	  torg;

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

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

	/* make sure pressure has been evaluated */
	strong = TotalPressure();

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

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

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

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

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

	/* map out te-ionization-cooling space before punching out. */
	if( called.lgTalk )
	{
		fprintf( ioQQQ, " map of heating, cooling, vs temp, follows.\n");
		fprintf( ioQQQ, 
			"     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)(MapPar.nMapStep));
		phycon.te = (float)tmin;
	}

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

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

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

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

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

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

		/* now make the space */
		MapPar.temap = (float *)malloc( sizeof(float)*(MapPar.nMapStep+4) );
		if( MapPar.temap == NULL )
		{ 
			printf( " not enough memory to allocate MapPar.temap in punt\n" );
			puts( "[Stop in punt]" );
			exit(1);
		}
		MapPar.cmap = (float *)malloc( sizeof(float)*(MapPar.nMapStep+4) );
		if( MapPar.cmap == NULL )
		{ 
			printf( " not enough memory to allocate MapPar.cmap in punt\n" );
			puts( "[Stop in punt]" );
			exit(1);
		}
		MapPar.hmap = (float *)malloc( sizeof(float)*(MapPar.nMapStep+4) );
		if( MapPar.hmap == NULL )
		{ 
			printf( " not enough memory to allocate MapPar.hmap in punt\n" );
			puts( "[Stop in punt]" );
			exit(1);
		}

	}

	cnegtk.lgCNegChk = FALSE;
	MapPar.nmap = 0;
	oldch = 0.;
	phycon.te *= (float)factor;

	while( phycon.te < tmax*0.999 && phycon.te > tmin*1.001 )
	{
		/* strong means nothing, but we must give value to something*/
		strong = TotalPressure();
		/* must reset upper and lower bounds for ionization distributions */
		nsset();
		/* this turns on constant reevaluation of everything */
		conv.lgSearch = TRUE;

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

		/* reevaluate gaunt factors and Boltzmann factors */
		boltgn();

		/* 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 inital ionization solution for this temperature */
		ConvEdenIoniz();
		MapPar.temap[MapPar.nmap] = phycon.te;
		MapPar.cmap[MapPar.nmap] = (float)tcool.ctot;
		MapPar.hmap[MapPar.nmap] = (float)heat.htot;

		lambda = 0;
		strong = 0.;

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

		cfrac = strong/tcool.ctot;
		strhet = 0.;
		/* these will be reset in following loop*/
		ih = -INT_MAX;
		jh = -INT_MAX;

		for( k=0; k < LIMELM; k++ )
		{
			for( j=0; j < LIMELM; j++ )
			{
				if( HeatingCom.heating[j][k] > strhet )
				{
					strhet = HeatingCom.heating[j][k];
					ih = k+1;
					jh = j+1;
				}
			}
		}

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

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

			else
			{
				lambda /= 10000;
				units = 'm';
			}

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

			/*fprintf( ioQQQ, 
				"%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(ioQQQ,"%s", PrintEfmt("%11.4e", phycon.te ) );
			fprintf(ioQQQ,"%s", PrintEfmt("%11.4e", heat.htot ) );
			fprintf(ioQQQ,"%4ld%4ld%6.3f",
			  ih, jh, 
			  hfrac);
			fprintf(ioQQQ,"%s", PrintEfmt("%11.4e", tcool.ctot ) );
			fprintf(ioQQQ," %4.4s %4ld%c%6.3f",
			  chLabel , 
			  lambda, 
			  units, 
			  cfrac );
			fprintf(ioQQQ,"%s", PrintEfmt(" %10.2e", heat.dHTotDT ) );
			fprintf(ioQQQ,"%s", PrintEfmt("%11.2e", tcool.dCooldT ) );
			fprintf(ioQQQ,"%s", PrintEfmt("%11.4e", phycon.eden ) );
			fprintf(ioQQQ,"%s", PrintEfmt("%11.4e", phycon.hden ) );
			fprintf(ioQQQ,"%6.2f%6.2f%6.2f%6.2f\n",
			  log10(MAX2(1e-9,h.hii/phycon.hden)), 
			  log10(MAX2(1e-9,he.hei/he.ahe)), 
			  log10(MAX2(1e-9,he.heii/he.ahe)), 
			  log10(MAX2(1e-9,he.heiii/he.ahe)) );
		}

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

		{
			enum {DEBUG=FALSE};
			if( DEBUG )
			{
				static int kount = 0;
				factor = 1.;
				phycon.te = 8674900.;
				++kount;
				if( kount >=100 )
				{
					fprintf(ioQQQ," exiting in punt\n");
					break;
				}
			}
		}
	}

	/* now check whether sharp inflections occurred, and also find the biggest jump
	 * in the heating and cooling */
	MapPar.lgMapOK = TRUE;
	 for( i=1; i< MapPar.nmap-2; ++i )
	 {
		float s1,s2,s3;/* the three slopes we will use */
		s1 = MapPar.cmap[i-2] - MapPar.cmap[i-1];
		s2 = MapPar.cmap[i-1] - MapPar.cmap[i];
		s3 = MapPar.cmap[i] - MapPar.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 */
			fprintf( ioQQQ,
				" cooling curve had double inflection at T=%.2e, warning being generated.\n",
				MapPar.temap[i] );
			MapPar.lgMapOK = FALSE;
		}

		s1 = MapPar.hmap[i-2] - MapPar.hmap[i-1];
		s2 = MapPar.hmap[i-1] - MapPar.hmap[i];
		s3 = MapPar.hmap[i]   - MapPar.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( ioQQQ,
				" heating curve had double inflection at T=%.2e, warning being generated.\n",
				MapPar.temap[i] );
			MapPar.lgMapOK = FALSE;
		}
	 }

	cnegtk.lgCNegChk = TRUE;
	phycon.te = (float)t1;

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