iso.h

Go to the documentation of this file.

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

#include "module.h"
#include "transition.h"
#include "container_classes.h"
 
class two_photon;
class freeBound;
 
extern long int max_num_levels;

#define KILL_BELOW_PLASMA(E_)           ( (rfield.lgPlasNu && ((E_)<rfield.plsfrq) ) ? 0.:1. )

#define N_(A_)  (iso_sp[ipISO][nelem].st[A_].n())
#define L_(A_)  (iso_sp[ipISO][nelem].st[A_].l())
#define S_(A_)  (iso_sp[ipISO][nelem].st[A_].S())
#define J_(A_)  (iso_sp[ipISO][nelem].st[A_].j())

 
const int ipH1s = 0;
const int ipH2s = 1;
const int ipH2p = 2;
const int ipH3s = 3;
const int ipH3p = 4;
const int ipH3d = 5;
const int ipH4s = 6;
const int ipH4p = 7;
const int ipH4d = 8;
const int ipH4f = 9;

 
/* level 1 */
const int ipHe1s1S = 0;
 
/* level 2 */
const int ipHe2s3S = 1;
const int ipHe2s1S = 2;
const int ipHe2p3P0 = 3;
const int ipHe2p3P1 = 4;
const int ipHe2p3P2 = 5;
const int ipHe2p1P = 6;
 
/* level 3 */
const int ipHe3s3S = 7;
const int ipHe3s1S = 8;
const int ipHe3p3P = 9;
const int ipHe3d3D = 10;
const int ipHe3d1D = 11;
const int ipHe3p1P = 12;

const int ipH_LIKE = 0;
const int ipHE_LIKE = 1;
const int ipLI_LIKE = 2;
const int ipBE_LIKE = 3;
const int ipB_LIKE = 4;
const int ipC_LIKE = 5;
const int ipN_LIKE = 6;
const int ipO_LIKE = 7;
const int ipF_LIKE = 8;
const int ipNE_LIKE = 9;
const int ipNA_LIKE = 10;
const int ipMG_LIKE = 11;
const int ipAL_LIKE = 12;
const int ipSI_LIKE = 13;
const int ipP_LIKE = 14;
const int ipS_LIKE = 15;
const int ipCL_LIKE = 16;
const int ipAR_LIKE = 17;
 
enum {
        ipSINGLET = 1, ipDOUBLET = 2, ipTRIPLET = 3, 
        ipMULTIPLET_END, ipMULTIPLET_BEGIN=ipSINGLET
};
 
const int IPRAD = 0;
const int IPCOLLIS = 1;
/*const int IPENERGY = 2;*/
 
/* following two macros used to define recombination coef arrays */
/* Max n desired in RRCoef file. */
/* Hydrogen and helium atoms will have precompiled recombination coefficients up to these maximum n. */
const int RREC_MAXN = 40;

inline int LIKE_RREC_MAXN(int nelem) { return ( nelem == ipHELIUM ) ? 40 : 20; }
 
const int N_ISO_TE_RECOMB = 41;

const int SumUpToThisN = 1000;
const int RECOMBMAGIC = 130216;
const int ENERGIESMAGIC = 20190102;
 
typedef uint64 QNPack;
 
inline QNPack QN2ind(long n, long l, long s, long g = -1)
{
        // Index routine for the iso sequences
        // Quantum physics guarantees 0 <= l < n and 0 <= j <= n, while s = 1, 2, or 3
        // In many cases g = 2*j+1 will be -1 to state that the level is not j-resolved
 
        // first check arguments
        ASSERT( n > 0 );
        l = max(l, -1);
        s = max(s, -1);
        g = max(g, -1);
 
        // make sure g value is resolved where possible, this avoids ambiguity in matches
        if( (l == 0 || s == 1) && g < 0 )
                g = max(2*l+1, s);
 
        // make sure that n, l, s, g can be extracted unambiguously from QNPack
        ASSERT( MAX4(n,l,s,g) < 0xffffL );
 
        uint64 nn = uint64(n&0xffff);
        uint64 ll = uint64(l&0xffff);
        uint64 ss = uint64(s&0xffff);
        uint64 gg = uint64(g&0xffff);
        return (nn<<48) | (ll<<32) | (ss<<16) | gg;
}
 
struct QNPair
{
        QNPack hi;
        QNPack lo;
        bool operator< (const QNPair& q2) const
        {
                return ( hi < q2.hi ) || (hi == q2.hi && lo < q2.lo);
        }
        QNPair(int nhi, int lhi, int shi, int ghi, int nlo, int llo, int slo, int glo)
        {
                hi = QN2ind(nhi, lhi, shi, ghi);
                lo = QN2ind(nlo, llo, slo, glo);
        }
};
 
inline QNPair QN2ind(int nhi, int lhi, int shi, int ghi, int nlo, int llo, int slo, int glo)
{
        return QNPair(nhi, lhi, shi, ghi, nlo, llo, slo, glo);
}
 
int getL(char l);

void iso_cascade( long ipISO, long nelem );

void iso_charge_transfer_update( long nelem );

void iso_collide( long ipISO, long nelem );

void iso_collisional_ionization( long ipISO, long nelem );

void iso_continuum_lower( long ipISO , long nelem );

void iso_cool( long ipISO , long nelem );

void iso_setRedisFun (long ipISO, long nelem, long ipLo, long ipHi);

void iso_setOpacity (long ipISO, long nelem, long ipLo, long ipHi);

void iso_create( void );

double iso_cross_section( double ERyd , double EthRyd, long n, long l, long S, long globalZ, long globalISO );

void iso_departure_coefficients( long ipISO, long nelem );

double iso_dielec_recomb_rate( long ipISO, long nelem, long ipLo );

void iso_error_generation( long ipISO, long nelem );

long iso_get_total_num_levels( long ipISO, long nmaxResolved, long numCollapsed );

void IonHydro( );

void iso_ionize_recombine( long ipISO , long nelem );

void iso_level( const long ipISO, const long nelem, double& renorm );

void iso_photo( long ipISO , long nelem );

void iso_prt_pops( long ipISO, long nelem, bool lgPrtDeparCoef );

void iso_put_error(long ipISO,long nelem,long ipHi,long ipLo,long whichData,realnum errorOpt,realnum errorPess);

void iso_put_error(long ipISO,long nelem,QNPack inHi,QNPack inLo,long whichData,realnum errorOpt,realnum errorPess);

void iso_radiative_recomb( long ipISO, long nelem );

void iso_radiative_recomb_effective( long ipISO, long nelem );

double iso_recomb_check( long ipISO, long nelem, long level, double temperature );

void iso_recomb_auxiliary_free();

void iso_recomb_alloc();

void iso_recomb_setup( long ipISO );

double iso_RRCoef_Te( long ipISO, long nelem, double temp, long n );

void iso_satellite_update( long nelem );
 
/* calculate radiative lifetime of an individual iso state 
\param ipISO
\param nelem
\param n
\param l
*/
double iso_state_lifetime( long ipISO, long nelem, long n, long l );

void iso_solve( long ipISO, long nelem, double &maxerr );

void iso_suprathermal( long ipISO, long nelem );

void iso_update_num_levels( long ipISO, long nelem );

void iso_update_rates( void );

long int iso_Max_Emitting_Level(long nelem, long ipISO, bool lgPrnIsoCollapsed);
 
void iso_init_energies();
 
double hydro_energy(long nelem, long n, long l, long s, long g);
 
void iso_collapsed_update( void );
 
void iso_set_ion_rates( long ipISO, long nelem);
 
void iso_init();
 
class t_isoCTRL : public module
{
public:
        void zero();
        void comment(t_warnings&) {}
 
        const char *chName() const
        {
                return "iso_ctrl";
        }
 
        bool lgPrintNumberOfLevels;
 
        const char *chISO[NISO];

        long int nLyman[NISO],
                nLyman_max[NISO],
                nLyman_alloc[NISO];

        long int nLymanHLike[LIMELM];

        double Resolution;

        bool lgColl_l_mixing[NISO];

        bool lgColl_excite[NISO];

        bool lgColl_ionize[NISO];
 
        bool lgLTE_levels[NISO];

        bool lgCollStrenThermAver;

        bool lgInd2nu_On;
 
        /* option to disable continuum lowering due to stark broadening, particle packing, etc. */
        bool lgContinuumLoweringEnabled[NISO];

        realnum stat_ion[NISO];

        bool lgDielRecom[NISO];

        bool lgNoMaser[NISO][LIMELM];

        realnum SmallA;

        int ipLyaRedist[NISO] , ipResoRedist[NISO] , ipSubRedist[NISO];

        int nLyaLevel[NISO];

        bool lgCompileRecomb[NISO];

        bool lgNoRecombInterp[NISO];

        bool lgCS_Vriens[NISO] ,
                lgCS_Lebedev[NISO],
                lgCS_Fujim[NISO],
                lgCS_vrgm[NISO],
                lgCS_None[NISO] ,
                lgCS_Seaton[NISO],
                lgCS_B72[NISO],
                lgCS_PSdeg[NISO],
                lgCS_Vrinceanu[NISO],
                lgCS_PS64[NISO],
                lgCS_PSClassic[NISO],
                lgCS_PSM20[NISO],
                lgCS_VOS12[NISO],
                lgCS_VOS12QM[NISO],
                lgCS_therm_ave[NISO],
                lgCS_VOS_thermal[NISO];
        int nCS_new[NISO];//vals are 0, 1, and 2

        bool lgCritDensLMix[NISO];

        bool lgRandErrGen[NISO];
 
        bool lgPessimisticErrors;
 
        bool lgKeepFS;
 
        bool lgTopoff[NISO];
 
        double RRC_TeUsed[NISO][LIMELM];
 
        t_isoCTRL()
        {
                chISO[ipH_LIKE] = "H-like ";
                chISO[ipHE_LIKE] = "He-like";
                /* This is 1/10 of Athena's resolution in eV. */
                Resolution = 0.25;
        }
};
 
extern t_isoCTRL iso_ctrl;
 
class extra_tr
{
public:
        double pestrk;
        double pestrk_up;
 
        /* NB NB NB ---  Error and ErrorFactor need one more slot than all the rest of these! */
        /* and the last dimension can just be hardwired to 3 */

        /* first dimension is upper level,
         * second is lower level,
         * third is for radiative, collisional, or energy errors.
         * MACROS are used for the last dimension: IPRAD, IPCOLLIS, and IPENERGY. */
        realnum Error[3];

        realnum ErrorFactor[3];

        double SigmaCascadeProb;
};
 
class t_iso_sp
{
        map<QNPack, long> QNPack2Index;
 
public:
        TransitionProxy trans( const long ipHi, const long ipLo ) 
        {
                return (*tr)[ ipTrans[ipHi][ipLo] ];
        }
        multi_arr<long,2> ipTrans;
        multi_arr<extra_tr,2> ex;
        multi_arr<double,2> CascadeProb;
        multi_arr<double,2> BranchRatio;
        vector<freeBound> fb;
        qList st;
        qList stJ05; /* stJ05 should only be used for H iso-sequence, otherwise need to change iso_assign_quantum_numbers() routine in is_create.cpp */
        qList stJ15;
        TransitionList* tr;

        long QN2Index(QNPack ind);
        long QN2Index(long n, long l, long s, long g = -1)
        {
                QNPack ind = QN2ind(n, l, s, g);
                return QN2Index(ind);
        }

        map<QNPack, double> Energy;

        double IonPot;

        double energy(long n, long l, long s, long g = -1) const
        {
                QNPack ind = QN2ind(n, l, s, g);
                auto p = Energy.find(ind);
                if( p != Energy.end() )
                        return p->second;
                else
                        return -1.;
        }

        double energy_ioniz(long n, long l, long s, long g = -1) const
        {
                double ERelToground = energy(n, l, s, g);
                if( ERelToground >= 0. && ERelToground < IonPot )
                        return IonPot - ERelToground;
                else
                        return -1.;
        }

        map<QNPair, double> CachedAs;

        double xIonSimple;

        bool lgPrtDepartCoef;

        bool lgPrtNCrit;

        bool lgPrtLevelPops;

        bool lgLevelsLowered;

        bool lgLevelsEverLowered;
 
        /* flag that says we must reevaluate everything about this ion */
        bool lgMustReeval;
 
        /* set true if "element ionization" forces rescaling of pops */
        bool lgPopsRescaled;

        long int nCollapsed_max;
        long int nCollapsed_local;

        long int numLevels_max;

        long int numLevels_local;

        long int numLevels_alloc;

        long int n_HighestResolved_max;
        long int n_HighestResolved_local;

        realnum CaseBCheck;

        double RadRec_caseB;

        double RadRec_effec;

        double RecomCollisFrac;

        bool lgPopLTE_OK;

        double FreeBnd_net_Cool_Rate;

        double coll_ion;

        double cRest_cool;

        double xLineTotCool;

        double dLTot;

        double RadRecCool;

        double cBal_cool;

        double cLyrest_cool;

        double cLya_cool;

        double RecomInducCool_Rate;

        char chTypeAtomUsed[10];

        bool lgErrGenDone;

        double qTot2S;
 
        void Reset();
        vector<two_photon> TwoNu;
 
        vector<double> HighestLevelOpacStack;

        bool lgPrtMatrix;
        bool lgImgMatrix;
};
 
extern t_iso_sp iso_sp[NISO][LIMELM];

void iso_renorm( long nelem, long ipISO, double& renorm );

void iso_multiplet_opacities( void );
 
/* iso_comment_tran_levels - prepare comment string for entry to the line stack.
 * The comment has the form 'H-like, 1 3, 1^2S - 2^2P', where '1 3' are the energy
 * level indices, 1 being ground.
 *
 * \param ipISO         iso-sequence index
 * \param nelem         element index
 * \param ipLo, ipHi    lower and upper level indices
 * \return              comment string
 */
string iso_comment_tran_levels( long ipISO, long nelem, long ipLo, long ipHi );
 
string extraLymanJ_comment_tran_levels( const TransitionProxy &t );
 
inline bool lgIsLymanLine(const TransitionProxy &t)
{
        long ipISO = t.Lo()->nelem() - t.Lo()->IonStg();
        long nelem = t.Lo()->nelem() - 1;
        return (ipISO == ipH_LIKE && t.Lo()->n() == 1 && (t.Hi()->n() > iso_sp[ipISO][nelem].n_HighestResolved_local || t.Hi()->l() == 1));
}
 
inline bool lgIsLymanLineResolved(const TransitionProxy &t, const TransitionProxy &tJ05, const TransitionProxy &tJ15)
{
        double Ediff = tJ15.Hi()->energy().get("eV") - tJ05.Hi()->energy().get("eV");
        long nelem = t.Lo()->nelem() - 1;
        return (lgIsLymanLine(t) && Ediff > iso_ctrl.Resolution && nelem > ipHELIUM);
}
 
inline bool lgIsLymanLineUnresolved(const TransitionProxy &t)
{
        return (lgIsLymanLine(t) && t.Hi()->g() == 6);
}
 
inline bool lgIsM1Line(const TransitionProxy &t)
{
        long ipISO = t.Lo()->nelem() - t.Lo()->IonStg();
        return (ipISO == ipH_LIKE && t.Hi()->n() == 2 && t.Lo()->n() == 1 && t.Hi()->l() == 0);
}
#endif /* ISO_H_ */