/*cdIonFrac get ionization fractions for a constituent */
/* tested 99 mar 12 */
#include <string.h>
#include "cddefines.h"
#include "elementnames.h"
#include "mean.h"
#include "caps.h"
#include "cddrive.h"

int cdIonFrac(
	/* 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 fractional ionization */
	double *fracin, 
	/* how to weight the average, must be "VOLUME" or "RADIUS" */
	char *chWeight ) 
	/* return value is logical flag saying whether 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 chCARD[81];

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

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

	/*caps(chWeight);*/

	if( strcmp(chCARD,"RADIUS") == 0 )
	{
		lgVol = FALSE;
	}

	else if( strcmp(chCARD,"VOLUME") == 0 )
	{
		lgVol = TRUE;
	}

	else
	{
		fprintf( ioQQQ, " CDIONF: chWeight=%6.6s makes no sense to me\n", 
		  chWeight );
		*fracin = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(FALSE);
	}

	/* first ensure that chLabel is all caps */
	strcpy( chCARD, chLabel );
	/* make sure chWeight is all caps */
	caps(chCARD);/* convert to caps */

	/* find which element this is */
	nelem = 1;
	while( nelem <= LIMELM && 
		strcmp(chCARD,elementnames.chElementNameShort[nelem-1]) != 0 )
	{
		nelem += 1;
	}
	/* 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, " cdIonFrac called with unknown element chLabel, =%4.4s\n", 
		  chLabel );
		
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(FALSE);
	}

	/* sanity check - does this ionization stage exist? 
	 * both IonStage and nelem are on physical scales, H atoms are 1 1 */

	/* 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 || 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 );
		*fracin = -1.;
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\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('i', nelem,&ip,aaa);
		*fracin = pow(10.,aaa[IonStage-1]);
	}
	else
	{
		RadMean('i', nelem,&ip,aaa);
		*fracin = pow(10.,aaa[IonStage-1]);
	}

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