Show
Ignore:
Timestamp:
05/10/08 09:03:02 (4 months ago)
Author:
peter
Message:

Merging all changes from mainline upto r2033, except r1902.

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • branches/c08_branch/source/iso_radiative_recomb.cpp

    r1852 r2034  
    4141        double alpha; 
    4242 
     43        DEBUG_ENTRY( "iso_radrecomb_from_cross_section()" ); 
     44 
    4345        if( ipISO == ipH_LIKE ) 
    4446                alpha =  hlike_radrecomb_from_cross_section( temp, nelem, ipLo); 
    … …  
    6163        double topoff, TotMinusExplicitResolved, 
    6264                TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.; 
    63         double RecTriplets; 
    64         static double RecExplictLevels, TotalRadRecomb, RecCollapsed; 
     65        double RecExplictLevels, TotalRadRecomb, RecCollapsed; 
    6566        static double TeUsed[NISO][LIMELM]={ 
    6667                {0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
    … …  
    7879        /* this will be the sum of all levels explicitly included in the model atom */ 
    7980        RecExplictLevels = 0.; 
    80         RecTriplets = 0.; 
    81  
    82         /* this is the number of resolved levels, so first collapsed level is ipFirstCollapsed */ 
    83         /* >>chng 06 aug 17, from _max to _local */ 
     81 
     82        /* number of resolved levels, so first collapsed level is [ipFirstCollapsed] */ 
    8483        ipFirstCollapsed = iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem];  
    8584 
    … …  
    9594                        if( !iso.lgNoRecombInterp[ipISO] ) 
    9695                        { 
    97                                 if( ipISO == ipH_LIKE && ipLevel <= 2 ) 
    98                                 { 
    99                                         RadRec = t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, ipLevel, phycon.te); 
    100                                 } 
    101                                 else 
    102                                 { 
    103                                         RadRec = iso_RRCoef_Te( ipISO, nelem , ipLevel ); 
    104                                 } 
     96                                RadRec = iso_RRCoef_Te( ipISO, nelem , ipLevel ); 
    10597                        } 
    10698                        else 
    … …  
    120112                        } 
    121113 
     114                        ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. ); 
     115 
    122116                        RecExplictLevels += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]; 
    123                         if( S_(ipLevel)==3 && ipISO == ipHE_LIKE ) 
    124                                 RecTriplets += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]; 
    125  
    126                         ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. ); 
    127117 
    128118                        if( iso.lgDielRecom[ipISO] ) 
    … …  
    133123                        } 
    134124                } 
    135  
    136                 /* do not include ground in this sum    */ 
    137                 RecExplictLevels -= iso.RadRecomb[ipISO][nelem][0][ipRecRad];    
    138125 
    139126                /**************************************************/ 
    … …  
    176163                        } 
    177164                } 
    178  
    179                 if( ipISO==ipHE_LIKE ) 
    180                         RecTriplets += 0.75 * RecCollapsed; 
    181165 
    182166                /* >>chng 06 aug 17, from numLevels_max to numLevels_local */ 
    … …  
    220204                } 
    221205 
     206                /* Get total recombination into all levels, including those not explicitly considered. */ 
    222207                if( iso.lgNoRecombInterp[ipISO] ) 
    223208                { 
    224                         /* Get total recombination into all levels, including those not explicitly considered. */ 
    225                         if( ipISO == ipH_LIKE ) 
    226                         { 
    227                                 /* following returns case A rec coef */ 
    228                                 TotalRadRecomb = t_ADfA::Inst().H_rad_rec(nelem+1,-1,phycon.te); 
    229                         } 
    230                         else if( ipISO == ipHE_LIKE ) 
    231                         { 
    232                                 /* We are not interpolating, must calculate total now, as sum of explicit 
    233                                  * resolved and collapsed levels... */ 
    234                                 TotalRadRecomb = RecCollapsed + RecExplictLevels + 
    235                                         iso.RadRecomb[ipISO][nelem][ipHe1s1S][ipRecRad]; 
    236  
    237                                 /* Plus additional levels out to the macro SumUpToThisN */ 
    238                                 for( ipLevel= N_(ipFirstCollapsed) + iso.nCollapsed_max[ipISO][nelem]; ipLevel<=SumUpToThisN; ipLevel++ ) 
    239                                 { 
    240                                         TotalRadRecomb += Recomb_Seaton59( nelem, phycon.te, ipLevel ); 
    241                                 } 
    242                         } 
    243                         else 
    244                         { 
    245                                 /* Only remove this if other iso sequences are added. */ 
    246                                 TotalInsanity(); 
     209                        /* We are not interpolating, must calculate total now, as sum of resolved and collapsed levels... */ 
     210                        TotalRadRecomb = RecCollapsed + RecExplictLevels; 
     211 
     212                        /* Plus additional levels out to a predefined limit... */ 
     213                        for( long nLo = N_(ipFirstCollapsed) + iso.nCollapsed_max[ipISO][nelem]; nLo < NHYDRO_MAX_LEVEL; nLo++ ) 
     214                        { 
     215                                TotalRadRecomb += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, nLo, phycon.te); 
     216                        } 
     217                        /* Plus a bunch more for good measure */ 
     218                        for( long nLo = NHYDRO_MAX_LEVEL; nLo<=SumUpToThisN; nLo++ ) 
     219                        { 
     220                                TotalRadRecomb += Recomb_Seaton59( nelem+1-ipISO, phycon.te, nLo ); 
    247221                        } 
    248222                } 
    … …  
    250224                { 
    251225                        /* We are interpolating, and total was calculated here in iso_recomb_setup */ 
    252                         TotalRadRecomb = iso_RRCoef_Te( ipISO, nelem , 
    253                                 iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] ); 
     226                        TotalRadRecomb = iso_RRCoef_Te( ipISO, nelem, iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] ); 
    254227                } 
    255228 
    … …  
    264237                } 
    265238 
    266                 if( ipISO == ipHE_LIKE ) 
    267                         RecTriplets += 0.75*(TotalRadRecomb - RecCollapsed - RecExplictLevels - 
    268                                         iso.RadRecomb[ipISO][nelem][ipHe1s1S][ipRecRad] ); 
    269  
    270239                /* this is case B recombination, sum without the ground included */ 
    271240                iso.RadRec_caseB[ipISO][nelem] = TotalRadRecomb - iso.RadRecomb[ipISO][nelem][0][ipRecRad]; 
    272241                ASSERT( iso.RadRec_caseB[ipISO][nelem] > 0.); 
    273  
    274                 /** \todo 2 RecTriplets can probably be removed entirely.  */ 
    275                 if( ipISO == ipHE_LIKE ) 
    276                         iso.RadRec_caseB[ipISO][nelem] = RecTriplets; 
    277  
    278                 /* this is used for radiative recombination cooling in CoolEvaluate, this does not include 
    279                  * optical depth effects - may be reevaluated in level solver with opt dep included */ 
    280                 /** \todo 1 why does he-like include total recomb, but h-like is only case b? */ 
    281                 if( ipISO == ipH_LIKE ) 
    282                         ionbal.RR_rate_coef_used[nelem][nelem-ipISO] = iso.RadRec_caseB[ipISO][nelem]; 
    283                 else 
    284                         ionbal.RR_rate_coef_used[nelem][nelem-ipISO] = TotalRadRecomb; 
    285242 
    286243                /**********************************************************************/ 
    … …  
    293250                         * to all levels explicitly included in the model atom the total  
    294251                         * recombination rate.  The difference is the amount of "topoff" we will need to do */ 
    295                         TotMinusExplicitResolved = TotalRadRecomb -  
    296                                 iso.RadRecomb[ipISO][nelem][0][ipRecRad] - RecExplictLevels; 
     252                        TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels; 
    297253 
    298254                        topoff = TotMinusExplicitResolved - RecCollapsed; 
    … …  
    337293        } 
    338294 
     295        /**************************************************************/ 
     296        /***  Stuff escape probabilities, and effective rad recomb  ***/ 
     297        /**************************************************************/ 
     298 
     299        /* total effective radiative recombination, initialize to zero */ 
     300        iso.RadRec_effec[ipISO][nelem] = 0.; 
     301 
     302        for( ipLevel=0; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel ) 
     303        { 
     304                /* option for case b conditions, kill ground state recombination */ 
     305                if( opac.lgCaseB && ipLevel==0 ) 
     306                { 
     307                        iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 1e-10; 
     308                        iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 1e-10; 
     309                } 
     310                else 
     311                { 
     312                        iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] =  
     313                                RT_recom_effic(iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]); 
     314 
     315                        /* net escape prob includes dest by background opacity */ 
     316                        iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] =  
     317                                MIN2( 1., 
     318                                iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] + 
     319                                (1.-iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc])* 
     320                                iso.ConOpacRatio[ipISO][nelem][ipLevel] ); 
     321                } 
     322 
     323                ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] >= 0. ); 
     324                ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] >= 0. ); 
     325                ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] <= 1. ); 
     326 
     327                /* sum of all effective rad rec */ 
     328                iso.RadRec_effec[ipISO][nelem] += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]* 
     329                  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]; 
     330        } 
     331 
     332        /* zero out escape probabilities of levels above numLevels_local */ 
     333        for( ipLevel=iso.numLevels_local[ipISO][nelem]; ipLevel<iso.numLevels_max[ipISO][nelem]; ++ipLevel ) 
     334        { 
     335                /* this is escape probability */ 
     336                iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 0.;  
     337                /* net escape prob includes dest by background opacity */ 
     338                iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 0.; 
     339        } 
     340 
     341        /* trace escape probabilities */ 
     342        if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) ) 
     343        { 
     344                fprintf( ioQQQ, "       iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem ); 
     345 
     346                /* print continuum escape probability */ 
     347                fprintf( ioQQQ, "       iso_radiative_recomb recomb effic" ); 
     348                for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
     349                { 
     350                        fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] )); 
     351                } 
     352                fprintf( ioQQQ, "\n" ); 
     353 
     354                /* net recombination efficiency factor, including background opacity*/ 
     355                fprintf( ioQQQ, "       iso_radiative_recomb recomb net effic" ); 
     356                for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
     357                { 
     358                        fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]) ); 
     359                } 
     360                fprintf( ioQQQ, "\n" ); 
     361 
     362                /* inward continuous optical depths */ 
     363                fprintf( ioQQQ, "       iso_radiative_recomb in optic dep" ); 
     364                for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
     365                {        
     366                        fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] )); 
     367                } 
     368                fprintf( ioQQQ, "\n" ); 
     369 
     370                /* outward continuous optical depths*/ 
     371                fprintf( ioQQQ, "       iso_radiative_recomb out op depth" ); 
     372                for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
     373                { 
     374                        fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] )); 
     375                } 
     376                fprintf( ioQQQ, "\n" ); 
     377 
     378                /* print radiative recombination coefficients */ 
     379                fprintf( ioQQQ, "       iso_radiative_recomb rad rec coef " ); 
     380                for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
     381                { 
     382                        fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]) ); 
     383                } 
     384                fprintf( ioQQQ, "\n" ); 
     385        } 
     386 
     387        if( trace.lgTrace ) 
     388        { 
     389                fprintf( ioQQQ, "     iso_radiative_recomb total rec coef" ); 
     390                fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_effec[ipISO][nelem] )); 
     391                fprintf( ioQQQ, " case A=" ); 
     392                fprintf( ioQQQ,PrintEfmt("%10.3e",  
     393                        iso.RadRec_caseB[ipISO][nelem] + iso.RadRecomb[ipISO][nelem][ipH1s][ipRecRad] ) ); 
     394                fprintf( ioQQQ, " case B="); 
     395                fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_caseB[ipISO][nelem] )); 
     396                fprintf( ioQQQ, "\n" ); 
     397        } 
     398 
     399        /****************************/ 
     400        /***  begin sanity check  ***/ 
     401        /****************************/ 
     402        { 
     403                bool lgOK = true; 
     404                for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
     405                { 
     406                        if( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] <= 0. ) 
     407                        { 
     408                                fprintf( ioQQQ,  
     409                                        " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",  
     410                                  ipISO, nelem, ipLevel, iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] , phycon.te); 
     411                                        lgOK = false; 
     412                        } 
     413                } 
     414                /* bail if we found problems */ 
     415                if( !lgOK ) 
     416                { 
     417                        ShowMe(); 
     418                        cdEXIT(EXIT_FAILURE); 
     419                } 
     420                /*end sanity check */ 
     421        } 
     422 
     423        /* confirm that we have good rec coef at bottom and top of atom/ion */ 
     424        ASSERT( iso.RadRecomb[ipISO][nelem][0][ipRecRad] > 0. ); 
     425        ASSERT( iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecRad] > 0. ); 
     426 
     427        /* set true when punch recombination coefficients command entered */ 
     428        if( punch.lgioRecom ) 
     429        { 
     430                /* this prints Z on physical, not C, scale */ 
     431                fprintf( punch.ioRecom, "%s %s %2li ",  
     432                        iso.chISO[ipISO], elementnames.chElementSym[nelem], nelem+1 ); 
     433                fprintf( punch.ioRecom,PrintEfmt("%9.2e", iso.RadRec_caseB[ipISO][nelem] )); 
     434                fprintf( punch.ioRecom, "\n" ); 
     435        } 
     436 
     437        return; 
     438} 
     439 
     440void iso_radiative_recomb_effective( long ipISO, long nelem ) 
     441{ 
     442        DEBUG_ENTRY( "iso_radiative_recomb_effective()" ); 
     443 
    339444        /* Find effective recombination coefficients */ 
    340445        for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 
    … …  
    379484        if( iso.lgRandErrGen[ipISO] ) 
    380485        { 
    381                 long ipLo, ipHi, ipHigher; 
    382  
    383486                dprintf( ioQQQ, "ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" ); 
    384487 
    385488                /* >>chng 06 aug 17, from numLevels_max to numLevels_local */ 
    386                 for( ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 
     489                for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 
    387490                { 
    388491                        iso.SigmaRadEffec[ipISO][nelem][ipHi] = 0.; 
    389492 
    390493                        /* >>chng 06 aug 17, from numLevels_max to numLevels_local */ 
    391                         for( ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ ) 
     494                        for( long ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ ) 
    392495                        { 
    393496                                ASSERT( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][ipHigher][IPRAD] >= 0. ); 
    … …  
    403506                        iso.SigmaRadEffec[ipISO][nelem][ipHi] = sqrt( iso.SigmaRadEffec[ipISO][nelem][ipHi] ); 
    404507 
    405                         for( ipLo = 0; ipLo < ipHi; ipLo++ ) 
     508                        for( long ipLo = 0; ipLo < ipHi; ipLo++ ) 
    406509                        { 
    407510                                if( (( L_(ipLo) == L_(ipHi) + 1 ) || ( L_(ipHi) == L_(ipLo) + 1 )) ) 
    … …  
    438541        } 
    439542 
    440         /**************************************************************/ 
    441         /***  Stuff escape probabilities, and effective rad recomb  ***/ 
    442         /**************************************************************/ 
    443  
    444         /* total effective radiative recombination, initialize to zero */ 
    445         iso.RadRec_effec[ipISO][nelem] = 0.; 
    446  
    447         for( ipLevel=0; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel ) 
    448         { 
    449                 /* option for case b conditions, kill ground state recombination */ 
    450                 if( opac.lgCaseB && ipLevel==0 ) 
    451                 { 
    452                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 1e-10; 
    453                 } 
    454                 else 
    455                 { 
    456                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] =  
    457                                 RT_recom_effic(iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]); 
    458                 } 
    459  
    460                 /* net escape prob includes dest by background opacity */ 
    461                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] =  
    462                         MIN2( 1., 
    463                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] + 
    464                         (1.-iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc])* 
    465                         iso.ConOpacRatio[ipISO][nelem][ipLevel] ); 
    466  
    467                 /** \todo 2 remove this entirely? */ 
    468                 if( ipISO == ipHE_LIKE && opac.lgCaseB && ipLevel==0 ) 
    469                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 1e-10; 
    470  
    471                 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] >= 0. ); 
    472                 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] >= 0. ); 
    473                 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] <= 1. ); 
    474  
    475                 /* sum of all effective rad rec */ 
    476                 iso.RadRec_effec[ipISO][nelem] += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]* 
    477                   iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]; 
    478         } 
    479  
    480         /* zero out escape probabilities of levels above numLevels_local */ 
    481         for( ipLevel=iso.numLevels_local[ipISO][nelem]; ipLevel<iso.numLevels_max[ipISO][nelem]; ++ipLevel ) 
    482         { 
    483                 /* this is escape probability */ 
    484                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 0.;  
    485                 /* net escape prob includes dest by background opacity */ 
    486                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 0.; 
    487         } 
    488  
    489         /* trace escape probabilities */ 
    490         if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) ) 
    491         { 
    492                 fprintf( ioQQQ, "       iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem ); 
    493  
    494                 /* print continuum escape probability */ 
    495                 fprintf( ioQQQ, "       iso_radiative_recomb recomb effic" ); 
    496                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
    497                 { 
    498                         fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] )); 
    499                 } 
    500                 fprintf( ioQQQ, "\n" ); 
    501  
    502                 /* net recombination efficiency factor, including background opacity*/ 
    503                 fprintf( ioQQQ, "       iso_radiative_recomb recomb net effic" ); 
    504                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
    505                 { 
    506                         fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]) ); 
    507                 } 
    508                 fprintf( ioQQQ, "\n" ); 
    509  
    510                 /* inward continuous optical depths */ 
    511                 fprintf( ioQQQ, "       iso_radiative_recomb in optic dep" ); 
    512                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
    513                 {        
    514                         if( ipISO==ipH_LIKE  && ipLevel==1 ) 
    515                                 continue; 
    516                         fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] )); 
    517                 } 
    518                 fprintf( ioQQQ, "\n" ); 
    519  
    520                 /* outward continuous optical depths*/ 
    521                 fprintf( ioQQQ, "       iso_radiative_recomb out op depth" ); 
    522                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
    523                 { 
    524                         if( ipISO==ipH_LIKE  && ipLevel==1 ) 
    525                                 continue; 
    526                         fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] )); 
    527                 } 
    528                 fprintf( ioQQQ, "\n" ); 
    529  
    530                 /* print radiative recombination coefficients */ 
    531                 fprintf( ioQQQ, "       iso_radiative_recomb rad rec coef " ); 
    532                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
    533                 { 
    534                         fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]) ); 
    535                 } 
    536                 fprintf( ioQQQ, "\n" ); 
    537         } 
    538  
    539         if( trace.lgTrace ) 
    540         { 
    541                 fprintf( ioQQQ, "     iso_radiative_recomb total rec coef" ); 
    542                 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_effec[ipISO][nelem] )); 
    543                 fprintf( ioQQQ, " case A=" ); 
    544                 fprintf( ioQQQ,PrintEfmt("%10.3e",  
    545                         iso.RadRec_caseB[ipISO][nelem] + iso.RadRecomb[ipISO][nelem][ipH1s][ipRecRad] ) ); 
    546                 fprintf( ioQQQ, " case B="); 
    547                 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_caseB[ipISO][nelem] )); 
    548                 fprintf( ioQQQ, "\n" ); 
    549         } 
    550  
    551         /****************************/ 
    552         /***  begin sanity check  ***/ 
    553         /****************************/ 
    554         { 
    555                 bool lgOK = true; 
    556                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 
    557                 { 
    558                         if( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] <= 0. ) 
    559                         { 
    560                                 fprintf( ioQQQ,  
    561                                         " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",  
    562                                   ipISO, nelem, ipLevel, iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] , phycon.te); 
    563                                         lgOK = false; 
    564                         } 
    565                 } 
    566                 /* bail if we found problems */ 
    567                 if( !lgOK ) 
    568                 { 
    569                         ShowMe(); 
    570                         cdEXIT(EXIT_FAILURE); 
    571                 } 
    572                 /*end sanity check */ 
    573         } 
    574  
    575         /* confirm that we have good rec coef at bottom and top of atom/ion */ 
    576         ASSERT( iso.RadRecomb[ipISO][nelem][0][ipRecRad] > 0. ); 
    577         ASSERT( iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecRad] > 0. ); 
    578  
    579         /* set true when punch recombination coefficients command entered */ 
    580         if( punch.lgioRecom ) 
    581         { 
    582                 if( nelem==ipISO ) 
    583                 { 
    584                         if( ipISO==ipH_LIKE ) 
    585                         { 
    586                                 fprintf( punch.ioRecom, "Hydrogenic  species\nH-like  "); 
    587                         } 
    588                         else if( ipISO==ipHE_LIKE ) 
    589                         { 
    590                                 fprintf( punch.ioRecom, "Helium-like species\nHe-like "); 
    591                         } 
    592                         else 
    593                                 /* Must change this if other iso sequences are added. */ 
    594                                 TotalInsanity(); 
    595                 } 
    596                 /* this prints Z on physical, not C, scale */ 
    597                 fprintf( punch.ioRecom, "%2li %s ",  
    598                         nelem+1 , elementnames.chElementSym[nelem] ); 
    599                 fprintf( punch.ioRecom,PrintEfmt("%9.2e", iso.RadRec_caseB[ipISO][nelem] )); 
    600                 fprintf( punch.ioRecom, "\n" ); 
    601         } 
    602  
    603543        return; 
    604544} 
    605  
    606545/*iso_RRCoef_Te evaluated radiative recombination coef at some temperature */ 
    607546double iso_RRCoef_Te( long ipISO, long nelem , long n ) 
    … …  
    659598} 
    660599 
     600/* malloc space needed for iso recombination tables */ 
     601void iso_recomb_malloc( void ) 
     602{ 
     603        DEBUG_ENTRY( "iso_recomb_malloc()" ); 
     604 
     605        NumLevRecomb = (long **)MALLOC(sizeof(long *)*(unsigned)NISO ); 
     606        TotalRecomb = (double ***)MALLOC(sizeof(double **)*(unsigned)NISO ); 
     607        RRCoef = (double ****)MALLOC(sizeof(double ***)*(unsigned)NISO ); 
     608 
     609        for( long ipISO=0; ipISO<NISO; ipISO++ ) 
     610        { 
     611                TotalRecomb[ipISO] = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 
     612                RRCoef[ipISO] = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM ); 
     613                /* The number of recombination coefficients to be read from file for each element.      */ 
     614                NumLevRecomb[ipISO] = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM ); 
     615 
     616                for( long nelem=ipISO; nelem < LIMELM; ++nelem ) 
     617                { 
     618                        long int MaxLevels, maxN; 
     619 
     620                        TotalRecomb[ipISO][nelem] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB ); 
     621 
     622                        if( nelem == ipISO ) 
     623                                maxN = RREC_MAXN; 
     624                        else 
     625                                maxN = LIKE_RREC_MAXN( nelem ); 
     626 
     627                        NumLevRecomb[ipISO][nelem] = iso_get_total_num_levels( ipISO, maxN, 0 ); 
     628 
     629                        if( nelem == ipISO || dense.lgElmtOn[nelem] ) 
     630                        { 
     631                                /* must always have at least NumLevRecomb[ipISO][nelem] levels since that is number  
     632                                * that will be read in from he rec data file, but possible to require more */ 
     633                                MaxLevels = MAX2( NumLevRecomb[ipISO][nelem] , iso.numLevels_max[ipISO][nelem] ); 
     634 
     635                                /* always define this */ 
     636                                /* >>chng 02 jan 24, RRCoef will be iso.numLevels_max[ipISO][nelem] levels, not iso.numLevels_max, 
     637                                * code will stop if more than this is requested */ 
     638                                RRCoef[ipISO][nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MaxLevels) ); 
     639 
     640                                for( long ipLo=0; ipLo < MaxLevels;++ipLo ) 
     641                                { 
     642                                        RRCoef[ipISO][nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB ); 
     643                                } 
     644                        } 
     645                } 
     646        } 
     647 
     648        for(long i = 0; i < N_ISO_TE_RECOMB; i++) 
     649        { 
     650                /* this is the vector of temperatures */ 
     651                TeRRCoef[i] = 0.25*(i); 
     652        } 
     653 
     654        /* >>chng 06 jun 06, NP, assert thrown at T == 1e10 K, just bump the  
     655         * high temperature end slightly.  */ 
     656        TeRRCoef[N_ISO_TE_RECOMB-1] += 0.01f; 
     657 
     658        return; 
     659} 
     660 
     661void iso_recomb_auxiliary_free( void ) 
     662{ 
     663        DEBUG_ENTRY( "iso_recomb_auxiliary_free()" ); 
     664 
     665        for( long i = 0; i< NISO; i++ ) 
     666        { 
     667                free( NumLevRecomb[i] ); 
     668        } 
     669        free( NumLevRecomb ); 
     670 
     671        return; 
     672} 
     673 
    661674void iso_recomb_setup( long ipISO ) 
    662675{ 
    … …  
    664677        long int i, i1, i2, i3, i4, i5; 
    665678        long int ipLo, nelem; 
    666         static bool lgFirstCall = true, lgSecondCall = false; 
    667679