/* cdColm get the column density for a constituent 
 * can be used outside code to examine predictions after
 * calculation, and is used within code while optimizing 
 * column densities */
#include <string.h>
#include "cddefines.h"
#include "elementnames.h"
#include "colden.h"
#include "hevmolec.h"
#include "mean.h"
#include "caps.h"
#include "cddrive.h"

int cdColm(

	/* 4-char + eol string that is first 
	 * 4 char of element name as spelled by cloudy */
	char *chLabel,	

	/* integer stage of ionization, 1 for atom, 0 for CO or H2 */
	long int ion,

	/* the theoretical column density derived by the code */
	/* return value indicates whether code could find the column density,
	 * TRUE if OK, FALSE if not OK */
	double *theocl )	
{
	long int nelem;
	/* use following to store local image of character strings */
	char chCARD[81];

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

	strcpy( chCARD, chLabel );
	/* convert element label to all caps */
	caps(chCARD);

	/* there are two molecules, "H2  " and "CO  " that are recognized by this sub.
	 * these are indicated by negative ion */
	if( ion < 0 )
	{
		fprintf( ioQQQ, " cdColm called with insane ion, =%li\n", 
		  ion );
#		ifdef DEBUG_FUN
		fputs( " <->cdColm()\n", debug_fp );
#		endif
		return(FALSE);
	}
	else if( ion == 0 )
	{
		/* this case molecular column density */
		/* want the molecular hydrogen column density */
		if( strcmp( chCARD , "H2  " )==0 )
		{
			*theocl = coldenCom.colden[IPCOLH2-1];
		}

		/* carbon monoxide column density */
		else if( strcmp( chCARD , "CO  " )==0 )
		{
			*theocl = hevmolec.hevcol[IPCO-1];
		}

		/* OH column density */
		else if( strcmp( chCARD , "OH  " )==0 )
		{
			*theocl = hevmolec.hevcol[IPOH-1];
		}

		/* clueless as to what was meant - bail */
		else
		{
			fprintf( ioQQQ, " cdColm called with unknown element chLabel, =%4.4s\n", 
			  chLabel );
#			ifdef DEBUG_FUN
			fputs( " <->cdColm()\n", debug_fp );
#			endif
			return(FALSE);
		}
	}
	else
	{
		/* this case, ionization stage of some element */
		/* find which element this is */
		nelem = 0;
		while( nelem < LIMELM && 
			strncmp(chCARD,elementnames.chElementNameShort[nelem],4) != 0 )
		{
			++nelem;
		}

		/* this is true if we have one of the first 30 elements in the label */
		if( nelem < LIMELM )
		{

			/* sanity check - does this ionization stage exist? */
			if( ion > nelem + 2 )
			{
				fprintf( ioQQQ, 
				  " cdColm asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 
				  ion, chLabel );
				return(FALSE);
			}

			/* the column density */
			*theocl = IonMeans.xIonMeans[0][ion][nelem];
		}
	}

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