Changeset 1942

Show
Ignore:
Timestamp:
04/12/08 09:30:35 (9 months ago)
Author:
rjrw
Message:

Merged from trunk r1917:1941

Location:
branches/newmole/source
Files:
46 modified

Legend:

Unmodified
Added
Removed

  • branches/newmole/source/assert_results.cpp

    r1918 r1942  
    1010#include "iso.h" 
    1111#include "called.h" 
     12#include "atmdat.h" 
    1213#include "hcmap.h" 
    1314#include "thermal.h" 
    … …  
    4748static char **chAssertLineLabel; 
    4849 
     50/* element to assert H-like case B */ 
     51static char chElemLabelCaseB[5]; 
     52static long int nelemCaseB , nISOCaseB; 
     53 
    4954/* this will be = for equal, < or > for limit */ 
    5055static char *chAssertLimit; 
    … …  
    344349 
    345350        /* assert line "LINE"  --  key is ine_ since linear option appears on some commands */ 
    346         else if( nMatch("INE ",input.chCARDCAPS ) ) 
     351        else if( nMatch(" LINE ",input.chCARDCAPS ) ) 
    347352        { 
    348353                if(  nMatch("LUMI",input.chCARDCAPS ) || nMatch("INTE",input.chCARDCAPS )) 
    … …  
    440445 
    441446        /* assert departure coefficients */ 
     447        else if( nMatch("CASE",input.chCARDCAPS ) ) 
     448        { 
     449                /* this is Case B for some element */ 
     450                strcpy( chAssertType[nAsserts] , "CB" ); 
     451                i = 5; 
     452                /* this is relative error */ 
     453                AssertError[nAsserts] = 
     454                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 
     455                if( lgEOL ) 
     456                        /* default error was set in define above */ 
     457                        AssertError[nAsserts] = DEF_ERROR; 
     458                AssertQuantity[nAsserts] = 0; 
     459                strcpy( chAssertLineLabel[nAsserts] , "Ca B" );  
     460                wavelength[nAsserts] = 0.; 
     461 
     462                /* assert case b for H or He */ 
     463                if( nMatch("H-LI",input.chCARDCAPS ) ) 
     464                { 
     465                        /* H-like - now get an element */ 
     466                        if( (nelem = GetElem( input.chCARDCAPS )) < 0 ) 
     467                        { 
     468                                /* no name found */ 
     469                                fprintf(ioQQQ, "assert H-like case B did not find an element on this line, sorry %s\n", 
     470                                        input.chCARDCAPS ); 
     471                                cdEXIT(EXIT_FAILURE); 
     472                        } 
     473                        if( nelem>7 ) 
     474                        { 
     475                                /* beyond reach of tables */ 
     476                                fprintf(ioQQQ, "assert H-like cannot do elements heavier than O, sorry %s\n", 
     477                                        input.chCARDCAPS ); 
     478                                cdEXIT(EXIT_FAILURE); 
     479                        } 
     480                        nelemCaseB = nelem; 
     481                        nISOCaseB = 0; 
     482                        strcpy( chElemLabelCaseB , elementnames.chElementSym[nelem] ); 
     483                        char chNumb[4]; 
     484                        sprintf( chNumb , "%2li" , nelem+1 ); 
     485                        strcat( chElemLabelCaseB , chNumb ); 
     486                } 
     487                else if( nMatch("HE-L",input.chCARDCAPS ) ) 
     488                { 
     489                        /* He-like - only helium itself */ 
     490                        nelemCaseB = ipHELIUM; 
     491                        nISOCaseB = 0; 
     492                        strcpy( chElemLabelCaseB , "He 1" ); 
     493                        /* beyond reach of tables */ 
     494                        fprintf(ioQQQ, "assert He-like case b not implemented yet %s\n", 
     495                                input.chCARDCAPS ); 
     496                        cdEXIT(EXIT_FAILURE); 
     497                } 
     498        } 
    442499        else if( nMatch("DEPA",input.chCARDCAPS ) ) 
    443500        { 
    … …  
    447504                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 
    448505                if( lgEOL ) 
    449                 { 
    450506                        NoNumb(input.chOrgCard); 
    451                 } 
    452507 
    453508                /* this is relative error, max departure from unity of any level or std */ 
    … …  
    455510                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 
    456511                if( lgEOL ) 
    457                 { 
    458512                        /* default error was set in define above */ 
    459513                        AssertError[nAsserts] = DEF_ERROR; 
    460                 } 
    461514 
    462515                /* H-like key means do one of the hydrogenic ions */ 
    … …  
    18281881                        /* uncertainty in difference is zero */ 
    18291882                        AssertError[i] = 0.; 
     1883                } 
     1884 
     1885                else if( strcmp( chAssertType[i] , "CB") == 0 ) 
     1886                { 
     1887                        int iCase = 1; 
     1888                        RelError[i] = 0.; 
     1889                        fprintf(ioQQQ,"        Species  nHi nLo  Wl  Computed  Asserted       error\n"); 
     1890                        long nHighestPrinted = StatesElem[nISOCaseB][nelemCaseB][iso.numPrintLevels[nISOCaseB][nelemCaseB]-1].n; 
     1891                        for( long int ipLo=1+iCase; ipLo< MIN2(24,nHighestPrinted-1); ++ipLo ) 
     1892                        { 
     1893                                for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi ) 
     1894                                { 
     1895                                        /* assert the line */ 
     1896                                        realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo]; 
     1897                                        double relint , absint,CBrelint , CBabsint; 
     1898                                        cdLine( "Ca B" , wl , &CBrelint , &CBabsint ); 
     1899                                        CBabsint = pow( 10., CBabsint ); 
     1900                                        double error; 
     1901                                        if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0) 
     1902                                        { 
     1903                                                absint = pow( 10., absint ); 
     1904                                                error = (CBabsint - absint)/MAX2(CBabsint , absint); 
     1905                                                fprintf(ioQQQ,"ChkAssert %s %3li %3li ",  
     1906                                                        chElemLabelCaseB , ipHi , ipLo ); 
     1907                                                prt_wl(ioQQQ, wl ); 
     1908                                                fprintf(ioQQQ," %.2e %.2e %10.3f",  
     1909                                                        absint , CBabsint , error ); 
     1910                                        } 
     1911                                        else 
     1912                                                TotalInsanity(); 
     1913                                        if( fabs(error) > AssertError[i]  ) 
     1914                                                fprintf(ioQQQ , " BOTCH \n"); 
     1915                                        else 
     1916                                                fprintf(ioQQQ , "\n"); 
     1917 
     1918                                        PredQuan[i] = 0; 
     1919                                        AssertQuantity[i] = 0.; 
     1920                                        RelError[i] = MAX2( RelError[i] , fabs(error) ); 
     1921                                } 
     1922                        } 
     1923                        fprintf(ioQQQ,"\n"); 
    18301924                } 
    18311925 
    … …  
    23922486 
    23932487                                { 
    2394                                         /*@-redef@*/ 
    23952488                                        enum {DEBUG_LOC=false}; 
    2396                                         /*@+redef@*/ 
    23972489 
    23982490                                        if( DEBUG_LOC ) 

  • branches/newmole/source/atmdat.h

    r1830 r1942  
    202202         * with "no collisional ionization" command */ 
    203203        bool lgCollIonOn; 
     204 
     205        /** wavelengths of Hummer & Storey case B lines for H - O */ 
     206        realnum WaveLengthCaseB[8][25][24]; 
    204207 
    205208        }       atmdat; 

  • branches/newmole/source/atmdat_chianti.cpp

    r1918 r1942  
    263263                strncpy(atmolStates[intNS][j].chLabel,Species[intNS].chptrSpName, 4); 
    264264                atmolStates[intNS][j].chLabel[4] = '\0'; 
    265  
    266                 /*Index for keep track of species index*/ 
    267                 atmolStates[intNS][j].AtmolIndex = intNS; 
    268265 
    269266                fstatwt = (realnum)atof(&chLine[SWS]); 

  • branches/newmole/source/atmdat_lamda.cpp

    r1918 r1942  
    141141                strncpy(atmolStates[intNS][nMolLevs].chLabel,Species[intNS].chptrSpName, 4); 
    142142                atmolStates[intNS][nMolLevs].chLabel[4] = '\0'; 
    143                 /*Keeping track of species index*/ 
    144                 atmolStates[intNS][nMolLevs].AtmolIndex = intNS; 
    145143 
    146144                long i = 1; 

  • branches/newmole/source/cddefines.h

    r1918 r1942  
    12791279        char chConfig[11]; 
    12801280 
    1281         /*Index to identify the quantum state with the species index*/ 
    1282         int AtmolIndex; 
    1283  
    12841281        species   *sp; 
    12851282 

  • branches/newmole/source/cont_createpointers.cpp

    r1918 r1942  
    607607                                        atmolEmis[i].tran->EnergyWN/PI4); 
    608608                atmolEmis[i].damp = -1000.0; 
    609                 //atmolEmis[i].damp = fdampXvel/(realnum)atmolEmis[i].tran->Hi->lifetime ; 
    610609                /*Put null terminated line label into chLab*/ 
    611610                strncpy(chLab,atmolEmis[i].tran->Hi->chLabel,4); 

  • branches/newmole/source/cpu.cpp

    r1830 r1942  
    208208 
    209209        /* if the environment variable was not set, the default set in path.h takes effect */ 
    210         string chSearchPathRaw = ( path != NULL ) ? string( path ) : string( LOCAL_DATA_PATH ); 
     210        string chSearchPathRaw = ( path != NULL ) ? string( path ) : string( CLOUDY_DATA_PATH ); 
    211211 
    212212#       ifdef _WIN32 
    … …  
    490490                        fprintf( ioQQQ, "the most likely reason is that the path has not been properly set.\n"); 
    491491                        fprintf( ioQQQ, "Please have a look at the file path.h in the source directory.\n"); 
    492                         fprintf( ioQQQ, "Check how the variable LOCAL_DATA_PATH is set - "); 
     492                        fprintf( ioQQQ, "Check how the variable CLOUDY_DATA_PATH is set - "); 
    493493                        fprintf( ioQQQ, "it should give the location of the data files I need.\n"); 
    494494                        fprintf( ioQQQ, "These are the files in the data download from the web site.\n"); 

  • branches/newmole/source/helike_cs.cpp

    r1918 r1942  
    11541154        global_s = s; 
    11551155        global_temp = temp; 
    1156         global_collider_charge = 1.; 
    11571156        global_Collider = Collider; 
     1157        global_collider_charge = ColliderCharge[Collider]; 
     1158        ASSERT( global_collider_charge > 0. ); 
    11581159 
    11591160        /* no need to do this for h-like */ 

  • branches/newmole/source/helike_einsta.cpp

    r1918 r1942  
    353353                        /* This is an extrapolation to the data given above.  The fit reproduces 
    354354                         * the above values to 10% or better. */ 
    355                         A= 7.22E-9*pow(nelem, 9.33); 
     355                        A= 7.22E-9*pow((double)nelem, 9.33); 
    356356                        iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f); 
    357357                } 

  • branches/newmole/source/helike_recom.cpp

    r1739 r1942  
    21112111        double lambda = TE1RYD * nelem * nelem / temp; 
    21122112        /* smallest x ever used here should be lowest Z, highest T, highest n... 
    2113          * using helium, logt = 10., and n = 100, we get xmin = 1.5789E-9.      */ 
     2113         * using helium, logt = 10., and n = 1000, we get xmin = 1.5789E-11.    */ 
    21142114        double x = lambda / n / n; 
    21152115        double AlphaN; 

  • branches/newmole/source/helike_recom.h

    r1739 r1942  
    44#ifndef _HELIKE_RECOM_H_ 
    55#define _HELIKE_RECOM_H_ 
    6  
    7  
    8 /* following two macros used to define recombination coef arrays */ 
    9 /* Max n desired in RRCoef file.        */ 
    10 /* >>chng 02 jan 28, from 10 to 20 so large HeI tests are fast */ 
    11 /** this is the number of levels used with the 
    12  * atom he-like levels large command */ 
    13 /* Helium itself will have precompiled recombination coefficients up to this maximum n. */ 
    14 #define HE_RREC_MAXN            40 
    15  
    16 /** Ions of the sequence will go up to this n.  */ 
    17 #define HE_LIKE_RREC_MAXN               20 
    18  
    19 #define N_HE_TE_RECOMB          41 
    20  
    21 /** This is the n to go up to when calculating total recombination.     Any change  
    22  * here will not be reflected in total recomb until "compile helike" is run     */ 
    23 #define SumUpToThisN    1000 
    24 /** the magic number for the table of recombination coefficients, YYMMDD */ 
    25 /* >>chng 02 mar 08, data structure the same but improved fits to He photo cs */ 
    26 #define RECOMBMAGIC             (51213) 
    276 
    287double DielecRecombRate( long ipISO, long nelem, long ipLo ); 

  • branches/newmole/source/hydro_bauman.cpp

    r1846 r1942  
    99/* In this version, quantities that would normally cause a 64-bit floating point processor      */ 
    1010/* to either underflow or overflow are evaluated using logs instead of floating point math.     */ 
    11 /* This allows us to use an upper principle quantum number `n' greater than  which the          */ 
     11/* This allows us to use an upper principal quantum number `n' greater than  which the          */ 
    1212/* other version begins to fail. The trade-off is, of course, lower accuracy                    */ 
    1313/* ( or is it precision ). We use LOG_10 for convenience.                                       */ 
    … …  
    4343/*    The recombination coefficient is obtained from the                                        */ 
    4444/*    photoionization cross-section (see Burgess 1965).                                         */ 
    45 /*    We have,                                                                                  */ 
     45/*    We have:                                                                                  */ 
    4646/*                                                                                              */ 
    4747/*                  -                           -               l'=l+1                          */ 
    … …  
    124124/*                                               r^2                                            */ 
    125125/*                                                                                              */ 
    126 /*      Similiarly, the recombination co-efficient satisfies                                    */ 
     126/*      Similiarly, the recombination coefficient satisfies                                    */ 
    127127/************************************************************************************************/ 
    128128 
    … …  
    144144/** \verbatim 
    145145************************* for LOG version of the file  *************************************** 
    146  In this version, quantites that would normal cause a 64-bit floating point processor         
     146 In this version, quantities that would normal cause a 64-bit floating point processor         
    147147 to underflowed or overflow on intermediate values (ones internal to the calculation)          
    148  are evaluated using logs. This allows us to use an upper principle quantum number `n'         
     148 are evaluated using logs. This allows us to use an upper principal quantum number `n'         
    149149 greater than 50 which is where the other version begins to fail. The trade-off is,            
    150150 of course, lower accuracy( or is it precision ) and perhaps speed.                            
    … …  
    352352 
    353353/** \verbatim 
    354      This routine, hri(), calculates the hydrogen radial intergral,   
     354     This routine, hri(), calculates the hydrogen radial integral,   
    355355      for the transition n,l --> n',l'                                 
    356356      It is, of course, dimensionless.                                 
    … …  
    373373 
    374374/**\verbatim 
    375   This routine, hri_log10(), calculates the hydrogen radial intergral,   
     375  This routine, hri_log10(), calculates the hydrogen radial integral,   
    376376  for the transition n,l --> n',l'                                       
    377377  It is, of course, dimensionless.                                       
    … …  
    539539 
    540540/************************************************************************/ 
    541 /* (4.0/3.0) * PI * FINE_STRUCTURE_CONSTANT * BORHRADIUS * BORHRADIUS   */ 
     541/* (4.0/3.0) * PI * FINE_STRUCTURE_CONSTANT * BOHRRADIUS * BOHRRADIUS   */ 
    542542/*                                                                      */ 
    543543/*                                                                      */ 
    … …  
    742742        /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd)    */ 
    743743        double K, 
    744         /* principle quantum number                             */ 
     744        /* principal quantum number                             */ 
    745745        long int n, 
    746746        /* angular  momentum quantum number                     */ 
    747747        long int l, 
    748         /* Temperary storage for intermediate                   */ 
     748        /* Temporary storage for intermediate                   */ 
    749749        /*  results of the recursive routine                    */ 
    750750        double *rcsvV 
    … …  
    754754 
    755755        long int lp = 0; /* l' */ 
    756         double sigma = 0.; /*      Sum in brocklehurst eq. 3.13 */ 
     756        double sigma = 0.; /*      Sum in Brocklehurst eq. 3.13 */ 
    757757 
    758758        ASSERT( n > 0 ); 
    … …  
    779779        /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd)    */ 
    780780        double K, 
    781         /* principle quantum number                             */ 
     781        /* principal quantum number                             */ 
    782782        long int n, 
    783783        /* angular  momentum quantum number                     */ 
    784784        long int l, 
    785         /* Temperary storage for intermediate                   */ 
     785        /* Temporary storage for intermediate                   */ 
    786786        mxq *rcsvV_mxq 
    787787        /*   results of the recursive routine                   */ 
    … …  
    791791 
    792792        long int lp = 0; /* l' */ 
    793         double sigma = 0.; /*      Sum in brocklehurst eq. 3.13 */ 
     793        double sigma = 0.; /*      Sum in Brocklehurst eq. 3.13 */ 
    794794 
    795795        ASSERT( n > 0 ); 
    … …  
    835835        long int l, 
    836836        long int lp, 
    837         /* Temperary storage for intermediate         */ 
     837        /* Temporary storage for intermediate         */ 
    838838        /*  results of the recursive routine          */ 
    839839        double *rcsvV 
    … …  
    913913        long int l, 
    914914        long int lp, 
    915         /* Temperary storage for intermediate         */ 
     915        /* Temporary storage for intermediate         */ 
    916916        /* results of the recursive routine           */ 
    917917        mxq *rcsvV_mxq 
    … …  
    10101010        long int l, 
    10111011        long int lp, 
    1012         /* Temperary storage for intermediate         */ 
     1012        /* Temporary storage for intermediate         */ 
    10131013        /* results of the recursive routine           */ 
    10141014        double *rcsvV 
    … …  
    11041104        long int l, 
    11051105        long int lp, 
    1106         /* Temperary storage for intermediate        */ 
     1106        /* Temporary storage for intermediate        */ 
    11071107        /* results of the recursive routine          */ 
    11081108        mxq *rcsvV_mxq 
    … …  
    11211121        double Ksqrd = K * K; 
    11221122 
    1123         mx GK_mx = {0.0,0}, null_mx = {0.0,0}; 
     1123        mx GK_mx = {0.0,0}; 
    11241124 
    11251125        /* l=l'-1 or l=l'+1 */ 
    … …  
    13061306/*                                                                                              */ 
    13071307/*    from the reference                                                                        */ 
    1308 /*    M. BrocKlehurst                                                                           */ 
     1308/*    M. Brocklehurst                                                                           */ 
    13091309/*    Mon. Note. R. astr. Soc. (1971) 153, 471-490                                              */ 
    13101310/*                                                                                              */ 
    … …  
    13221322        long int l, 
    13231323        long int lp, 
    1324         /* Temperary storage for intermediate        */ 
     1324        /* Temporary storage for intermediate        */ 
    13251325        /*  results of the recursive routine         */ 
    13261326        double *rcsvV, 
    … …  
    14731473        long int l, 
    14741474        long int lp, 
    1475         /* Temperary storage for intermediate       */ 
     1475        /* Temporary storage for intermediate       */ 
    14761476        /* results of the recursive routine         */ 
    14771477        mxq *rcsvV_mxq, 
    … …  
    16821682/*                                                                                              */ 
    16831683/*    from the reference                                                                        */