/* 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 */
/*HeatPunch punch contributors to local heating, with punch heat command, called by punch_do */
#include "cddefines.h"
#include "thermal.h"
#include "radius.h"
#include "conv.h"
#include "lines_service.h"
#include "dense.h"
#include "taulines.h"
#include "elementnames.h"
#include "dynamics.h"
#include "punch.h"

/* limit for number of heat agents that are saved */
/* limit to number to print */
#	ifdef IPRINT
#undef IPRINT
#	endif
#define IPRINT 9

/*HeatPunch punch contributors to local heating, with punch heat command, called by punch_do */
void HeatPunch(FILE* io)
{
	char **chLabel, 
	  chLbl[11];
	int lgHeatLine,
		lgFail;
	long int i, 
	  ipStrong, 
	  ipnt, 
	  *ipOrdered,
	  *ipsave, 
	  j, 
	  *jpsave, 
	  k, 
	  level;
	double CS, 
	  ColHeat, 
	  EscP, 
	  Pump, 
	  Strong, 
	  TauIn,
		coolok,
		heatok;
	float *save;

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

	save = (float *) CALLOC(LIMELM*LIMELM, sizeof(float));
	ipsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
	jpsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
	ipOrdered = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
	chLabel = (char **) CALLOC(LIMELM*LIMELM, sizeof(char *));

	for( i=0; i<LIMELM*LIMELM; ++i )
	{
		ipsave[i] = INT_MIN;
		jpsave[i] = INT_MIN;
		save[i] = -FLT_MAX;
		chLabel[i] = (char *) CALLOC(10, sizeof(char));
	}

	coolok = thermal.ctot;
	heatok = thermal.htot;
	if (dynamics.Cool > dynamics.Heat) 
	{
		coolok -= dynamics.Heat;
		heatok -= dynamics.Heat;
	} 
	else
	{
		coolok -= dynamics.Cool;
		heatok -= dynamics.Cool;
	}

	ipnt = 0;

	/* heat sources are saved in a 2d square array */
	/* WeakHeatCool set with 'set weakheatcool' command
	 * default is 0.05 */
	for( i=0; i < LIMELM; i++ )
	{
		for( j=0; j < LIMELM; j++ )
		{
			if( thermal.heating[i][j]/heatok > 0. )
			{
				ipsave[ipnt] = i;
				jpsave[ipnt] = j;
				save[ipnt] = (float)(thermal.heating[i][j]/heatok);
				ipnt++;
			}
		}
	}

	/* now check for possible line heating not counted in 1,23
	 * there should not be any significant heating source heree
	 * since they would not be counted in deriv correctly */
	for( i=0; i < thermal.ncltot; i++ )
	{
		if( thermal.heatnt[i]/heatok > punch.WeakHeatCool )
		{
			float awl;
			awl = thermal.collam[i];
			/* force to save wavelength convention as printout */
			if( awl > 100000 )
				awl /= 10000;
			fprintf( io, " Negative coolant was %s %.2f %.2e\n", 
			  thermal.chClntLab[i], awl, thermal.heatnt[i]/heatok );
		}
	}

	if( !conv.lgConvTemp )
	{
		fprintf( io, " >>>>  Temperature notconverged.\n" );
	}
	else if( !conv.lgConvEden )
	{
		fprintf( io, " >>>>  Electron density not convergenced.\n" );
	}
	else if( !conv.lgConvIoniz )
	{
		fprintf( io, " >>>>  Ionization not converged.\n" );
	}
	else if( !conv.lgConvPres )
	{
		fprintf( io, " >>>>  Pressure not converged.\n" );
	}

	/* this is mainly to keep the compiler from flagging possible paths with j not being set */
	i = INT_MIN;
	j = INT_MIN;
	/* following loop tries to identify strongest agents and turn to labels */
	for( k=0; k < ipnt; k++ )
	{
		/* generate labels that make sense in printout 
		 * if not identified with a specific agent, will print indices as [i][j],
		 * heating is thermal.heating[i][j] */
		i = ipsave[k];
		j = jpsave[k];
		if( i >= j )
		{
			if( dense.xIonDense[i][j] == 0.)
				fprintf(ioQQQ,"hit it\n");
			ASSERT( dense.xIonDense[i][j] > 0.);
			/* this is case of photoionization of atom or ion */
			strcpy( chLabel[k], elementnames.chElementSym[i] );
			strcat( chLabel[k], elementnames.chIonStage[j]  );
		}
		/* notice that in test i and j are swapped from order in heating array */
		else if( i == 0 && j == 1 )
		{
			/* photoionization from all excited states of Hydrogenic species,
			 * heating[0][1] */
			strcpy( chLabel[k], "Hn=2" );
		}
		else if( i == 0 && j == 3 )
		{
			/* collisional ionization of H from ground, heating[0][3] */
			strcpy( chLabel[k], "Hion" );
		}
		else if( i == 0 && j == 8 )
		{
			/* thermal.heating[0][8] is heating due to collisions within X of H2,
			 * code var hmi.HeatH2Dexc_used, hmi.HeatH2Dexc_BigH2, hmi.HeatH2Dexc_TH85 */
			strcpy( chLabel[k], "H2vH" );
		}
		else if( i == 0 && j == 17 )
		{
			/* thermal.heating[0][17] is heating due to photodissociation heating of X within H2,
			 * code var hmi.HeatH2Dish_used */
			strcpy( chLabel[k], "H2dH" );
		}
		else if( i == 0 && j == 9 )
		{
			/* CO dissociation, co.CODissHeat, heating[0][9] */
			strcpy( chLabel[k], "COds" );
		}
		else if( i == 0 && j == 20 )
		{
			/* extra heat thermal.heating[0][20]*/
			strcpy( chLabel[k], "extH" );
		}
		else if( i == 0 && j == 11 )
		{
			/* free free heating, heating[0][11] */
			strcpy( chLabel[k], "H FF" );
		}
		else if( i == 0 && j == 12 )
		{
			/* heating line, that was supposed to cool, heating[0][12] */
			strcpy( chLabel[k], "Clin" );
		}
		else if( i == 0 && j == 13 )
		{
			/* grain photoionization, heating[0][13] */
			strcpy( chLabel[k], "GrnP" );
		}
		else if( i == 0 && j == 14 )
		{
			/* grain collisions, heating[0][14] */
			strcpy( chLabel[k], "GrnC" );
		}
		else if( i == 0 && j == 15 )
		{
			/* H- heating, heating[0][15] */
			strcpy( chLabel[k], "H-  " );
		}
		else if( i == 0 && j == 16 )
		{
			/* H2+ heating, heating[0][16] */
			strcpy( chLabel[k], "H2+ " );
		}
		else if( i == 0 && j == 19 )
		{
			/* Compton heating, heating[0][19] */
			strcpy( chLabel[k], "Comp" );
		}
		else if( i == 0 && j == 22 )
		{
			/* line heating, heating[0][22] */
			strcpy( chLabel[k], "line" );
		}
		else if( i == 0 && j == 23 )
		{
			/* iso-sequence line heating - all elements together, heating[0][23] */
			strcpy( chLabel[k], "Hlin" );
		}
		else if( i == 0 && j == 24 )
		{
			/* charge transfer heating, heating[0][24] */
			strcpy( chLabel[k], "ChaT" );
		}
		else if( i == 1 && j == 3 )
		{
			/* helium triplet line heating, heating[1][3] */
			strcpy( chLabel[k], "He3l" );
		}
		else if( i == 1 && j == 5 )
		{
			/* advective heating, heating[1][5] */
			strcpy( chLabel[k], "adve" );
		}
		else if( i == 1 && j == 6 )
		{
			/* cosmic ray heating thermal.heating[1][6]*/
			strcpy( chLabel[k], "CR H" );
		}
		else if( i == 25 && j == 27 )
		{
			/* Fe 2 line heating, heating[25][27] */
			strcpy( chLabel[k], "Fe 2" );
		}
		else
		{
			sprintf( chLabel[k], "[%ld][%ld]" , i , j );
		}
	}

	/* now sort by decreasing importance */
	/*spsort netlib routine to sort array returning sorted indices */
	spsort(
		  /* input array to be sorted */
		  save, 
		  /* number of values in x */
		  ipnt, 
		  /* permutation output array */
		  ipOrdered, 
		  /* flag saying what to do - 1 sorts into increasing order, not changing
		   * the original routine */
		  -1, 
		  /* error condition, should be 0 */
		  &lgFail);

	/* begin the print out with zone number, total heating and cooling */
	fprintf( io, "%.5e\t%.4e\t%.4e", 
		radius.depth_mid_zone, 
		heatok, 
		coolok );

	/* we only want to print the IPRINT strongest of the group */
	ipnt = MIN2( ipnt , IPRINT );

	for( k=0; k < ipnt; k++ )
	{
		int ip = ipOrdered[k];
		i = ipsave[ip];
		j = jpsave[ip];
		ASSERT( i<LIMELM && j<LIMELM );
		/*lint -e662 out of bounds array index */
		if (k > 4 && thermal.heating[i][j]/heatok < punch.WeakHeatCool)
			break;
		/*lint +e662 out of bounds array index */
		fprintf( io, "\t%s\t%.7f ", 
			chLabel[ip], save[ip] );
	}
	fprintf( io, " \n" );

	/* a negative pointer in the heating array is probably a problem,
	 * indicating that some line has become a heat source */
	lgHeatLine = FALSE;

	/* check if any lines were major heat sources */
	for( i=0; i < ipnt; i++ )
	{
		/* heating[22][0] is line heating - identify line if important */
		if( ipsave[i] == 0 && jpsave[i] == 22 )
			lgHeatLine = TRUE;
	}

	if( lgHeatLine )
	{
		/* a line was a major heat source - identify it */
		FndLineHt(&level,&ipStrong,&Strong);
		if( Strong/heatok > 0.005 )
		{
			if( level == 1 )
			{
				strcpy( chLbl, chLineLbl(&TauLines[ipStrong] ) );
				TauIn = TauLines[ipStrong].TauIn;
				Pump = TauLines[ipStrong].pump;
				EscP = TauLines[ipStrong].Pesc;
				CS = TauLines[ipStrong].cs;
				/* ratio of line to total heating */
				ColHeat = TauLines[ipStrong].heat/
				  heatok;
			}
			else if( level == 2 )
			{
				strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]) );
				TauIn = TauLine2[ipStrong].TauIn;
				Pump = TauLine2[ipStrong].pump;
				EscP = TauLine2[ipStrong].Pesc;
				CS = TauLine2[ipStrong].cs;
				ColHeat = TauLine2[ipStrong].heat/
				  heatok;
			}
			else if( level == 3 )
				/* C12O16 */
			{
				strcpy( chLbl, chLineLbl(&C12O16Rotate[ipStrong]) );
				TauIn = C12O16Rotate[ipStrong].TauIn;
				Pump = C12O16Rotate[ipStrong].pump;
				EscP = C12O16Rotate[ipStrong].Pesc;
				CS = C12O16Rotate[ipStrong].cs;
				ColHeat = C12O16Rotate[ipStrong].heat/
				  heatok;
			}
			else if( level == 4 )
				/* C13O16 */
			{
				strcpy( chLbl, chLineLbl(&C13O16Rotate[ipStrong]) );
				TauIn = C13O16Rotate[ipStrong].TauIn;
				Pump = C13O16Rotate[ipStrong].pump;
				EscP = C13O16Rotate[ipStrong].Pesc;
				CS = C13O16Rotate[ipStrong].cs;
				ColHeat = C13O16Rotate[ipStrong].heat/
				  heatok;
			}
			else if( level == 5 )
				/* hyperfine levels */
			{
				strcpy( chLbl, chLineLbl(&HFLines[ipStrong]) );
				TauIn = HFLines[ipStrong].TauIn;
				Pump = HFLines[ipStrong].pump;
				EscP = HFLines[ipStrong].Pesc;
				CS = HFLines[ipStrong].cs;
				ColHeat = HFLines[ipStrong].heat/
				  heatok;
			}
			else
			{
				fprintf( ioQQQ, " HeatPunch insane level%4ld\n", 
				  level );
				puts( "[Stop in HeatPunch]" );
				cdEXIT(EXIT_FAILURE);
			}
			fprintf( io, "  LHeat lv%2ld %10.10s TIn%10.2e Pmp%9.1e EscP%9.1e CS%9.1e Hlin/tot%10.2e\n", 
			  level, chLbl, TauIn, Pump, EscP, CS, ColHeat );
		}
	}
	for( i=0; i<LIMELM*LIMELM; ++i )
	{
		free(chLabel[i]);
	}

	free(chLabel);
	free(ipOrdered);
	free(jpsave);
	free(ipsave);
	free(save);

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

#undef IPRINT

