/*ofit compute cross sections for all shells of atomic oxygen */
#include "cddefines.h"
#include "ofit.h"

void ofit(double e, 
	  float *otot, 
	  float opart[])
{
	long int i, 
	  _r;

	double q, 
	  x;

	static double y[7][5];
	static double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.};
	static long l[7]={1,1,1,0,1,1,0};
	static int _aini = 1;

	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {8.915,3995.,3.242,10.44,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[0][i-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {11.31,1498.,5.27,7.319,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[1][i-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {10.5,1.059e05,1.263,13.04,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[2][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {19.49,48.47,8.806,5.983,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[3][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {50.,4.244e04,0.1913,7.012,4.454e-02};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[4][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {110.5,0.1588,148.3,-3.38,3.589e-02};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[5][i-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {177.4,32.37,381.2,1.083,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[6][i-1] = _itmp6[_r++];
			}
		}
		_aini = 0;
	}

#	ifdef DEBUG_FUN
	fputs( "<+>ofit()\n", debug_fp );
#	endif
	/*compute cross sections for all shells of atomic oxygen
	 * Photoionization of OI
	 * Input parameter:   e - photon energy, eV
	 * Output parameters: otot - total photoionization cross section, Mb
	 *  opart(1) - 2p-shell photoionization, goes to 4So
	 *  opart(2) - 2p-shell photoionization, goes to 2Do
	 *  opart(3) - 2p-shell photoionization, goes to 2Po
	 *  opart(4) - 2s-shell photoionization
	 *  opart(5) - double photoionization, goes to O++
	 *  opart(6) - triple photoionization, goes to O+++
	 *  opart(7) - 1s-shell photoionization */

	*otot = 0.0;

	for( i=0; i < 7; i++ )
	{
		opart[i] = 0.0;
	}

	for( i=0; i < 7; i++ )
	{
		if( e >= eth[i] )
		{
			q = 5.5 - 0.5*y[i][3] + l[i];

			x = e/y[i][0];

			opart[i] = (float)(y[i][1]*(POW2(x - 1.0) + POW2(y[i][4]))/
			  pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3]));

			*otot += opart[i];
		}
	}


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