cloudy trunk
|
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 /*ion_recomb generate recombination coefficients for any species */ 00004 /*ion_recombAGN generate recombination coefficients for AGN table */ 00005 #include "cddefines.h" 00006 #include "phycon.h" 00007 #include "heavy.h" 00008 #include "hmi.h" 00009 #include "grainvar.h" 00010 #include "dense.h" 00011 #include "conv.h" 00012 #include "thermal.h" 00013 #include "iso.h" 00014 #include "abund.h" 00015 #include "punch.h" 00016 #include "elementnames.h" 00017 #include "atmdat.h" 00018 #include "ionbal.h" 00019 00020 /*ion_recomb generate recombination coefficients for any species */ 00021 void ion_recomb( 00022 /* this is debug flag */ 00023 bool lgPrintIt, 00024 const double *dicoef, 00025 const double *dite, 00026 const double ditcrt[], 00027 const double aa[], 00028 const double bb[], 00029 const double cc[], 00030 const double dd[], 00031 const double ff[], 00032 /* nelem is the atomic number on the C scale, 0 for H */ 00033 long int nelem/*, 00034 double tlow[]*/) 00035 { 00036 #define DICOEF(I_,J_) (*(dicoef+(I_)*(nelem+1)+(J_))) 00037 #define DITE(I_,J_) (*(dite+(I_)*(nelem+1)+(J_))) 00038 long int i, 00039 ion, 00040 limit; 00041 double 00042 fac2, 00043 factor, 00044 DielRecomRateCoef_HiT[LIMELM] , 00045 DielRecomRateCoef_LowT[LIMELM] , 00046 ChargeTransfer[LIMELM] , 00047 t4m1, 00048 tefac; 00049 00050 /* these are used for adding noise to rec coef */ 00051 static double RecNoise[LIMELM][LIMELM]; 00052 static bool lgNoiseNeedEval=true; 00053 00054 DEBUG_ENTRY( "ion_recomb(false,)" ); 00055 00056 /* set up ionization balance matrix, C(I,1)=destruction, 2=creation 00057 * heating rates saved in array B(I) in same scratch block 00058 * factor is for Aldrovandi+Pequignot fit, FAC2 is for Nuss+Storey 00059 * fit for dielectronic recombination 00060 * GrnIonRec is rate ions recombine on grain surface, normally zero; 00061 * set in hmole, already has factor of hydrogen density 00062 * */ 00063 00064 /* routine only used for Li on up */ 00065 ASSERT( nelem < LIMELM); 00066 ASSERT( nelem > 1 ); 00067 00068 /* check that range of ionization is correct */ 00069 ASSERT( dense.IonLow[nelem] >= 0 ); 00070 ASSERT( dense.IonLow[nelem] <= nelem+1 ); 00071 00072 atmdat.nsbig = MAX2(dense.IonHigh[nelem]+1,atmdat.nsbig); 00073 t4m1 = 1e4/phycon.te; 00074 fac2 = 1e-14*phycon.sqrte; 00075 00076 /* option to put noise into rec coefficient - 00077 * one time initialization of noise - this is set with command 00078 * set dielectronic recombination kludge noise */ 00079 if( lgNoiseNeedEval ) 00080 { 00081 int n; 00082 if( ionbal.guess_noise !=0. ) 00083 { 00084 for( n=ipHYDROGEN; n<LIMELM; ++n ) 00085 { 00086 for( ion=0; ion<=n; ++ion ) 00087 { 00088 /* log normal noise with dispersion entered on command line */ 00089 /* NB the seed for rand was set when the command was parsed */ 00090 RecNoise[n][ion] = pow(10., RandGauss( 0. , ionbal.guess_noise ) ); 00091 } 00092 } 00093 } 00094 else 00095 { 00096 for( n=ipHYDROGEN; n<LIMELM; ++n ) 00097 { 00098 for( ion=0; ion<=n; ++ion ) 00099 { 00100 RecNoise[n][ion] = 1.; 00101 } 00102 } 00103 } 00104 lgNoiseNeedEval = false; 00105 } 00106 00107 /* this routine only does simple two-level species, 00108 * loop over ions will be <= limit, IonHigh is -1 since very 00109 * highest stage of ionization is not recombined into. 00110 * for Li, will do only atom, since ions are H and He like, 00111 * so limit is zero */ 00112 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1); 00113 ASSERT( limit >= -1 ); 00114 00115 /* zero-out loop comes before main loop since there are off-diagonal 00116 * elements in the main ionization loop, due to multi-electron processes */ 00117 /* >>chng 00 dec 07, limit changed to identical to ion_solver */ 00118 for( ion=0; ion <= limit; ion++ ) 00119 { 00120 ionbal.RateRecomTot[nelem][ion] = 0.; 00121 ChargeTransfer[ion] = 0.; 00122 DielRecomRateCoef_LowT[ion] = 0.; 00123 DielRecomRateCoef_HiT[ion] = 0.; 00124 ionbal.RR_rate_coef_used[nelem][ion] = 0.; 00125 ionbal.DR_rate_coef_used[nelem][ion] = 0.; 00126 } 00127 for( ion=limit+1; ion < LIMELM; ion++ ) 00128 { 00129 /* >>chng 01 dec 18, do not set this to FLT_MAX since it clobbers what 00130 * had been set in h-like and he-like routines - that would only affect 00131 * the printout */ 00132 ChargeTransfer[ion] = -FLT_MAX; 00133 DielRecomRateCoef_LowT[ion] = -FLT_MAX; 00134 DielRecomRateCoef_HiT[ion] = -FLT_MAX; 00135 } 00136 00137 DielRecomRateCoef_HiT[nelem] = 0.; 00138 DielRecomRateCoef_HiT[nelem-1] = 0.; 00139 00140 DielRecomRateCoef_LowT[nelem] = 0.; 00141 DielRecomRateCoef_LowT[nelem-1] = 0.; 00142 00143 /* these are counted elsewhere and must not be added here */ 00144 Heavy.xLyaHeavy[nelem][nelem] = 0.; 00145 Heavy.xLyaHeavy[nelem][nelem-1] = 0.; 00146 00147 /* IonLow is 0 for the atom, limit chosen to NOT include iso sequences */ 00148 for( ion=dense.IonLow[nelem]; ion <= limit; ion++ ) 00149 { 00150 /* number of bound electrons of the ion after recombination, 00151 * for an atom (ion=0) this is 00152 * equal to nelem+1, the element on the physical scale, since nelem is 00153 * on the C scale, being 5 for carbon */ 00154 long n_bnd_elec_after_recomb = nelem+1 - ion; 00155 00156 /* >>chng 02 nov 06 add charge transfer here rather than in ion_solver */ 00157 /* charge transfer recombination of this species by ionizing hydrogen and helium */ 00158 ChargeTransfer[ion] = 00159 /* He0 + ion charge transfer recombination */ 00160 atmdat.HeCharExcRecTo[nelem][ion]* 00161 /* following product is density [cm-3] of ground state of He0 */ 00162 StatesElem[ipHE_LIKE][ipHELIUM][ipHe1s1S].Pop*dense.xIonDense[ipHELIUM][1] + 00163 /* H0 + ion charge transfer recombination */ 00164 atmdat.HCharExcRecTo[nelem][ion]* 00165 /* following product is density [cm-3] of ground state of H0 */ 00166 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]; 00167 00168 /*>>chng 04 feb 20, add this, had always been in for destruction for H- */ 00169 /* charge transfer recombination of first ion to neutral, by reaction with H- 00170 * the ion==0 is right, the first array element is the */ 00171 if( ion==0 && nelem>ipHELIUM && atmdat.lgCTOn ) 00172 ChargeTransfer[ion] += hmi.hmin_ct_firstions * hmi.Hmolec[ipMHm]; 00173 00174 /* >>chng 06 feb 01, add option to use Badnell RR data rather than Verner */ 00175 if( ionbal.lgRR_recom_Badnell_use && ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] ) 00176 { 00177 ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Badnell_rate_coef[nelem][ion]; 00178 } 00179 else 00180 { 00181 ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Verner_rate_coef[nelem][ion]; 00182 } 00183 00184 /* Burgess or high-T dielectronic recombination */ 00185 DielRecomRateCoef_HiT[ion] = 0.; 00186 /* >>chng 03 oct 29, do fe here rather than after loop */ 00187 if( nelem==ipIRON ) 00188 { 00189 /* implements dn > 0 DR from Arnaud & Raymond 1992 */ 00190 DielRecomRateCoef_HiT[ion] = atmdat_dielrec_fe(ion+1,phycon.te); 00191 } 00192 else if( phycon.te > (ditcrt[ion]*0.1) ) 00193 { 00194 DielRecomRateCoef_HiT[ion] = ionbal.DielSupprs[0][ion]/phycon.te32* 00195 DICOEF(0,ion)*exp(-DITE(0,ion)/phycon.te)* 00196 (1. + DICOEF(1,ion)* 00197 sexp(DITE(1,ion)/phycon.te)); 00198 } 00199 00200 /* begin dn = 0 dielectronic recombination 00201 * do not include it for rec from 00202 * a closed shell, n_bnd_elec_after_recomb-1 is number of bound electrons in parent ion */ 00203 DielRecomRateCoef_LowT[ion] = 0.; 00204 if( ((n_bnd_elec_after_recomb-1) != 2) && 00205 ((n_bnd_elec_after_recomb-1) != 10) && 00206 ((n_bnd_elec_after_recomb-1) != 18) ) 00207 { 00208 tefac = ff[ion]*t4m1; 00209 /* do not do iron here since all dn = 0 DR either Badnell or a guess */ 00210 if( ff[ion] != 0. && nelem!=ipIRON ) 00211 { 00212 /* >>chng 06 feb 14, O+3 ff[ion] is very negative, as a result exp goes to +inf 00213 * at very low temperatures. This is a error in the Nussbaumer & Storey fits to DR. 00214 * do not use them is tefac = ff[ion] / t4 is very negative */ 00215 if( tefac > -5. ) 00216 { 00217 factor = (((aa[ion]*t4m1+bb[ion])*t4m1+cc[ion])*t4m1+dd[ion])* sexp(tefac); 00218 DielRecomRateCoef_LowT[ion] = ionbal.DielSupprs[1][ion]*fac2* 00219 MAX2(0.,factor ); 00220 } 00221 else 00222 { 00223 DielRecomRateCoef_LowT[ion] = 0.; 00224 } 00225 } 00226 else if( ionbal.lg_guess_coef ) 00227 { 00228 /* first guesses are those based on Nussbaumer & Storey and are from 00229 * >>refer ion DR Ali, B., Blum, R. D., Bumgardner, T. E., Cranmer, S. R., 00230 * >>refercon Ferland, G. J., Haefner, R. I., & Tiede, G. P. 1991, PASP, 103, 1182*/ 00231 if( (ff[ion] == 0.) && (ion <= 3) ) 00232 { 00233 /* these are guessed dielectronic rec rates, taken from 00234 * >>refer all dielectronic Ali, B., Blum, R. D., Bumgardner, T. E., Cranmer, S. R., Ferland, G. J., 00235 * >>refercon Haefner, R. I., & Tiede, G. P. 1991, PASP, 103, 1182 */ 00236 static double cludge[4]={3e-13,3e-12,1.5e-11,2.5e-11}; 00237 00238 /* >>chng 01 jun 30, make GuessDiel array */ 00239 DielRecomRateCoef_LowT[ion] = ionbal.DielSupprs[1][ion]*cludge[ion]* 00240 /* this is optional scale factor on set dielectronic rec command 00241 * >>chng 06 feb 07, add Boltzmann factor to induce behavior 00242 * like in Nussbaumer & Storey */ 00243 ionbal.GuessDiel[ion] * sexp( 1e3 / phycon.te ); 00244 } 00245 /* this implements the low-T kludge dielectronic command */ 00246 else 00247 { 00248 /* assume the recombination rate is the coefficient scanned off the steve option, times the charge 00249 * before recombination raised to a power */ 00250 double fitcoef[3][3] = 00251 { 00252 /* L- shell, Z=17-23 */ 00253 {-52.5073,+5.19385,-0.126099} , 00254 /* M-shell (Z=9-15) */ 00255 {-10.9679,1.66397,-0.0605965} , 00256 /* N shell */ 00257 {-3.95599,1.61884,-0.194540} 00258 }; 00259 00260 /* these implement the guesses in 00261 * Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ, 604, 556 */ 00262 if( nelem==ipIRON ) 00263 { 00264 int nshell; 00265 if( (n_bnd_elec_after_recomb>=4) && (n_bnd_elec_after_recomb<=11) ) 00266 nshell = 0; 00267 else if( (n_bnd_elec_after_recomb>=12) && (n_bnd_elec_after_recomb<=19 )) 00268 nshell = 1; 00269 else 00270 nshell = 2; 00271 /* n_bnd_elec_after_recomb is the number of bound electrons */ 00272 DielRecomRateCoef_LowT[ion] = fitcoef[nshell][0] + 00273 fitcoef[nshell][1]*(ion+1) + 00274 fitcoef[nshell][2]*POW2(ion+1.); 00275 DielRecomRateCoef_LowT[ion] = 1e-10*pow(10.,DielRecomRateCoef_LowT[ion]); 00276 00277 } 00278 else 00279 { 00280 /* this is guess for all other elements presented in 00281 * >>refer all DR Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ, 604, 556 */ 00282 DielRecomRateCoef_LowT[ion] = 3e-12*pow(10.,(double)(ion+1)*0.1); 00283 } 00284 /* >>chng 06 feb 02, add option to use mean of Badnell DR in place 00285 * of these hacks */ 00286 if( ionbal.lg_use_DR_Badnell_rate_coef_mean_ion ) 00287 DielRecomRateCoef_LowT[ion] = ionbal.DR_Badnell_rate_coef_mean_ion[ion]; 00288 } 00289 /* include optional noise here 00290 * >>chng 06 feb 07, move noise down to here so that use for both 00291 * guesses of DR rates */ 00292 DielRecomRateCoef_LowT[ion] *= RecNoise[nelem][ion]; 00293 } 00294 } 00295 /* >>chng 05 dec 19, add option to use Badnell numbers */ 00296 /* this is total old DR rates - may not use it */ 00297 ionbal.DR_old_rate_coef[nelem][ion] = DielRecomRateCoef_HiT[ion] + DielRecomRateCoef_LowT[ion]; 00298 00299 /* set total DR rate - either Badnell if it exists or on with set dielectronic recombination badnell */ 00300 if( ionbal.lgDR_recom_Badnell_use && ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] ) 00301 { 00302 ionbal.DR_rate_coef_used[nelem][ion] = ionbal.DR_Badnell_rate_coef[nelem][ion]; 00303 } 00304 else 00305 { 00306 ionbal.DR_rate_coef_used[nelem][ion] = ionbal.DR_old_rate_coef[nelem][ion]; 00307 } 00308 00309 /* sum of recombination rates [units s-1] for radiative, three body, charge transfer */ 00310 ionbal.RateRecomTot[nelem][ion] = 00311 dense.eden* ( 00312 ionbal.RR_rate_coef_used[nelem][ion] + 00313 ionbal.DR_rate_coef_used[nelem][ion] + 00314 ionbal.CotaRate[ion] ) + 00315 ChargeTransfer[ion]; 00316 00317 /* >>chng 01 jun 30, FRAC_LINE was 0.1, not 1, did not include anything except 00318 * radiative recombination, the radrec term */ 00319 # define FRAC_LINE 1. 00320 /* was 0.1 */ 00321 /*Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden*radrec*FRAC_LINE );*/ 00322 Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden* 00323 (ionbal.RR_rate_coef_used[nelem][ion]+DielRecomRateCoef_LowT[ion]+DielRecomRateCoef_HiT[ion])*FRAC_LINE ); 00324 } 00325 00326 /* option to punch rec coefficients */ 00327 if( punch.lgioRecom || lgPrintIt ) 00328 { 00329 /* >>chng 04 feb 22, make option to print ions for single element */ 00330 FILE *ioOut; 00331 if( lgPrintIt ) 00332 ioOut = ioQQQ; 00333 else 00334 ioOut = punch.ioRecom; 00335 00336 /* print name of element */ 00337 fprintf( ioOut, 00338 " %s recombination coefficients fnzone:%.2f \tte\t%.4e\tne\t%.4e\n", 00339 elementnames.chElementName[nelem] , fnzone , phycon.te , dense.eden ); 00340 00341 /*limit = MIN2(11,dense.IonHigh[nelem]);*/ 00342 /* >>chng 05 sep 24, just print one long line - need info */ 00343 limit = dense.IonHigh[nelem]; 00344 for( i=0; i < limit; i++ ) 00345 { 00346 fprintf( ioOut, "%10.2e", ionbal.RR_rate_coef_used[nelem][i] ); 00347 } 00348 fprintf( ioOut, " radiative used vs Z\n" ); 00349 00350 for( i=0; i < limit; i++ ) 00351 { 00352 fprintf( ioOut, "%10.2e", ionbal.RR_Verner_rate_coef[nelem][i] ); 00353 } 00354 fprintf( ioOut, " old Verner vs Z\n" ); 00355 00356 for( i=0; i < limit; i++ ) 00357 { 00358 fprintf( ioOut, "%10.2e", ionbal.RR_Badnell_rate_coef[nelem][i] ); 00359 } 00360 fprintf( ioOut, " new Badnell vs Z\n" ); 00361 00362 for( i=0; i < limit; i++ ) 00363 { 00364 /* >>chng 06 jan 19, from div by eden to div by H0 - want units of cm3 s-1 but 00365 * no single collider does this so not possible to get rate coefficient easily 00366 * H0 is more appropriate than electron density */ 00367 fprintf( ioOut, "%10.2e", ChargeTransfer[i]/SDIV(dense.xIonDense[ipHYDROGEN][0]) ); 00368 } 00369 fprintf( ioOut, " CT/n(H0)\n" ); 00370 00371 for( i=0; i < limit; i++ ) 00372 { 00373 fprintf( ioOut, "%10.2e", ionbal.CotaRate[ion] ); 00374 } 00375 fprintf( ioOut, " 3body vs Z /ne\n" ); 00376 00377 /* note different upper limit - this routine does grain rec for all ions */ 00378 for( i=0; i < dense.IonHigh[nelem]; i++ ) 00379 { 00380 fprintf( ioOut, "%10.2e", gv.GrainChTrRate[nelem][i+1][i]/dense.eden ); 00381 } 00382 fprintf( ioOut, " Grain vs Z /ne\n" ); 00383 00384 for( i=0; i < limit; i++ ) 00385 { 00386 fprintf( ioOut, "%10.2e", DielRecomRateCoef_HiT[i] ); 00387 } 00388 fprintf( ioOut, " Burgess vs Z\n" ); 00389 00390 for( i=0; i < limit; i++ ) 00391 { 00392 fprintf( ioOut, "%10.2e", ionbal.DR_old_rate_coef[nelem][i] ); 00393 } 00394 fprintf( ioOut, " old Nussbaumer Storey DR vs Z\n" ); 00395 00396 for( i=0; i < limit; i++ ) 00397 { 00398 fprintf( ioOut, "%10.2e", ionbal.DR_Badnell_rate_coef[nelem][i] ); 00399 } 00400 fprintf( ioOut, " new Badnell DR vs Z\n" ); 00401 00402 for( i=0; i < limit; i++ ) 00403 { 00404 fprintf( ioOut, "%10.2e", DielRecomRateCoef_LowT[i] ); 00405 } 00406 fprintf( ioOut, " low T DR used vs Z\n" ); 00407 00408 /* total recombination rate, with density included - this goes into the matrix */ 00409 for( i=0; i < limit; i++ ) 00410 { 00411 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] ); 00412 } 00413 fprintf( ioOut, 00414 " total rec rate (with density) for %s\n", 00415 elementnames.chElementSym[nelem] ); 00416 for( i=0; i < limit; i++ ) 00417 { 00418 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i]/dense.eden ); 00419 } 00420 fprintf( ioOut, 00421 " total rec rate / ne for %s\n\n", 00422 elementnames.chElementSym[nelem] ); 00423 00424 /* spill over to next line for many stages of ionization */ 00425 if( dense.IonHigh[nelem] > 11 ) 00426 { 00427 limit = MIN2(29,dense.IonHigh[nelem]); 00428 fprintf( ioOut, " R " ); 00429 for( i=11; i < limit; i++ ) 00430 { 00431 fprintf( ioOut, "%10.2e", dense.eden*ionbal.CotaRate[ion] ); 00432 } 00433 fprintf( ioOut, "\n" ); 00434 00435 fprintf( ioOut, " B " ); 00436 for( i=11; i < limit; i++ ) 00437 { 00438 fprintf( ioOut, "%10.2e", DielRecomRateCoef_HiT[i] ); 00439 } 00440 fprintf( ioOut, "\n" ); 00441 00442 fprintf( ioOut, " NS" ); 00443 for( i=11; i < limit; i++ ) 00444 { 00445 fprintf( ioOut, "%10.2e", DielRecomRateCoef_LowT[i] ); 00446 } 00447 fprintf( ioOut, "\n" ); 00448 00449 fprintf( ioOut, " " ); 00450 for( i=11; i < limit; i++ ) 00451 { 00452 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] ); 00453 } 00454 fprintf( ioOut, "\n\n" ); 00455 } 00456 } 00457 00458 /* >>chng 02 nov 09, from -2 to -NISO */ 00459 /*limit = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/ 00460 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1); 00461 for( i=dense.IonLow[nelem]; i <= limit; i++ ) 00462 { 00463 ASSERT( Heavy.xLyaHeavy[nelem][i] > 0. ); 00464 ASSERT( ionbal.RateRecomTot[nelem][i] > 0. ); 00465 } 00466 return; 00467 # undef DITE 00468 # undef DICOEF 00469 } 00470 00471 /*ion_recombAGN generate recombination coefficients for AGN table */ 00472 void ion_recombAGN( FILE * io ) 00473 { 00474 # define N1LIM 3 00475 # define N2LIM 4 00476 double te1[N1LIM]={ 5000., 10000., 20000.}; 00477 double te2[N2LIM]={ 20000.,50000.,100000.,1e6}; 00478 /* this is boundary between two tables */ 00479 double BreakEnergy = 100./13.0; 00480 long int nelem, ion , i; 00481 /* this will hold element symbol + ionization */ 00482 char chString[100], 00483 chOutput[100]; 00484 /* save temp here */ 00485 double TempSave = phycon.te; 00486 /* save ne here */ 00487 double EdenSave = dense.eden; 00488 00489 DEBUG_ENTRY( "ion_recomb(false,)" ); 00490 00491 dense.eden = 1.; 00492 /*atmdat_readin();*/ 00493 00494 /* first put header on file */ 00495 fprintf(io,"X+i\\Te"); 00496 for( i=0; i<N1LIM; ++i ) 00497 { 00498 phycon.te = te1[i]; 00499 fprintf(io,"\t%.0f K",phycon.te); 00500 } 00501 fprintf(io,"\n"); 00502 00503 /* now do loop over temp, but add elements */ 00504 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem ) 00505 { 00506 /* this list of elements included in the AGN tables is defined in zeroabun.c */ 00507 if( abund.lgAGN[nelem] ) 00508 { 00509 for( ion=0; ion<=nelem; ++ion ) 00510 { 00511 ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 ); 00512 00513 if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy ) 00514 break; 00515 00516 /* print chemical symbol */ 00517 sprintf(chOutput,"%s", 00518 elementnames.chElementSym[nelem]); 00519 /* some elements have only one letter - this avoids leaving a space */ 00520 if( chOutput[1]==' ' ) 00521 chOutput[1] = chOutput[2]; 00522 /* now ionization stage */ 00523 if( ion==0 ) 00524 { 00525 sprintf(chString,"0 "); 00526 } 00527 else if( ion==1 ) 00528 { 00529 sprintf(chString,"+ "); 00530 } 00531 else 00532 { 00533 sprintf(chString,"+%li ",ion); 00534 } 00535 strcat( chOutput , chString ); 00536 fprintf(io,"%5s",chOutput ); 00537 00538 for( i=0; i<N1LIM; ++i ) 00539 { 00540 TempChange(te1[i] , false); 00541 dense.IonLow[nelem] = 0; 00542 dense.IonHigh[nelem] = nelem+1; 00543 if( ConvBase(0) ) 00544 fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n"); 00545 fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]); 00546 } 00547 fprintf(io,"\n"); 00548 } 00549 fprintf(io,"\n"); 00550 } 00551 } 00552 00553 /* second put header on file */ 00554 fprintf(io,"X+i\\Te"); 00555 for( i=0; i<N2LIM; ++i ) 00556 { 00557 TempChange(te2[i] , false); 00558 fprintf(io,"\t%.0f K",phycon.te); 00559 } 00560 fprintf(io,"\n"); 00561 00562 /* now do same loop over temp, but add elements */ 00563 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00564 { 00565 /* this list of elements included in the AGN tables is defined in zeroabun.c */ 00566 if( abund.lgAGN[nelem] ) 00567 { 00568 for( ion=0; ion<=nelem; ++ion ) 00569 { 00570 ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 ); 00571 00572 if( Heavy.Valence_IP_Ryd[nelem][ion] <= BreakEnergy ) 00573 continue; 00574 00575 /* print chemical symbol */ 00576 fprintf(io,"%s", 00577 elementnames.chElementSym[nelem]); 00578 /* now ionization stage */ 00579 if( ion==0 ) 00580 { 00581 fprintf(io,"0 "); 00582 } 00583 else if( ion==1 ) 00584 { 00585 fprintf(io,"+ "); 00586 } 00587 else 00588 { 00589 fprintf(io,"+%li",ion); 00590 } 00591 00592 for( i=0; i<N2LIM; ++i ) 00593 { 00594 TempChange(te2[i] , false); 00595 dense.IonLow[nelem] = 0; 00596 dense.IonHigh[nelem] = nelem+1; 00597 if( ConvBase(0) ) 00598 fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n"); 00599 fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]); 00600 } 00601 fprintf(io,"\n"); 00602 } 00603 fprintf(io,"\n"); 00604 } 00605 } 00606 00607 TempChange(TempSave , true); 00608 dense.eden = EdenSave; 00609 return; 00610 }