/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*IonBoron ionization balance for boron */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	6
#define	NELEM	5
#include "dense.h"
#include "trace.h"
#include "ionbal.h"

void IonBoron(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double ditcrt[NDIM - 1]={1.2e4,1.2e4,1.1e4,4.4e5,1e20};
	static double aa[NDIM - 1]={1.8267,2.3196,0.,0.,0.};
	static double bb[NDIM - 1]={4.1012,10.7328,0.,0.,0.};
	static double cc[NDIM - 1]={4.8443,6.8830,0.,0.,0.};
	static double dd[NDIM - 1]={.2261,-0.1824,0.,0.,0.};
	static double ff[NDIM - 1]={0.5960,0.4101,0.1,0.1,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {2.54e-3,6.15e-3,1.62e-3,4.78e-2,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {4.42e-2,5.88e-2,0.343,0.362,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {1.57e5,1.41e5,8.19e4,3.44e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {3.74e5,1.41e5,1.59e5,5.87e5,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* boron nelem=5
	 * data are for carbon
	 *
	 * real CollidRate(nelem,2)
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */
	/* DATA GRDEFF/0.10,0.10,0.10,0.053,0.10/
	 * GRDEFF is fraction of recombinations to ground state, used for
	 * outward diffuse fields
	 *
	 * rec from +3, +4 from Arnaud et al Ast Ap Sup 60 425. (1985)
	 * rec from fully ionized uses Seaton '79 in ionrat */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	if( !dense.lgElmtOn[NELEM-1] )
	{
		
#		ifdef DEBUG_FUN
		fputs( " <->IonBoron()\n", debug_fp );
#		endif
		return;
	}

	/* zero out ionization balance arrays */
	ion_zero(NELEM-1);

	ion_photo(NELEM-1,FALSE);

	/* find collisional ionization rates */
	ion_collis(NELEM-1);

	/* get recombination coefficients */
	ion_recomb(FALSE,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
	  dd,ff,NELEM-1);

	/* solve for ionization balance */
	ion_solver(NELEM-1,FALSE);

	if( trace.lgTrace && trace.lgHeavyBug )
	{
		fprintf( ioQQQ, "     Boroni returns; frac=" );
		for( i=1; i <= (NELEM + 1); i++ )
		{
			fprintf( ioQQQ, "%10.3e", dense.xIonDense[NELEM-1][i-1]/
			  dense.gas_phase[NELEM-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

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

