/*IonMagne ionization balance for magnesium */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	13
#define	NELEM	12
#include "opacpoint.h"
#include "mgexc.h"
#include "trace.h"
#include "elmton.h"
#include "he.h"
#include "timed.h"
#include "charexc.h"
#include "photrate.h"
#include "nh.h"
#include "zonecnt.h"
#include "ionrange.h"
#include "ionfracs.h"
#include "ionzer.h"
#include "collidionize.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "gammas.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonMagne(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={6.86e-10,3.04e-9,3.02e-9,5.73e-9,1.02e-8,
	  1.38e-8,1.86e-8,3.21e-8,7.84e-8,2.34e-7,4.15e-7};
	static double pl[NDIM - 2]={-0.855,-0.838,-0.734,-0.718,-0.716,
	  -0.695,-0.691,-0.711,-0.765,-0.829,-0.821};
	static double tlow[2]={2.01e-10,-0.681};
	static double ditcrt[NDIM - 1]={4.0e3,7.4e4,6.6e4,5.5e4,4.4e4,4.5e4,
	  4.5e4,5.0e5,3.4e4,2.4e6,4.0e6,1e20};
	static double aa[NDIM - 1]={1.2044,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.};
	static double bb[NDIM - 1]={-4.6836,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.};
	static double cc[NDIM - 1]={7.6620,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.};
	static double dd[NDIM - 1]={-0.5930,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.};
	static double ff[NDIM - 1]={1.6260,0.1,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {4.49e-4,1.95e-3,5.12e-3,7.74e-3,
		  1.17e-2,3.69e-2,3.63e-2,4.15e-2,8.86e-3,.252,.144,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {.021,.074,.323,.636,.807,.351,.548,
		  .233,.318,.315,.291,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {5.01e4,6.06e5,4.69e5,3.74e5,3.28e5,
		  4.80e5,3.88e5,3.39e5,2.11e5,1.40e7,1.50e7,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {2.81e4,1.44e6,7.55e5,7.88e5,1.02e6,
		  9.73e5,7.38e5,3.82e5,1.54e6,2.64e6,3.09e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* magnesium nelem=12
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */

	/* rec from +9, +10, +11 from Arnauld et al '85 */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* mg i from nuss+storey, who say that =>Mgii is very small */

	/* magnesium, atomic number 12 */
	if( !elmton.lgElmtOn[11] )
	{
		mgexc.xMg2Max = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->IonMagne()\n", debug_fp );
#		endif
		return;
	}

	ionzer(NELEM);

	PhotoIonize(NELEM-1,FALSE);
	/* debugging printout for shell photo rates - 0 for atom */
	/*GammaPrtRate( ioQQQ , 0 , NELEM-1 );*/

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

	/* get recombination coefficients */
	MakeRecomb(rec,pl,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
	  dd,ff,NELEM-1,tlow);

	/* charge exchange ionization, Allan et al. 1988, MNRAS 235, 1245.
	 * >>chng 96 Oct 1, moved all H ct to bidiag
	 * CTHeavy(1,1) = CTHeavy(1,1) + HII * 7.5E-19 * TE*TE*TE10
	 * charge exchange, rates from Butler and Dalgarno 1980B
	 * MG+2 => MG+
	 * CTHeavy(2,2) = CTHeavy(2,2) + 9E-14 * HI
	 * T = MIN( SQRTE , 200. )
	 * MG+3 => MG+2
	 * >>chng 96 Oct 1, moved all H ct to bidiag
	 * CTHeavy(3,2) = CTHeavy(3,2) + 6.5E-11 * T*HI + HEI*7E-10 */
	CharExc.CTHeavy[1][2] += (float)(he.hei*7e-10);
	/* MG+4 => MG+3
	 * >>chng 96 Oct 1, moved all H ct to bidiag
	 * CTHeavy(4,2) = CTHeavy(4,2) + 6.5E-11 * T * HI + HEI*3E-9 */
	CharExc.CTHeavy[1][3] += (float)(he.hei*3e-9);

	if( IonRange.IonLow[NELEM-1] <= 1 )
	{
		/* can only do this one time
		 * photoionization from excited upper state of 2798 */
		mgexc.rmg2l = (float)GammaK(OpacPoint.ipmgex,nh.ipHn[0][0],OpacPoint.ipOpMgEx,1.);
		PhotRate.PhotoRate[0][3][1][NELEM-1] += mgexc.rmg2l*mgexc.popmg2;
		if( ZoneCnt.nzone <= 1 )
		{
			mgexc.xMg2Max = 0.;
		}
		else if( PhotRate.PhotoRate[0][3][1][NELEM-1] > 1e-30 )
		{
			/* remember max rel photo rate */
			mgexc.xMg2Max = (float)(MAX2(mgexc.xMg2Max,mgexc.rmg2l*mgexc.popmg2/
			  PhotRate.PhotoRate[0][3][1][NELEM-1]));
		}
	}
	else
	{
		mgexc.xMg2Max = 0.;
	}

	if( timed.itime == -1 )
	{
		theavy(NELEM);
		
#		ifdef DEBUG_FUN
		fputs( " <->IonMagne()\n", debug_fp );
#		endif
		return;
	}

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

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

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