/* 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 */
/*CoolPunch punch coolants */
#include "cddefines.h"
#include "thermal.h"
#include "dynamics.h"
#include "converge.h"
#include "phycon.h"
#include "punch.h"

/* this is limit to number of coolants to print out */
#define  IPRINT 100

/*CoolPunch punch coolants */
void CoolPunch(FILE * io)
{
	long int i, 
	  ip, 
	  is, 
		*index;

	int
		lgFail;

	double 	  cset, 
		coolok,
		heatok;
	float
		*csav,
		*sgnsav;

#	ifdef DEBUG_FUN
	fputs( "<+>CoolPunch()\n", debug_fp );
#	endif
	index = (long int *) CALLOC((size_t)thermal.ncltot,sizeof(long int));
	csav = (float *) CALLOC((size_t)thermal.ncltot,sizeof(float));
	sgnsav = (float *) CALLOC((size_t)thermal.ncltot,sizeof(float));

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

	/* cset will be weakest cooling to consider
	 * WeakHeatCool set with 'set weakheatcool' command
	 * default is 0.05 */
	cset = coolok*punch.WeakHeatCool;

	/* first find all strong lines, both + and - sign */
	ip = thermal.ncltot;

	for( i=0; i < ip; i++ )
	{
		csav[i] = (float)(MAX2(thermal.cooling[i],thermal.heatnt[i])/
			coolok);
		
		/* save sign to remember if heating or cooling line */
		if( thermal.heatnt[i] == 0. )
		{
			sgnsav[i] = 1.;
		}
		else
		{
			sgnsav[i] = -1.;
		}
	}

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

	/* warn if tcovergence failure occurred */
	if( !conv.lgConvTemp )
	{
		fprintf( io, " >>>>  Temperature not converged.\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" );
	}

	/* start the printout with total heating and total cooling */
	fprintf( io, "%.4e\t%.4e\t%.4e", 
		phycon.te, 
		heatok,
		coolok );

	/* print only up to IPRINT, which is defined above */
	ip = MIN2( ip , IPRINT );

	/* now print the coolants 
	 * keep sign of coolant, for strong negative cooling 
	 * order is ion, wavelength, fraction of total */
	for( is=0; is < ip; is++ )
	{
		if (is > 4 && (thermal.cooling[index[is]] < cset && thermal.heatnt[index[is]] < cset))
			break;
		fprintf( io, "\t%s %.1f\t%.7f", 
			thermal.chClntLab[index[is]], 
			thermal.collam[index[is]], 
			sign(csav[index[is]],sgnsav[index[is]]) );
	}
	fprintf( io, " \n" );

	/* finished, now free space */
	free(sgnsav);
	free(csav);
	free(index);

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

