Changeset 2346 for branches/newmole

Show
Ignore:
Timestamp:
08/19/08 10:14:48 (5 months ago)
Author:
rjrw
Message:

Merged from trunk r2137:2345

Location:
branches/newmole/source
Files:
119 modified

Legend:

Unmodified
Added
Removed

  • branches/newmole/source/assert_results.cpp

    r2139 r2346  
    544544                        strcpy( chAssertLineLabel[nAsserts] , "+Col"); 
    545545                } 
     546                else 
     547                { 
     548                        /*no option found */ 
     549                        fprintf( ioQQQ,  
     550                                " I could not identify an iso-sequence on this Case A/B command.\n"); 
     551                        fprintf( ioQQQ, " Sorry.\n" ); 
     552                        cdEXIT(EXIT_FAILURE); 
     553                } 
    546554        } 
    547555        else if( nMatch("DEPA",input.chCARDCAPS ) ) 
    … …  
    19421950                        strcat( chElemLabelCaseB , chNumb ); 
    19431951 
    1944                         int iCase = 1; 
     1952                        /* sets lowest quantum number index */ 
     1953                        int iCase; 
    19451954                        if( strcmp( "Ca A", chAssertLineLabel[i] ) == 0 ) 
    19461955                                iCase = 0; 
    … …  
    19571966                        { 
    19581967                                fprintf(ioASSERT,"                 Species  nHi nLo  Wl  Computed  Asserted       error\n"); 
    1959                                 /* limit of 10 is because that is all we printed and saved in prt_lines_hydro */ 
     1968                                /* limit of 10 is because that is all we printed and saved in prt_lines_hydro  
     1969                                 * wavelengths will come out of atmdat.WaveLengthCaseB - first index is 
     1970                                 * nelem on C scale, H is 0, second two are configurations of line on 
     1971                                 * physics scale, so Ha is 3-2 */ 
    19601972                                for( long int ipLo=1+iCase; ipLo< MIN2(10,nHighestPrinted-1); ++ipLo ) 
    19611973                                { 
    … …  
    25892601 
    25902602                /* write start of title and version number of code */ 
    2591                 fprintf( ioASSERT, "=============Results of asserts: Cloudy %s ", version::Inst().chVersion ); 
     2603                fprintf( ioASSERT, "=============Results of asserts: Cloudy %s ", t_version::Inst().chVersion ); 
    25922604 
    25932605                /* usually print date and time info - do not if "no times" command entered,  
    … …  
    27842796                if( prt.lgPrintTime ) 
    27852797                { 
    2786                         fprintf( ioQQQ, " %s\n\n", version::Inst().chInfo ); 
     2798                        fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo ); 
    27872799                } 
    27882800        } 

  • branches/newmole/source/atmdat.h

    r2023 r2346  
    203203        bool lgCollIonOn; 
    204204 
    205         /** wavelengths of Hummer & Storey case B lines for H - O */ 
     205        /** wavelengths of Hummer & Storey case B lines for H - O  
     206         * first dimension is atomic number of C scale, H is 0 
     207         * next two are upper and lower configurations on physics  
     208         * scale - Lya is 2-1, Lyb is 3-1, Ha is 3-2, etc */ 
    206209        realnum WaveLengthCaseB[8][25][24]; 
    207210 

  • branches/newmole/source/atmdat_chianti.cpp

    r2079 r2346  
    103103        { 
    104104                lgProtonData = true; 
     105                fclose( atmolProColDATA ); 
     106                atmolProColDATA = NULL; 
    105107        } 
    106108        else 

  • branches/newmole/source/atmdat_lines_setup.cpp

    r2023 r2346  
    11401140                /* check that energy is positive*/ 
    11411141                ASSERT( HFLines[i].EnergyWN > 0 ); 
     1142                ASSERT( HFLines[i].Emis->Aul>0 ); 
     1143                ASSERT( HFLines[i].Emis->damp>0 ); 
    11421144 
    11431145                /* HFLines[i].Emis->gf this is gf if positive, A if negative */ 

  • branches/newmole/source/atmdat_readin.cpp

    r2139 r2346  
    240240        struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 
    241241 
    242         struc.AccelTot = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 
     242        struc.AccelTotalOutward = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 
    243243 
    244244        struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 

  • branches/newmole/source/atom_feii.cpp

    r2139 r2346  
    17621762 
    17631763        /* get the normalization line */ 
    1764         if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. ) 
    1765         { 
    1766                 renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]; 
     1764        if( LineSv[LineSave.ipNormWavL].SumLine[LineSave.lgLineEmergent] > 0. ) 
     1765        { 
     1766                renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].SumLine[LineSave.lgLineEmergent]; 
    17671767        } 
    17681768        else 
    … …  
    30653065                                        Fe2LevN[ipHi][ipLo].Emis->PopOpc > smallfloat ) 
    30663066                                { 
    3067                                         RadPres1 =  RadiationPressure( &Fe2LevN[ipHi][ipLo] ); 
     3067                                        RadPres1 =  PressureRadiationLine( &Fe2LevN[ipHi][ipLo], 1.0 ); 
    30683068 
    30693069#                                       ifdef DEBUGFE 

  • branches/newmole/source/atom_hyperfine.cpp

    r2139 r2346  
    1010/*H21cm_proton - evaluate proton spin changing H atom collision rate, */ 
    1111#include "cddefines.h" 
     12#include "conv.h" 
    1213#include "lines_service.h" 
    1314#include "phycon.h" 
    … …  
    4546                rate12 , rate21 , pump12 , pump21 , coll12 , coll21, 
    4647                texc , occnu_lya_23 , occnu_lya_13,occnu_lya_24,occnu_lya_14, texc1, texc2; 
    47         a31 = 2.08e8;   /* einstein co-efficient for transition 1p1/2 to 0s1/2 */ 
    48         a32 = 4.16e8;   /* einstein co-efficient for transition 1p1/2 to 1s1/2 */ 
    49         a41 = 4.16e8;   /* einstein co-efficient for transition 1p3/2 to 0s1/2 */ 
    50         a42 = 2.08e8;   /* einstein co-efficient for transition 1p3/2 to 1s1/2 */ 
     48 
     49        PopTot = StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]; 
     50        /* population can be zero in certain tests where H is turned off, 
     51         * also if initial solver does not see any obvious source of ionization 
     52         * also possible to set H0 density to zero with element ionization command, 
     53         * as is done in func_set_ion test case */ 
     54        if( PopTot <0 ) 
     55                TotalInsanity(); 
     56        else if( PopTot == 0 ) 
     57        { 
     58                /*return after zeroing local variables */ 
     59                HFLines[0].Hi->Pop = 0.; 
     60                HFLines[0].Lo->Pop = 0.; 
     61                HFLines[0].Emis->PopOpc = 0.; 
     62                HFLines[0].Emis->phots = 0.; 
     63                HFLines[0].Emis->xIntensity = 0.; 
     64                HFLines[0].Emis->ColOvTot = 0.; 
     65                hyperfine.Tspin21cm = 0.; 
     66                return; 
     67        } 
     68 
     69        a31 = 2.08e8;   /* Einstein co-efficient for transition 1p1/2 to 0s1/2 */ 
     70        a32 = 4.16e8;   /* Einstein co-efficient for transition 1p1/2 to 1s1/2 */ 
     71        a41 = 4.16e8;   /* Einstein co-efficient for transition 1p3/2 to 0s1/2 */ 
     72        a42 = 2.08e8;   /* Einstein co-efficient for transition 1p3/2 to 1s1/2 */ 
    5173        /* These A values are determined from eqn. 17.64 of "The theory of Atomic structure 
    5274         * and Spectra" by R. D. Cowan  
    5375         * A hyperfine level has degeneracy Gf=(2F + 1) 
    54          * a2p1s = 6.24e8;  einstein co-efficient for transition 2p to 1s */ 
    55         a21 = 2.85e-15; /* einstein co-efficient for transition 1s1/2 to 0s1/2 */ 
     76         * a2p1s = 6.24e8;  Einstein co-efficient for transition 2p to 1s */ 
     77        a21 = 2.85e-15; /* Einstein co-efficient for transition 1s1/2 to 0s1/2 */ 
    5678 
    5779        /* above is spontaneous rate - the net rate is this times escape and destruction 
    5880         * probabilities */ 
    5981        a21 *= (HFLines[0].Emis->Pdest + HFLines[0].Emis->Pesc + HFLines[0].Emis->Pelec_esc); 
    60  
    61         /* >>chng 04 dec 01, add hyperfine.lgLya_pump_21cm, option to turn off Lya pump 
    62          * of 21 cm, with no 21cm lya pump */ 
     82        ASSERT( a21>0. ); 
     83 
     84        /* hyperfine.lgLya_pump_21cm is option to turn off Lya pump 
     85         * of 21 cm, with no 21cm lya pump command - note that this 
     86         * can be negative if Lya mases - can occur during search phase */ 
    6387        occnu_lya = OccupationNumberLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) * 
    6488                hyperfine.lgLya_pump_21cm; 
    65  
    66         /* >>chng 05 apr 20, GS, the lya occupation number for the hyperfine levels 0S1/2 and 1S1/2 are different*/ 
     89        if( occnu_lya<0 ) 
     90        { 
     91                static bool lgCommentDone = false; 
     92                /* Lya is masing - could be due to very bad solution in search phase */ 
     93                if( !conv.lgSearch && !lgCommentDone ) 
     94                { 
     95                        fprintf(ioQQQ, 
     96                        "NOTE Lya masing will invert 21 cm, occupation number set zero\n"); 
     97                        lgCommentDone = true; 
     98                } 
     99                occnu_lya = 0.; 
     100        } 
     101 
     102        /* Lya occupation number for the hyperfine levels 0S1/2 and 1S1/2 are different*/ 
    67103        texc = TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ); 
    68         /* >>chng 05 apr 21, GS, Energy difference between 2p1/2 and 2p3/2 is taken from NSRDS */ 
     104        /* Energy difference between 2p1/2 and 2p3/2 taken from NSRDS */ 
    69105        if( texc > 0. ) 
    70106        { 
    71                 /* convert to boltz factor, which will applied to occupation number of hiher energy transition */ 
     107                /* convert to Boltzmann factor, which will applied to occupation  
     108                 * number of higher energy transition */ 
    72109                texc1 = sexp(0.068/texc); 
    73110                texc2 = sexp(((82259.272-82258.907)*T1CM)/texc); 
    … …  
    95132         * were multiplied by collider density when evaluated in CoolEvaluate */ 
    96133        /* ContBoltz is Boltzmann factor for wavelength of line */ 
     134        ASSERT( HFLines[0].Coll.col_str>0. ); 
    97135        coll12 = HFLines[0].Coll.col_str*dense.cdsqte/HFLines[0].Lo->g*rfield.ContBoltz[HFLines[0].ipCont-1]; 
    98136        coll21 = HFLines[0].Coll.col_str*dense.cdsqte/HFLines[0].Hi->g; 
    … …  
    128166        x = rate12 / SDIV(a21 + rate21); 
    129167 
    130         PopTot = StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]; 
    131168        /* the Transitions term is the total population of 1s */ 
    132169        HFLines[0].Hi->Pop = (x/(1.+x))* PopTot; 
    133170        HFLines[0].Lo->Pop = (1./(1.+x))* PopTot; 
     171        ASSERT( HFLines[0].Hi->Pop >0. ); 
     172        ASSERT( HFLines[0].Lo->Pop >0. ); 
    134173 
    135174        /* the population with correction for stimulated emission */ 
    136         /* >>chng 04 dec 03, rewrite for precision in limit where a21 is very small */ 
    137         /*HFLines[0].Emis->PopOpc = HFLines[0].Hi->Pop*(a21/SDIV(rate12) + (3.*rate21-rate12)/SDIV(3.*rate12) );*/ 
    138         /* >>chng 04 dec 14, GS, error in defin of x above */ 
    139         /* HFLines[0].Emis->PopOpc = HFLines[0].Lo->Pop*(3*(rate21+a21)- rate12)/SDIV(3*(a21+ rate21));*/ 
    140175        HFLines[0].Emis->PopOpc = HFLines[0].Lo->Pop*((3*rate21- rate12) + 3*a21)/SDIV(3*(a21+ rate21)); 
    141176 
    142         /*fprintf(ioQQQ,"DEBUG 21cm\t%e", HFLines[0].Emis->PopOpc ); 
    143         HFLines[0].Emis->PopOpc = HFLines[0].Lo->Pop - HFLines[0].Hi->Pop*(1./3.);  
    144         fprintf(ioQQQ,"\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",  
    145                 HFLines[0].Emis->PopOpc , a21,rate12,rate12,HFLines[0].Lo->Pop,HFLines[0].Hi->Pop,PopTot,x);*/ 
     177        /* number of escaping line photons, used elsewhere for outward beam */ 
     178        HFLines[0].Emis->phots = a21*HFLines[0].Hi->Pop; 
     179        ASSERT( HFLines[0].Emis->phots >= 0. ); 
     180        /* intensity of line */ 
     181        HFLines[0].Emis->xIntensity = HFLines[0].Emis->phots*HFLines[0].EnergyErg; 
     182 
     183        /* ratio of collisional to total (collisional + pumped) excitation */ 
     184        HFLines[0].Emis->ColOvTot = coll12 / rate12; 
    146185 
    147186        /* finally save the spin temperature */ 
    … …  
    487526                        wavelength = (realnum)help[2]; 
    488527                        HFLines[j].Emis->Aul = (realnum)Aul; 
     528                        HFLines[j].Emis->damp = 1e-20f; 
    489529                        HFLines[j].Hi->g = (realnum) (2*(spin + .5) + 1); 
    490530                        HFLines[j].Lo->g = (realnum) (2*(spin - .5) + 1); 
    … …  
    501541 
    502542        } 
     543        ASSERT( j==nHFLines ); 
    503544 
    504545        /* closing the file */ 

  • branches/newmole/source/atom_level2.cpp

    r2139 r2346  
    2222          nelem; 
    2323        double AbunxIon,  
    24           Enr12,  
    25           Enr21,  
    2624          a21,  
    2725          boltz,  
    … …  
    177175         * so that negative cooling does not occur */ 
    178176 
    179         Enr12 = plower*col12*t->EnergyErg; 
    180         Enr21 = pfs2*col21*t->EnergyErg; 
     177        //double Enr12 = plower*col12*t->EnergyErg; 
     178        //double Enr21 = pfs2*col21*t->EnergyErg; 
    181179 
    182180        /* energy exchange due to this level 

  • branches/newmole/source/cddefines.h

    r2139 r2346  
    7373template<bool> struct StaticAssertFailed; 
    7474template<> struct StaticAssertFailed<true> {}; 
    75 #define STATIC_ASSERT(x) (StaticAssertFailed< (x) == true >()) 
     75#define STATIC_ASSERT(x) ((void)StaticAssertFailed< (x) == true >()) 
    7676 
    7777typedef float sys_float; 

  • branches/newmole/source/cddrive.cpp

    r1918 r2346  
    359359void cdVersion(char chString[]) 
    360360{ 
    361         strcpy( chString , version::Inst().chVersion ); 
     361        strcpy( chString , t_version::Inst().chVersion ); 
    362362        return; 
    363363} 
    … …  
    374374void cdDate(char chString[]) 
    375375{ 
    376         strcpy( chString , version::Inst().chDate ); 
     376        strcpy( chString , t_version::Inst().chDate ); 
    377377        return; 
    378378} 
    … …  
    993993        /* gas pressure */ 
    994994        double GasPressure[],                            
    995         /* radiation pressure */ 
     995        /* line radiation pressure */ 
    996996        double RadiationPressure[]) 
    997997{ 
    … …  
    10181018                double *PresTotal,  /* total pressure, all forms, for the last computed zone*/ 
    10191019                double *PresGas,    /* gas pressure */ 
    1020                 double *PresRad)    /* radiation pressure */ 
     1020                double *PresRad)    /* line radiation pressure */ 
    10211021{ 
    10221022 
    … …  
    10251025        *PresGas = pressure.PresGasCurr; 
    10261026        *PresRad = pressure.pres_radiation_lines_curr; 
    1027         *PresTotal = pressure.PresGasCurr + pressure.pres_radiation_lines_curr; 
     1027        *PresTotal = pressure.PresTotlCurr; 
    10281028        return; 
    10291029} 
    … …  
    13771377 
    13781378                                /* does the normalization line have a positive intensity*/ 
    1379                                 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. ) 
     1379                                if( LineSv[LineSave.ipNormWavL].SumLine[LineSave.lgLineEmergent] > 0. ) 
    13801380                                { 
    1381                                         *relint = LineSv[ipobs].sumlin[LineSave.lgLineEmergent]/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]* 
     1381                                        *relint = LineSv[ipobs].SumLine[LineSave.lgLineEmergent]/LineSv[LineSave.ipNormWavL].SumLine[LineSave.lgLineEmergent]* 
    13821382                                          LineSave.ScaleNormLine; 
    13831383                                } 
    … …  
    13881388 
    13891389                                /* return log of current line intensity if it is positive */ 
    1390                                 if( LineSv[ipobs].sumlin[LineSave.lgLineEmergent] > 0. ) 
     1390                                if( LineSv[ipobs].SumLine[LineSave.lgLineEmergent] > 0. ) 
    13911391                                { 
    1392                                         *absint = log10(LineSv[ipobs].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten; 
     1392                                        *absint = log10(LineSv[ipobs].SumLine[LineSave.lgLineEmergent]) + radius.Conv2PrtInten; 
    13931393                                } 
    13941394                                else 
    … …  
    14521452 
    14531453        /* does the normalization line have a positive intensity*/ 
    1454         if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. ) 
    1455         { 
    1456                 *relint = LineSv[ipLine].sumlin[LineSave.lgLineEmergent]/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]* 
     1454        if( LineSv[LineSave.ipNormWavL].SumLine[LineSave.lgLineEmergent] > 0. ) 
     1455        { 
     1456                *relint = LineSv[ipLine].SumLine[LineSave.lgLineEmergent]/LineSv[LineSave.ipNormWavL].SumLine[LineSave.lgLineEmergent]* 
    14571457                  LineSave.ScaleNormLine; 
    14581458        } 
    … …  
    14631463 
    14641464        /* return log of current line intensity if it is positive */ 
    1465         if( LineSv[ipLine].sumlin[LineSave.lgLineEmergent] > 0. ) 
    1466         { 
    1467                 *absint = log10(LineSv[ipLine].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten; 
     1465        if( LineSv[ipLine].SumLine[LineSave.lgLineEmergent] > 0. ) 
     1466        {