root/branches/newmole/source/atom_pop5.cpp

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  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 "dense.h"
#include "thirdparty.h"
#include "atoms.h"

/*atom_pop5 do five level atom population and cooling */
void atom_pop5(
        /* vector giving statistical weights on the five levels */
        double g[], 
        /* vector giving the excitation energy differences of the 5 levels.  The energies
         * are the energy in wavenumbers between adjacent levels.  So ex[0] is the energy
         * 1-2, ex[1] is the energy between 2-3, etc */
        double ex[], 
        /* the collision strengths for the levels */
  double cs12, 
  double cs13, 
  double cs14, 
  double cs15, 
  double cs23, 
  double cs24, 
  double cs25, 
  double cs34, 
  double cs35, 
  double cs45, 
  /* the transition probabilities between the various levels */
  double a21, 
  double a31, 
  double a41, 
  double a51, 
  double a32, 
  double a42, 
  double a52, 
  double a43, 
  double a53, 
  double a54, 
  /* the destroyed level populations (units cm-3) for the five levels */
  double p[], 
  /* the total density of this ion */
  realnum 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];

        DEBUG_ENTRY( "atom_pop5()" );

        /* 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.;
                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.;
                return;
        }

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

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

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

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

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

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

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

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

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

        c54 = dense.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;

        /* solve matrix */
        getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner);
        getrs_wrapper('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" );
                cdEXIT(EXIT_FAILURE);
        }

        /* now put results back into z so rest of code treats 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];
        return;
}