/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*atom_pop5 do five level atom population and cooling */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "lapack.h"
#include "atoms.h"

void atom_pop5(double g[], 
  double ex[], 
  double cs12, 
  double cs13, 
  double cs14, 
  double cs15, 
  double cs23, 
  double cs24, 
  double cs25, 
  double cs34, 
  double cs35, 
  double cs45, 
  double a21, 
  double a31, 
  double a41, 
  double a51, 
  double a32, 
  double a42, 
  double a52, 
  double a43, 
  double a53, 
  double a54, 
  double p[], 
  float abund)
{
	int32 ipiv[5], ner;

	long int i, j;

	double c12, 
	  c13, 
	  c14, 
	  c15, 
	  c21, 
	  c23, 
	  c24, 
	  c25, 
	  c31, 
	  c32, 
	  c34, 
	  c35, 
	  c41, 
	  c42, 
	  c43, 
	  c45, 
	  c51, 
	  c52, 
	  c53, 
	  c54, 
	  e12, 
	  e13, 
	  e14, 
	  e15, 
	  e23, 
	  e24, 
	  e25, 
	  e34, 
	  e35, 
	  e45, 
	  tf;

	double amat[5][5], 
	  bvec[5], 
	  dmax, 
	  zz[6][6];

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

	/* following used for linpac matrix inversion */

	/* tf = 1.438 / te */
	tf = T1CM/phycon.te;

	/* quit if no species present */
	if( abund <= 0. )
	{
		p[0] = 0.;
		p[1] = 0.;
		p[2] = 0.;
		p[3] = 0.;
		p[4] = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->atom_pop5()\n", debug_fp );
#		endif
		return;
	}

	/* get some Boltzmann factors */
	e12 = sexp(ex[0]*tf);
	e23 = sexp(ex[1]*tf);
	e34 = sexp(ex[2]*tf);
	e45 = sexp(ex[3]*tf);
	e13 = e12*e23;
	e14 = e13*e34;
	e15 = e14*e45;
	e24 = e23*e34;
	e25 = e24*e45;
	e35 = e34*e45;

	/* quit it highest level Boltzmann factor too large */
	if( e15 == 0. )
	{
		p[0] = 0.;
		p[1] = 0.;
		p[2] = 0.;
		p[3] = 0.;
		p[4] = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->atom_pop5()\n", debug_fp );
#		endif
		return;
	}

	/* get collision rates,
	 * gcdsqte IS 8.629e-6 / sqrte * eden */
	c21 = phycon.cdsqte*cs12/g[1];
	c12 = c21*g[1]/g[0]*e12;

	c31 = phycon.cdsqte*cs13/g[2];
	c13 = c31*g[2]/g[0]*e13;

	c41 = phycon.cdsqte*cs14/g[3];
	c14 = c41*g[3]/g[0]*e14;

	c51 = phycon.cdsqte*cs15/g[4];
	c15 = c51*g[4]/g[0]*e15;

	c32 = phycon.cdsqte*cs23/g[2];
	c23 = c32*g[2]/g[1]*e23;

	c42 = phycon.cdsqte*cs24/g[3];
	c24 = c42*g[3]/g[1]*e24;

	c52 = phycon.cdsqte*cs25/g[4];
	c25 = c52*g[4]/g[1]*e25;

	c43 = phycon.cdsqte*cs34/g[3];
	c34 = c43*g[3]/g[2]*e34;

	c53 = phycon.cdsqte*cs35/g[4];
	c35 = c53*g[4]/g[2]*e35;

	c54 = phycon.cdsqte*cs45/g[4];
	c45 = c54*g[4]/g[3]*e45;

	/* rate equations equal zero */
	for( i=0; i < 5; i++ )
	{
		zz[i][4] = 1.0;
		zz[5][i] = 0.;
	}
	zz[5][4] = abund;

	/* level one balance equation */
	zz[0][0] = c12 + c13 + c14 + c15;
	zz[1][0] = -a21 - c21;
	zz[2][0] = -a31 - c31;
	zz[3][0] = -a41 - c41;
	zz[4][0] = -a51 - c51;

	/* level two balance equation */
	zz[0][1] = -c12;
	zz[1][1] = c21 + a21 + c23 + c24 + c25;
	zz[2][1] = -c32 - a32;
	zz[3][1] = -c42 - a42;
	zz[4][1] = -c52 - a52;

	/* level three balance equation */
	zz[0][2] = -c13;
	zz[1][2] = -c23;
	zz[2][2] = a31 + a32 + c31 + c32 + c34 + c35;
	zz[3][2] = -c43 - a43;
	zz[4][2] = -c53 - a53;

	/* level four balance equation */
	zz[0][3] = -c14;
	zz[1][3] = -c24;
	zz[2][3] = -c34;
	zz[3][3] = a41 + c41 + a42 + c42 + a43 + c43 + c45;
	zz[4][3] = -c54 - a54;

	/* divide both sides of equation by largest number to stop overflow */
	dmax = -1e0;

	for( i=0; i < 6; i++ )
	{
		for( j=0; j < 5; j++ )
		{
			dmax = MAX2(zz[i][j],dmax);
		}
	}

	for( i=0; i < 6; i++ )
	{
		for( j=0; j < 5; j++ )
		{
			zz[i][j] /= dmax;
		}
	}

	/* this one may be more robust */
	for( j=0; j < 5; j++ )
	{
		for( i=0; i < 5; i++ )
		{
			amat[i][j] = zz[i][j];
		}
		bvec[j] = zz[5][j];
	}

	ner = 0;

	getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner);
	getrs_wrapper('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner);

	/*DGETRF(5,5,(double*)amat,5,ipiv,&ner);*/
	/*DGETRS('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner);*/

	if( ner != 0 )
	{
		fprintf( ioQQQ, " atom_pop5: dgetrs finds singular or ill-conditioned matrix\n" );
		puts( "[Stop in atom_pop5]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now put results back into z so rest of code treates only
		* one case - as if matin1 had been used */
	for( i=0; i < 5; i++ )
	{
		zz[5][i] = bvec[i];
	}

	/* TODO	2	p(5) was very slightly negative (1e-40) for SII in dqher.in */
	p[1] = MAX2(0.e0,zz[5][1]);
	p[2] = MAX2(0.e0,zz[5][2]);
	p[3] = MAX2(0.e0,zz[5][3]);
	p[4] = MAX2(0.e0,zz[5][4]);
	p[0] = abund - p[1] - p[2] - p[3] - p[4];

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

