/* 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 */
/*SaveHeat save contributors to local heating, with save 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 "phycon.h"
#include "elementnames.h"
#include "dynamics.h"
#include "save.h"

/* limit for number of heat agents that are saved */
/* limit to number to print */
static const int IPRINT= 9;

/*SaveHeat save contributors to local heating, with save heat command, 
 * called by punch_do */
void SaveHeat(FILE* io)
{
	char **chLabel, 
	  chLbl[11];
	bool lgHeatLine;
	int nFail;
	long int i, 
	  ipnt, 
	  *ipOrdered,
	  *ipsave, 
	  j, 
	  *jpsave, 
	  k; 
	double CS, 
	  ColHeat, 
	  EscP, 
	  Pump, 
	  TauIn,
		cool_total,
		heat_total;
	realnum *SaveVal;

	DEBUG_ENTRY( "SaveHeat()" );

	SaveVal = (realnum *) CALLOC(LIMELM*LIMELM, sizeof(realnum));
	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;
		SaveVal[i] = -FLT_MAX;
		chLabel[i] = (char *) CALLOC(10, sizeof(char));
	}

	cool_total = thermal.ctot;
	heat_total = thermal.htot;

	/* >>chng 06 mar 17, comment out following block and replace with this 
	 * removing dynamics heating & cooling and report only physical
	 * heating and cooling 
	 * NB the heating and cooling as punched no longer need be
	 * equal for a converged model */
	cool_total -= dynamics.Cool();
	heat_total -= dynamics.Heat();
#	if 0
	if(dynamics.Cool() > dynamics.Heat()) 
	{
		cool_total -= dynamics.Heat();
		heat_total -= dynamics.Heat();
	} 
	else
	{
		cool_total -= dynamics.Cool();
		heat_total -= dynamics.Cool();
	}
#	endif

	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]/SDIV(heat_total) > SMALLFLOAT )
			{
				ipsave[ipnt] = i;
				jpsave[ipnt] = j;
				SaveVal[ipnt] = (realnum)(thermal.heating[i][j]/SDIV(heat_total));
				ipnt++;
			}
		}
	}

	/* now check for possible line heating not counted in 1,23
	 * there should not be any significant heating source here
	 * since they would not be counted in derivative correctly */
	for( i=0; i < thermal.ncltot; i++ )
	{
		if( thermal.heatnt[i]/SDIV(heat_total) > save.WeakHeatCool )
		{
			realnum 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]/SDIV(heat_total) );
		}
	}

	if( !conv.lgConvTemp )
	{
		fprintf( io, "#>>>>  Temperature not converged.\n" );
	}
	else if( !conv.lgConvEden )
	{
		fprintf( io, "#>>>>  Electron density not converged.\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];
		/* i >= j indicates agent is one of the first LIMELM elements */
		if( i >= j )
		{
			if( dense.xIonDense[i][j] == 0. && thermal.heating[i][j]>SMALLFLOAT )
				fprintf(ioQQQ,"DISASTER assert about to be thrown - search for hit it\n");
			/*fprintf(ioQQQ,"DEBUG %li %li %.2e %.2e\n", i , j , 
				dense.xIonDense[i][j],
				thermal.heating[i][j]);*/
			ASSERT( dense.xIonDense[i][j] > 0. || thermal.heating[i][j]<SMALLFLOAT );
			/* 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 all iso-seq from all levels, 
			 * heating[0][3] */
			strcpy( chLabel[k], "Hion" );
		}
		else if( i == 0 && j == 7 )
		{
			/* UTA ionization heating[0][7] */
			strcpy( chLabel[k], " UTA" );
		}
		else if( i == 0 && j == 8 )
		{
			/* thermal.heating[0][8] is heating due to collisions within 
			 * X of H2, code var hmi.HeatH2Dexc_used */
			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 == 21 )
		{
			/* pair heating thermal.heating[0][21]*/
			strcpy( chLabel[k], "pair" );
		}
		else if( i == 0 && j == 11 )
		{
			/* free free heating, heating[0][11] */
			strcpy( chLabel[k], "H FF" );
		}
		else if( i == 0 && j == 12 )
		{
			/* heating coolant (not line), physical cooling process, often a bug, heating[0][12] */
			strcpy( chLabel[k], "Hcol" );
		}
		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 == 18 )
		{
			/* H2 photoionization heating, heating[0][18] */
			strcpy( chLabel[k], "H2ph" );
		}
		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 */
		  SaveVal, 
		  /* 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 */
		  &nFail);

	/*>>chng 06 jun 06, change start of save to give same info as cooling 
	 * as per comment by Yumihiko Tsuzuki */
	/* begin the print out with zone number, total heating and cooling */
	fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e", 
		radius.depth_mid_zone, 
		phycon.te, 
		heat_total, 
		cool_total );

	/* 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 );
		if(k > 4 && thermal.heating[i][j]/SDIV(heat_total) < save.WeakHeatCool )
			break;
		fprintf( io, "\t%s\t%.7f ", 
			chLabel[ip], SaveVal[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 )
	{
		long level = -1;
		/* a line was a major heat source - identify it */
		TransitionProxy t = FndLineHt(&level);
		if( t.Coll().heat()/SDIV(heat_total) > 0.005 )
		{
			ASSERT( t.associated() );
			strcpy( chLbl, chLineLbl(t) );
			TauIn = t.Emis().TauIn();
			Pump = t.Emis().pump();
			EscP = t.Emis().Pesc();
			CS = t.Coll().col_str();
			/* ratio of line to total heating */
			ColHeat = t.Coll().heat()/SDIV(heat_total);
			
			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(SaveVal);
	return;
}