root/branches/newmole/source/atom_feii.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 */
/*FeII_Colden maintain H2 column densities within X */
/*FeIILevelPops  main FeII routine, called by CoolIron to evaluate iron cooling */
/*FeIICreate read in needed data from file 
 * convert form of FeII data, as read in from file within routine FeIICreate
 * into physical form.  called by atmdat_readin  */
/*FeIIPunPop - punch level populations */
/*AssertFeIIDep called by assert FeII depart coef command */
/*FeIIPrint print FeII information */
/*FeIICollRatesBoltzmann evaluate collision strenths, 
 * both interpolating on r-mat and creating g-bar 
 * find Boltzmann factors, evaluate collisional rate coefficients */
/*FeIIPrint print output from large FeII atom, called by prtzone */
/*FeIISumBand sum up large FeII emission over certain bands, called in lineset4 */
/*FeII_RT_TauInc called once per zone in RT_tau_inc to increment large FeII atom line optical depths */
/*FeII_RT_tau_reset reset optical depths for large FeII atom, called by update after each iteration */
/*FeIIPoint called by ContCreatePointers to create pointers for lines in large FeII atom */
/*FeIIAccel called by rt_line_driving to compute radiative acceleration due to FeII lines */
/*FeIIAddLines save accumulated FeII intensities called by lineset4 */
/*FeIIPunchLines punch FeII lines at end of calculation, if punch verner set, called by punch_do*/
/*FeIIPunchOpticalDepth punch FeII line optical depth at end of calculation, called by punch_do*/
/*FeIIPunchLevels punch FeII levels and energies */
/*FeII_LineZero zero out storage for large FeII atom, called by tauout */
/*FeIIIntenZero zero out intensity of FeII atom */
/*FeIIZero initialize some variables, called by zero one time before commands parsed */
/*FeIIReset reset some variables, called by zero */
/*FeIIPunData punch line data */ 
/*FeIIPunDepart punch some departure coef for large atom, 
 * set with punch FeII departure command*/
/*FeIIPun1Depart send the departure coef for physical level nPUN to unit ioPUN */
/*FeIIContCreate create FeII continuum bins to add lines into ncell cells between wavelengths lambda low and high,
 * returns number of cells used */
/*FeIIBandsCreate returns number of FeII bands */
/*FeII_RT_Out do outward rates for FeII, called by RT_diffuse */
/*FeII_OTS do ots rates for FeII, called by RT_OTS */
/*FeII_RT_Make called by RT_line_all, does large FeII atom radiative transfer */
/*FeIILyaPump find rate of Lya excitation of the FeII atom */
/*FeIIOvrLap handle overlapping FeII lines */
/*ParseAtomFeII parse the atom FeII command */
/*FeIIPunchLineStuff include FeII lines in punched optical depths, etc, called from PunchLineStuff */
#include "cddefines.h"
#include "cddrive.h"
#include "thermal.h"
#include "physconst.h"
#include "doppvel.h"
#include "taulines.h"
#include "dense.h"
#include "rfield.h"
#include "radius.h"
#include "lines_service.h"
#include "ipoint.h"
#include "thirdparty.h"
#include "hydrogenic.h"
#include "lines.h"
#include "rt.h"
#include "trace.h"
#include "punch.h"
#include "phycon.h"
#include "atomfeii.h"
#include "iso.h"
#include "pressure.h"

/* FeIIOvrLap handle overlapping FeII lines */
STATIC void FeIIOvrLap(void);

/* add FeII lines into ncell cells between wavelengths lambda low and high */
STATIC void FeIIContCreate(double xLamLow , double xLamHigh , long int ncell );

/* read in the FeII Bands file, and sets nFeIIBands, the number of bands,
 * if argument is "" then use default file with bands, 
 * if filename within "" then use it instead,
 * return value is 0 if success, 1 if failed */
STATIC int FeIIBandsCreate( const char chFile[] );

/* this will be the smallest collision strength we will permit with the gbar
 * collision strengths, or for the data that have zeroes */
/* >>chng 99 jul 24, this was 1e-9 in the old fortran version and in baldwin et al. 96,
 * and has been changed several times, and affects results.  this is the smallest
 * non-zero collision strength in the r-matrix calculations */
realnum CS2SMALL = (realnum)1e-5;
/* routines used only within this file */

/*FeIICollRatesBoltzmann evaluate collision strenths, 
 * both interpolating on r-mat and creating g-bar 
 * find Boltzmann factors, evaluate collisional rate coefficients */
STATIC void FeIICollRatesBoltzmann(void);
/* find rate of Lya excitation of the FeII atom */
STATIC void FeIILyaPump(void);

/*extern realnum Fe2LevN[NFE2LEVN][NFE2LEVN][NTA];*/
/*extern realnum Fe2LevN[ipHi][ipLo].t[NTA];*/
/*realnum ***Fe2LevN;*/
/* >>chng 06 mar 01, boost to global namespace */
/*transition **Fe2LevN;*/
/* flag for the collision strength */
int **ncs1;

/* all following variables have file scope
#define NFEIILINES      68635 */

/* this is size of nnPradDat array */
#define NPRADDAT 159

/* band wavelength, lower and upper bounds, in vacuum Angstroms */
/* FeII_Bands[n][3], where n is the number of bands in fe2Bands.dat
 * these bands are defined in fe2Bands.dat and read in at startup
 * of calculation */
realnum **FeII_Bands; 

/* continuum wavelengths, lower and upper bounds, in vacuum Angstroms,
 * third is integrated intensity */
/* FeII_Cont[n][3], where n is the number of cells needed
 * these bands are defined in FeIIContCreate */
realnum **FeII_Cont; 

/* this is the number of bands read in from FeII_bands.ini */
long int nFeIIBands;

/* number of bands in continuum array */
long int nFeIIConBins;

/* the dim of this vector this needs one extra since there are 20 numbers per line,
 * highest not ever used for anything */
/*long int nnPradDat[NPRADDAT+1];*/
static long int *nnPradDat;

/* malloced in feiidata */
/* realnum sPradDat[NPRADDAT][NPRADDAT][8];*/
/* realnum sPradDat[ipHi][ipLo].t[8];*/
static realnum ***sPradDat;

/* array used to integrate FeII line intensities over model,
 * Fe2SavN[upper][lower],
 *static realnum Fe2SavN[NFE2LEVN][NFE2LEVN];*/
static double **Fe2SavN;

/* save effective transition rates */
static double **Fe2A;

/* induced transition rates */
static double **Fe2LPump, **Fe2CPump;

/* energies read in from fe2energies.dat data file */
realnum *Fe2Energies;

/* collision rates */
static realnum **Fe2Coll;

/* Fe2DepCoef[NFE2LEVN];
realnum cli[NFEIILINES], 
  cfe[NFEIILINES];*/
static double 
        /* departure coefficients */
        *Fe2DepCoef , 
        /* level populations */
        *Fe2LevelPop ,
        /* column densities */
        *Fe2ColDen ,
        /* this will become array of Boltzmann factors */
        *FeIIBoltzmann;
        /*FeIIBoltzmann[NFE2LEVN] ,*/

static double EnerLyaProf1, 
  EnerLyaProf4, 
  PhotOccNumLyaCenter;
static double 
                /* the yVector - will become level populations after matrix inversion */
                *yVector,
          /* this is used to call matrix routines */
          /*xMatrix[NFE2LEVN][NFE2LEVN] ,*/
          **xMatrix , 
          /* this will become the very large 1-D array that 
           * is passed to the matrix inversion routine*/
          *amat;


/*FeII_Colden maintain H2 column densities within X */
void FeII_Colden( const char *chLabel )
{
        long int n;

        DEBUG_ENTRY( "FeII_Colden()" );

        /* >>chng 05 nov 29, FeII always on, always want to evaluate this */
        /* nothing to do if no FeII 
        if( !FeII. lgFeIION )
                return;*/

        if( strcmp(chLabel,"ZERO") == 0 )
        {
                /* zero out column density */
                for( n=0; n < FeII.nFeIILevel; ++n )
                {
                        /* space for the rotation quantum number */
                        Fe2ColDen[n] = 0.;
                }
        }

        else if( strcmp(chLabel,"ADD ") == 0 )
        {
                /*  add together column densities */
                for( n=0; n < FeII.nFeIILevel; ++n )
                {
                        /* state-specific FeII column density */
                        Fe2ColDen[n] += Fe2LevelPop[n]*radius.drad_x_fillfac;
                }
        }

        /* check for the print option, which we can't do, else  we have a problem */
        else if( strcmp(chLabel,"PRIN") != 0 )
        {
                fprintf( ioQQQ, " FeII_Colden does not understand the label %s\n", 
                  chLabel );
                cdEXIT(EXIT_FAILURE);
        }

        return;
}

/*
 *=====================================================================
 */
/* FeIICreate read in FeII data from files on disk.  called by atmdat_readin 
 * but only if FeII. lgFeIION is true, set with atom FeII verner command */
void FeIICreate(void)
{
        FILE *ioDATA;

        char chLine[FILENAME_PATH_LENGTH_2];

        long int i, 
          ipHi ,
          ipLo,
          lo,
          ihi,
          k, 
          m1, 
          m2, 
          m3;

        DEBUG_ENTRY( "FeIICreate()" );

        if( lgFeIIMalloc )
        {
                /* we have already been called one time, just bail out */

                return;
        }

        /* now set flag so never done again - this is set false when initi
         * when this is true it is no longer possible to change the number of levels
         * in the model atom with the atom FeII levels command */
        lgFeIIMalloc = true;

        /* remember how many levels this was, so that in future calculations
         * we can reset the atom to full value */
        FeII.nFeIILevelAlloc = FeII.nFeIILevel;

        /* set up array to save FeII transition probabilities */
        Fe2A = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc );

        /* second dimension, lower level, for line save array */
        for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ++ipHi )
        {
                Fe2A[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevelAlloc);
        }

        /* set up array to save FeII pumping rates */
        Fe2CPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc );

        /* set up array to save FeII pumping rates */
        Fe2LPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc );

        /* second dimension, lower level, for line save array */
        for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ++ipHi )
        {
                Fe2CPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevelAlloc);

                Fe2LPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevelAlloc);
        }

        /* set up array to save FeII collision rates */
        Fe2Energies = (realnum *)MALLOC(sizeof(realnum)*(unsigned long)FeII.nFeIILevelAlloc );

        /* set up array to save FeII collision rates */
        Fe2Coll = (realnum **)MALLOC(sizeof(realnum *)*(unsigned long)FeII.nFeIILevelAlloc );

        /* second dimension, lower level, for line save array */
        for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ++ipHi )
        {
                Fe2Coll[ipHi]=(realnum*)MALLOC(sizeof(realnum )*(unsigned long)FeII.nFeIILevelAlloc);
        }

        /* MALLOC space for the 2D matrix array */
        xMatrix = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc );

        /* now do the second dimension */
        for( i=0; i<FeII.nFeIILevelAlloc; ++i )
        {
                xMatrix[i] = (double *)MALLOC(sizeof(double)*(unsigned long)FeII.nFeIILevelAlloc );
        }
        /* MALLOC space for the  1-yVector array */
        amat=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevelAlloc*FeII.nFeIILevelAlloc) ));

        /* MALLOC space for the  1-yVector array */
        yVector=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevelAlloc) ));

        /* set up array to save FeII line intensities */
        Fe2SavN = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc );

        /* second dimension, lower level, for line save array */
        for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ++ipHi )
        {
                Fe2SavN[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)ipHi);
        }

        /* now MALLOC space for energy level table*/
        nnPradDat = (long*)MALLOC( (NPRADDAT+1)*sizeof(long) );

        /*Fe2DepCoef[NFE2LEVN];*/
        Fe2DepCoef = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) );

        /*Fe2LevelPop[NFE2LEVN];*/
        Fe2LevelPop = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) );

        /*Fe2ColDen[NFE2LEVN];*/
        Fe2ColDen = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) );

        /*FeIIBoltzmann[NFE2LEVN];*/
        FeIIBoltzmann = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) );


        /* MALLOC the realnum sPradDat[NPRADDAT][NPRADDAT][8] array */
        /* MALLOC the realnum sPradDat[ipHi][ipLo].t[8] array */
        sPradDat = ((realnum ***)MALLOC(NPRADDAT*sizeof(realnum **)));

        /* >>chng 00 dec 06, changed lower limit of loop to 1, Tru64 alpha's will not allocate 0 bytes!, PvH */
        sPradDat[0] = NULL;
        for(ipHi=1; ipHi < NPRADDAT; ipHi++) 
        {
                /* >>chng 00 dec 06, changed sizeof(realnum) to sizeof(realnum*), PvH */
                sPradDat[ipHi] = (realnum **)MALLOC((unsigned long)ipHi*sizeof(realnum *));

                /* now make space for the second dimension */
                for( ipLo=0; ipLo< ipHi; ipLo++ )
                { 
                        sPradDat[ipHi][ipLo] = (realnum *)MALLOC(8*sizeof(realnum ));
                }
        }

        /* now set junk to make sure reset before used */
        for(ipHi=0; ipHi < NPRADDAT; ipHi++) 
        {
                for( ipLo=0; ipLo< ipHi; ipLo++ )
                { 
                        for( k=0; k<8; ++k )
                        {
                                sPradDat[ipHi][ipLo][k] = -FLT_MAX;
                        }
                }
        }

        /* now create the main emission line array and a helper array for the cs flag  */
        Fe2LevN=(transition**)MALLOC(sizeof(transition *)*(unsigned long)FeII.nFeIILevelAlloc);
        ncs1=(int**)MALLOC(sizeof(int *)*(unsigned long)FeII.nFeIILevelAlloc);

        for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ++ipHi )
        {
                Fe2LevN[ipHi]=(transition*)MALLOC(sizeof(transition)*(unsigned long)ipHi);

                ncs1[ipHi]=(int*)MALLOC(sizeof(int)*(unsigned long)ipHi);
        }

#if     0
        /* now that its created, set to junk */
        for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ )
        {
                for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
                {
                        TransitionJunk( &Fe2LevN[ipHi][ipLo] );
                }
        }
#endif

        /* now assign state and Emis pointers */
        Fe2LevN[1][0].Lo = AddState2Stack();
        /* now that its created, set to junk */
        for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
        {
                /* add this upper level */
                Fe2LevN[ipHi][0].Hi = AddState2Stack();
                for( ipLo=0; ipLo < ipHi; ipLo++ )
                {
                        if( ipLo == 0 )
                        {
                                Fe2LevN[ipHi][ipLo].Lo = Fe2LevN[1][0].Lo;
                        }
                        else
                        {
                                /* lower state is same as a previous upper state. */
                                Fe2LevN[ipHi][ipLo].Lo = Fe2LevN[ipLo][0].Hi;
                        }

                        Fe2LevN[ipHi][ipLo].Hi = Fe2LevN[ipHi][0].Hi;
                        Fe2LevN[ipHi][ipLo].Emis = AddLine2Stack( true );
                }
        }

        if( trace.lgTrace )
        {
                fprintf( ioQQQ," FeIICreate opening fe2nn.dat:");
        }

        ioDATA = open_data( "fe2nn.dat", "r" );

        ASSERT( ioDATA !=NULL );
        /* read in the fe2nn.dat file - this gives the Zheng and Pradhan number of level
         * for every cloudy level.  So nnPradDat[1] is the index in the cloudy level
         * counting for level 1 of Zheng & Pradan
         * note that the order of some levels is different, the nnPradDat file goes 
         * 21 23 22 - also that many levels are missing, the file goes 95 99 94 93 116
         */
        for( i=0; i < 8; i++ )
        {
                if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                {
                        fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* get these integers from fe2nn.dat */
                k = 20*i;
                /* NPRADDAT is size of nnPradDat array, 159+1, make sure we do not exceed it */
                ASSERT( k+19 < NPRADDAT+1 );
                sscanf( chLine ,
                        "%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld",
                        &nnPradDat[k+0], &nnPradDat[k+1],  &nnPradDat[k+2], &nnPradDat[k+3], &nnPradDat[k+4],
                        &nnPradDat[k+5], &nnPradDat[k+6],  &nnPradDat[k+7], &nnPradDat[k+8], &nnPradDat[k+9],
                        &nnPradDat[k+10],&nnPradDat[k+11], &nnPradDat[k+12],&nnPradDat[k+13],&nnPradDat[k+14],
                        &nnPradDat[k+15],&nnPradDat[k+16], &nnPradDat[k+17],&nnPradDat[k+18],&nnPradDat[k+19]
                        );
#               if !defined(NDEBUG)
                for( m1=0; m1<20; ++m1 )
                {
                        ASSERT( nnPradDat[k+m1] >= 0 && nnPradDat[k+m1] <= NFE2LEVN );
                }
#               endif
        }
        fclose(ioDATA);

        /* now get energies */
        if( trace.lgTrace )
        {
                fprintf( ioQQQ," FeIICreate opening fe2energies.dat:");
        }

        ioDATA = open_data( "fe2energies.dat", "r" );

        /* file now open, read the data */
        for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
        {
                /* keep reading until have non-comment line, one that does not
                 * start with # */
                chLine[0] = '#';
                while( chLine[0] == '#' )
                {
                        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                        {
                                fprintf( ioQQQ, " fe2energies.dat error reading data\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                /* first and last number on this line */
                double help;
                sscanf( chLine, "%lf", &help );
                Fe2Energies[ipHi] = (realnum)help;
        }
        fclose(ioDATA);

        /* transition probabilities */

        if( trace.lgTrace )
                fprintf( ioQQQ," FeIICreate opening fe2rad.dat:");

        ioDATA = open_data( "fe2rad.dat", "r" );

        /* get the first line, this is a version number */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " fe2rad.dat error reading data\n" );
                cdEXIT(EXIT_FAILURE);
        }
        /* scan off three integers */
        sscanf( chLine ,"%ld%ld%ld",&lo, &ihi, &m1);
        const int nYrFe2Rad=8, nMoFe2Rad=8, nDyFe2Rad=24;
        if( lo!=nYrFe2Rad || ihi!=nMoFe2Rad || m1!=nDyFe2Rad )
        {
                fprintf( ioQQQ, "DISASTER fe2rad.dat has the wrong magic numbers, expected "
                        "%2i %2i %2i and got %2li %2li %2li\n" ,
                        nYrFe2Rad, nMoFe2Rad, nDyFe2Rad,
                        lo, ihi, m1);
                cdEXIT(EXIT_FAILURE);
        }

        /* file now open, read the data */
        for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
        {
                for( ipLo=0; ipLo < ipHi; ipLo++ )
                {
                        /* following double since used in sscanf */
                        double save[2];
                        /* keep reading until have non-comment line, one that does not
                         * start with # */
                        chLine[0] = '#';
                        while( chLine[0] == '#' )
                        {
                                if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                                {
                                        fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }

                        /* first and last number on this line */
                        sscanf( chLine ,
                                "%ld%ld%ld%ld%lf%lf%ld",
                                &lo, &ihi, &m1, &m2 ,
                                &save[0], &save[1] , &m3);
                        /* the indices ihi and lo within this array were 
                         * read in from the line above */
                        Fe2LevN[ihi-1][lo-1].Lo->g = (realnum)m1;
                        Fe2LevN[ihi-1][lo-1].Hi->g = (realnum)m2;
                        Fe2LevN[ihi-1][lo-1].Emis->Aul = (realnum)save[0];
                        /*>>chng 06 apr 10, option to use table of energies */
#                       define  USE_OLD true
                        if( USE_OLD )
                                Fe2LevN[ihi-1][lo-1].EnergyWN = (realnum)save[1];
                        else
                                Fe2LevN[ihi-1][lo-1].EnergyWN = Fe2Energies[ihi-1]-Fe2Energies[lo-1];

                        /* Aul == -1 indicates intercombination line with no real Aul */
                        if( fp_equal( Fe2LevN[ihi-1][lo-1].Emis->Aul , (realnum)-1.0f ) )
                        {
                                /* these are made-up intercombination lines, set gf to 1e-5 */
                                Fe2LevN[ihi-1][lo-1].Emis->gf = 1e-5f;
                                
                                /* get corresponding A */
                                Fe2LevN[ihi-1][lo-1].Emis->Aul = Fe2LevN[ihi-1][lo-1].Emis->gf*(realnum)TRANS_PROB_CONST*
                                        POW2(Fe2LevN[ihi-1][lo-1].EnergyWN)*Fe2LevN[ihi-1][lo-1].Lo->g/Fe2LevN[ihi-1][lo-1].Hi->g;
                        }

                        /* the last column of fe2rad.dat, and is 1, 2, or 3.  
                         * 1 signifies that transition is permitted,
                         * 2 is semi-forbidden
                         * 3 forbidden, within lowest 63 levels are forbidden, first permitted
                         * transition is from 64 */
                        ncs1[ihi-1][lo-1] = (int)m3;
                }
        }
        fclose(ioDATA);

        /* now read collision data in fe2col.dat 
         * These are from the following sources
         >>refer        fe2     cs      Zhang, H.L., & Pradhan, A.K., 1995, A&A 293, 953 
         >>refer        fe2     cs      Bautista, M., (private communication), 
         >>refer        fe2     cs      Mewe, R., 1972, A&AS 20, 215 (the g-bar approximation)
         */

        if( trace.lgTrace )
        {
                fprintf( ioQQQ," FeIICreate opening fe2col.dat:");
        }

        ioDATA = open_data( "fe2col.dat", "r" );

        ASSERT( ioDATA !=NULL);
        for( ipHi=1; ipHi<NPRADDAT; ++ipHi )
        {
                for( ipLo=0; ipLo<ipHi; ++ipLo )
                {
                        /* double since used in sscanf */
                        double save[8];
                        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
                        {
                                fprintf( ioQQQ, " fe2col.dat error reading data\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        sscanf( chLine ,
                                "%ld%ld%lf%lf%lf%lf%lf%lf%lf%lf",
                                &m1, &m2,
                                &save[0], &save[1] , &save[2],&save[3], &save[4] , &save[5],
                                &save[6], &save[7] 
                                );
                        for( k=0; k<8; ++k )
                        {
                                /* the max is here because there are some zeroes in the data file.
                                 * this is unphysical but is part of their distribution.  As a result
                                 * don't let the cs fall below 0.01 */
                                sPradDat[m2-1][m1-1][k] = max(CS2SMALL , (realnum)save[k] );
                        }
                }
        }

        fclose( ioDATA );

        /*generate needed opacity data for the large FeII atom */

        /* this routine is called only one time when cloudy is init
         * for the very first time.  It converts the FeII data stored
         * in the large FeII arrays into the array storage needed by cloudy
         * for its line optical depth arrays
         */

        /* convert large FeII line arrays into standard heavy el ar */
        for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ )
        {
                for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
                {
                        /* pull information out of existing data arrays */

                        ASSERT( Fe2LevN[ipHi][ipLo].EnergyWN > 0. );
                        ASSERT( Fe2LevN[ipHi][ipLo].Emis->Aul> 0. );

                        /* now put into standard internal line format
                        Fe2LevN[ipHi][ipLo].WLAng = (realnum)(1.e8/Fe2LevN[ipHi][ipLo].EnergyWN); */
                        /* >>chng 03 oct 28, above neglected index of refraction of air -
                         * convert to below */
                        Fe2LevN[ipHi][ipLo].WLAng = 
                                (realnum)(1.0e8/
                        Fe2LevN[ipHi][ipLo].EnergyWN/
                        RefIndex( Fe2LevN[ipHi][ipLo].EnergyWN ));

                        /* generate gf from A */
                        Fe2LevN[ipHi][ipLo].Emis->gf = 
                                (realnum)(Fe2LevN[ipHi][ipLo].Emis->Aul*Fe2LevN[ipHi][ipLo].Hi->g/
                                TRANS_PROB_CONST/POW2(Fe2LevN[ipHi][ipLo].EnergyWN));

                        Fe2LevN[ipHi][ipLo].Hi->IonStg = 2;
                        Fe2LevN[ipHi][ipLo].Hi->nelem = 26;
                        /* which redistribution function??  
                         * For resonance line use K2 (-1),
                         * for subordinate lines use CRD with wings */
                        /* >>chng 01 mar 09, all had been 1, inc redis with wings */
                        /* to reset this, so that the code works as it did pre 01 mar 29,
                         * use command 
                         * atom FeII redistribution resonance pdr 
                         * atom FeII redistribution subordinate pdr */
                        if( ipLo == 0 )
                        {
                                Fe2LevN[ipHi][ipLo].Emis->iRedisFun = FeII.ipRedisFcnResonance;
                        }
                        else
                        {
                                /* >>chng 01 feb 27, had been -1, crd with core only,
                                 * change to crd with wings as per discussion with Ivan Hubeny */
                                Fe2LevN[ipHi][ipLo].Emis->iRedisFun = FeII.ipRedisFcnSubordinate;
                        }
                        Fe2LevN[ipHi][ipLo].Emis->phots = 0.;
                        Fe2LevN[ipHi][ipLo].Emis->ots = 0.;
                        Fe2LevN[ipHi][ipLo].Emis->FracInwd = 1.;
                }
        }

        /* finally get FeII bands, this sets  */
        if( FeIIBandsCreate("") )
        {
                fprintf( ioQQQ," failed to open FeII bands file \n");
                cdEXIT(EXIT_FAILURE);
        }

        /* now establish the FeII continuum, these are set to
         * 1000, 7000, and 1000 in FeIIZero in this file, and
         * are reset with the atom FeII continuum command */
        FeIIContCreate( FeII.fe2con_wl1 , FeII.fe2con_wl2 , FeII.nfe2con );

        /* this must be true */
        ASSERT( FeII.nFeIILevelAlloc == FeII.nFeIILevel );

        return;
}

/*
 *=====================================================================
 */
/***********************************************************************
 *** version of FeIILevelPops.f with overlap in procces 05.28.97 ooooooo
 ******change in common block *te* sqrte 05.28.97
 *******texc is fixed 03.28.97
 *********this version has work on overlap 
 *********this version has # of zones (ZoneCnt) 07.03.97
 *********taux - optical depth depends on iter correction 03.03.97
 *********calculate cooling (Fe2_large_cool) and output fecool from Cloudy 01.29.97god
 *********and cooling derivative (ddT_Fe2_large_cool)
 ************ this program for 371 level iron model 12/14/1995
 ************ this program for 371 level iron model 1/11/1996
 ************ this program for 371 level iron model 3/21/1996
 ************ this program without La 3/27/1996
 ************ this program for 371 level iron model 4/9/1996
 ************ includes:FeIICollRatesBoltzmann;cooling;overlapping in lines */
void FeIILevelPops( void )
{
        long int  i, 
                ipHi ,
                ipLo ,
                n;
        /* used in test for non-positive level populations */
        bool lgPopNeg;

        double 
          EnergyWN,
          pop ,
          sum;

        int32 info;
        int32 ipiv[NFE2LEVN];

        DEBUG_ENTRY( "FeIILevelPops()" );

        if( trace.lgTrace )
        {
                fprintf( ioQQQ,"   FeIILevelPops fe2 pops called\n");
        }

        /* FeII.lgSimulate was set true with simulate flag on atom FeII command,
         * for bebugging without actually calling the routine */
        if( FeII.lgSimulate )
        {

                return;
        }

        /* zero out some arrays */
        for( n=0; n<FeII.nFeIILevel; ++n)
        {
                for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo )
                {
                        Fe2CPump[ipLo][n] = 0.;
                        Fe2LPump[ipLo][n] = 0.;
                        Fe2A[ipLo][n] = 0.;
                        xMatrix[ipLo][n] = 0.;
                }
        }

        /* make the g-bar collision strengths and do linear interpolation on r-matrix data.
         * this also sets Boltzmann factors for all levels,
         * sets values of FeColl used below, but only if temp has changed */
        FeIICollRatesBoltzmann();

        /* pumping and spontantous decays */
        for( n=0; n<FeII.nFeIILevel; ++n)
        {
                for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
                {
                        /* continuum pumping rate from n to upper ipHi */
                        Fe2CPump[n][ipHi] = Fe2LevN[ipHi][n].Emis->pump;

                        /* continuum pumping rate from ipHi to lower n */
                        Fe2CPump[ipHi][n] = Fe2LevN[ipHi][n].Emis->pump*
                                Fe2LevN[ipHi][n].Hi->g/Fe2LevN[ipHi][n].Lo->g;

                        /* spontaneous decays */
                        Fe2A[ipHi][n] = Fe2LevN[ipHi][n].Emis->Aul*(Fe2LevN[ipHi][n].Emis->Pesc + Fe2LevN[ipHi][n].Emis->Pelec_esc +
                          Fe2LevN[ipHi][n].Emis->Pdest);
                }
        }

        /* now do pumping of atom by Lya */
        FeIILyaPump();

        /* ************************************************************************** 
         *
         * final rates into matrix 
         *
         ***************************************************************************/

        /* fill in xMatrix with matrix elements */
        for( n=0; n<FeII.nFeIILevel; ++n)
        {
                /* all processes leaving level n going down*/
                for( ipLo=0; ipLo<n; ++ipLo )
                {
                        xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipLo] + Fe2LPump[n][ipLo]+ Fe2A[n][ipLo] + 
                                Fe2Coll[n][ipLo]*dense.eden;
                }
                /* all processes leaving level n going up*/
                for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
                {
                        xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipHi] + Fe2LPump[n][ipHi] + Fe2Coll[n][ipHi]*dense.eden;
                }
                /* all processes entering level n from below*/
                for( ipLo=0; ipLo<n; ++ipLo )
                {
                        xMatrix[ipLo][n] = xMatrix[ipLo][n] - Fe2CPump[ipLo][n] - Fe2LPump[ipLo][n] - Fe2Coll[ipLo][n]*dense.eden;
                }
                /* all processes entering level n from above*/
                for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
                {
                        xMatrix[ipHi][n] = xMatrix[ipHi][n] - Fe2CPump[ipHi][n] - Fe2LPump[ipHi][n] - Fe2Coll[ipHi][n]*dense.eden - 
                                Fe2A[ipHi][n];
                }

                /* the top row of the matrix is the sum of level populations */
                xMatrix[n][0] = 1.0;
        }

        {
                /* option to print out entire matrix */
                enum {DEBUG_LOC=false};
                if( DEBUG_LOC )
                {
                        /* print the matrices */
                        for( n=0; n<FeII.nFeIILevel; ++n)
                        {
                                /*fprintf(ioQQQ,"\n");*/
                                /* now print the matrix*/
                                for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo )
                                {
                                        fprintf(ioQQQ," %.1e", xMatrix[n][ipLo]);
                                }
                                fprintf(ioQQQ,"\n");
                        }
                }
        }

        /* define the Y Vector.  The oth element is the sum of all level populations
         * adding up to the total population.  The remaining elements are the level
         * balance equations adding up to zero */
        yVector[0] = 1.0;
        for( n=1; n < FeII.nFeIILevel; n++ )
        {
                yVector[n] = 0.0;
        }

        /* create the 1-yVector array that will save vector,
         * this is the macro trick */
#       ifdef AMAT
#               undef AMAT
#       endif
#       define AMAT(I_,J_)      (*(amat+(I_)*FeII.nFeIILevel+(J_)))

        /* copy current contents of xMatrix array over to special amat array*/
        for( ipHi=0; ipHi < FeII.nFeIILevel; ipHi++ )
        {
                for( i=0; i < FeII.nFeIILevel; i++ )
                {
                        AMAT(i,ipHi) = xMatrix[i][ipHi];
                }
        }

        info = 0;

        /* do the linear algebra to find the level populations */
        getrf_wrapper(FeII.nFeIILevel, FeII.nFeIILevel, amat, FeII.nFeIILevel, ipiv, &info);
        getrs_wrapper('N', FeII.nFeIILevel, 1, amat, FeII.nFeIILevel, ipiv, yVector, FeII.nFeIILevel, &info);

        if( info != 0 )
        {
                fprintf( ioQQQ, "DISASTER FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* yVector now contains the level populations */

        /* this better be false after this loop - if not then non-positive level pops */
        lgPopNeg = false;
        /* copy all level pops over to Fe2LevelPop */
        for( ipLo=0; ipLo < FeII.nFeIILevel; ipLo++ )
        {
                if(yVector[ipLo] < 0. )
                {
                        lgPopNeg = true;
                        fprintf(ioQQQ,"PROBLEM FeIILevelPops finds non-positive level population, level is %ld pop is %g\n",
                                ipLo , yVector[ipLo] );
                }
                /* this is now correct level population, cm^-3 */
                Fe2LevelPop[ipLo] = yVector[ipLo] * dense.xIonDense[ipIRON][1];
        }
        if( lgPopNeg )
        {
                /* this is here to use the lgPopNeg value for something, if not here then
                 * lint and some compilers will note that var is never used */
                fprintf(ioQQQ , "PROBLEM FeIILevelPops exits with negative level populations.\n");
        }

        /* >>chng 06 jun 24, make sure remainder of populations up through max
         * limit are zero - this makes safe the case where the number
         * of levels actually computed has been trimmed down from largest
         * possible number of levels, for instance, in cool gas */
        for( ipLo=FeII.nFeIILevel; ipLo < FeII.nFeIILevelAlloc; ++ipLo )
        {
                Fe2LevelPop[ipLo] = 0.;
        }

        /* now set line opacities, intensities, and level populations 
         * >>chng 06 jun ipLo should go up to FeII.nFeIILevel-1 since this
         * is the largest lower level with non-zero population */
        for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
        {
                /* >>chng 06 jun 24, upper level should go to limit
                 * of all that were allocated */
                /*for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ )*/
                for( ipHi=ipLo+1; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
                {
                        /* >>chng 06 jun 24, in all of these the product 
                         * yVector[ipHi]*dense.xIonDense[ipIRON][1] has been replaced
                         * with Fe2LevelPop[ipLo] - should always have been this way,
                         * and saves a mult */
                        Fe2LevN[ipHi][ipLo].Emis->PopOpc = (Fe2LevelPop[ipLo] - 
                                Fe2LevelPop[ipHi]*Fe2LevN[ipHi][ipLo].Lo->g/F