/* 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 */
/*PrtMeanIon print mean ionization fractions or temperatures for all elements */
#include "cddefines.h"
#include "geometry.h"
#include "dense.h"
#include "mean.h"
#include "elementnames.h"
#include "prt.h"

void PrtMeanIon(
	/* this is either 't' or 'i', for mean ionization or temperature */
	char chType ,
	/* TRUE include density, FALSE do not */
	int lgDensity ,
	/* this is stream used for io, is stdout when called by final,
	 * is punch unit when punch output generated */
	 FILE *ioMEAN )
{
	long int i, 
	  limit, 
	  n, 
	  nelem;
	static int lgPrtLots=FALSE;

	float aa[64];

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

	/* print mean ionization or temperature?? */
	/* chType must be i for ionization, or t for temperature */
	ASSERT( chType == 'i' || chType == 't' );

	/*  print mean ionization or temperature for the computed model
	 * the means are relative to the total atoms and ions and does not include molecules or grains */
	
	/* only print sphereical if not plane parallel */
	/* >>>chng 99 may 01, test included one for sphere being set,
	 * so no printout when not set.  now print even when sphere not
	 * set but geo is not plane parallel */
	if( /*sphere.lgSphere &&*/ (!geometry.lgGeoPP) )
	{

		/* get means for hydrogen and helium, MeanIonVolume is in mean.c */
		MeanIonVolume(chType,ipHYDROGEN,&n,aa,lgDensity);
		/*MeanIonVolume(chType,2,&n,c ,lgDensity);*/

		/* now print hydrogen, special since part of title goes here */
		fprintf( ioMEAN, "\n Hydrogen  " );
		for( i=0; i < 3; i++ )
		{
			fprintf( ioMEAN, "%7.3f", aa[i] );
		}
		fprintf(ioMEAN," (H2)");
		if( chType=='i' && lgDensity )
		{
			fprintf( ioMEAN, 
				"         Log10 Mean Ionisation (over volume*electron density)\n" );
		}
		else if( chType=='i' )
		{
			fprintf( ioMEAN, 
				"                 Log10 Mean Ionisation (over volume)\n" );
		}
		else if( chType=='t' && lgDensity )
		{
			fprintf( ioMEAN, 
				"          Log10 Mean Temperature (over volume*electron density)\n" );
		}
		else if( chType=='t' )
		{
			fprintf( ioMEAN, 
				"                  Log10 Mean Temperature (over volume)\n" );
		}
		else
		{
			fprintf(ioQQQ,"PrtMeanIon called with insane job\n");
		}

		/* ionization fractions for remaining elements */
		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				MeanIonVolume(chType,nelem,&n,aa,lgDensity);
				limit = MIN2(17,n);
				fprintf( ioMEAN, " %10.10s", elementnames.chElementName[nelem]  );

				for( i=0; i < limit; i++ )
				{
					fprintf( ioMEAN, "%7.3f", aa[i] );
				}
				fprintf( ioMEAN, "\n" );

				if( n > 17 )
				{
					lgPrtLots = TRUE;
					fprintf( ioMEAN, "           " );
					for( i=17; i < n; i++ )
					{
						fprintf( ioMEAN, "%7.3f", aa[i] );
					}
					fprintf( ioMEAN, "\n" );
				}
			}
		}
	}

	fprintf( ioMEAN, "\n         " );
	for( i=1; i <= 17; i++ )
	{
		fprintf( ioMEAN, "%7ld", i );
	}
	fprintf( ioMEAN, "\n" );
	if( lgPrtLots )
	{
		fprintf( ioMEAN, "         " );
		for( i=18; i <= LIMELM; i++ )
		{
			fprintf( ioMEAN, "%7ld", i );
		}
		fprintf( ioMEAN, "\n" );
	}

	/* get hydrogen ionization fractions */
	MeanIonRadius(chType,ipHYDROGEN,&n,aa,lgDensity);

	/* get helium ionization fractions 
	MeanIonRadius(chType,ipHELIUM,&n,c,lgDensity);*/

	fprintf( ioMEAN, "\n Hydrogen  " );

	for( i=0; i < 3; i++ )
	{
		/* print hydrogen fractions */
		fprintf( ioMEAN, "%7.3f", aa[i] );
	}
	fprintf(ioMEAN," (H2)");
	if( chType=='i' && lgDensity )
	{
		fprintf( ioMEAN, 
			"         Log10 Mean Ionisation (over radius*electron density)\n" );
	}
	else if( chType=='i' )
	{
		fprintf( ioMEAN, 
			"                 Log10 Mean Ionisation (over radius)\n" );
	}
	else if( chType=='t' && lgDensity )
	{
		fprintf( ioMEAN, 
			"         Log10 Mean Temperature (over radius*electron density)\n" );
	}
	else if( chType=='t' )
	{
		fprintf( ioMEAN, 
			"                 Log10 Mean Temperature (over radius)\n" );
	}
	else
	{
		fprintf(ioQQQ,"PrtMeanIon called with insane job\n");
	}

	/* now the remaining elements, he on up */
	for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			MeanIonRadius(chType,nelem,&n,aa,lgDensity);
			limit = MIN2(17,n);
			fprintf( ioMEAN, " %10.10s", 
				elementnames.chElementName[nelem]  );
			for( i=0; i < limit; i++ )
			{
				fprintf( ioMEAN, "%7.3f", aa[i] );
			}
			fprintf( ioMEAN, "\n" );

			if( n > 17 )
			{
				lgPrtLots = TRUE;
				fprintf( ioMEAN, "           " );

				for( i=17; i < n; i++ )
				{
					fprintf( ioMEAN, "%7.3f", aa[i] );
				}
				fprintf( ioMEAN, "\n" );
			}
		}
	}

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