save_linedata.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 */
/*SaveLineData punches selected line data for all lines transferred in code */
/*Save1LineData save data for one line */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "phycon.h"
#include "physconst.h"
#include "elementnames.h"
#include "hydrogenic.h"
#include "lines_service.h"
#include "dense.h"
#include "atomfeii.h"
#include "conv.h"
#include "lines.h"
#include "atmdat.h"
#include "prt.h"
#include "h2.h"
#include "thermal.h"
#include "cooling.h"
#include "save.h"
#include "mole.h"
 
NORETURN void SaveLineData(FILE * ioPUN)
{
 
        long int i, 
          j, 
          limit ,
          nelem ,
          ipHi , 
          ipLo;
 
        const long nskip=2; /* number of emission lines per line of output */
        double tot;
        bool lgElemOff=false;
 
        DEBUG_ENTRY( "SaveLineData()" );
 
        /* routine punches out (on unit ioPUN) line data
         * for all recombination lines, and all transitions that are transferred */
 
        /* say what is happening so we know why we stopped */
        fprintf( ioQQQ, " saving line data, then stopping\n" );
 
        /* first check that all lines are turned on */
        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
        {
                if( !dense.lgElmtOn[nelem] )
                {
                        fprintf(ioQQQ," WARNING - I am saving line data but element %s is turned off.\n",
                        elementnames.chElementName[nelem]);
                        lgElemOff = true;
                }
        }
        if( lgElemOff )
        {
                fprintf(ioQQQ,"Some elements are turned off and save line data requested.\n");
                fprintf(ioQQQ,"Code is now designed to do save line data only with all elements on.\n");
                fprintf(ioQQQ,"Please try again with all elements on.\n");
                fprintf(ioQQQ,"Please try again with all elements on.\n");
                cdEXIT(EXIT_FAILURE);
        }
 
        /* evaluate rec coefficient at constant temperature if this is set, else
         * use 10,000K */
        double TeNew;
        if( thermal.lgTemperatureConstant )
        {
                TeNew = thermal.ConstTemp;
        }
        else
        {
                TeNew = 1e4;
        }
        TempChange(TeNew , false);
 
        /* this is set of Dima's recombination lines */
        t_ADfA::Inst().rec_lines(phycon.te,LineSave.RecCoefCNO);
        fprintf( ioPUN, "\n       Recombination lines of C, N, O\n" );
        fprintf( ioPUN, "#Ion\tWL(A)\tCoef\tIon\tWL(A)\tCoef\n" );
        for( i=0; i<471; i+=nskip)
        {
                /* nskip is set to 2 above */
                limit = MIN2(471,i+nskip);
                fprintf( ioPUN, "    " );
                for( j=i; j < limit; j++ )
                {
                        fprintf( ioPUN, "%2.2s%2ld\t%6ld\t%8.3f\t",
                                elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][j])-1], 
                          (long)(LineSave.RecCoefCNO[0][j]-LineSave.RecCoefCNO[1][j]+1.01), 
                          (long)(LineSave.RecCoefCNO[2][j]+0.5), 
                          log10(SDIV(LineSave.RecCoefCNO[3][j]) ) );
                }
                fprintf( ioPUN, "    \n" );
        }
        fprintf( ioPUN, "\n\n" );
 
        dense.SetGasPhaseDensity( ipHYDROGEN, 1. );
        dense.EdenHCorr = 1.;
        EdenChange( 1. );
 
        /* want very small neutral fractions so get mostly e- cs */
        dense.xIonDense[ipHYDROGEN][1] = 1.e-5f;
        findspecieslocal("H2")->den = 0.;
        dense.xIonDense[ipHYDROGEN][1] = 1.;
        for( i=1; i <= nLevel1; i++ )
        {
                TauLines[i].Lo()->Pop() = 1.;
        }
 
        for( i=0; i < nWindLine; i++ )
        {
                TauLine2[i].Lo()->Pop() = 1.;
        }
 
        for( i=0; i < nUTA; i++ )
        {
                (*UTALines[i].Lo()).Pop() = 1.;
        }
 
        for( i=0; i < LIMELM; i++ )
        {
                for( j=0; j < LIMELM+1; j++ )
                {
                        dense.xIonDense[i][j] = 1.;
                }
        }
 
        /* evaluate cooling, this forces evaluation of collision strengths */
        CoolEvaluate(&tot);
 
        fprintf( ioPUN, "       Level 1 transferred lines\n" );
        bool lgPrint = true;
        for( i=1; i <= nLevel1; i++ )
        {
                /* chLineLbl generates a 1 char string from the line transfer array info -
                 * "Ne 2  128" the string is null terminated,
                 * in following printout the first 4 char is used first, followed by
                 * an integer, followed by the rest of the array*/
                Save1LineData( TauLines[i] , ioPUN , true , lgPrint );
        }
 
        fprintf( ioPUN, "\n\n\n end level 1, start level 2\n" );
        lgPrint = true;
        for( i=0; i < nWindLine; i++ )
        {
                if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
                {
                        Save1LineData( TauLine2[i] , ioPUN , true , lgPrint );
                }
        }
 
        fprintf( ioPUN, "\n\n\n  end level 2, start inner shell UTA\n" );
        lgPrint = true;
        for( i=0; i < nUTA; i++ )
        {
                Save1LineData( UTALines[i] , ioPUN , true , lgPrint);
        }
 
        fprintf( ioPUN, "\n\n\n  end inner shell, start H-like iso seq\n" );
        /* h-like iso sequence */
        /* the hydrogen like iso-sequence */
        lgPrint = true;
        for( nelem=0; nelem < LIMELM; nelem++ )
        {
                iso_collide( ipH_LIKE, nelem );
                if( nelem < 2 || dense.lgElmtOn[nelem] )
                {
                        /* arrays are dim'd iso_sp[ipH_LIKE][nelem].numLevels_max+1 */
                        /* keep this limit to iso.numLevels_max, instead of _local.  */
                        for( ipLo=ipH1s; ipLo < iso_sp[ipH_LIKE][nelem].numLevels_max-1; ipLo++ )
                        {
                                for( ipHi=ipLo+1; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_max; ipHi++ )
                                {
                                        Save1LineData( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo) , ioPUN , false,lgPrint );
                                }
                        }
                }
        }
 
        fprintf( ioPUN, "\n\n\n  end H-like iso seq, start He-like iso seq\n" );
        lgPrint = true;
        for( nelem=1; nelem < LIMELM; nelem++ )
        {
                if( nelem < 2 || dense.lgElmtOn[nelem] )
                {
                        /* arrays are dim'd iso_sp[ipH_LIKE][nelem].numLevels_max+1  */
                        for( ipLo=ipHe1s1S; ipLo < iso_sp[ipHE_LIKE][nelem].numLevels_max-1; ipLo++ )
                        {
                                for( ipHi=ipLo+1; ipHi < iso_sp[ipHE_LIKE][nelem].numLevels_max; ipHi++ )
                                {
                                        Save1LineData( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo) , ioPUN , false,lgPrint );
                                }
                        }
                }
        }
 
        fprintf( ioPUN, "\n\n\n   end he-like iso seq, start hyperfine structure lines\n" );
        /* fine structure lines */
        lgPrint = true;
        for( i=0; i < nHFLines; i++ )
        {
                Save1LineData( HFLines[i] , ioPUN , true , lgPrint);
        }
 
        fprintf( ioPUN, "\n\n\n end hyperfine, start database lines\n" );
        /* Databases: Atoms & Molecules*/
        lgPrint = true;
        for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
        {
                for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
                          em != dBaseTrans[ipSpecies].Emis().end(); ++em)
                {
                        Save1LineData( (*em).Tran() , ioPUN , true , lgPrint);
                }
        }
 
        fprintf( ioPUN, "\n\n\n  end database, start satellite lines\n" );
        lgPrint = true;
        for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
        {
                for( long nelem = ipISO; nelem < LIMELM; nelem++ )
                {
                        if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] ) 
                        {
                                for( i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
                                {
                                        Save1LineData( SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][i]],
                                                        ioPUN , true , lgPrint );
                                }
                        }
                }
        }
 
        /* want very small ionized fractions so get mostly H2 cs */
        dense.SetGasPhaseDensity( ipHYDROGEN, 1e-6f );
        dense.EdenHCorr = 1e-6f;
        findspecieslocal("H2")->den = 1.;
        findspecieslocal("H2*")->den = 1.;
        dense.xIonDense[ipHYDROGEN][1] = 1e-6;
        EdenChange( 1e-6 );
 
        /* H2 molecule */
        fprintf( ioPUN, "\n\n\n" );
        fprintf( ioPUN, "       end satellite, start H2 lines\n" );
 
        /* ioPUN unit, and option to print all possible lines - false indicates
         * save only significant lines */
        for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
        {
                bool lgPopsConverged;
                double old_val, new_val;
                (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
                (*diatom)->H2_Punch_line_data( ioPUN, false );
        }
 
        /* save FeII data if atom is currently enabled */
        fprintf( ioPUN, "\n\n\n" );
        fprintf( ioPUN, "       end H2, start FeII lines\n" );
 
        /* ioPUN unit, and option to print all possible lines - false indicates
         * save only significant lines */
        FeIIPunData( ioPUN , false );
 
        /* ChkMonitorstring is searched for by one of the scripts in the nightly run
         * this run will terminate with no asserts but that is the correct behavior */
        fprintf( ioQQQ , "\n The code is left in a disturbed state after creating the SAVE LINE DATA file.\n"
                        " No calculation is actually performed, only the SAVE LINE DATA file is produced.\n"
                        " Remove the SAVE LINE DATA command to do the calculation.\n\n ChkMonitorend is ok.\n" );
 
        /* stop when done, we have done some serious damage to the code */
        cdEXIT(EXIT_SUCCESS);
}
 
/*Save1LineData save data for one line */
void Save1LineData( const TransitionProxy& t , FILE * ioPUN  , 
        /* flag saying whether to give collision strength too - in multi level atoms
         * it will be not valid without a great deal more work */
         bool lgCS_2 ,
         /* print header line? */
         bool &lgPrintHead )
{
 
        double CritDen;
        /* these are used for parts of the line label */
        char chLbl[11];
 
        DEBUG_ENTRY( "Save1LineData()" );
 
        if( lgPrintHead )
                fprintf( ioPUN, "#Ion\tWL\tgl\tgu\tgf\tA\tCS\tn(crt)\tdamp\n" );
        lgPrintHead = false;
 
        if( t.ipCont() <= 0 )
        {
                // not a real line, just give \n
                return;
        }
 
        /*iWL = iWavLen( t , &chUnits , &chShift );*/
        /* ion label, like C  1 */
        chIonLbl(chLbl , t );
        fprintf(ioPUN,"%s\t", chLbl );
 
        /* this is the second piece of the line label, pick up string after start */
 
        /* the wavelength */
        if( strcmp( save.chConPunEnr[save.ipConPun], "labl" )== 0 )
        {
                prt_wl( ioPUN , t.WLAng() );
        }
        else
        {
                /* this converts energy in Rydbergs into any of the other units */
                fprintf( ioPUN , "%.5e", AnuUnit((realnum)(t.EnergyRyd())) );
        }
 
        fprintf( ioPUN, "\t%3ld\t%3ld",
          /* lower and upper stat weights */
                                (long)((*t.Lo()).g()), 
                                (long)((*t.Hi()).g()) );
 
        /* oscillator strength */
        fprintf( ioPUN,PrintEfmt("\t%9.2e",  t.Emis().gf()));
 
        /* Einstein A for transition */
        fprintf( ioPUN,PrintEfmt("\t%9.2e",  t.Emis().Aul()));
 
        /* next collision strengths, use different formats depending on size 
         * of collision strength */
        if( t.Coll().col_str() > 100. )
        {
                fprintf( ioPUN, "\t%7.1f", t.Coll().col_str() );
        }
        else if( t.Coll().col_str() > 10. )
        {
                fprintf( ioPUN, "\t%7.2f", t.Coll().col_str() );
        }
        else if( t.Coll().col_str() > 1. )
        {
                fprintf( ioPUN, "\t%7.3f", t.Coll().col_str() );
        }
        else if( t.Coll().col_str() > .01 )
        {
                fprintf( ioPUN, "\t%7.4f", t.Coll().col_str() );
        }
        else if( t.Coll().col_str() > 0.0 )
        {
                fprintf( ioPUN, "\t%.3e", t.Coll().col_str() );
        }
        else
        {
                fprintf( ioPUN, "\t%7.4f", 0. );
        }
 
        /* now print critical density but only if cs is positive 
         * >>chng 06 mar 24, add flag lgCS_2 - in multi-level systems do not want
         * to save cs since not computed properly */
        if( lgCS_2 && t.Coll().col_str()> 0. )
        {
                CritDen = t.Emis().Aul() * (*t.Hi()).g()*phycon.sqrte / (t.Coll().col_str()*COLL_CONST);
        }
        else
        {
                CritDen = 0.;
        }
        fprintf( ioPUN, "\t%.3e",CritDen );
 
        // damping constant for current conditions
        fprintf( ioPUN,PrintEfmt("\t%9.2e",  t.Emis().damp() ));
 
        fprintf( ioPUN, "\n" );
 
        return;
}