/* 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 */
/*molcol generate and print molecular column densities */
#include "cddefines.h"
#include "radius.h"
#include "dense.h"
#include "doppvel.h"
#include "colden.h"
#include "h2.h"
#include "co.h"
#include "molcol.h"

void molcol(char *chLabel,
	/* file for printout */
	FILE *ioMEAN )
{
	long int i;

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

	if( strcmp(chLabel,"PRIN") == 0 )
	{
		/* total hydrogen column density, all forms */
		fprintf( ioMEAN, "\n                                                     Log10 Column density (cm^-2)\n");
		fprintf( ioMEAN, "   Htot :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_HTOT])));
		fprintf( ioMEAN, "   HII  :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCHII])));
		fprintf( ioMEAN, "   HI   :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCHI])));
		fprintf( ioMEAN, "   H-   :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCHMIN])));
		fprintf( ioMEAN, "   H2   :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOLH2])));
		fprintf( ioMEAN, "   H2+  :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCH2PLS])));
		fprintf( ioMEAN, "   HeH+ :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCHEHP] )));
		fprintf( ioMEAN, "   H3+  :");
		fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCH3PLS] )));
		fprintf( ioMEAN, "\n");
	}

	/* call large H2 and CO column density routines which will do their jobs */
	H2_Colden( chLabel);
	CO_colden( chLabel);

	if( strcmp(chLabel,"ZERO") == 0 )
	{
		/*  zero out the column densities */
		for( i=0; i < NUM_HEAVY_MOLEC; i++ )
		{
			co.hevcol[i] = 0.;
		}
	}

	else if( strcmp(chLabel,"ADD ") == 0 )
	{
		/*  add together column densities */
		for( i=0; i < NUM_COMOLE_CALC; i++ )
		{
			co.hevcol[i] += co.hevmol[i]*(float)radius.drad_x_fillfac;
		}
	}

	else if( strcmp(chLabel,"PRIN") == 0 )
	{
		/*  print the molecular column densities
		 * want to print all the molecules, not including the atoms/ions
		 * that are part of the co solver.  use first to print them all */
		/*for( i=0; i < NUM_COMOLE_CALC; i++ )*/
		for( i=0; i < NUM_HEAVY_MOLEC; i++ )
		{
			/* print 7 column densities per line */
			if( i!=0 && i%8==0 )
				fprintf( ioMEAN, "\n" );
			fprintf( ioMEAN, "   %s:", co.chLab[i] );
			fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,co.hevcol[i] )));
		}
		fprintf( ioMEAN, "\n" );
	}

	else
	{
		fprintf( ioMEAN, " molcol does not understand the label %4.4s\n", 
		  chLabel );
		puts( "[Stop in molcol]" );
		cdEXIT(EXIT_FAILURE);
	}

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

}

