/*IonCoppe do ionization balance for copper */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	30
#define	NELEM	29
#include "elmton.h"
#include "timed.h"
#include "collidionize.h"
#include "ionzer.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonCoppe(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={2.27e-10,6.31e-10,8.83e-10,1.01e-09,
	  3.71e-09,4.84e-09,1.09e-08,2.48e-08,9.68e-08,1.72e-07,2.61e-07,
	  3.95e-07,4.70e-07,5.70e-07,6.98e-07,8.32e-07,9.76e-07,1.19e-06,
	  1.16e-06,1.14e-06,1.17e-06,1.24e-06,1.35e-06,1.38e-06,6.64e-07,
	  4.67e-07,8.98e-07,8.98e-07};
	static double pl[NDIM - 2]={-0.7,-0.7,-0.7,-0.7,-0.746,-0.682,-0.699,
	  -0.728,-0.759,-0.79,-0.81,-0.829,-0.828,-0.834,-0.836,-0.84,
	  -0.846,-0.85,-0.836,-0.824,-0.816,-0.811,-0.808,-0.8,-0.718,
	  -0.677,-0.732,-0.732};
	static double tlow[2]={1.81e-10,-0.681};
	static double ditcrt[NDIM - 1]={6e3,2e4,4e4,5e4,7e4,8e4,8e4,3e4,
	  3e4,3e4,3e4,9e4,4e4,5e4,3e4,9e5,2e5,2e5,2e5,2e5,1e5,7e4,4e4,
	  6e6,6e6,1e20,1e20,1e20,1e20};
	static double aa[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	static double bb[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	static double cc[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	static double dd[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	static double ff[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {1.41e-3,5.2e-3,.0138,.023,.0419,
		  .0683,.122,.300,.15,.697,.709,.644,.525,.446,.363,.302,.102,
		  .27,.0467,.0835,.0996,.199,.24,.115,.0316,.803,.575,.575,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {.469,.357,.281,.128,.0417,.0558,
		  .0346,0.,1.9,.277,.135,.134,.192,.332,.337,.121,.0514,.183,
		  7.56,4.55,4.87,2.19,1.15,1.23,.132,.289,.286,.286,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {9.82e4,2.01e5,3.05e5,4.20e5,5.56e5,
		  6.72e5,7.93e5,9.00e5,1.00e6,7.81e5,7.64e5,7.44e5,6.65e5,
		  5.97e5,5.24e5,4.96e5,4.46e5,8.46e6,1.36e6,1.23e6,1.06e6,
		  1.25e6,1.23e6,3.32e5,6.45e5,6.65e7,6.81e7,6.81e7,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {1.01e5,1.91e5,2.32e5,3.18e5,4.55e5,
		  5.51e5,5.28e5,0.,5.50e5,8.87e5,1.80e6,1.25e6,1.89e6,8.84e5,
		  1.29e6,6.24e5,1.59e6,8.01e6,9.32e6,9.45e6,9.45e6,8.01e6,
		  7.57e6,2.64e6,1.93e6,1.19e7,9.08e6,9.08e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* copper nelem=29
	 * rec data not for real
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */

	/* GRDEFF is fraction of recombinations to ground state, used for
	 * DATA GRDEFF/28*0.2/
	 *
	 * all rates from Shull and Van Steenberg apj sup 48, 95. */

	/* Pequignot and Aldrovandi Ast Ap 161, 169. */


	/* copper */
	if( !elmton.lgElmtOn[NELEM-1] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonCoppe()\n", debug_fp );
#	endif
		return;
	}

	ionzer(NELEM);

	PhotoIonize(NELEM-1,FALSE);

	/* 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);

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

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

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