Changeset 2139

Show
Ignore:
Timestamp:
06/29/08 09:04:54 (5 months ago)
Author:
rjrw
Message:

Merged from trunk r2124:2137

Location:
branches/newmole/source
Files:
65 modified

Legend:

Unmodified
Added
Removed

  • branches/newmole/source/abundances.cpp

    r1780 r2139  
    272272        /* if stop temp set below default then we are going into cold and possibly  
    273273         * molecular gas - check some parameters in this case */ 
    274         if( called.lgTalk && (StopCalc.tend < phycon.TEMP_STOP_DEFAULT ||  
     274        if( called.lgTalk && (StopCalc.TempLoStopZone < phycon.TEMP_STOP_DEFAULT ||  
    275275                /* thermal.ConstTemp def is zero, set pos when used */ 
    276276                (thermal.ConstTemp > 0. && thermal.ConstTemp < phycon.TEMP_STOP_DEFAULT ) ) ) 
    … …  
    278278 
    279279                /* print warning if temperature set below default but C > O */ 
    280                 if( dense.gas_phase[ipCARBON]/SDIV( dense.gas_phase[ipOXYGEN]) >= 1. ) 
     280                if( dense.lgElmtOn[ipOXYGEN] && dense.gas_phase[ipCARBON]/SDIV( dense.gas_phase[ipOXYGEN]) >= 1. ) 
    281281                { 
    282282                        fprintf( ioQQQ, "\n >>> \n" 

  • branches/newmole/source/assert_results.cpp

    r2023 r2139  
    19431943 
    19441944                        int iCase = 1; 
     1945                        if( strcmp( "Ca A", chAssertLineLabel[i] ) == 0 ) 
     1946                                iCase = 0; 
     1947                        else if( strcmp( "Ca B", chAssertLineLabel[i] ) == 0 ) 
     1948                                iCase = 1; 
     1949                        else if( strcmp( "+Col", chAssertLineLabel[i] ) == 0 ) 
     1950                                iCase = 1; 
     1951                        else 
     1952                                TotalInsanity(); 
     1953 
    19451954                        RelError[i] = 0.; 
    19461955                        long nHighestPrinted = StatesElem[nISOCaseB][nelemCaseB][iso.numPrintLevels[nISOCaseB][nelemCaseB]-1].n; 

  • branches/newmole/source/atmdat_readin.cpp

    r1830 r2139  
    538538 
    539539                        /* this is new in c94.01 and returns nothing (0.) for most lines */ 
    540                         TauLines[i].Coll.cs = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 
     540                        TauLines[i].Coll.col_str = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 
    541541 
    542542                        /* finally increment i */ 
    … …  
    956956 
    957957                        /*ASSERT( UTALines[i].Emis->Aul / sum_spon_auto <= 1. && UTALines[i].Aul / sum_spon_auto >= 0. );*/ 
    958                         UTALines[i].Coll.cs = (realnum)MIN2(0., UTALines[i].Emis->Aul / sum_spon_auto - 1. ); 
     958                        UTALines[i].Coll.col_str = (realnum)MIN2(0., UTALines[i].Emis->Aul / sum_spon_auto - 1. ); 
    959959 
    960960                        /* option to print UTAs */ 
    … …  
    10281028                                double frac_ioniz = Aauto/(Aul + Aauto); 
    10291029                                ASSERT( frac_ioniz >=0. &&  frac_ioniz <=1. ); 
    1030                                 UTALines[i].Coll.cs = -(realnum)frac_ioniz; 
     1030                                UTALines[i].Coll.col_str = -(realnum)frac_ioniz; 
    10311031 
    10321032                                /* save gf scanned from line */ 
    … …  
    11301130                                /* >>chng 03 feb 18, cs to negative of branching ratio */ 
    11311131                                ASSERT( frac_relax >=0. &&  frac_relax <=1. ); 
    1132                                 UTALines[i].Coll.cs = (realnum)(frac_relax-1.); 
     1132                                UTALines[i].Coll.col_str = (realnum)(frac_relax-1.); 
    11331133#                               if 0 
    11341134                                fprintf(ioQQQ,"data set %li line %li %2s%2li\t%.4f\t%.3e\t%.3e\n", 
    … …  
    12141214                        /* >>chng 03 feb 18, cs to negative of br */ 
    12151215                        ASSERT( frac_relax >=0. &&  frac_relax <=1. ); 
    1216                         UTALines[i].Coll.cs = (realnum)(frac_relax-1.); 
     1216                        UTALines[i].Coll.col_str = (realnum)(frac_relax-1.); 
    12171217 
    12181218                        /* finally increment i */ 

  • branches/newmole/source/atom_feii.cpp

    r2023 r2139  
    11431143                          POW2(y + 1.0)) + cg*y); 
    11441144 
    1145                         Fe2LevN[ipHi][ipLo].Coll.cs = 23.861f*1e5f*gb* 
     1145                        Fe2LevN[ipHi][ipLo].Coll.col_str = 23.861f*1e5f*gb* 
    11461146                          Fe2LevN[ipHi][ipLo].Emis->Aul*Fe2LevN[ipHi][ipLo].Hi->g/ 
    11471147                          POW3(Fe2LevN[ipHi][ipLo].EnergyWN); 
    11481148 
    11491149                        /* confirm that collision strength is positive */ 
    1150                         ASSERT( Fe2LevN[ipHi][ipLo].Coll.cs > 0.); 
     1150                        ASSERT( Fe2LevN[ipHi][ipLo].Coll.col_str > 0.); 
    11511151 
    11521152                        /* g-bar cs becomes unphysically small for forbidden transitions - 
    11531153                         * this sets a lower limit to the final cs - CS2SMALL is defined above */ 
    1154                         Fe2LevN[ipHi][ipLo].Coll.cs = MAX2( CS2SMALL, Fe2LevN[ipHi][ipLo].Coll.cs); 
     1154                        Fe2LevN[ipHi][ipLo].Coll.col_str = MAX2( CS2SMALL, Fe2LevN[ipHi][ipLo].Coll.col_str); 
    11551155                        /* this was the logic used in the old fortran version, 
    11561156                         * and reproduces the results in Baldwin et al '96 
    … …  
    12351235                                /*fprintf(ioQQQ,"DEBUG\t%.2e", Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs);*/ 
    12361236                                /* do linear interpolation */ 
    1237                                 Fe2LevN[ipHiFe2-1][ipLoFe2-1].Coll.cs =  
     1237                                Fe2LevN[ipHiFe2-1][ipLoFe2-1].Coll.col_str =  
    12381238                                        sPradDat[ipHi][ipLo][ipTemp] * FracLowTe +  
    12391239                                        sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe; 
    … …  
    12411241 
    12421242                                /* confirm that this is positive */ 
    1243                                 ASSERT( Fe2LevN[ipHiFe2-1][ipLoFe2-1].Coll.cs > 0. ); 
     1243                                ASSERT( Fe2LevN[ipHiFe2-1][ipLoFe2-1].Coll.col_str > 0. ); 
    12441244                        } 
    12451245                } 
    … …  
    12941294                        { 
    12951295                                fprintf(ioQQQ,"%e %e\n",  
    1296                                         Fe2LevN[51][ipLo].Coll.cs,Fe2LevN[52][ipLo].Coll.cs); 
     1296                                        Fe2LevN[51][ipLo].Coll.col_str,Fe2LevN[52][ipLo].Coll.col_str); 
    12971297                        } 
    12981298                        cdEXIT(EXIT_FAILURE); 
    … …  
    13081308                         * note that it needs eden to become rate 
    13091309                         * units cm3 s-1 */ 
    1310                         Fe2Coll[ipHi][ipLo] = (realnum)(COLL_CONST/phycon.sqrte*Fe2LevN[ipHi][ipLo].Coll.cs/ 
     1310                        Fe2Coll[ipHi][ipLo] = (realnum)(COLL_CONST/phycon.sqrte*Fe2LevN[ipHi][ipLo].Coll.col_str/ 
    13111311                          Fe2LevN[ipHi][ipLo].Hi->g); 
    13121312                        /*fprintf(ioQQQ,"DEBUG fe2coll %li %li %.2e", ipHi , ipLo , Fe2Coll[ipHi][ipLo] );*/ 

  • branches/newmole/source/atom_hyperfine.cpp

    r1830 r2139  
    9595         * were multiplied by collider density when evaluated in CoolEvaluate */ 
    9696        /* ContBoltz is Boltzmann factor for wavelength of line */ 
    97         coll12 = HFLines[0].Coll.cs*dense.cdsqte/HFLines[0].Lo->g*rfield.ContBoltz[HFLines[0].ipCont-1]; 
    98         coll21 = HFLines[0].Coll.cs*dense.cdsqte/HFLines[0].Hi->g; 
     97        coll12 = HFLines[0].Coll.col_str*dense.cdsqte/HFLines[0].Lo->g*rfield.ContBoltz[HFLines[0].ipCont-1]; 
     98        coll21 = HFLines[0].Coll.col_str*dense.cdsqte/HFLines[0].Hi->g; 
    9999 
    100100        /* set up rate (s-1) equations 

  • branches/newmole/source/atom_level2.cpp

    r1830 r2139  
    6161        /* collision strength for this transition, omega is zero for hydrogenic 
    6262         * species which are done in special hydro routines */ 
    63         omega = t->Coll.cs; 
     63        omega = t->Coll.col_str; 
    6464 
    6565        /* ROUGH check whether upper level has significant population,*/ 

  • branches/newmole/source/atom_level3.cpp

    r1830 r2139  
    292292 
    293293        /* collision strengths */ 
    294         o10 = t10->Coll.cs; 
    295         o21 = t21->Coll.cs; 
    296         o20 = t20->Coll.cs; 
     294        o10 = t10->Coll.col_str; 
     295        o21 = t21->Coll.col_str; 
     296        o20 = t20->Coll.col_str; 
    297297 
    298298        /* begin sanity checks, check statistic weights,  

  • branches/newmole/source/atom_oi.cpp

    r1918 r2139  
    318318        cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01; 
    319319        PutCS(cs,&TauLines[ipTO1025]); 
    320         c21 = dense.cdsqte*TauLines[ipT1304].Coll.cs/g[1]; 
    321         c51 = dense.cdsqte*TauLines[ipTO1025].Coll.cs/g[4]; 
     320        c21 = dense.cdsqte*TauLines[ipT1304].Coll.col_str/g[1]; 
     321        c51 = dense.cdsqte*TauLines[ipTO1025].Coll.col_str/g[4]; 
    322322 
    323323        /* all following are g-bar approx, g-bar = 0.2 */ 
    324324        c31 = dense.cdsqte*1.0/g[2]; 
    325325        PutCS(0.27,&TauLines[ipT1039]); 
    326         c41 = dense.cdsqte*TauLines[ipT1039].Coll.cs/g[3]; 
     326        c41 = dense.cdsqte*TauLines[ipT1039].Coll.col_str/g[3]; 
    327327        c61 = dense.cdsqte*1./g[5]; 
    328328 

  • branches/newmole/source/atom_seq_beryllium.cpp

    r2023 r2139  
    9191 
    9292        /* set the cs before if below, since we must reset to line cs in all cases */ 
    93         cs1u = t->Coll.cs; 
     93        cs1u = t->Coll.col_str; 
    9494        /* >>chng 01 sep 09, reset cs coming in, to propoer cs for the ^3P_1 level only 
    9595         * so that critical density is printed correctly with punch lines data command */ 
    96         t->Coll.cs /= 3.f; 
     96        t->Coll.col_str /= 3.f; 
    9797 
    9898        /* low density cutoff to keep matrix happy */ 

  • branches/newmole/source/atom_seq_boron.cpp

    r1830 r2139  
    175175        } 
    176176 
    177         ASSERT( t10->Coll.cs > 0.); 
    178         ASSERT( t20->Coll.cs > 0.); 
    179         ASSERT( t30->Coll.cs > 0.); 
    180         ASSERT( t21->Coll.cs > 0.); 
    181         ASSERT( t31->Coll.cs > 0.); 
    182         ASSERT( t41->Coll.cs > 0.); 
     177        ASSERT( t10->Coll.col_str > 0.); 
     178        ASSERT( t20->Coll.col_str > 0.); 
     179        ASSERT( t30->Coll.col_str > 0.); 
     180        ASSERT( t21->Coll.col_str > 0.); 
     181        ASSERT( t31->Coll.col_str > 0.); 
     182        ASSERT( t41->Coll.col_str > 0.); 
    183183        ASSERT( cs40>0.); 
    184184        ASSERT( cs32>0.); 
    … …  
    232232        AulEscp[1][0] = t10->Emis->Aul*(t10->Emis->Pesc + t10->Emis->Pelec_esc); 
    233233        AulDest[1][0] = t10->Emis->Aul*t10->Emis->Pdest; 
    234         col_str[1][0] = t10->Coll.cs; 
     234        col_str[1][0] = t10->Coll.col_str; 
    235235        AulPump[0][1] = t10->Emis->pump; 
    236236 
    … …  
    240240        AulEscp[2][0] = t20->Emis->Aul*(t20->Emis->Pesc + t20->Emis->Pelec_esc); 
    241241        AulDest[2][0] = t20->Emis->Aul*t20->Emis->Pdest; 
    242         col_str[2][0] = t20->Coll.cs; 
     242        col_str[2][0] = t20->Coll.col_str; 
    243243        AulPump[0][2] = t20->Emis->pump; 
    244244 
    245245        AulEscp[3][0] = t30->Emis->Aul*(t30->Emis->Pesc + t30->Emis->Pelec_esc); 
    246246        AulDest[3][0] = t30->Emis->Aul*t30->Emis->Pdest; 
    247         col_str[3][0] = t30->Coll.cs; 
     247        col_str[3][0] = t30->Coll.col_str; 
    248248        AulPump[0][3] = t30->Emis->pump; 
    249249 
    … …  
    255255        AulEscp[2][1] = t21->Emis->Aul*(t21->Emis->Pesc + t21->Emis->Pelec_esc); 
    256256        AulDest[2][1] = t21->Emis->Aul*t21->Emis->Pdest; 
    257         col_str[2][1] = t21->Coll.cs; 
     257        col_str[2][1] = t21->Coll.col_str; 
    258258        AulPump[1][2] = t21->Emis->pump; 
    259259 
    260260        AulEscp[3][1] = t31->Emis->Aul*(t31->Emis->Pesc + t31->Emis->Pelec_esc); 
    261261        AulDest[3][1] = t31->Emis->Aul*t31->Emis->Pdest; 
    262         col_str[3][1] = t31->Coll.cs; 
     262        col_str[3][1] = t31->Coll.col_str; 
    263263        AulPump[1][3] = t31->Emis->pump; 
    264264 
    265265        AulEscp[4][1] = t41->Emis->Aul*(t41->Emis->Pesc + t41->Emis->Pelec_esc); 
    266266        AulDest[4][1] = t41->Emis->Aul*t41->Emis->Pdest; 
    267         col_str[4][1] = t41->Coll.cs; 
     267        col_str[4][1] = t41->Coll.col_str; 
    268268        AulPump[1][4] = t41->Emis->pump; 
    269269 

  • branches/newmole/source/cddefines.h

    r2112 r2139  
    12681268 
    12691269        /** [dimensionless] collision strength of rates for transition */ 
    1270         realnum cs; 
     1270        realnum col_str; 
    12711271 
    12721272        /** [dimensionless] collision strength of rates for individual colliders */ 
    1273         realnum csi[ipNCOLLIDER]; 
     1273        realnum col_stri[ipNCOLLIDER]; 
    12741274 
    12751275        /** cooling and heating due to collisional excitation [erg s-1 cm-3] */ 

  • branches/newmole/source/cont_createmesh.cpp

    r2125 r2139  
    462462         *  
    463463         * rfield.fine_opac_nelem is the most massive (hence sharpest line) 
    464          * we will worry about.  By default this is irion but can be changed 
     464         * we will worry about.  By default this is iron but can be changed 
    465465         * with SET FINE CONTINUUM command  
    466466         *  

  • branches/newmole/source/cont_createpointers.cpp

    r1942 r2139  
    632632        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 
    633633        { 
    634                 /* do remaining part of the he iso sequence */ 
     634                /* do remaining part of the iso sequences */ 
    635635                for( nelem=2; nelem < LIMELM; nelem++ ) 
    636636                { 
    … …  
    639639                                /* generate label for this ion */ 
    640640                                sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO ); 
    641                                 /* Lya itself is the only transition below n=3 - we explicitly do not 
    642                                  * want to generate pointers for 2s-1s or 2p-2s */ 
    643                                 /* array index for continuum edges for levels in He-like ions  */ 
     641                                /* array index for continuum edges */ 
    644642                                iso.ipIsoLevNIonCon[ipISO][nelem][0] =  
    645643                                        ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab); 
    … …  
    647645                                for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 
    648646                                { 
    649                                         /* array index for continuum edges for levels in He-like ions  */ 
     647                                        /* array index for continuum edges */ 
    650648                                        iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab); 
    651649 
    652                                         /* define all he-like line pointers */ 
     650                                        /* define all line pointers */ 
    653651                                        for( ipLo=0; ipLo < ipHi; ipLo++ ) 
    654652                                        { 

  • branches/newmole/source/conv_base.cpp

    r2125 r2139  
    487487                                                /* cs is actually the negative of the branching ratio for autoionization,  
    488488                                                * rateone is inverse lifetime of level against autoionization */ 
    489                                                 double rateone = UTALines[i].Emis->pump * -UTALines[i].Coll.cs; 
     489                                                double rateone = UTALines[i].Emis->pump * -UTALines[i].Coll.col_str; 
    490490                                                ionbal.UTA_ionize_rate[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] += rateone; 
    491491                                                /* heating rate in erg atom-1 s-1 */ 
    … …  
    609609                } 
    610610 
     611                fixit();/* following test is very coarse rel error of 0.2 - should be much 
     612                                 * smaller error.  Due to source/sink terms in ionization balance 
     613                                 * equaitons */ 
    611614                for(long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 
    612615                { 
    … …  
    618621                                                        iso.xIonSimple[ipISO][nelem]<1e-10) ) 
    619622                                { 
    620                                         const realnum err_tol = 1e-2f; 
    621                                         realnum abund = 0.; 
     623                                        /* sim blr_hizqso has largest errors in following */ 
     624                                        double abund = 0.; 
    622625                                        for (long n=0; n<iso.numLevels_local[ipISO][nelem]; ++n) 
    623626                                        { 
    624                                                 abund += (realnum)StatesElem[ipISO][nelem][n].Pop; 
     627                                                abund += StatesElem[ipISO][nelem][n].Pop; 
    625628                                        } 
    626                                         if( conv.lgConvIoniz && iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 
    627                                                         ! fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 
    628                                                                                                                  err_tol*dense.xIonDense[nelem][nelem-ipISO])) 
     629                                        double rel_err = fabs(abund*dense.xIonDense[nelem][nelem+1-ipISO]-dense.xIonDense[nelem][nelem-ipISO])/ 
     630                                                SDIV(dense.xIonDense[nelem][nelem-ipISO]); 
     631 
     632                                        fixit();/* the following test should pass cleanly, rm the 0 && to  
     633                                                         * activate.  when it passes rm this block and  
     634                                                         * change err_tol in following block to a small number */ 
     635                                        if(0 &  conv.lgConvIoniz &&  
     636                                                iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 
     637                                                dense.xIonDense[nelem][nelem+1-ipISO]>SMALLFLOAT && 
     638                                                rel_err > 0.02f )  
    629639                                        { 
    630                                                 //fprintf(ioQQQ,"Inconsistent iso %ld %ld %ld %g %g %g\n",loop_ion,nelem,ipISO, 
    631                                                 //                              abund*dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO])-1 
    632                                                 //                              ,abund*dense.xIonDense[nelem][nelem+1-ipISO],(dense.xIonDense[nelem][nelem-ipISO]) 
    633                                                 //      ); 
     640                                                /* remove this block when error is made small - this block 
     641                                                 * is here */ 
     642                                                fprintf(ioQQQ,"PROBLEM ISO/ion do not agree ipISO=%li " 
     643                                                        "nelem=%li iso soln=%.2e ion soln=%.2e rel_err %.2e \n", 
     644                                                        ipISO , nelem , 
     645                                                        abund*dense.xIonDense[nelem][nelem+1-ipISO], 
     646                                                        dense.xIonDense[nelem][nelem-ipISO] , 
     647                                                        rel_err); 
     648                                        } 
     649 
     650                                        const realnum err_tol = 0.4f; 
     651                                        if( conv.lgConvIoniz &&  
     652                                                iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 
     653                                                dense.xIonDense[nelem][nelem+1-ipISO]>SMALLFLOAT && 
     654                                                rel_err > err_tol )  
     655                                        { 
     656                                                //fprintf(ioQQQ,"Inconsistent iso %ld %ld %ld %g\n",loop_ion,nelem,ipISO, 
     657                                                //                              abund*dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO])-1); 
    634658                                                conv.lgConvIoniz = false; 
    635659                                                sprintf( conv.chConvIoniz , "Iso!=ion" ); 
    … …  
    856880 
    857881        ASSERT(lgElemsConserved()); 
     882 
    858883 
    859884        /* update some counters that keep track of how many times this routine 

  • branches/newmole/source/conv_fail.cpp

    r1780 r2139  
    119119                if( called.lgTalk ) 
    120120                { 
    121                         fprintf( ioQQQ, " PROBLEM  ConvFail %li, %s ionization not converged iteration %li zone %li fnzone %.2f \n",