ion_collis.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ion_collis fill in collisional ionization rates, and resulting cooling */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "rfield.h"
#include "heavy.h"
#include "dense.h"
#include "atmdat.h"
#include "ionbal.h"

void ion_collis(
        /* element number on c scale, H is 0 */
        long int nelem )
{
        long int ion,
                limit;
        double DimaRate, crate;

        DEBUG_ENTRY( "ion_collis()" );

        /* compute collisional ionization rate */

        /* CollidRate[nelem][ion][0] is collisional ionization rate, s-1
         * CollidRate[nelem][ion][1] is collisional ionization cooling, erg/s
         */

        /* zero out rates below lowest ionization stage we will consider */
        for( ion=0; ion < (dense.IonLow[nelem] - 1); ion++ )
        {
                ionbal.CollIonRate_Ground[nelem][ion][0] = 0.;
                ionbal.CollIonRate_Ground[nelem][ion][1] = 0.;
        }

        /* chng logic changed to be precisely same as ion_solver */
        /* >>chng 02 nov 08, change 2 to NISO */
        /*limit = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/
        limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
        ASSERT( limit < LIMELM );

        for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
        {
                //Get the rate coefficients using either Dima or Hybrid
                DimaRate = t_ADfA::Inst().coll_ion_wrapper( nelem, ion , phycon.te );

                crate = DimaRate*dense.EdenHCorr;

                /* total collisional ionization rate 
                 * with only thermal suprathermal electrons */
                ionbal.CollIonRate_Ground[nelem][ion][0] = crate;

                /* cooling due to collisional ionization, which only includes thermal */
                ionbal.CollIonRate_Ground[nelem][ion][1] = (crate*
                        rfield.anu[Heavy.ipHeavy[nelem][ion]-1]* EN1RYD);
        }

        for( ion=dense.IonHigh[nelem]; ion <= nelem; ion++ )
        {
                ionbal.CollIonRate_Ground[nelem][ion][0] = 0.;
                ionbal.CollIonRate_Ground[nelem][ion][1] = 0.;
        }

        /* check not rates are negative - in release mode this loop will optimize out */
        for( ion=0; ion <= nelem; ion++ )
        {
                /* there can be no negative rates */
                ASSERT( ionbal.CollIonRate_Ground[nelem][ion][0] >= 0. );
        }
        return;
}