cloudy: /home66/gary/public_html/cloudy/c08_branch/source/iso_radiative

00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*iso_radiative_recomb find state specific creation and destruction rates for iso sequences */
00004 /*iso_RRCoef_Te - evaluate radiative recombination coef at some temperature */
00005 /*iso_recomb_check - called by SanityCheck to confirm that recombination coef are ok,
00006  * return value is relative error between new calculation of recom, and interp value */
00007 
00008 #include "cddefines.h"
00009 #include "atmdat.h"
00010 #include "conv.h"
00011 #include "dense.h"
00012 #include "elementnames.h"
00013 #include "helike.h"
00014 #include "helike_recom.h"
00015 #include "hydrogenic.h"
00016 #include "ionbal.h"
00017 #include "iso.h"
00018 #include "opacity.h"
00019 #include "phycon.h"
00020 #include "physconst.h" 
00021 #include "prt.h"
00022 #include "punch.h"
00023 #include "thermal.h"
00024 #include "thirdparty.h"
00025 #include "trace.h"
00026 #include "rt.h"
00027 
00028 /* this will save log of radiative recombination rate coefficients at N_ISO_TE_RECOMB temperatures.
00029  * there will be NumLevRecomb[ipISO][nelem] levels saved in RRCoef[nelem][level][temp] */
00030 static double ****RRCoef/*[LIMELM][NumLevRecomb[ipISO][nelem]][N_ISO_TE_RECOMB]*/;
00031 static long **NumLevRecomb;
00032 static double ***TotalRecomb;   /*[ipISO][nelem][i]*/
00033 
00034 /* the array of logs of temperatures at which RRCoef was defined */
00035 static double TeRRCoef[N_ISO_TE_RECOMB];
00036 
00037 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements );
00038 
00039 double iso_radrecomb_from_cross_section(long ipISO, double temp, long nelem, long ipLo)
00040 {
00041         double alpha;
00042 
00043         DEBUG_ENTRY( "iso_radrecomb_from_cross_section()" );
00044 
00045         if( ipISO == ipH_LIKE )
00046                 alpha =  hlike_radrecomb_from_cross_section( temp, nelem, ipLo);
00047         else if( ipISO == ipHE_LIKE )
00048                 alpha =  helike_radrecomb_from_cross_section( temp, nelem, ipLo);
00049         else
00050                 TotalInsanity();
00051 
00052         return alpha;
00053 }
00054 
00055 /*=======================================================*/
00056 /* iso_radiative_recomb get rad recomb rate coefficients for iso sequences */
00057 void iso_radiative_recomb( 
00058                                                   long ipISO,
00059                                                   /* nelem on the c scale, He is 1 */
00060                                                   long nelem )
00061 {
00062         long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel; 
00063         double topoff, TotMinusExplicitResolved,
00064                 TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.;
00065         double RecExplictLevels, TotalRadRecomb, RecCollapsed;
00066         static double TeUsed[NISO][LIMELM]={
00067                 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00068                  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00069                  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00070                 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00071                  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00072                  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
00073 
00074         DEBUG_ENTRY( "iso_radiative_recomb()" );
00075 
00076         /* evaluate recombination escape probability for all levels */
00077 
00078         /* define radiative recombination rates for all levels */ 
00079         /* this will be the sum of all levels explicitly included in the model atom */
00080         RecExplictLevels = 0.;
00081 
00082         /* number of resolved levels, so first collapsed level is [ipFirstCollapsed] */
00083         ipFirstCollapsed = iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem]; 
00084 
00085         if( (fabs(1.-TeUsed[ipISO][nelem]/phycon.te)> 0.001) || hydro.lgReevalRecom || !conv.nTotalIoniz )
00086         {
00087                 TeUsed[ipISO][nelem] = phycon.te;
00088 
00089                 for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel )
00090                 {
00091                         /* this is radiative recombination rate coefficient */
00092                         double RadRec;
00093 
00094                         if( !iso.lgNoRecombInterp[ipISO] )
00095                         {
00096                                 RadRec = iso_RRCoef_Te( ipISO, nelem , ipLevel );
00097                         }
00098                         else
00099                         {
00100                                 RadRec = iso_radrecomb_from_cross_section( ipISO, phycon.te, nelem, ipLevel);
00101                         }
00102 
00103                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] = RadRec;
00104 
00105                         if( iso.lgRandErrGen[ipISO] )
00106                         {
00107                                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *=
00108                                         /* this has to be from iso.numLevels_max instead of iso.numLevels_local because
00109                                          * the error factors for rrc are always stored at iso.numLevels_max, regardless of
00110                                          * continuum lowering effects. */
00111                                         iso.ErrorFactor[ipISO][nelem][ iso.numLevels_max[ipISO][nelem] ][ipLevel][IPRAD];
00112                         }
00113 
00114                         ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. );
00115 
00116                         RecExplictLevels += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
00117 
00118                         if( iso.lgDielRecom[ipISO] )
00119                         {
00120                                 /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */
00121                                 iso.DielecRecomb[ipISO][nelem][ipLevel] = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
00122                                 Total_DR_Added  += iso.DielecRecomb[ipISO][nelem][ipLevel];
00123                         }
00124                 }
00125 
00126                 /**************************************************/
00127                 /***  Add on recombination to collapsed levels  ***/
00128                 /**************************************************/
00129                 RecCollapsed = 0.;
00130                 for( ipLevel=ipFirstCollapsed; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel )
00131                 {
00132                         /* use hydrogenic for collapsed levels */
00133                         double RadRec = t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, N_(ipLevel), phycon.te);
00134 
00135                         /* this is radiative recombination rate coefficient */
00136                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] = RadRec;
00137 
00138                         /* This kills recombination into the collapsed level so that the forced
00139                          * statistical weighting can be bypassed */
00140                         if( !iso.lgTopoff[ipISO] )
00141                         {
00142                                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *= 1E-10;
00143                         }
00144 
00145                         if( iso.lgRandErrGen[ipISO] )
00146                         {
00147                                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *= 
00148                                         /* this has to be from iso.numLevels_max instead of iso.numLevels_local because
00149                                          * the error factors for rrc are always stored at iso.numLevels_max, regardless of
00150                                          * continuum lowering effects. */
00151                                         iso.ErrorFactor[ipISO][nelem][ iso.numLevels_max[ipISO][nelem] ][ipLevel][IPRAD];
00152                         }
00153 
00154                         RecCollapsed += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
00155 
00156                         ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. );
00157 
00158                         if( iso.lgDielRecom[ipISO] )
00159                         {
00160                                 /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */
00161                                 iso.DielecRecomb[ipISO][nelem][ipLevel] = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
00162                                 Total_DR_Added  += iso.DielecRecomb[ipISO][nelem][ipLevel];
00163                         }
00164                 }
00165 
00166                 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
00167                 for( ipLevel = 0; ipLevel<iso.numLevels_local[ipISO][nelem]; ipLevel++ )
00168                 {
00169                         if( N_(ipLevel) == ThisN )
00170                         {
00171                                 TotRRThisN += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
00172                         }
00173                         else
00174                         {
00175                                 ASSERT( (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) );
00176                                 LastN = ThisN;
00177                                 ThisN = N_(ipLevel);
00178                                 TotRRLastN = TotRRThisN;
00179                                 TotRRThisN = iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
00180 
00181                                 {
00182                                         /* Print the sum of recombination coefficients per n at current temp.   */
00183                                         /*@-redef@*/
00184                                         enum {DEBUG_LOC=false};
00185                                         /*@+redef@*/
00186                                         static long RUNONCE = false;
00187 
00188                                         if( !RUNONCE && DEBUG_LOC )
00189                                         {
00190                                                 static long FIRSTTIME = true;
00191 
00192                                                 if( FIRSTTIME )
00193                                                 {
00194                                                         fprintf( ioQQQ,"Sum of Radiative Recombination at current iso, nelem, temp = %li %li %.2f\n", 
00195                                                                 ipISO, nelem, phycon.te);
00196                                                         FIRSTTIME= false;
00197                                                 }
00198 
00199                                                 fprintf( ioQQQ,"%li\t%.2e\n",LastN,TotRRLastN );
00200                                         }
00201                                         RUNONCE = true;
00202                                 }
00203                         }
00204                 }
00205 
00206                 /* Get total recombination into all levels, including those not explicitly considered. */
00207                 if( iso.lgNoRecombInterp[ipISO] )
00208                 {
00209                         /* We are not interpolating, must calculate total now, as sum of resolved and collapsed levels... */
00210                         TotalRadRecomb = RecCollapsed + RecExplictLevels;
00211 
00212                         /* Plus additional levels out to a predefined limit... */
00213                         for( long nLo = N_(ipFirstCollapsed) + iso.nCollapsed_max[ipISO][nelem]; nLo < NHYDRO_MAX_LEVEL; nLo++ )
00214                         {
00215                                 TotalRadRecomb += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, nLo, phycon.te);
00216                         }
00217                         /* Plus a bunch more for good measure */
00218                         for( long nLo = NHYDRO_MAX_LEVEL; nLo<=SumUpToThisN; nLo++ )
00219                         {
00220                                 TotalRadRecomb += Recomb_Seaton59( nelem+1-ipISO, phycon.te, nLo );
00221                         }
00222                 }
00223                 else
00224                 {
00225                         /* We are interpolating, and total was calculated here in iso_recomb_setup */
00226                         TotalRadRecomb = iso_RRCoef_Te( ipISO, nelem, iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] );
00227                 }
00228 
00229                 /* If generating random error, apply one to total recombination */
00230                 if( iso.lgRandErrGen[ipISO] )
00231                 {
00232                         /* this has to be from iso.numLevels_max instead of iso.numLevels_local because
00233                          * the error factors for rrc are always stored at iso.numLevels_max, regardless of
00234                          * continuum lowering effects. */
00235                         TotalRadRecomb *= 
00236                                 iso.ErrorFactor[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][iso.numLevels_max[ipISO][nelem]][IPRAD];
00237                 }
00238 
00239                 /* this is case B recombination, sum without the ground included */
00240                 iso.RadRec_caseB[ipISO][nelem] = TotalRadRecomb - iso.RadRecomb[ipISO][nelem][0][ipRecRad];
00241                 ASSERT( iso.RadRec_caseB[ipISO][nelem] > 0.);
00242 
00243                 /**********************************************************************/
00244                 /***  Add topoff (excess) recombination to top level.  This is only     ***/
00245                 /***  done if atom is not full size.                                ***/
00246                 /**********************************************************************/
00247                 if( !iso.lgLevelsLowered[ipISO][nelem] )
00248                 {
00249                         /* at this point we have RecExplictLevels, the sum of radiative recombination 
00250                          * to all levels explicitly included in the model atom the total 
00251                          * recombination rate.  The difference is the amount of "topoff" we will need to do */
00252                         TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels;
00253 
00254                         topoff = TotMinusExplicitResolved - RecCollapsed;
00255 
00256                         /* the t_ADfA::Inst().rad_rec fits are too small at high temperatures, so this atom is
00257                          * better than the topoff.  Only a problem for helium itself, at high temperatures.
00258                          * complain if the negative topoff is not for this case */
00259                         if( topoff < 0. && (nelem!=ipHELIUM || phycon.te < 1e5 ) &&
00260                                 fabs( topoff/TotalRadRecomb ) > 0.01 )
00261                         {
00262                                 fprintf(ioQQQ," PROBLEM  negative RR topoff for %li, rel error was %.2e, temp was %.2f\n",  
00263                                         nelem, topoff/TotalRadRecomb, phycon.te );
00264                         }
00265 
00266                         if( !iso.lgTopoff[ipISO] )
00267                                 topoff *= 1E-20;
00268 
00269                         /* We always have at least one collapsed level if continuum is not lowered.  Put topoff there.  */
00270                         iso.RadRecomb[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1][ipRecRad] += MAX2(0., topoff );
00271 
00272                         /* check for negative DR topoff, but only if Total_DR_Added is not negligible compared with TotalRadRecomb */
00273                         if( Total_DR_Added > TotalRadRecomb/100. )
00274                         {
00275                                 if( Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] > 1.02 )
00276                                 {
00277                                         fprintf(ioQQQ," PROBLEM  negative DR topoff for %li, rel error was %.2e, temp was %.2f\n",  
00278                                                 nelem,
00279                                                 Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - 1.0,
00280                                                 phycon.te );
00281                                 }
00282                         }
00283 
00284                         if( iso.lgDielRecom[ipISO] )
00285                         {
00286                                 /* \todo 2 suppress this total rate for continuum lowering using factors from Jordan (1969). */
00287                                 /* put extra DR in top level */
00288                                 iso.DielecRecomb[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] = 
00289                                         ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - Total_DR_Added;
00290                         }
00291                 }
00292 
00293         }
00294 
00295         /**************************************************************/
00296         /***  Stuff escape probabilities, and effective rad recomb  ***/
00297         /**************************************************************/
00298 
00299         /* total effective radiative recombination, initialize to zero */
00300         iso.RadRec_effec[ipISO][nelem] = 0.;
00301 
00302         for( ipLevel=0; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel )
00303         {
00304                 /* option for case b conditions, kill ground state recombination */
00305                 if( opac.lgCaseB && ipLevel==0 )
00306                 {
00307                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 1e-10;
00308                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 1e-10;
00309                 }
00310                 else
00311                 {
00312                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 
00313                                 RT_recom_effic(iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]);
00314 
00315                         /* net escape prob includes dest by background opacity */
00316                         iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 
00317                                 MIN2( 1.,
00318                                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] +
00319                                 (1.-iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc])*
00320                                 iso.ConOpacRatio[ipISO][nelem][ipLevel] );
00321                 }
00322 
00323                 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] >= 0. );
00324                 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] >= 0. );
00325                 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] <= 1. );
00326 
00327                 /* sum of all effective rad rec */
00328                 iso.RadRec_effec[ipISO][nelem] += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]*
00329                   iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc];
00330         }
00331 
00332         /* zero out escape probabilities of levels above numLevels_local */
00333         for( ipLevel=iso.numLevels_local[ipISO][nelem]; ipLevel<iso.numLevels_max[ipISO][nelem]; ++ipLevel )
00334         {
00335                 /* this is escape probability */
00336                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 0.; 
00337                 /* net escape prob includes dest by background opacity */
00338                 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 0.;
00339         }
00340 
00341         /* trace escape probabilities */
00342         if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
00343         {
00344                 fprintf( ioQQQ, "       iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem );
00345 
00346                 /* print continuum escape probability */
00347                 fprintf( ioQQQ, "       iso_radiative_recomb recomb effic" );
00348                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
00349                 {
00350                         fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] ));
00351                 }
00352                 fprintf( ioQQQ, "\n" );
00353 
00354                 /* net recombination efficiency factor, including background opacity*/
00355                 fprintf( ioQQQ, "       iso_radiative_recomb recomb net effic" );
00356                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
00357                 {
00358                         fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]) );
00359                 }
00360                 fprintf( ioQQQ, "\n" );
00361 
00362                 /* inward continuous optical depths */
00363                 fprintf( ioQQQ, "       iso_radiative_recomb in optic dep" );
00364                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
00365                 {       
00366                         fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] ));
00367                 }
00368                 fprintf( ioQQQ, "\n" );
00369 
00370                 /* outward continuous optical depths*/
00371                 fprintf( ioQQQ, "       iso_radiative_recomb out op depth" );
00372                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
00373                 {
00374                         fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] ));
00375                 }
00376                 fprintf( ioQQQ, "\n" );
00377 
00378                 /* print radiative recombination coefficients */
00379                 fprintf( ioQQQ, "       iso_radiative_recomb rad rec coef " );
00380                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
00381                 {
00382                         fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]) );
00383                 }
00384                 fprintf( ioQQQ, "\n" );
00385         }
00386 
00387         if( trace.lgTrace )
00388         {
00389                 fprintf( ioQQQ, "     iso_radiative_recomb total rec coef" );
00390                 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_effec[ipISO][nelem] ));
00391                 fprintf( ioQQQ, " case A=" );
00392                 fprintf( ioQQQ,PrintEfmt("%10.3e", 
00393                         iso.RadRec_caseB[ipISO][nelem] + iso.RadRecomb[ipISO][nelem][ipH1s][ipRecRad] ) );
00394                 fprintf( ioQQQ, " case B=");
00395                 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_caseB[ipISO][nelem] ));
00396                 fprintf( ioQQQ, "\n" );
00397         }
00398 
00399         /****************************/
00400         /***  begin sanity check  ***/
00401         /****************************/
00402         {
00403                 bool lgOK = true;
00404                 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
00405                 {
00406                         if( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] <= 0. )
00407                         {
00408                                 fprintf( ioQQQ, 
00409                                         " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n", 
00410                                   ipISO, nelem, ipLevel, iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] , phycon.te);
00411                                         lgOK = false;
00412                         }
00413                 }
00414                 /* bail if we found problems */
00415                 if( !lgOK )
00416                 {
00417                         ShowMe();
00418                         cdEXIT(EXIT_FAILURE);
00419                 }
00420                 /*end sanity check */
00421         }
00422 
00423         /* confirm that we have good rec coef at bottom and top of atom/ion */
00424         ASSERT( iso.RadRecomb[ipISO][nelem][0][ipRecRad] > 0. );
00425         ASSERT( iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecRad] > 0. );
00426 
00427         /* set true when punch recombination coefficients command entered */
00428         if( punch.lgioRecom )
00429         {
00430                 /* this prints Z on physical, not C, scale */
00431                 fprintf( punch.ioRecom, "%s %s %2li ", 
00432                         iso.chISO[ipISO], elementnames.chElementSym[nelem], nelem+1 );
00433                 fprintf( punch.ioRecom,PrintEfmt("%9.2e", iso.RadRec_caseB[ipISO][nelem] ));
00434                 fprintf( punch.ioRecom, "\n" );
00435         }
00436 
00437         return;
00438 }
00439 
00440 void iso_radiative_recomb_effective( long ipISO, long nelem )
00441 {
00442         DEBUG_ENTRY( "iso_radiative_recomb_effective()" );
00443 
00444         /* Find effective recombination coefficients */
00445         for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
00446         {
00447                 iso.RadEffec[ipISO][nelem][ipHi] = 0.;
00448 
00449                 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
00450                 for( long ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ )
00451                 {
00452                         ASSERT( iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] >= 0. );
00453                         ASSERT( iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad] >= 0. );
00454 
00455                         iso.RadEffec[ipISO][nelem][ipHi] += iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] *
00456                                 iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad];
00457                 }
00458         }
00459 
00460         /**************************************************************/
00461         /***  option to print effective recombination coefficients  ***/
00462         /**************************************************************/
00463         {
00464                 enum {DEBUG_LOC=false};
00465 
00466                 if( DEBUG_LOC )
00467                 {
00468                         const int maxPrt=10;
00469 
00470                         fprintf( ioQQQ,"Effective recombination, ipISO=%li, nelem=%li, Te = %e\n", ipISO, nelem, phycon.te );
00471                         fprintf( ioQQQ, "N\tL\tS\tRadEffec\tLifetime\n" );
00472 
00473                         for( long i=0; i<maxPrt; i++ )
00474                         {
00475                                 fprintf( ioQQQ, "%li\t%li\t%li\t%e\t%e\n", N_(i), L_(i), S_(i),
00476                                         iso.RadEffec[ipISO][nelem][i],
00477                                         MAX2( 0., StatesElem[ipISO][nelem][i].lifetime ) );
00478                         }
00479                         fprintf( ioQQQ, "\n" );
00480                 }
00481         }
00482 
00483         /* If we have the variable set, find errors in rad rates */
00484         if( iso.lgRandErrGen[ipISO] )
00485         {
00486                 dprintf( ioQQQ, "ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" );
00487 
00488                 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
00489                 for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
00490                 {
00491                         iso.SigmaRadEffec[ipISO][nelem][ipHi] = 0.;
00492 
00493                         /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
00494                         for( long ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ )
00495                         {
00496                                 ASSERT( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][ipHigher][IPRAD] >= 0. );
00497                                 ASSERT( iso.SigmaCascadeProb[ipISO][nelem][ipHigher][ipHi] >= 0. );
00498 
00499                                 /* iso.RadRecomb has to appear here because iso.Error is only relative error */ 
00500                                 iso.SigmaRadEffec[ipISO][nelem][ipHi] += pow( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][ipHigher][IPRAD] *
00501                                         iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] * iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad], 2.) + 
00502                                         pow( iso.SigmaCascadeProb[ipISO][nelem][ipHigher][ipHi] * iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad], 2.);
00503                         }
00504 
00505                         ASSERT( iso.SigmaRadEffec[ipISO][nelem][ipHi] >= 0. );
00506                         iso.SigmaRadEffec[ipISO][nelem][ipHi] = sqrt( iso.SigmaRadEffec[ipISO][nelem][ipHi] );
00507 
00508                         for( long ipLo = 0; ipLo < ipHi; ipLo++ )
00509                         {
00510                                 if( (( L_(ipLo) == L_(ipHi) + 1 ) || ( L_(ipHi) == L_(ipLo) + 1 )) )
00511                                 {       
00512                                         double EnergyInRydbergs = iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi];
00513                                         realnum wavelength = (realnum)(RYDLAM/MAX2(1E-8,EnergyInRydbergs));
00514                                         double emissivity = iso.RadEffec[ipISO][nelem][ipHi] * iso.BranchRatio[ipISO][nelem][ipHi][ipLo] * EN1RYD * EnergyInRydbergs;
00515                                         double sigma_emiss = 0., SigmaBranchRatio = 0.;
00516 
00517                                         if( ( emissivity > 2.E-29 ) && ( wavelength < 1.E6 ) && (N_(ipHi)<=5) )
00518                                         {
00519                                                 SigmaBranchRatio = iso.BranchRatio[ipISO][nelem][ipHi][ipLo] * sqrt( 
00520                                                         pow( (double)iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD], 2. ) +
00521                                                         pow( iso.SigmaAtot[ipISO][nelem][ipHi]*StatesElem[ipISO][nelem][ipHi].lifetime, 2.) );
00522 
00523                                                 sigma_emiss =  EN1RYD * EnergyInRydbergs * sqrt( 
00524                                                         pow( (double)iso.SigmaRadEffec[ipISO][nelem][ipHi] * iso.BranchRatio[ipISO][nelem][ipHi][ipLo], 2.) +
00525                                                         pow( SigmaBranchRatio * iso.RadEffec[ipISO][nelem][ipHi], 2.) );
00526 
00527                                                 /* \todo 2 make this a trace option */
00528                                                 dprintf( ioQQQ, "%li\t%li\t", ipHi, ipLo );
00529                                                 prt_wl( ioQQQ, wavelength );
00530                                                 fprintf( ioQQQ, "\t%e\t%e\t%e\t%e\t%e\t%e\n", 
00531                                                         emissivity,
00532                                                         sigma_emiss,
00533                                                         iso.RadEffec[ipISO][nelem][ipHi],
00534                                                         iso.SigmaRadEffec[ipISO][nelem][ipHi],
00535                                                         iso.BranchRatio[ipISO][nelem][ipHi][ipLo],
00536                                                         SigmaBranchRatio);
00537                                         }
00538                                 }
00539                         }
00540                 }
00541         }
00542 
00543         return;
00544 }
00545 /*iso_RRCoef_Te evaluated radiative recombination coef at some temperature */
00546 double iso_RRCoef_Te( long ipISO, long nelem , long n )
00547 {
00548         double rate = 0.;
00549 
00550         DEBUG_ENTRY( "iso_RRCoef_Te()" );
00551 
00552         ASSERT( !iso.lgNoRecombInterp[ipISO] );
00553 
00554         /* if n is equal to the number of levels, return the total recomb, else recomb for given level. */
00555         if( n == iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] )
00556         {
00557                 rate = TempInterp( TeRRCoef, TotalRecomb[ipISO][nelem], N_ISO_TE_RECOMB );
00558         }
00559         else
00560         {
00561                 rate = TempInterp( TeRRCoef, RRCoef[ipISO][nelem][n], N_ISO_TE_RECOMB );
00562         }
00563 
00564         /* that was the log, now make linear */
00565         rate = pow( 10. , rate );
00566 
00567         return rate;
00568 }
00569 
00570 /*iso_recomb_check called by SanityCheck to confirm that recombination coef are ok
00571  * return value is relative error between new calculation of recom, and interp value */
00572 double iso_recomb_check( long ipISO, long nelem, long level, double temperature )
00573 {
00574         double RecombRelError ,
00575                 RecombInterp,
00576                 RecombCalc,
00577                 SaveTemp;
00578 
00579         DEBUG_ENTRY( "iso_recomb_check()" );
00580 
00581         /* first set temp to desired value */
00582         SaveTemp = phycon.te;
00583         TempChange(temperature);
00584 
00585         /* actually calculate the recombination coefficient from the Milne relation,
00586          * normally only done due to compile he-like command */
00587         RecombCalc = iso_radrecomb_from_cross_section( ipISO, temperature , nelem , level );
00588 
00589         /* interpolate the recombination coefficient, this is the usual method */
00590         RecombInterp = iso_RRCoef_Te( ipISO, nelem , level );
00591 
00592         /* reset temp */
00593         TempChange(SaveTemp);
00594 
00595         RecombRelError = ( RecombInterp - RecombCalc ) / MAX2( RecombInterp , RecombCalc );
00596 
00597         return RecombRelError;
00598 }
00599 
00600 /* malloc space needed for iso recombination tables */
00601 void iso_recomb_malloc( void )
00602 {
00603         DEBUG_ENTRY( "iso_recomb_malloc()" );
00604 
00605         NumLevRecomb = (long **)MALLOC(sizeof(long *)*(unsigned)NISO );
00606         TotalRecomb = (double ***)MALLOC(sizeof(double **)*(unsigned)NISO );
00607         RRCoef = (double ****)MALLOC(sizeof(double ***)*(unsigned)NISO );
00608 
00609         for( long ipISO=0; ipISO<NISO; ipISO++ )
00610         {
00611                 TotalRecomb[ipISO] = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
00612                 RRCoef[ipISO] = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
00613                 /* The number of recombination coefficients to be read from file for each element.      */
00614                 NumLevRecomb[ipISO] = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM );
00615 
00616                 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00617                 {
00618                         long int MaxLevels, maxN;
00619 
00620                         TotalRecomb[ipISO][nelem] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
00621 
00622                         if( nelem == ipISO )
00623                                 maxN = RREC_MAXN;
00624                         else
00625                                 maxN = LIKE_RREC_MAXN( nelem );
00626 
00627                         NumLevRecomb[ipISO][nelem] = iso_get_total_num_levels( ipISO, maxN, 0 );
00628 
00629                         if( nelem == ipISO || dense.lgElmtOn[nelem] )
00630                         {
00631                                 /* must always have at least NumLevRecomb[ipISO][nelem] levels since that is number 
00632                                 * that will be read in from he rec data file, but possible to require more */
00633                                 MaxLevels = MAX2( NumLevRecomb[ipISO][nelem] , iso.numLevels_max[ipISO][nelem] );
00634 
00635                                 /* always define this */
00636                                 /* >>chng 02 jan 24, RRCoef will be iso.numLevels_max[ipISO][nelem] levels, not iso.numLevels_max,
00637                                 * code will stop if more than this is requested */
00638                                 RRCoef[ipISO][nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MaxLevels) );
00639 
00640                                 for( long ipLo=0; ipLo < MaxLevels;++ipLo )
00641                                 {
00642                                         RRCoef[ipISO][nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
00643                                 }
00644                         }
00645                 }
00646         }
00647 
00648         for(long i = 0; i < N_ISO_TE_RECOMB; i++)
00649         {
00650                 /* this is the vector of temperatures */
00651                 TeRRCoef[i] = 0.25*(i);
00652         }
00653 
00654         /* >>chng 06 jun 06, NP, assert thrown at T == 1e10 K, just bump the 
00655          * high temperature end slightly.  */
00656         TeRRCoef[N_ISO_TE_RECOMB-1] += 0.01f;
00657 
00658         return;
00659 }
00660 
00661 void iso_recomb_auxiliary_free( void )
00662 {
00663         DEBUG_ENTRY( "iso_recomb_auxiliary_free()" );
00664 
00665         for( long i = 0; i< NISO; i++ )
00666         {
00667                 free( NumLevRecomb[i] );
00668         }
00669         free( NumLevRecomb );
00670 
00671         return;
00672 }
00673 
00674 void iso_recomb_setup( long ipISO )
00675 {
00676         double RadRecombReturn;
00677         long int i, i1, i2, i3, i4, i5;
00678         long int ipLo, nelem;
00679 
00680         const int chLine_LENGTH = 1000;
00681         char chLine[chLine_LENGTH];
00682         /* this must be longer than data path string, set in path.h/cpu.cpp */
00683         const char* chFilename[NISO] = 
00684                 { "h_iso_recomb.dat", "he_iso_recomb.dat" };
00685 
00686         FILE *ioDATA;
00687         bool lgEOL;
00688 
00689         DEBUG_ENTRY( "iso_recomb_setup()" );
00690 
00691         /* if we are compiling the recombination data file, we must interpolate in temperature */
00692         if( iso.lgCompileRecomb[ipISO] )
00693         {
00694                 iso.lgNoRecombInterp[ipISO] = false;
00695         }
00696 
00697         if( !iso.lgNoRecombInterp[ipISO] )
00698         {
00699                 /******************************************************************/
00701                 /******************************************************************/
00702                 /* This flag says we are not compiling the data file    */
00703                 if( !iso.lgCompileRecomb[ipISO] )
00704                 {
00705                         if( trace.lgTrace )
00706                                 fprintf( ioQQQ," iso_recomb_setup opening %s:", chFilename[ipISO] );
00707 
00708                         /* Now try to read from file...*/
00709                         try
00710                         {
00711                                 ioDATA = open_data( chFilename[ipISO], "r" );
00712                         }
00713                         catch( cloudy_exit )
00714                         {
00715                                 fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
00716                                 fprintf( ioQQQ, " but the code runs much faster if you compile %s!\n", chFilename[ipISO]);
00717                                 ioDATA = NULL;
00718                         }
00719                         if( ioDATA == NULL )
00720                         {
00721                                 /* Do on the fly computation of R.R. Coef's instead.    */
00722                                 for( nelem = ipISO; nelem < LIMELM; nelem++ )
00723                                 {
00724                                         if( dense.lgElmtOn[nelem] )
00725                                         {
00726                                                 /* Zero out the recombination sum array.        */
00727                                                 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00728                                                 {
00729                                                         TotalRecomb[ipISO][nelem][i] = 0.;
00730                                                 }
00731 
00732                                                 /* NumLevRecomb[ipISO][nelem] corresponds to n = 40 for H and He and 20 for ions, at present    */
00733                                                 /* There is no need to fill in values for collapsed levels, because we do not need to
00734                                                 * interpolate for a given temperature, just calculate it directly with a hydrogenic routine.    */
00735                                                 for( ipLo=0; ipLo < iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipLo++ )
00736                                                 {
00737                                                         /* loop over temperatures to produce array of recombination coefficients        */
00738                                                         for(i = 0; i < N_ISO_TE_RECOMB; i++)
00739                                                         {
00740                                                                 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
00741                                                                 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
00742                                                                 TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
00743                                                                 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
00744                                                         }
00745                                                 }
00746                                                 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00747                                                 {
00748                                                         for( i1 = iso.n_HighestResolved_max[ipISO][nelem]+1; i1< NHYDRO_MAX_LEVEL; i1++ )
00749                                                         {
00750                                                                 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
00751                                                         }
00752                                                         for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
00753                                                         {
00754                                                                 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
00755                                                         }
00756                                                         TotalRecomb[ipISO][nelem][i] = log10( TotalRecomb[ipISO][nelem][i] );
00757                                                 }
00758                                         }
00759                                 }
00760                         }
00761                         /* Data file is present and readable...begin read.      */
00762                         else 
00763                         {
00764                                 /* check that magic number is ok */
00765                                 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00766                                 {
00767                                         fprintf( ioQQQ, " iso_recomb_setup could not read first line of %s.\n", chFilename[ipISO]);
00768                                         cdEXIT(EXIT_FAILURE);
00769                                 }
00770                                 i = 1;
00771                                 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00772                                 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00773                                 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00774                                 i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00775                                 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
00776                                 {
00777                                         fprintf( ioQQQ, 
00778                                                 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
00779                                         fprintf( ioQQQ, 
00780                                                 " iso_recomb_setup: I expected to find the numbers  %i %li %li %i and got %li %li %li %li instead.\n" ,
00781                                                 RECOMBMAGIC ,
00782                                                 NumLevRecomb[ipISO][ipISO],
00783                                                 NumLevRecomb[ipISO][ipISO+1],
00784                                                 N_ISO_TE_RECOMB,
00785                                                 i1 , i2 , i3, i4 );
00786                                         fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00787                                         fprintf( ioQQQ, 
00788                                                 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
00789                                         cdEXIT(EXIT_FAILURE);
00790                                 }
00791 
00792                                 i5 = 1;
00793                                 /* now read in the data */
00794                                 for( nelem = ipISO; nelem < LIMELM; nelem++ )
00795                                 {
00796                                         for( ipLo=0; ipLo <= NumLevRecomb[ipISO][nelem]; ipLo++ )
00797                                         {
00798                                                 i5++;
00799                                                 /* get next line image */
00800                                                 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00801                                                 {
00802                                                         fprintf( ioQQQ, " iso_recomb_setup could not read line %li of %s.\n", i5,
00803                                                                 chFilename[ipISO] );
00804                                                         cdEXIT(EXIT_FAILURE);
00805                                                 }
00806                                                 /* each line starts with element and level number */
00807                                                 i3 = 1;
00808                                                 i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
00809                                                 i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
00810                                                 /* check that these numbers are correct */
00811                                                 if( i1!=nelem || i2!=ipLo )
00812                                                 {
00813                                                         fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
00814                                                         fprintf( ioQQQ, 
00815                                                                 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
00816                                                         cdEXIT(EXIT_FAILURE);
00817                                                 }
00818 
00819                                                 /* loop over temperatures to produce array of recombination coefficients        */
00820                                                 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00821                                                 {
00822                                                         double ThisCoef = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
00823 
00824                                                         if( nelem == ipISO || dense.lgElmtOn[nelem] )
00825                                                         {
00826                                                                 /* The last line for each element is the total recombination for each temp.     */
00827                                                                 if( ipLo == NumLevRecomb[ipISO][nelem] )
00828                                                                 {
00829                                                                         TotalRecomb[ipISO][nelem][i] = ThisCoef;
00830                                                                 }
00831                                                                 else
00832                                                                         RRCoef[ipISO][nelem][ipLo][i] = ThisCoef;
00833                                                         }
00834 
00835                                                         if( lgEOL )
00836                                                         {
00837                                                                 fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
00838                                                                 fprintf( ioQQQ, 
00839                                                                         " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
00840                                                                 cdEXIT(EXIT_FAILURE);
00841                                                         }
00842                                                 }
00843                                         }
00844 
00845                                         /* following loop only executed if we need more levels than are
00846                                                 * stored in the recom coef data set
00847                                                 * do not do collapsed levels since will use H-like recom there */
00848                                         if( nelem == ipISO || dense.lgElmtOn[nelem] ) 
00849                                         {
00850                                                 for( ipLo=NumLevRecomb[ipISO][nelem]; ipLo < iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipLo++ )
00851                                                 {
00852                                                                 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00853                                                                 {
00854                                                                         /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
00855                                                                         RRCoef[ipISO][nelem][ipLo][i] = log10(iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo));
00856                                                                 }
00857                                                         }
00858                                                 }
00859                                 }
00860 
00861                                 /* check that ending magic number is ok */
00862                                 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00863                                 {
00864                                         fprintf( ioQQQ, " iso_recomb_setup could not read last line of %s.\n", chFilename[ipISO]);
00865                                         cdEXIT(EXIT_FAILURE);
00866                                 }
00867                                 i = 1;
00868                                 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00869                                 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00870                                 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00871                                 i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00872 
00873                                 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
00874                                 {
00875                                         fprintf( ioQQQ, 
00876                                                 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
00877                                         fprintf( ioQQQ, 
00878                                                 " iso_recomb_setup: I expected to find the numbers  %i %li %li %i and got %li %li %li %li instead.\n" ,
00879                                                 RECOMBMAGIC ,
00880                                                 NumLevRecomb[ipISO][ipISO],
00881                                                 NumLevRecomb[ipISO][ipISO+1],
00882                                                 N_ISO_TE_RECOMB,
00883                                                 i1 , i2 , i3, i4 );
00884                                         fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00885                                         fprintf( ioQQQ, 
00886                                                 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
00887                                         cdEXIT(EXIT_FAILURE);
00888                                 }
00889 
00890                                 /* close the data file */
00891                                 fclose( ioDATA );
00892                         }
00893                 }
00894                 /* We are compiling the he_iso_recomb.dat data file.    */
00895                 else if( iso.lgCompileRecomb[ipISO] )
00896                 {
00897                         /* option to create table of recombination coefficients,
00898                                 * executed with the compile he-like command */
00899                         FILE *ioRECOMB;
00900 
00901                         ASSERT( !iso.lgNoRecombInterp[ipISO] );
00902 
00903                         /*RECOMBMAGIC the magic number for the table of recombination coefficients */
00904                         /*NumLevRecomb[ipISO][nelem] the number of levels that will be done */
00905                         /* create recombination coefficients  */
00906                         ioRECOMB = open_data( chFilename[ipISO], "w", AS_LOCAL_ONLY );
00907                         fprintf(ioRECOMB,"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient H-LIke [or HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
00908                                 RECOMBMAGIC ,
00909                                 NumLevRecomb[ipISO][ipISO],
00910                                 NumLevRecomb[ipISO][ipISO+1],
00911                                 N_ISO_TE_RECOMB,
00912                                 iso.chISO[ipISO],
00913                                 NumLevRecomb[ipISO][ipISO],
00914                                 elementnames.chElementSym[ipISO],
00915                                 NumLevRecomb[ipISO][ipISO+1],
00916                                 N_ISO_TE_RECOMB );
00917 
00918                         for( nelem = ipISO; nelem < LIMELM; nelem++ )
00919                         {
00920                                 /* this must pass since compile xx-like command reset numLevels to the macro */
00921                                 ASSERT( NumLevRecomb[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] );
00922 
00923                                 /* Zero out the recombination sum array.        */
00924                                 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00925                                 {
00926                                         TotalRecomb[ipISO][nelem][i] = 0.;
00927                                 }
00928 
00929                                 for( ipLo=ipHe1s1S; ipLo < NumLevRecomb[ipISO][nelem]; ipLo++ )
00930                                 {
00931                                         fprintf(ioRECOMB, "%li\t%li", nelem, ipLo );
00932                                         /* loop over temperatures to produce array of recombination coefficients        */
00933                                         for(i = 0; i < N_ISO_TE_RECOMB; i++)
00934                                         {
00935                                                 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
00936                                                 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
00937                                                 TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
00938                                                 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
00939                                                 fprintf(ioRECOMB, "\t%f", RRCoef[ipISO][nelem][ipLo][i] );
00940                                         }
00941                                         fprintf(ioRECOMB, "\n" );
00942                                 }
00943 
00944                                 /* Store one additional line in XX_iso_recomb.dat that gives the total recombination,
00945                                  * as computed by the sum so far, plus levels up to NHYDRO_MAX_LEVEL using Verner's fits,
00946                                  * plus levels up to SumUpToThisN using Seaton 59, for each element and each temperature.       */
00947                                 fprintf(ioRECOMB, "%li\t%li", nelem, NumLevRecomb[ipISO][nelem] );
00948                                 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00949                                 {
00950                                         for( i1 = ( (nelem == ipISO) ? (RREC_MAXN + 1) : (LIKE_RREC_MAXN( nelem ) + 1) ); i1< NHYDRO_MAX_LEVEL; i1++ )
00951                                         {
00952                                                 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
00953                                         }
00954                                         for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
00955                                         {
00956                                                 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
00957                                         }
00958                                         fprintf(ioRECOMB, "\t%f", log10( TotalRecomb[ipISO][nelem][i] ) );
00959                                 }
00960                                 fprintf(ioRECOMB, "\n" );
00961                         }
00962                         /* end the file with the same information */
00963                         fprintf(ioRECOMB,"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient [H-LIke/HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
00964                                 RECOMBMAGIC ,
00965                                 NumLevRecomb[ipISO][ipISO],
00966                                 NumLevRecomb[ipISO][ipISO+1],
00967                                 N_ISO_TE_RECOMB,
00968                                 iso.chISO[ipISO],
00969                                 NumLevRecomb[ipISO][ipISO],
00970                                 elementnames.chElementSym[ipISO],
00971                                 NumLevRecomb[ipISO][ipISO+1],
00972                                 N_ISO_TE_RECOMB );
00973 
00974                         fclose( ioRECOMB );
00975 
00976                         fprintf( ioQQQ, "iso_recomb_setup: compilation complete, %s created.\n", chFilename[ipISO] );
00977                         fprintf( ioQQQ, "The compilation is completed successfully.\n");
00978                         cdEXIT(EXIT_SUCCESS);
00979                 }
00980         }
00981 
00982         return;
00983 }
00984 
00985 double iso_dielec_recomb_rate( long ipISO, long nelem, long ipLo )
00986 {
00987         double rate;
00988         long ipTe, i;
00989         double TeDRCoef[NUM_DR_TEMPS];
00990         const double Te_over_Z1_Squared[NUM_DR_TEMPS] = {
00991                 1.00000,        1.30103,        1.69897,        2.00000,        2.30103,        2.69897,        3.00000,
00992                 3.30103,        3.69897,        4.00000,        4.30103,        4.69897,        5.00000,        5.30103,
00993                 5.69897,        6.00000,        6.30103,        6.69897,        7.00000 };
00994 
00995         DEBUG_ENTRY( "iso_dielec_recomb_rate()" );
00996 
00997         /* currently only two iso sequences and only he-like is applicable. */
00998         ASSERT( ipISO == ipHE_LIKE );
00999         ASSERT( ipLo >= 0 );
01000 
01001         /* temperature grid is nelem^2 * constant temperature grid above. */
01002         for( i=0; i<NUM_DR_TEMPS; i++ )
01003         {
01004                 TeDRCoef[i] = Te_over_Z1_Squared[i] + 2. * log10( (double) nelem );
01005         }
01006 
01007         if( ipLo == ipHe1s1S )
01008         {
01009                 rate = 0.;
01010         }
01011         else if( ipLo<iso.numLevels_max[ipISO][nelem] )
01012         {
01013                 if( phycon.alogte <= TeDRCoef[0] )
01014                 {
01015                         /* Take lowest tabulated value for low temperature end. */
01016                         rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][0];
01017                 }
01018                 else if( phycon.alogte >= TeDRCoef[NUM_DR_TEMPS-1] )
01019                 {
01020                         /* use T^-1.5 extrapolation at high temperatures. */
01021                         rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][NUM_DR_TEMPS-1] * 
01022                                 pow( 10., 1.5* (TeDRCoef[NUM_DR_TEMPS-1] - phycon.alogte ) ) ;
01023                 }
01024                 else
01025                 {
01026                         /* find temperature in tabulated values.  */
01027                         ipTe = hunt_bisect( TeDRCoef, NUM_DR_TEMPS, phycon.alogte );                    
01028 
01029                         ASSERT( (ipTe >=0) && (ipTe < NUM_DR_TEMPS-1)  );
01030 
01031                         if( iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1] == 0. )
01032                                 rate = 0.;
01033                         else if( iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe] == 0. )
01034                                 rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1];
01035                         else
01036                         {
01037                                 /* interpolate between tabulated points */
01038                                 rate = log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe]) + 
01039                                         (phycon.alogte-TeDRCoef[ipTe])*
01040                                         (log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1])-log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe]))/
01041                                         (TeDRCoef[ipTe+1]-TeDRCoef[ipTe]);
01042 
01043                                 rate = pow( 10., rate );
01044                         }
01045                 }
01046         }
01047         else 
01048         {
01049                 rate = 0.;
01050         }
01051 
01052         ASSERT( rate >= 0. && rate < 1.0e-12 );
01053 
01054         return rate*iso.lgDielRecom[ipISO];
01055 }
01056 
01057 /* TempInterp - interpolate on an array */
01059 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements )
01060 {
01061         static long int ipTe=-1;
01062         double rate = 0.;
01063         long i0;
01064 
01065         DEBUG_ENTRY( "TempInterp()" );
01066 
01067         ASSERT( fabs( 1. - (double)phycon.alogte/log10((double)phycon.te) ) < 0.0001 );
01068 
01069         if( ipTe < 0 )
01070         {
01071                 /* te totally unknown */
01072                 if( ( phycon.alogte < TempArray[0] ) || 
01073                         ( phycon.alogte > TempArray[NumElements-1] ) )
01074                 {
01075                         fprintf(ioQQQ," TempInterp called with te out of bounds \n");
01076                         cdEXIT(EXIT_FAILURE);
01077                 }
01078                 ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte );                    
01079         }
01080         else if( phycon.alogte < TempArray[ipTe] )
01081         {
01082                 /* temp is too low, must also lower ipTe */
01083                 ASSERT( phycon.alogte > TempArray[0] );
01084                 /* decrement ipTe until it is correct */
01085                 while( ( phycon.alogte < TempArray[ipTe] ) && ipTe > 0)
01086                         --ipTe;
01087         }
01088         else if( phycon.alogte > TempArray[ipTe+1] )
01089         {
01090                 /* temp is too high */
01091                 ASSERT( phycon.alogte <= TempArray[NumElements-1] );
01092                 /* increment ipTe until it is correct */
01093                 while( ( phycon.alogte > TempArray[ipTe+1] ) && ipTe < NumElements-1)
01094                         ++ipTe;
01095         }
01096 
01097         ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
01098 
01099         /* ipTe should now be valid */
01100         ASSERT( ( phycon.alogte >= TempArray[ipTe] )
01101                 && ( phycon.alogte <= TempArray[ipTe+1] ) && ( ipTe < NumElements-1 ) );
01102 
01103         if( ValueArray[ipTe+1] == 0. && ValueArray[ipTe] == 0. )
01104         {
01105                 rate = 0.;
01106         }
01107         else
01108         {
01109                 /* Do a four-point interpolation */
01110                 const int ORDER = 3; /* order of the fitting polynomial */
01111                 i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
01112                 rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte );
01113         }
01114 
01115         return rate;
01116 }
/home66/gary/public_html/cloudy/c08_branch/source/iso_radiative_recomb.cpp
