/*PunHeat punch contributors to local heating, with punch heat command */
#include <string.h>
#include <float.h>
#include <limits.h>
#include "cddefines.h"
#include "coolants.h"
#include "heat.h"
#include "tcool.h"
#include "zonecnt.h"
#include "heating.h"
#include "converge.h"
#include "taulines.h"
#include "elementnames.h"
#include "wkhtcl.h"
#include "fndlineht.h"
#include "chlinelbl.h"
#include "punheat.h"
/* limit for number of heat agents that are saved */
#	ifdef IPNSAVE
#undef IPNSAVE
#	endif
#define	IPNSAVE	30
/* limit to number to print */
#	ifdef IPRINT
#undef IPRINT
#	endif
#define IPRINT 9

void PunHeat(FILE* io)
{
	char chLabel[IPNSAVE][5], 
	  chLbl[11];
	int lgHeatLine;
	long int i, 
	  ipStrong, 
	  ipnt, 
	  ipsave[IPNSAVE], 
	  j, 
	  jpsave[IPNSAVE], 
	  k, 
	  level;
	double CS, 
	  ColHeat, 
	  EscP, 
	  Pump, 
	  Strong, 
	  TauIn, 
	  save[IPNSAVE];

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

	for( i=0; i<IPNSAVE ; ++i )
	{
		ipsave[i] = INT_MIN;
		jpsave[i] = INT_MIN;
		save[i] = -DBL_MAX;
	}

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

	/* 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 < Coolants.ncltot; i++ )
	{
		if( Coolants.heatnt[i]/heat.htot > WkHtCl.WeakHeatCool )
		{
			j = Coolants.collam[i];
			/* force to save wavelength convention as printout */
			if( j > 100000 )
				j /= 10000;
			fprintf( io, " Negative coolant was %4.4s%6ld%10.2e\n", 
			  chCoolants.chClntLab[i], j, Coolants.heatnt[i]/heat.htot );
		}
	}

	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" );
	}


	/* following loop tries to identify strongest agents and turn to labels */
	for( k=0; k < ipnt; k++ )
	{
		/* generate labels that make sense in printout 
		 * heating is [j][i] */
		i = ipsave[k];
		j = jpsave[k];
		if( i >= j )
		{
			/* this is case of photoionization of atom or ion */
			strcpy( chLabel[k], elementnames.chElementSym[i] );
			strcat( chLabel[k], elementnames.chIonStage[j]  );
		}
		else if( i == 0 && j == 1 )
		{
			/* photoionization from all excited states of Hydrogenic species,
			 * heating[1][0] */
			strcpy( chLabel[k], "Hn=2" );
		}
		else if( i == 0 && j == 3 )
		{
			/* collisional ionization of H from ground */
			strcpy( chLabel[k], "Hion" );
		}
		else if( i == 0 && j == 9 )
		{
			/* CO dissociation */
			strcpy( chLabel[k], "COds" );
		}
		else if( i == 0 && j == 11 )
		{
			/* free free heating */
			strcpy( chLabel[k], "H FF" );
		}
		else if( i == 0 && j == 12 )
		{
			/* heating line, that was supposed to cool */
			strcpy( chLabel[k], "Clin" );
		}
		else if( i == 0 && j == 13 )
		{
			/* grain photoionization */
			strcpy( chLabel[k], "dust" );
		}
		else if( i == 0 && j == 15 )
		{
			/* H- heating */
			strcpy( chLabel[k], "H-  " );
		}
		else if( i == 0 && j == 17 )
		{
			/* H2 dissociation */
			strcpy( chLabel[k], "H2ds" );
		}
		else if( i == 0 && j == 18 )
		{
			/* Compton recoil of bound electrons */
			strcpy( chLabel[k], "CBnd" );
		}
		else if( i == 0 && j == 19 )
		{
			/* Compton heating */
			strcpy( chLabel[k], "Comp" );
		}
		else if( i == 0 && j == 22 )
		{
			/* line heating */
			strcpy( chLabel[k], "line" );
		}
		else if( i == 0 && j == 23 )
		{
			/* hydrogen line heating */
			strcpy( chLabel[k], "Hlin" );
		}
		else if( i == 1 && j == 3 )
		{
			/* helium triplet line heating */
			strcpy( chLabel[k], "He3l" );
		}
		else
		{
			sprintf( chLabel[k], "%2ld%2ld", i, j );
		}
	}

	/* begin the print out with zone number, total heating and cooling */
	fprintf( io, "%4ld%10.2e%10.2e", ZoneCnt.nzone, heat.htot, 
	  tcool.ctot );

	fprintf( io, " " );

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

	for( i=1; i <= ipnt; i++ )
	{
		fprintf( io, "%4.4s%6.3f ", chLabel[i-1], save[i-1] );
	}
	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(0,22) 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/heat.htot > 0.005 )
		{
			if( level == 1 )
			{
				strcpy( chLbl, chLineLbl(&TauLines[ipStrong-1] ) );
				TauIn = TauLines[ipStrong-1].TauIn;
				Pump = TauLines[ipStrong-1].pump;
				EscP = TauLines[ipStrong-1].Pesc;
				CS = TauLines[ipStrong-1].cs;
				/* ratio of line to total heating */
				ColHeat = TauLines[ipStrong-1].heat/
				  heat.htot;
			}
			else if( level == 2 )
			{
				strcpy( chLbl, chLineLbl(&TauLine2[ipStrong-1]) );
				TauIn = TauLine2[ipStrong-1].TauIn;
				Pump = TauLine2[ipStrong-1].pump;
				EscP = TauLine2[ipStrong-1].Pesc;
				CS = TauLine2[ipStrong-1].cs;
				ColHeat = TauLine2[ipStrong-1].heat/
				  heat.htot;
			}
			else
			{
				fprintf( ioQQQ, " PunHeat insane level%4ld\n", 
				  level );
				puts( "[Stop in punheat]" );
				exit(1);
			}
			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 );
		}
	}


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