/*cdTemp get mean electron temperature for any element */
#include <string.h>
#include "cddefines.h"
#include "elementnames.h"
#include "mean.h"
#include "caps.h"
#include "cddrive.h"

int cdTemp(
	/* four char string, null terminzed, giving the element name */
	char *chLabel, 
	/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
	long int IonStage, 
	/* will be temperature */
	double *TeMean, 
	/* how to weight the average, must be "VOLUME" or "RADIUS" */
	char *chWeight ) 
	/* return value is logical flag saying whether element was found, 
	 * FALSE if element not found and there are problems; 
	 * TRUE if things are ok and element was found*/
{

	int lgVol;
	long int ip, 
		ipaaa, /* used as pointer within aaa vector*/
		nelem;
	float aaa[LIMELM + 1];
	/* use following to store local image of character strings */
	char chWGHT[81] , chELEM[81];

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

	/* make sure chWeight is all caps */
	strcpy( chWGHT, chWeight );
	caps(chWGHT);/* convert to caps */

	/* ensure that chLabel is all caps */
	strcpy( chELEM, chLabel );
	caps(chELEM);/* convert to caps */

	/* special case, element is 21cm mean */
	if( strcmp(chELEM,"21CM") == 0 )
	{
		
		/* this is mean of inverse temperature */
		*TeMean = IonMeans.HarMeanTemp[0]/IonMeans.HarMeanTemp[1];
		if( *TeMean > SMALLFLOAT )
		{
			*TeMean = 1./ *TeMean;
		}
		else
		{
			*TeMean = 0.;
		}
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		/* say things are ok */
		return(TRUE);
	}

	/* now see if volume or radius weighting */
	if( strcmp(chWGHT,"RADIUS") == 0 )
	{
		lgVol = FALSE;
	}
	else if( strcmp(chWGHT,"VOLUME") == 0 )
	{
		lgVol = TRUE;
	}
	else
	{
		fprintf( ioQQQ, " cdTemp: chWeight=%6.6s makes no sense to me\n", 
		  chWeight );
		*TeMean = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(FALSE);
	}

	/* find which element this is */
	nelem = 0;
	while( nelem < LIMELM && 
		strcmp(chELEM,elementnames.chElementNameShort[nelem]) != 0 )
	{
		++nelem;
	}
	/* after this loop nelem is atomic number of element, H is 1 */

	/* if element not recognized and above loop falls through, nelem is LIMELM+1 
	 * nelem counter is on physical scale, H = 1 Zn = 30 */
	if( nelem >= LIMELM )
	{
		fprintf( ioQQQ, " cdTemp called with unknown element chLabel, =%4.4s\n", 
		  chLabel );
		
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(FALSE);
	}

	/* sanity check - does this ionization stage exist? 
	 * both IonStage on physical scale, nelem on c */

	/* ipaaa will be used as pointer within the aaa array that contains average values,
	 * done this way to prevent lint from falsing in acess to aaa array */
	ipaaa = IonStage - 1;

	/*if( IonStage > nelem + 1 || IonStage < 1 || IonStage > LIMELM+1 )*/
	if( ipaaa > nelem+1 || ipaaa < 0 || ipaaa > LIMELM )
	{
		fprintf( ioQQQ, " CDIONF asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 
		  IonStage, chLabel );
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(FALSE);
	}

	/* get either volume or radius average, aaa is filled in from 0 */
	if( lgVol )
	{
		/* aaa is dim'd LIMELM+1 so largest arg is LIMELM */
		VolMean('t', nelem+1,&ip,aaa);
		*TeMean = pow(10.,aaa[ipaaa]);
	}
	else
	{
		RadMean('t', nelem+1,&ip,aaa);
		*TeMean = pow(10.,aaa[ipaaa]);
	}

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