/*IonOxyge derive ionization balance for oxygen */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	9
#define	NELEM	8
#include "physconst.h"
#include "opacpoint.h"
#include "o3exc.h"
#include "heating.h"
#include "elmton.h"
#include "hevmolec.h"
#include "nh.h"
#include "trace.h"
#include "he.h"
#include "augro3.h"
#include "pmp2s.h"
#include "s3727.h"
#include "i3727.h"
#include "timed.h"
#include "rfield.h"
#include "h.h"
#include "photrate.h"
#include "coatom.h"
#include "ionfracs.h"
#include "zonecnt.h"
#include "charexc.h"
#include "ionrange.h"
#include "collidionize.h"
#include "ionzer.h"
#include "popexc.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "gammas.h"
#include "ionheavy.h"

void IonOxyge(void)
{
	long int i, 
	  iup, 
	  _r;
	double aeff, 
	  save;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={1.60e-10,7.67e-10,2.35e-09,4.59e-09,
	  1.73e-08,3.04e-8,1.56e-7};
	static double pl[NDIM - 2]={-0.678,-0.646,-0.666,-0.670,-0.759,
	  -0.774,-0.834};
	static double tlow[2]={0.,0.};
	static double ditcrt[NDIM - 1]={2.7e4,2.2e4,2.4e4,2.5e4,1.6e4,1.0e6,
	  1.5e6,1e20};
	static double aa[NDIM - 1]={0.,-0.0036,0.,0.0061,-2.8425,0.,0.,
	  0.};
	static double bb[NDIM - 1]={0.0238,0.7519,21.8790,0.2269,0.2283,
	  0.,0.,0.};
	static double cc[NDIM - 1]={0.0659,1.5252,16.2730,32.1419,40.4072,
	  0.,0.,0.};
	static double dd[NDIM - 1]={0.0349,-0.0838,-0.7020,1.9939,-3.4956,
	  0.,0.,0.};
	static double ff[NDIM - 1]={0.5334,0.2769,1.1899,-0.0646,1.7558,
	  0.,0.,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {1.11e-3,5.07e-3,1.48e-2,1.84e-2,
		  4.13e-3,1.06e-1,6.23e-2,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {.0925,.181,.305,.1,.162,.34,.304,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {1.75e5,1.98e5,2.41e5,2.12e5,1.25e5,
		  6.25e6,7.01e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {1.45e5,3.35e5,2.83e5,2.83e5,2.27e5,
		  1.12e6,1.47e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* oxygen nelem=8
	 *
	 *
	 * highest three from Arnaud et al '85 */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* oxygen, atomic number 8 */
	if( !elmton.lgElmtOn[7] )
	{
		o3exc.poiii2Max = 0.;
		o3exc.poiii3Max = 0.;
		o3exc.r4363Max = 0.;
		o3exc.r5007Max = 0.;
		o3exc.poiii2 = 0.;
		pmp2s.p1666 = 0.;
		coatom.oatmic = 0.;
		augro3.AugerO3 = 0.;
		pmp2s.p1401 = 0.;
		s3727Com.s3727 = 0.;
		s3727Com.s7325 = 0.;
		HeatingCom.heating[9][7] = 0.;
		o3exc.poimax = 0.;
		
#	ifdef DEBUG_FUN
		fputs( " <->IonOxyge()\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);

	/* O-H charge exchange
	 *
	 * O+2 => O+ Butler Heil and Dalgarno 1980
	 * CTHeavy(2,2) = CTHeavy(2,2) + 2.8E-11*TE10*TE10*TE10*TE03*TE03*HI
	 * T = MIN( SQRTE,200. ) * HI
	 * O+3 => O+2 Bulter Heil and Dalgarno 1980
	 * CTHeavy(3,2) = CTHeavy(3,2) + 8.6E-11 * T + HEI*1E-9 */
	CharExc.CTHeavy[1][2] += (float)(he.hei*1e-9);
	/* write(6,*) 'o3', 8.6e-11*t ,HCharExcRec(3,8)*hi
	 * O+4 => O+3 Butler and Dalgarno 1980B
	 * CTHeavy(4,2) = CTHeavy(4,2) + 2.6E-12 * t + HEI*6E-10 */
	CharExc.CTHeavy[1][3] += (float)(he.hei*6e-10);

	/* photoexcitation of O III 1666 and O IV 1401
	 * 666 error! this will be zero in current form of phfit
	 * set 2s**2 rate to rate for O V */
	pmp2s.p1666 = PhotRate.PhotoRate[0][1][3][7];

	pmp2s.p1401 = PhotRate.PhotoRate[0][1][2][7];

	/* photoionization from O++ 1D
	 *
	 * estimate gamma function by assuming no frequency dependence
	 * betwen 1D and O++3P edge */
	o3exc.d5007r = (float)(GammaK(OpacPoint.ipo3exc[0],OpacPoint.ipo3exc[1],
	  OpacPoint.ipo3exc[2] , 1. ));

	o3exc.d4363 = (float)(GammaK(OpacPoint.ipo3exc3[0],OpacPoint.ipo3exc3[1],
	  OpacPoint.ipo3exc3[2] , 1. ));

	o3exc.d6300 = (float)(GammaK(OpacPoint.ipo1exc[0],OpacPoint.ipo1exc[1],
	  OpacPoint.ipo1exc[2] , 1. ));

	/* A21 = 0.0263 */
	aeff = 0.0263 + o3exc.d5007r;

	/* 1. as last arg makes this the relative population */
	o3exc.poiii2 = (float)(popexc(2.5,9.,5.,aeff,2.88e4,1.)/aeff);

	/* photoionization from excited states */
	if( ZoneCnt.nzone > 0 )
	{
		PhotRate.PhotoRate[0][2][0][7] = PhotRate.PhotoRate[0][2][0][7]*
		  (1. - o3exc.poiexc) + o3exc.d6300*o3exc.poiexc;

		PhotRate.PhotoRate[0][2][2][7] = PhotRate.PhotoRate[0][2][2][7]*
		  (1. - o3exc.poiii2 - o3exc.poiii3) + o3exc.d5007r*o3exc.poiii2 + 
		  o3exc.d4363*o3exc.poiii3;

		if( PhotRate.PhotoRate[0][2][2][7] > 1e-30 && IonRange.IonLow[NELEM-1] <= 3 )
		{
			if( (o3exc.d5007r*o3exc.poiii2 + o3exc.d4363*o3exc.poiii3)/
			  PhotRate.PhotoRate[0][2][2][7] > (o3exc.r4363Max + 
			  o3exc.r5007Max) )
			{
				o3exc.poiii2Max = (float)(o3exc.d5007r*o3exc.poiii2/PhotRate.PhotoRate[0][2][2][7]);
				o3exc.poiii3Max = (float)(o3exc.d4363*o3exc.poiii3/PhotRate.PhotoRate[0][2][2][7]);
			}
			o3exc.r4363Max = (float)(MAX2(o3exc.r4363Max,o3exc.d4363));
			o3exc.r5007Max = (float)(MAX2(o3exc.r5007Max,o3exc.d5007r));
		}
		if( IonRange.IonLow[NELEM-1] <= 1 && (PhotRate.PhotoRate[0][2][0][7] + 
		  CharExc.HCharExcIon[7][0]*h.hii) > 1e-30 )
		{
			o3exc.poimax = (float)(MAX2(o3exc.poimax,o3exc.d6300*o3exc.poiexc/
			  (PhotRate.PhotoRate[0][2][0][7]+CharExc.HCharExcIon[7][0]* h.hii)));
		}
	}
	else
	{
		o3exc.poiii2Max = 0.;
		o3exc.poiii3Max = 0.;
		o3exc.r4363Max = 0.;
		o3exc.r5007Max = 0.;
		o3exc.poimax = 0.;
	}

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

	/* save photodistruction rate for 3727 creation */
	if( IonRange.IonLow[NELEM-1] == 1 && i3727.i2d < rfield.nflux )
	{
		s3727Com.s3727 = (float)(GammaK(i3727.i2d,i3727.i2p,OpacPoint.iopo2d , 1. ));

		iup = MIN2(nh.ipHn[1][0],rfield.nflux);
		s3727Com.s7325 = (float)(GammaK(i3727.i2d,iup,OpacPoint.iopo2d , 1. ));

		s3727Com.s7325 -= s3727Com.s3727;
		s3727Com.s3727 = s3727Com.s3727 + s3727Com.s7325;

		/* ratio of cross sections */
		s3727Com.s7325 *= 0.66f;
	}
	else
	{
		s3727Com.s3727 = 0.;
		s3727Com.s7325 = 0.;
	}

	augro3.AugerO3 = (float)PhotRate.PhotoRate[0][0][0][7];

	/* following accounts for part of oxygen that is in CO, since BiDiag ony
	 * wants atoms/ions, but [0][7] includes molecules too.  Must temperariyly
	 * convert [0][7] to atoms/molecules only */
	coatom.oatmic = (float)(xIonFracs[0][7]*
		MAX2(1e-7 , 1.-hevmolec.hevmol[IPCO-1]/xIonFracs[0][7] ) );

	save = xIonFracs[0][7];

	xIonFracs[0][7] = coatom.oatmic;

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

	xIonFracs[0][7] = (float)save;

	/* charge exchange onto hydrogen */
	CharExc.CTHionAgnt[2] = xIonFracs[2][7]*CharExc.HCharExcRec[7][0];

	/* write(72,'(i4,1p,10e10.2)')nzone,xIonFracs(8,2),HCharExcRec(1,8),
	 *  1 CTHionAgnt(3), xIonFracs(8,1) , HCharExcIon(1,8) ,
	 *  1  xIonFracs(8,1) * HCharExcIon(1,8) */
	CharExc.CTHion[0][0] += CharExc.CTHionAgnt[2];
	CharExc.CTHrec[0][0] += xIonFracs[1][7]*CharExc.HCharExcIon[7][0];

	CharExc.chrener += 
		(float)(-(h.hi*xIonFracs[2][7]*CharExc.HCharExcRec[7][0] - 
	  h.hii*xIonFracs[1][7]*CharExc.HCharExcIon[7][0])*230./TE1RYD);

	/* 1666 ratio corrected for phot crs at 50ev */
	pmp2s.p1666 *= xIonFracs[2][7]*0.3;
	pmp2s.p1401 *= xIonFracs[3][7]*0.43;
	s3727Com.s3727 *= xIonFracs[1][7];
	s3727Com.s7325 *= xIonFracs[1][7];
	augro3.AugerO3 *= xIonFracs[1][7];

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     IonOxyge returns; frac=" );
		for( i=1; i <= 9; i++ )
		{
			fprintf( ioQQQ, " %10.3e", xIonFracs[i][7]/
			  xIonFracs[0][7] );
		}
		fprintf( ioQQQ, "\n" );
	}

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