ion_recomb.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_recomb generate recombination coefficients for any species */
/*ion_recombAGN generate recombination coefficients for AGN table */
#include "cddefines.h"
#include "phycon.h"
#include "heavy.h"
#include "mole.h"
#include "hmi.h"
#include "grainvar.h"
#include "dense.h"
#include "conv.h"
#include "thermal.h"
#include "iso.h"
#include "abund.h"
#include "save.h"
#include "elementnames.h"
#include "atmdat.h"
#include "ionbal.h"
 
/*ion_recomb generate recombination coefficients for any species */
void ion_recomb(
  /* this is debug flag */
  bool lgPrintIt,
  /* nelem is the atomic number on the C scale, 0 for H */
  long int nelem/*, 
  double tlow[]*/)
{
        long int i, 
          ion, 
          limit;
 
        DEBUG_ENTRY( "ion_recomb(false,)" );
 
        ASSERT( nelem < LIMELM);
        ASSERT( nelem > 1 );
 
        /* check that range of ionization is correct */
        ASSERT( dense.IonLow[nelem] >= 0 );
        ASSERT( dense.IonLow[nelem] <= nelem+1 );
 
        atmdat.nsbig = MAX2(dense.IonHigh[nelem]+1,atmdat.nsbig);
 
        /* this routine only does simple two-level species, 
         * loop over ions will be <= limit, IonHigh is -1 since very
         * highest stage of ionization is not recombined into.  
         * for Li, will do only atom, since ions are H and He like,
         * so limit is zero */
 
        /* zero-out loop comes before main loop since there are off-diagonal
         * elements in the main ionization loop, due to multi-electron processes */
        /* >>chng 00 dec 07, limit changed to identical to ion_solver */
        for( ion=0; ion <= dense.IonHigh[nelem]-1; ion++ )
        {
                ionbal.RateRecomTot[nelem][ion] = 0.;
        }
 
        for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
        {
                if( ((dense.IonHigh[nelem] >= nelem - ipISO) &&
                          (dense.IonLow[nelem] <= nelem - ipISO)) )
                {
                        ionbal.RateRecomTot[nelem][nelem-ipISO] = ionbal.RateRecomIso[nelem][ipISO];
                }
        }
 
        limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
        ASSERT( limit >= -1 );
 
        /* these are counted elsewhere and must not be added here */
        Heavy.xLyaHeavy[nelem][nelem] = 0.;
        Heavy.xLyaHeavy[nelem][nelem-1] = 0.;
 
        /* IonLow is 0 for the atom, limit chosen to NOT include iso sequences  */
        for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
        {
                /* number of bound electrons of the ion after recombination,
                 * for an atom (ion=0) this is
                 * equal to nelem+1, the element on the physical scale, since nelem is
                 * on the C scale, being 5 for carbon */
                ASSERT(ionbal.DR_Badnell_rate_coef[nelem][ion] >= 0);
                ASSERT(ionbal.RR_rate_coef_used[nelem][ion] >= 0);
 
                /* sum of recombination rates [units s-1] for radiative, three body, charge transfer */
                ionbal.RateRecomTot[nelem][ion] = 
                        dense.eden* (
                        ionbal.RR_rate_coef_used[nelem][ion] + 
                        ionbal.DR_Badnell_rate_coef[nelem][ion] +
                        ionbal.CotaRate[ion] );
 
                /* >>chng 01 jun 30, FRAC_LINE was 0.1, not 1, did not include anything except
                 * radiative recombination, the radrec term */
#               define FRAC_LINE 1.
                /* was 0.1 */
                /*Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden*radrec*FRAC_LINE );*/
                Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden*
                        (ionbal.RR_rate_coef_used[nelem][ion]+ionbal.DR_Badnell_rate_coef[nelem][ion])*FRAC_LINE );
        }
 
        /* option to save rec coefficients */
        if( save.lgioRecom || lgPrintIt )
        {
                /* >>chng 04 feb 22, make option to print ions for single element */
                FILE *ioOut;
                if( lgPrintIt )
                        ioOut = ioQQQ;
                else
                        ioOut = save.ioRecom;
 
                /* print name of element */
                fprintf( ioOut, 
                        " %s recombination coefficients fnzone:%.2f \tte\t%.4e\tne\t%.4e\n", 
                        elementnames.chElementName[nelem] , fnzone , phycon.te , dense.eden );
 
                /*limit = MIN2(11,dense.IonHigh[nelem]);*/
                /* >>chng 05 sep 24, just print one long line - need info */
                limit = dense.IonHigh[nelem];
                // give ion stage
                for( i=0; i<limit; ++i )
                        fprintf( ioOut, "%10ld",i+1);
                fprintf( ioOut, "\n");
                
                for( i=0; i < limit; i++ )
                {
                        fprintf( ioOut, "%10.2e", ionbal.RR_rate_coef_used[nelem][i] );
                }
                fprintf( ioOut, " radiative used vs Z\n" );
 
                for( i=0; i < limit; i++ )
                {
                        fprintf( ioOut, "%10.2e", ionbal.RR_Verner_rate_coef[nelem][i] );
                }
                fprintf( ioOut, " old Verner vs Z\n" );
 
                for( i=0; i < limit; i++ )
                {
                        fprintf( ioOut, "%10.2e", ionbal.RR_Badnell_rate_coef[nelem][i] );
                }
                fprintf( ioOut, " new Badnell vs Z\n" );
 
                for( i=0; i < limit; i++ )
                {
                        /* >>chng 06 jan 19, from div by eden to div by H0 - want units of cm3 s-1 but
                         * no single collider does this so not possible to get rate coefficient easily
                         * H0 is more appropriate than electron density */
                        fprintf( ioOut, "%10.2e", ionbal.CX_recomb_rate_used[nelem][i]/SDIV(dense.xIonDense[ipHYDROGEN][0]) );
                }
                fprintf( ioOut, " CT/n(H0)\n" );
 
                for( i=0; i < limit; i++ )
                {
                        fprintf( ioOut, "%10.2e", ionbal.CotaRate[ion] );
                }
                fprintf( ioOut, " 3body vs Z /ne\n" );
 
                /* note different upper limit - this routine does grain rec for all ions */
                for( i=0; i < dense.IonHigh[nelem]; i++ )
                {
                        fprintf( ioOut, "%10.2e", gv.GrainChTrRate[nelem][i+1][i]/dense.eden );
                }
                fprintf( ioOut, " Grain vs Z /ne\n" );
                fprintf( ioOut, " old Nussbaumer Storey DR vs Z\n" );
 
                for( i=0; i < limit; i++ )
                {
                        fprintf( ioOut, "%10.2e", ionbal.DR_Badnell_rate_coef[nelem][i] );
                }
                fprintf( ioOut, " new Badnell DR vs Z\n" );
 
                /* total recombination rate, with density included - this goes into the matrix */
                for( i=0; i < limit; i++ )
                {
                        fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
                }
                fprintf( ioOut, 
                        " total rec rate (with density) for %s\n", 
                        elementnames.chElementSym[nelem] );
                for( i=0; i < limit; i++ )
                {
                        fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i]/dense.eden );
                }
                fprintf( ioOut, 
                        " total rec rate / ne for %s\n\n", 
                        elementnames.chElementSym[nelem] );
 
                /* spill over to next line for many stages of ionization */
                if( dense.IonHigh[nelem] > 11 )
                {
                        limit = MIN2(29,dense.IonHigh[nelem]);
                        fprintf( ioOut, " R " );
                        for( i=11; i < limit; i++ )
                        {
                                fprintf( ioOut, "%10.2e", dense.eden*ionbal.CotaRate[ion] );
                        }
                        fprintf( ioOut, "\n" );
 
                        fprintf( ioOut, "   " );
                        for( i=11; i < limit; i++ )
                        {
                                fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
                        }
                        fprintf( ioOut, "\n\n" );
                }
        }
 
        /* >>chng 02 nov 09, from -2 to -NISO */
        /*limit = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/
        limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
        for( i=dense.IonLow[nelem]; i <= limit; i++ )
        {
                ASSERT( Heavy.xLyaHeavy[nelem][i] > 0. );
                ASSERT( ionbal.RateRecomTot[nelem][i] > 0. );
        }
        return;
}
 
/*ion_recombAGN generate recombination coefficients for AGN table */
void ion_recombAGN( FILE * io )
{
#       define N1LIM 3
#       define N2LIM 4
        double te1[N1LIM]={ 5000., 10000., 20000.};
        double te2[N2LIM]={ 20000.,50000.,100000.,1e6};
        /* this is boundary between two tables */
        double BreakEnergy = 100./13.0;
        long int nelem, ion , i;
        /* this will hold element symbol + ionization */
        char chString[100],
                chOutput[100];
        /* save temp here       */
        double TempSave = phycon.te;
        /* save ne here */
        double EdenSave = dense.eden;
 
        DEBUG_ENTRY( "ion_recomb(false,)" );
 
        EdenChange( 1. );
        /*atmdat_readin();*/
 
        /* first put header on file */
        fprintf(io,"X+i\\Te");
        for( i=0; i<N1LIM; ++i )
        {
                phycon.te = te1[i];
                fprintf(io,"\t%.0f K",phycon.te);
        }
        fprintf(io,"\n");
 
        /* now do loop over temp, but add elements */
        for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
        {
                /* this list of elements included in the AGN tables is defined in zeroabun.c */
                if( abund.lgAGN[nelem] )
                {
                        for( ion=0; ion<=nelem; ++ion )
                        {
                                ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
 
                                if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
                                        break;
 
                                /* print chemical symbol */
                                sprintf(chOutput,"%s", 
                                        elementnames.chElementSym[nelem]);
                                /* some elements have only one letter - this avoids leaving a space */
                                if( chOutput[1]==' ' )
                                        chOutput[1] = chOutput[2];
                                /* now ionization stage */
                                if( ion==0 )
                                {
                                        sprintf(chString,"0 ");
                                }
                                else if( ion==1 )
                                {
                                        sprintf(chString,"+ ");
                                }
                                else
                                {
                                        sprintf(chString,"+%li ",ion);
                                }
                                strcat( chOutput , chString );
                                fprintf(io,"%5s",chOutput );
 
                                for( i=0; i<N1LIM; ++i )
                                {
                                        TempChange(te1[i] , false);
                                        dense.IonLow[nelem] = 0;
                                        dense.IonHigh[nelem] = nelem+1;
                                        if( ConvBase(0) )
                                                fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
                                        fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
                                }
                                fprintf(io,"\n");
                        }
                        fprintf(io,"\n");
                }
        }
 
        /* second put header on file */
        fprintf(io,"X+i\\Te");
        for( i=0; i<N2LIM; ++i )
        {
                TempChange(te2[i] , false);
                fprintf(io,"\t%.0f K",phycon.te);
        }
        fprintf(io,"\n");
 
        /* now do same loop over temp, but add elements */
        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
        {
                /* this list of elements included in the AGN tables is defined in zeroabun.c */
                if( abund.lgAGN[nelem] )
                {
                        for( ion=0; ion<=nelem; ++ion )
                        {
                                ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
 
                                if( Heavy.Valence_IP_Ryd[nelem][ion] <= BreakEnergy )
                                        continue;
 
                                /* print chemical symbol */
                                fprintf(io,"%s", 
                                        elementnames.chElementSym[nelem]);
                                /* now ionization stage */
                                if( ion==0 )
                                {
                                        fprintf(io,"0 ");
                                }
                                else if( ion==1 )
                                {
                                        fprintf(io,"+ ");
                                }
                                else
                                {
                                        fprintf(io,"+%li",ion);
                                }
 
                                for( i=0; i<N2LIM; ++i )
                                {
                                        TempChange(te2[i] , false);
                                        dense.IonLow[nelem] = 0;
                                        dense.IonHigh[nelem] = nelem+1;
                                        if( ConvBase(0) )
                                                fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
                                        fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
                                }
                                fprintf(io,"\n");
                        }
                        fprintf(io,"\n");
                }
        }
 
        TempChange(TempSave , true);
        EdenChange( EdenSave );
        return;
}