atmdat.h

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2023 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */

#ifndef ATMDAT_H_
#define ATMDAT_H_

#include "container_classes.h"
#include "module.h"

/*This structure is specifically used to hold the collision data in the format given in the LEIDEN Database
The data is available as collision rate coefficients(cm3 s-1) over different temperatures*/
class CollRateCoeffArray
{
public:
        /*Array of temps*/
        vector<double> temps;
        /*Matrix of collision rates(temp,up,lo)*/
        multi_arr<double,3> collrates;
        
}  ;

/*This structure is specifically used to hold the collision data in the format given in the CHIANTI Database
The data is available as spline fits to the Maxwellian averaged collision strengths */
class CollSplinesArray
{
public:
        /*Matrix of spline fits(hi,lo,spline index)*
         *The first five columns gives the no of spline pts,transition type,gf value,delta E
         *& Scaling parameter ,in the specified order*/
        /*The transition type basically tells how the temperature and collision
        strengths have been scaled*/
        vector<double> collspline;
        vector<double> SplineSecDer;

        long nSplinePts = 0L;
        long intTranType = 0L;
        double EnergyDiff;
        double ScalingParam;
        CollSplinesArray() : EnergyDiff(0.), ScalingParam(0.) {}
};

/*This structure is specifically used to hold the collision data in the format given in the STOUT Database
The data are available as collision strengths and rates over different temperatures*/
struct StoutColls
{
private:
        long m_offset, m_ntemps;
public:
        /*Number of temps*/
        long ntemps()
        {
                return m_ntemps;
        }
        /*Reset values to safe defaults*/
        void junk()
        {
                m_offset = -1;
                m_ntemps = -1;
                lgIsRate = false;
        }
        // Set range of tabulated data corresponding to this level & collider
        void setslice(long offset, long ntemps)
        {
                m_offset = offset;
                m_ntemps = ntemps;
        }
        long offset() const
        {
                return m_offset;
        }
        /*Is this a deexcitation rate or collision strength*/
        bool lgIsRate;

};

class StoutCollArray
{
        multi_arr<StoutColls,2> m_a;
        vector<double> m_temps,m_collstrs;
        static long ilev(long ipHi, long ipLo)
        {
                // Can use triangular arrangement, as ipHi > ipLo
                return (ipHi*(ipHi-1))/2+ipLo;
        }
public:
        // Allocates background data structures
        void alloc(long nHi, long , long nCollider)
        {
                // Size of array is first level which *isn't* required
                m_a.alloc(ilev(nHi+1,0),nCollider);
                m_temps.resize(0);
                m_collstrs.resize(0);
        }
        // Initializes per-collision data to safe values
        void junk(long ipHi, long ipLo, long ipCollider)
        {
                m_a[ilev(ipHi,ipLo)][ipCollider].junk();
        }
        bool& lgIsRate(long ipHi, long ipLo, long ipCollider)
        {
                return m_a[ilev(ipHi,ipLo)][ipCollider].lgIsRate;
        }
        // Reserves space for the number of table points for this level and
        // collider
        void setpoints(long ipHi, long ipLo, long ipCollider,long npoints)
        {
                long offset = m_temps.size();
                long length = offset+npoints;
                // testing shows that 1200 is a reasonable initial guess for the size
                m_temps.resize(max(length,1200));
                m_collstrs.resize(max(length,1200));
                m_a[ilev(ipHi,ipLo)][ipCollider].setslice(offset,npoints);
        }
        // Returns number of table points
        long ntemps(long ipHi, long ipLo, long ipCollider)
        {
                return m_a[ilev(ipHi,ipLo)][ipCollider].ntemps();
        }
        // Returns pointer to temperature ordinates
        double* temps(long ipHi, long ipLo, long ipCollider)
        {
                return &m_temps[m_a[ilev(ipHi,ipLo)][ipCollider].offset()];
        }
        // Returns pointer to collision data
        double* collstrs(long ipHi, long ipLo, long ipCollider)
        {
                return &m_collstrs[m_a[ilev(ipHi,ipLo)][ipCollider].offset()];
        }
};

void LoadIsotopes ( );


double atmdat_2phot_shapefunction( double EbyE2nu, long ipISO, long nelem );

void atmdat_readin(void);

void atmdat_STOUT_readin( long intNS, const string& chFileName );

void atmdat_CHIANTI_readin( long intNS, const string& chFileName );

void atmdat_LAMDA_readin( long intNS, const string& chFileName );


void atmdat_outer_shell(
  long int iz, 
  long int in,
  long int *imax,
  long int *ig0, 
  long int *ig1);

void ChargTranEval( void );

double ChargTranSumHeat(void);

/*ChargTranPun save charge transfer rate coefficients */
void ChargTranPun( FILE* ipPnunit , char* chSave );

double CHIANTI_Upsilon(long, long, long, long,double);

double atmdat_dielrec_fe(long int ion, double t);

void atmdat_H_phot_cs(void);

void atmdat_3body(void);

double atmdat_HS_caseB( 
        long int iHi, long int iLo, long int iZ, double TempIn,
        double DenIn,
        char chCase 
        );

// arrays for Hummer & Storey 98 He1 cross sections and energies
extern multi_arr<double,4> HS_He1_Xsectn;
extern multi_arr<double,4> HS_He1_Energy;

// arrays for TOPbase Helike cross sections and energies
extern multi_arr<vector<double>,4> OP_Helike_Xsectn;
extern multi_arr<vector<double>,4> OP_Helike_Energy;
extern multi_arr<long,4> OP_Helike_NumPts;

/* these are the vectors that store the original Hummer and Storey case B
 * line data for H and He - the declaration for the interpolator follows */
#define NHSDIM 15 
#define NLINEHS 300  
#define HS_NZ 8 
#define NHCSTE  8 
#define NUM_HS98_DATA_POINTS 811

struct t_atmdat : public module {
        const char *chName() const
        {
                return "atmdat";
        }
        void zero();
        void comment(t_warnings&) {}
        enum {NCX=2};

        double EIonPot[LIMELM][LIMELM];

        inline double getIonPot( long int nelem, long int ion )
        {
                return EIonPot[nelem][ion];
        }

        double CharExcIonOf[NCX][LIMELM][LIMELM+1], //(cm3 s-1)
          CharExcRecTo[NCX][LIMELM][LIMELM+1],          //(cm3 s-1)
          HCharHeatMax, 
          HCharCoolMax, 
          HCharHeatOn;

        /* rate coefficient (cm3 s-1) for N+(3P) + H+ -> N(2D) + H+ charge transfer*/
        double HCharExcRecTo_N0_2D;

        double CharExcIonTotal[NCX],
                CharExcRecTotal[NCX];

        double HIonFrac;

        double HIonFracMax;

        double HCTAlex;

        bool lgCTOn;

        double Density[2][HS_NZ][NHSDIM], 
                ElecTemp[2][HS_NZ][NHSDIM],
                Emiss[2][HS_NZ][NHSDIM][NHSDIM][NLINEHS];

        long int nDensity[2][HS_NZ] , ntemp[2][HS_NZ];
        long int ncut[2][HS_NZ];

        bool lgHCaseBOK[2][HS_NZ];

        long int nsbig;

        bool lgCollIonOn;

        realnum WaveLengthCaseB[8][25][24];

        vector<realnum> CaseBWlHeI;

        const long nDefaultPhotoLevelsFe;
        const long nDefaultPhotoLevels;
        const long nDefaultCollLevelsFe;
        const long nDefaultCollLevels;
        const long nDefaultMolLevels;
        // Set the constant member variables to the default values.

        bool lgChiantiOn;
        bool lgChiantiLvl2Hybrid;
        bool lgChiantiPrint;
        bool lgChiantiExp;
        char chCloudyChiantiFile[FILENAME_PATH_LENGTH];
        long nChiantiMaxLevelsFe;
        long nChiantiMaxLevels;
        bool lgChiantiLevelsSet;

        bool lgLamdaOn;
        bool lgLamdaPrint;
        char chLamdaFile[FILENAME_PATH_LENGTH];
        long nLamdaMaxLevels;
        bool lgLamdaLevelsSet;

        bool lgStoutOn;
        bool lgStoutLvl2Hybrid;
        bool lgStoutPrint;
        char chStoutFile[FILENAME_PATH_LENGTH];
        long nStoutMaxLevelsFe;
        long nStoutMaxLevels;
        bool lgStoutLevelsSet;
        bool lgCalpgmOn;
        bool lgGbarOn;

        double collstrDefault;

        bool lgInnerShellLine_on;
        bool lgInnerShell_Kisielius;
        bool lgUTAprint;

        static const double aulThreshold;

        /* type and enum for determining what collisional ionization rate coefficient data to use.
        * Dima: Cloudy original data from Voronov97
        * Hybrid: Dima version scaled by the ratio of Dere07 to Dima */
        typedef enum { DIMA, HYBRID } CollIonRC;
        CollIonRC CIRCData;

        string chVersion;

        char chdBaseSources[LIMELM][LIMELM+1][10];
        bool lgdBaseSourceExists[LIMELM][LIMELM+1];

        multi_arr<long,2> ipSpecIon;

        t_atmdat() : nDefaultPhotoLevelsFe(25), nDefaultPhotoLevels(15), nDefaultCollLevelsFe(100),
                        nDefaultCollLevels(50), nDefaultMolLevels(70)
        {
                //set all sources to blank and false.
                //Sources will be added from masterlist files in species.cpp
                for( int nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
                {
                        for( int ion=0; ion < nelem+2; ++ion )
                        {
                                strcpy(chdBaseSources[nelem][ion]," ");
                                lgdBaseSourceExists[nelem][ion] = false;
                        }
                }

                ipSpecIon.alloc(LIMELM,LIMELM);
                ipSpecIon = -1;
        }
};

extern t_atmdat atmdat;

class Funct 
{
public:
        virtual void operator()( long&, long&, const char*, long&) = 0;
        virtual ~Funct() = 0;
};
 
inline Funct::~Funct() {}

typedef Funct* FunctPtr;

class FunctLAMDA : public Funct 
{
public:
        explicit FunctLAMDA(void) { }
        virtual void operator()( long& ipHi, long& ipLo, const char* chLine, long& i) 
        { 
                DEBUG_ENTRY( "FunctLAMDA()" );
                bool lgEOL;
                long index = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
                if( index <= 0 )
                {
                        fprintf( ioQQQ, " PROBLEM with LAMDA readin.\n");
                        cdEXIT( EXIT_FAILURE );
                }
                ipHi = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL ) - 1;
                ipLo = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL ) - 1;
                return;
        }
private:
};

#include "h2_priv.h"

class FunctDiatoms : public Funct 
{
public:
        explicit FunctDiatoms(const diatomics& diatom) : diatom_(diatom) { }
        virtual void operator()( long& ipHi, long& ipLo, const char* chLine, long& i) 
        { 
                diatom_.GetIndices( ipHi, ipLo, chLine, i );
        }
private:
        const diatomics& diatom_;
};

void ReadCollisionRateTable( CollRateCoeffArray& coll_rate_table, FILE* io, FunctPtr GetIndices, long nMolLevs, long nTemps = -1, long nTrans = -1 );

double InterpCollRate( const CollRateCoeffArray& rate_table, const long& ipHi, const long& ipLo, const double& ftemp);


/* Form a comment for a transition obtained from an external database.
 * The comment is the configuration of the lower and upper level */
inline string db_comment_tran_levels(const string& clo = "", const string& chi = "")
{
        stringstream Comment;
        if( clo.length() > 0 && chi.length() > 0 )
                Comment << clo << " -- " << chi;
        return Comment.str();
}

#endif /* ATMDAT_H_ */