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 /*CoolIron compute iron cooling */ 00004 /*Fe_10_11_13_cs evaluate collision strength for Fe */ 00005 /*fe14cs compute collision strengths for forbidden transitions */ 00006 /*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */ 00007 /*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */ 00008 /*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */ 00009 /*Fe2_cooling compute cooling due to FeII emission */ 00010 /*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */ 00011 /*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */ 00012 #include "cddefines.h" 00013 #include "physconst.h" 00014 #include "dense.h" 00015 #include "coolheavy.h" 00016 #include "taulines.h" 00017 #include "phycon.h" 00018 #include "iso.h" 00019 #include "conv.h" 00020 #include "trace.h" 00021 #include "hydrogenic.h" 00022 #include "ligbar.h" 00023 #include "cooling.h" 00024 #include "thermal.h" 00025 #include "lines_service.h" 00026 #include "atoms.h" 00027 #include "atomfeii.h" 00028 #include "fe.h" 00029 #define NLFE2 6 00030 00031 /*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */ 00032 STATIC void Fe11Lev5(void); 00033 00034 /*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */ 00035 STATIC void Fe13Lev5(void); 00036 00037 /*fe14cs compute collision strengths for forbidden transitions */ 00038 STATIC void fe14cs(double te1, 00039 double *csfe14); 00040 00041 /*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */ 00042 STATIC void Fe7Lev8(void); 00043 00044 /*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */ 00045 STATIC void Fe3Lev14(void); 00046 00047 /*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */ 00048 STATIC void Fe4Lev12(void); 00049 00050 /*Fe_10_11_13_cs evaluate collision strength for Fe */ 00051 STATIC double Fe_10_11_13_cs( 00052 /* ion, one of 10, 11, or 13 on physics scale */ 00053 int ion, 00054 /* initial and final index of levels, with lowest energy 1, highest of 5 00055 * on physics scale */ 00056 int init, 00057 int final ) 00058 { 00059 # define N 10 00060 static double Fe10cs[6][6][2]; 00061 static double Fe11cs[6][6][2]; 00062 static double Fe13cs[6][6][2]; 00063 int i, j; 00064 double cs; 00065 int index = 0; 00066 double temp_max, temp_min = 4; 00067 double temp_log = phycon.alogte; 00068 static bool lgFirstTime = true; 00069 00070 DEBUG_ENTRY( "Fe_10_11_13_cs()" ); 00071 00072 if( lgFirstTime ) 00073 { 00074 /* fit coefficients for collision strengths */ 00075 double aFe10[N] = {10.859,-1.1541,11.593,22.333,-0.4283,7.5663,3.087,1.0937,0.8261,59.678}; 00076 double bFe10[N] = {-1.4804,0.4956,-2.1096,-4.1288,0.1929,-1.3525,-0.5531,-0.1748,-0.1286,-11.081}; 00077 double aFe11[N] = {5.7269,1.2885,4.0877,0.4571,1.2911,2.2339,0.3621,0.7972,0.2225,1.1021}; 00078 double bFe11[N] = {-0.7559,-0.1671,-0.5678,-0.0653,-0.1589,-0.2924,-0.0506,-0.1038,-0.0302,-0.1062}; 00079 double aFe13[N] = {2.9102,1.8682,-0.353,0.0622,14.229,-4.3845,0.0375,-6.9222,0.688,-0.0609}; 00080 double bFe13[N] = {-0.4158,-0.242,0.1417,0.0023,-2.0643,1.2573,0.0286,2.0919,-0.083,0.1487}; 00081 /* do one-time initialization 00082 * assigning array initially to zero */ 00083 for(i=0; i<6; i++) 00084 { 00085 for(j=0; j<6; j++) 00086 { 00087 set_NaN( Fe10cs[i][j], 2L ); 00088 set_NaN( Fe11cs[i][j], 2L ); 00089 set_NaN( Fe13cs[i][j], 2L ); 00090 } 00091 } 00092 00093 /* reading coefficients into 3D array */ 00094 for(i=1; i<6; i++) 00095 { 00096 for(j=i+1; j<6; j++) 00097 { 00098 Fe10cs[i][j][0] = aFe10[index]; 00099 Fe10cs[i][j][1] = bFe10[index]; 00100 Fe11cs[i][j][0] = aFe11[index]; 00101 Fe11cs[i][j][1] = bFe11[index]; 00102 Fe13cs[i][j][0] = aFe13[index]; 00103 Fe13cs[i][j][1] = bFe13[index]; 00104 index++; 00105 } 00106 } 00107 lgFirstTime = false; 00108 } 00109 00110 /* Invalid entries returns '-1':the initial indices are smaller than 00111 * the final indices */ 00112 if(init >= final) 00113 { 00114 cs = -1; 00115 } 00116 /* Invalid returns '-1': the indices are greater than 5 or smaller than 0 */ 00117 else if(init < 1 || init > 4 || final < 2 || final > 5) 00118 { 00119 cs = -1; 00120 } 00121 else 00122 { 00123 /* cs = a + b*log(T) 00124 * if temp is out of range, return boundary values */ 00125 if(ion == 10) 00126 { 00127 temp_max = 5; 00128 temp_log = MAX2(temp_log, temp_min); 00129 temp_log = MIN2(temp_log, temp_max); 00130 cs = Fe10cs[init][final][0] + Fe10cs[init][final][1]*temp_log; 00131 } 00132 else if(ion == 11) 00133 { 00134 temp_max = 6.7; 00135 temp_log = MAX2(temp_log, temp_min); 00136 temp_log = MIN2(temp_log, temp_max); 00137 cs = Fe11cs[init][final][0] + Fe11cs[init][final][1]*temp_log; 00138 } 00139 else if(ion ==13) 00140 { 00141 temp_max = 5; 00142 temp_log = MAX2(temp_log, temp_min); 00143 temp_log = MIN2(temp_log, temp_max); 00144 cs = Fe13cs[init][final][0] + Fe13cs[init][final][1]*temp_log; 00145 } 00146 else 00147 /* this can't happen */ 00148 TotalInsanity(); 00149 } 00150 00151 return cs; 00152 00153 # undef N 00154 } 00155 00156 /*Fe2_cooling compute cooling due to FeII emission */ 00157 STATIC void Fe2_cooling( void ) 00158 { 00159 long int i , j; 00160 int nNegPop; 00161 00162 static double **AulPump, 00163 **CollRate, 00164 **AulEscp, 00165 **col_str , 00166 **AulDest, 00167 *depart, 00168 *pops, 00169 *destroy, 00170 *create; 00171 00172 static bool lgFirst=true; 00173 bool lgZeroPop; 00174 00175 /* stat weights for Fred's 6 level model FeII atom */ 00176 static double gFe2[NLFE2]={1.,1.,1.,1.,1.,1.}; 00177 /* excitation energies (Kelvin) for Fred's 6 level model FeII atom */ 00178 static double ex[NLFE2]={0.,3.32e4,5.68e4,6.95e4,1.15e5,1.31e5}; 00179 00180 /* these are used to only evaluated FeII one time per temperature, zone, 00181 * and abundance */ 00182 static double TUsed = 0.; 00183 static realnum AbunUsed = 0.; 00184 /* remember which zone last evaluated with */ 00185 static long int nZUsed=-1, 00186 /* make sure at least two calls per zone */ 00187 nCall=0; 00188 00189 DEBUG_ENTRY( "Fe2_cooling()" ); 00190 00191 /* return if nothing to do */ 00192 if( dense.xIonDense[ipIRON][1] == 0. ) 00193 { 00194 /* zero abundance, do nothing */ 00195 /* heating cooling and derivative from large atom */ 00196 FeII.Fe2_large_cool = 0.; 00197 FeII.Fe2_large_heat = 0.; 00198 FeII.ddT_Fe2_large_cool = 0.; 00199 00200 /* cooling and derivative from simple UV atom */ 00201 FeII.Fe2_UVsimp_cool = 0.; 00202 FeII.ddT_Fe2_UVsimp_cool = 0.; 00203 00204 /* now zero out intensities of all FeII lines and level populations */ 00205 FeIIIntenZero(); 00206 return; 00207 } 00208 00209 /* this can be the large 371 level FeII atom 00210 * >>chng 05 dec 04, always call this but with default of 16 levels only 00211 * >>chng 00 jan 06, total rewrite mostly done 00212 * >>chng 97 jan 17, evaluate large FeII atom cooling every time temp changes 00213 * >>chng 97 jan 31, added check on zone since must reevaluate even const temp 00214 * >>chng 99 may 21, reevaluate when zone or temperature changes, but not when 00215 * abundance changes, since we can simply rescale cooling */ 00216 00217 /* totally reevaluate large atom if new zone, or cooling is significant 00218 * and temperature changed, we are in search phase, 00219 * lgSlow option set true with atom FeII slow, forces constant 00220 * evaluation of atom */ 00221 if( FeII.lgSlow || 00222 nzone < 1 || 00223 nzone != nZUsed || 00224 /* on new call, nCall is now set at previous zone's number of calls. 00225 * it is set to zero below, so on second call, nCall is 0. On 00226 * third call nCall is 1. Check for <1 is to insure at least two calls */ 00227 nCall < 1 || 00228 /* check whether things have changed on later calls */ 00229 ( ! fp_equal( phycon.te, TUsed ) && fabs(FeII.Fe2_large_cool/thermal.ctot)> 0.002 && 00230 fabs(dense.xIonDense[ipIRON][1]-AbunUsed)/SDIV(AbunUsed)> 0.002 ) || 00231 ( ! fp_equal( phycon.te, TUsed ) && fabs(FeII.Fe2_large_cool/thermal.ctot)> 0.01) ) 00232 { 00233 00234 if( nZUsed == nzone ) 00235 { 00236 /* not first call, increment, check above to make sure at least 00237 * two evaluations */ 00238 ++nCall; 00239 } 00240 else 00241 { 00242 /* first call this zone set nCall to zero*/ 00243 nCall = 0; 00244 } 00245 00246 /* option to trace convergence and FeII calls */ 00247 if( trace.nTrConvg>=5 ) 00248 { 00249 fprintf( ioQQQ, " CoolIron5 calling FeIILevelPops since "); 00250 if( ! fp_equal( phycon.te, TUsed ) ) 00251 { 00252 fprintf( ioQQQ, 00253 "temperature changed, old new are %g %g, nCall %li ", 00254 TUsed, phycon.te , nCall); 00255 } 00256 else if( nzone != nZUsed ) 00257 { 00258 fprintf( ioQQQ, 00259 "new zone, nCall %li ", nCall ); 00260 } 00261 else if( FeII.lgSlow ) 00262 { 00263 fprintf( ioQQQ, 00264 "FeII.lgSlow set %li", nCall ); 00265 } 00266 else if( conv.lgSearch ) 00267 { 00268 fprintf( ioQQQ, 00269 " in search phase %li", nCall ); 00270 } 00271 else if( nCall < 2 ) 00272 { 00273 fprintf( ioQQQ, 00274 "not second nCall %li " , nCall ); 00275 } 00276 else if( ! fp_equal( phycon.te, TUsed ) && FeII.Fe2_large_cool/thermal.ctot > 0.001 ) 00277 { 00278 fprintf( ioQQQ, 00279 "temp or cooling changed, new are %g %g nCall %li ", 00280 phycon.te, FeII.Fe2_large_cool, nCall ); 00281 } 00282 else 00283 { 00284 fprintf(ioQQQ, "????"); 00285 } 00286 fprintf(ioQQQ, "\n"); 00287 } 00288 00289 /* remember parameters for current conditions */ 00290 TUsed = phycon.te; 00291 AbunUsed = dense.xIonDense[ipIRON][1]; 00292 nZUsed = nzone; 00293 00294 /* this print turned on with atom FeII print command */ 00295 if( FeII.lgPrint ) 00296 { 00297 fprintf(ioQQQ, 00298 " FeIILevelPops called zone %4li te %5f abun %10e c(fe/tot):%6f nCall %li\n", 00299 nzone,phycon.te,AbunUsed,FeII.Fe2_large_cool/thermal.ctot,nCall); 00300 } 00301 00302 /* this solves the multi-level problem, 00303 * sets FeII.Fe2_large_cool, 00304 FeII.Fe2_large_heat, & 00305 FeII.ddT_Fe2_large_cool 00306 but does nothing with them */ 00307 FeIILevelPops(); 00308 { 00309 /*@-redef@*/ 00310 enum{DEBUG_LOC=false}; 00311 /*@+redef@*/ 00312 if( DEBUG_LOC && iteration>1 && nzone>=4 ) 00313 { 00314 fprintf(ioQQQ,"DEBUG1\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 00315 nzone, 00316 phycon.te, 00317 dense.gas_phase[ipHYDROGEN], 00318 dense.eden, 00319 FeII.Fe2_large_cool , 00320 FeII.ddT_Fe2_large_cool , 00321 FeII.Fe2_large_cool/dense.eden/dense.gas_phase[ipHYDROGEN] , 00322 thermal.ctot ); 00323 } 00324 } 00325 00326 if( trace.nTrConvg>=5 || FeII.lgPrint) 00327 { 00328 /* spacing needed to get proper trace convergence output */ 00329 fprintf( ioQQQ, " FeIILevelPops5 returned cool=%.2e heat=%.2e derivative=%.2e\n ", 00330 FeII.Fe2_large_cool,FeII.Fe2_large_heat ,FeII.ddT_Fe2_large_cool); 00331 } 00332 00333 } 00334 else if( ! fp_equal( dense.xIonDense[ipIRON][1], AbunUsed ) ) 00335 { 00336 realnum ratio; 00337 /* this branch, same zone and temperature, but small change in abundance, so just 00338 * rescale cooling and derivative by this change. assumption is that very small changes 00339 * in abundance occurs as ots rates damp out */ 00340 if( trace.nTrConvg>=5 ) 00341 { 00342 fprintf( ioQQQ, 00343 " CoolIron rescaling FeIILevelPops since small change, CFe2=%.2e CTOT=%.2e\n", 00344 FeII.Fe2_large_cool,thermal.ctot); 00345 } 00346 ratio = dense.xIonDense[ipIRON][1]/AbunUsed; 00347 FeII.Fe2_large_cool *= ratio; 00348 FeII.ddT_Fe2_large_cool *= ratio; 00349 FeII.Fe2_large_heat *= ratio; 00350 AbunUsed = dense.xIonDense[ipIRON][1]; 00351 } 00352 else 00353 { 00354 /* this is case where temp is unchanged, so heating and cooling same too */ 00355 if( trace.nTrConvg>=5 ) 00356 { 00357 fprintf( ioQQQ, " CoolIron NOT calling FeIILevelPops\n"); 00358 } 00359 } 00360 00361 /* evaluate some strong lines that would have been evaluated by the 16 level atom */ 00362 FeIIFillLow16(); 00363 00364 /* now update total cooling and its derivative 00365 * all paths flow through here */ 00366 /* FeII.Fecool = FeII.Fe2_large_cool; */ 00367 CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_large_cool)); 00368 00369 /* add negative cooling to heating stack */ 00370 thermal.heating[25][27] = MAX2(0.,FeII.Fe2_large_heat); 00371 00372 /* counts as heating derivative if negative cooling */ 00373 if( FeII.Fe2_large_cool > 0. ) 00374 { 00375 /* >>chng 01 mar 16, add factor of 3 due to conv problems after changing damper */ 00376 thermal.dCooldT += 3.*FeII.ddT_Fe2_large_cool; 00377 } 00378 00379 if( trace.lgTrace && trace.lgCoolTr ) 00380 { 00381 fprintf( ioQQQ, " Large FeII returns te, cooling, dc=%11.3e%11.3e%11.3e\n", 00382 phycon.te, FeII.Fe2_large_cool, FeII.ddT_Fe2_large_cool ); 00383 } 00384 00385 /* >>chng 05 nov 29, still do simple UV atom if only ground term is done */ 00386 if( !FeII.lgFeIILargeOn ) 00387 { 00388 00389 /* following treatment of Fe II follows 00390 * >>refer fe2 model Wills, B.J., Wills, D., Netzer, H. 1985, ApJ, 288, 143 00391 * all elements are used, and must be set to zero if zero */ 00392 00393 /* set up space for simple model of UV FeII emission */ 00394 if( lgFirst ) 00395 { 00396 /* will never do this again in this core load */ 00397 lgFirst = false; 00398 /* allocate the 1D arrays*/ 00399 pops = (double *)MALLOC( sizeof(double)*(NLFE2) ); 00400 create = (double *)MALLOC( sizeof(double)*(NLFE2) ); 00401 destroy = (double *)MALLOC( sizeof(double)*(NLFE2) ); 00402 depart = (double *)MALLOC( sizeof(double)*(NLFE2) ); 00403 /* create space for the 2D arrays */ 00404 AulPump = ((double **)MALLOC((NLFE2)*sizeof(double *))); 00405 CollRate = ((double **)MALLOC((NLFE2)*sizeof(double *))); 00406 AulDest = ((double **)MALLOC((NLFE2)*sizeof(double *))); 00407 AulEscp = ((double **)MALLOC((NLFE2)*sizeof(double *))); 00408 col_str = ((double **)MALLOC((NLFE2)*sizeof(double *))); 00409 00410 for( i=0; i<(NLFE2); ++i ) 00411 { 00412 AulPump[i] = ((double *)MALLOC((NLFE2)*sizeof(double ))); 00413 CollRate[i] = ((double *)MALLOC((NLFE2)*sizeof(double ))); 00414 AulDest[i] = ((double *)MALLOC((NLFE2)*sizeof(double ))); 00415 AulEscp[i] = ((double *)MALLOC((NLFE2)*sizeof(double ))); 00416 col_str[i] = ((double *)MALLOC((NLFE2)*sizeof(double ))); 00417 } 00418 } 00419 00420 /*zero out all arrays, then check that upper diagonal remains zero below */ 00421 for( i=0; i < NLFE2; i++ ) 00422 { 00423 create[i] = 0.; 00424 destroy[i] = 0.; 00425 for( j=0; j < NLFE2; j++ ) 00426 { 00427 /*data[j][i] = 0.;*/ 00428 col_str[j][i] = 0.; 00429 AulEscp[j][i] = 0.; 00430 AulDest[j][i] = 0.; 00431 AulPump[j][i] = 0.; 00432 } 00433 } 00434 00435 /* now put in real data for lines */ 00436 AulEscp[1][0] = 1.; 00437 AulEscp[2][0] = ( TauLines[ipTuv3].Emis->Pesc + TauLines[ipTuv3].Emis->Pelec_esc)*TauLines[ipTuv3].Emis->Aul; 00438 AulDest[2][0] = TauLines[ipTuv3].Emis->Pdest*TauLines[ipTuv3].Emis->Aul; 00439 AulPump[0][2] = TauLines[ipTuv3].Emis->pump; 00440 00441 AulEscp[5][0] = (TauLines[ipTFe16].Emis->Pesc + TauLines[ipTFe16].Emis->Pelec_esc)*TauLines[ipTFe16].Emis->Aul; 00442 AulDest[5][0] = TauLines[ipTFe16].Emis->Pdest*TauLines[ipTFe16].Emis->Aul; 00443 /* continuum pumping of n=6 */ 00444 AulPump[0][5] = TauLines[ipTFe16].Emis->pump; 00445 /* Ly-alpha pumping */ 00446 00447 double PumpLyaFeII = StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*dense.xIonDense[ipHYDROGEN][1]* 00448 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul* 00449 hydro.dstfe2lya/SDIV(dense.xIonDense[ipIRON][1]); 00450 00451 if( iteration == 1 ) 00452 PumpLyaFeII = 0.; 00453 00454 AulPump[0][5] += PumpLyaFeII; 00455 00456 AulEscp[2][1] = (TauLines[ipTr48].Emis->Pesc + TauLines[ipTr48].Emis->Pelec_esc)*TauLines[ipTr48].Emis->Aul; 00457 AulDest[2][1] = TauLines[ipTr48].Emis->Pdest*TauLines[ipTr48].Emis->Aul; 00458 AulPump[1][2] = TauLines[ipTr48].Emis->pump; 00459 00460 AulEscp[5][1] = (TauLines[ipTFe26].Emis->Pesc + TauLines[ipTFe26].Emis->Pelec_esc)*TauLines[ipTFe26].Emis->Aul; 00461 AulDest[5][1] = TauLines[ipTFe26].Emis->Pdest*TauLines[ipTFe26].Emis->Aul; 00462 AulPump[1][5] = TauLines[ipTFe26].Emis->pump; 00463 00464 AulEscp[3][2] = (TauLines[ipTFe34].Emis->Pesc + TauLines[ipTFe34].Emis->Pelec_esc)*TauLines[ipTFe34].Emis->Aul; 00465 AulDest[3][2] = TauLines[ipTFe34].Emis->Pdest*TauLines[ipTFe34].Emis->Aul; 00466 AulPump[2][3] = TauLines[ipTFe34].Emis->pump; 00467 00468 AulEscp[4][2] = (TauLines[ipTFe35].Emis->Pesc + TauLines[ipTFe35].Emis->Pelec_esc)*TauLines[ipTFe35].Emis->Aul; 00469 AulDest[4][2] = TauLines[ipTFe35].Emis->Pdest*TauLines[ipTFe35].Emis->Aul; 00470 AulPump[2][4] = TauLines[ipTFe35].Emis->pump; 00471 00472 AulEscp[5][3] = (TauLines[ipTFe46].Emis->Pesc + TauLines[ipTFe46].Emis->Pelec_esc)*TauLines[ipTFe46].Emis->Aul; 00473 AulDest[5][3] = TauLines[ipTFe46].Emis->Pdest*TauLines[ipTFe46].Emis->Aul; 00474 AulPump[3][5] = TauLines[ipTFe46].Emis->pump; 00475 00476 AulEscp[5][4] = (TauLines[ipTFe56].Emis->Pesc + TauLines[ipTFe56].Emis->Pelec_esc)*TauLines[ipTFe56].Emis->Aul; 00477 AulDest[5][4] = TauLines[ipTFe56].Emis->Pdest*TauLines[ipTFe56].Emis->Aul; 00478 AulPump[4][5] = TauLines[ipTFe56].Emis->pump; 00479 00480 /* these are collision strengths */ 00481 col_str[1][0] = 1.; 00482 col_str[2][0] = 12.; 00483 col_str[3][0] = 1.; 00484 col_str[4][0] = 1.; 00485 col_str[5][0] = 12.; 00486 col_str[2][1] = 6.; 00487 col_str[3][1] = 1.; 00488 col_str[4][1] = 1.; 00489 col_str[5][1] = 12.; 00490 col_str[3][2] = 6.; 00491 col_str[4][2] = 12.; 00492 col_str[5][2] = 1.; 00493 col_str[4][3] = 1.; 00494 col_str[5][3] = 12.; 00495 col_str[5][4] = 6.; 00496 00497 /*void atom_levelN(long,long,realnum,double[],double[],double[],double*, 00498 double*,double*,long*,realnum*,realnum*,STRING,int*);*/ 00499 atom_levelN(NLFE2, 00500 dense.xIonDense[ipIRON][1], 00501 gFe2, 00502 ex, 00503 'K', 00504 pops, 00505 depart, 00506 &AulEscp , 00507 &col_str, 00508 &AulDest, 00509 &AulPump, 00510 &CollRate, 00511 create, 00512 destroy, 00513 false,/* say atom_levelN should evaluate coll rates from cs */ 00514 /*&&ipdest,*/ 00515 &FeII.Fe2_UVsimp_cool, 00516 &FeII.ddT_Fe2_UVsimp_cool, 00517 "FeII", 00518 &nNegPop, 00519 &lgZeroPop, 00520 false ); 00521 00522 /* nNegPop positive if negative pops occurred, negative if too cold */ 00523 if( nNegPop>0 /*negative if too cold - that is not negative and is OK */ ) 00524 { 00525 fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population for simple UV FeII.\n"); 00526 } 00527 00528 /* add heating - cooling by this process */; 00529 CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_UVsimp_cool)); 00530 thermal.heating[25][27] = MAX2(0.,-FeII.Fe2_UVsimp_cool); 00531 thermal.dCooldT += FeII.ddT_Fe2_UVsimp_cool; 00532 00533 /* LIMLEVELN is the dim of the PopLevels vector */ 00534 ASSERT( NLFE2 <= LIMLEVELN ); 00535 for( i=0; i<NLFE2; ++i) 00536 { 00537 atoms.PopLevels[i] = pops[i]; 00538 atoms.DepLTELevels[i] = depart[i]; 00539 } 00540 00541 TauLines[ipTuv3].Lo->Pop = pops[0]; 00542 TauLines[ipTuv3].Hi->Pop = pops[2]; 00543 TauLines[ipTuv3].Emis->PopOpc = (pops[0] - pops[2]); 00544 TauLines[ipTuv3].Emis->phots = pops[2]*AulEscp[2][0]; 00545 TauLines[ipTuv3].Emis->xIntensity = 00546 TauLines[ipTuv3].Emis->phots*TauLines[ipTuv3].EnergyErg; 00547 00548 TauLines[ipTr48].Lo->Pop = pops[1]; 00549 TauLines[ipTr48].Hi->Pop = pops[2]; 00550 TauLines[ipTr48].Emis->PopOpc = (pops[1] - pops[2]); 00551 TauLines[ipTr48].Emis->phots = pops[2]*AulEscp[2][1]; 00552 TauLines[ipTr48].Emis->xIntensity = 00553 TauLines[ipTr48].Emis->phots*TauLines[ipTr48].EnergyErg; 00554 00555 FeII.for7 = pops[1]*AulEscp[1][0]*4.65e-12; 00556 00557 TauLines[ipTFe16].Lo->Pop = pops[0]; 00558 TauLines[ipTFe16].Hi->Pop = pops[5]; 00559 /* Lyman alpha optical depths are not known on first iteration, 00560 * inward optical depths used, so line trapping overestimated, 00561 * this can cause artificial maser in FeII - prevent by not 00562 * including stimulated emission correction on first iteration */ 00563 if( iteration == 1 ) 00564 { 00565 /* prevent maser due to pumping by Ly-alpha */ 00566 TauLines[ipTFe16].Emis->PopOpc = pops[0]; 00567 } 00568 else 00569 { 00570 TauLines[ipTFe16].Emis->PopOpc = (pops[0] - pops[5]); 00571 } 00572 TauLines[ipTFe16].Emis->phots = pops[5]*AulEscp[5][0]; 00573 TauLines[ipTFe16].Emis->xIntensity = 00574 TauLines[ipTFe16].Emis->phots*TauLines[ipTFe16].EnergyErg; 00575 00576 TauLines[ipTFe26].Lo->Pop = pops[1]; 00577 TauLines[ipTFe26].Hi->Pop = pops[5]; 00578 if( iteration == 1 ) 00579 { 00580 /* prevent maser due to pumping by Ly-alpha */ 00581 TauLines[ipTFe26].Emis->PopOpc = pops[1]; 00582 } 00583 else 00584 { 00585 TauLines[ipTFe26].Emis->PopOpc = (pops[1] - pops[5]); 00586 } 00587 TauLines[ipTFe26].Emis->phots = pops[5]*AulEscp[5][1]; 00588 TauLines[ipTFe26].Emis->xIntensity = 00589 TauLines[ipTFe26].Emis->phots*TauLines[ipTFe26].EnergyErg; 00590 00591 TauLines[ipTFe34].Lo->Pop = pops[2]; 00592 TauLines[ipTFe34].Hi->Pop = pops[3]; 00593 if( iteration == 1 ) 00594 { 00595 /* prevent maser due to pumping by Ly-alpha */ 00596 /* this line is indirectly pumped after decays from 6 to 4 */ 00597 TauLines[ipTFe34].Emis->PopOpc = pops[2]; 00598 } 00599 else 00600 { 00601 TauLines[ipTFe34].Emis->PopOpc = (pops[2] - pops[3]); 00602 } 00603 TauLines[ipTFe34].Emis->phots = pops[3]*AulEscp[3][2]; 00604 TauLines[ipTFe34].Emis->xIntensity = 00605 TauLines[ipTFe34].Emis->phots*TauLines[ipTFe34].EnergyErg; 00606 00607 TauLines[ipTFe35].Lo->Pop = pops[2]; 00608 TauLines[ipTFe35].Hi->Pop = pops[4]; 00609 TauLines[ipTFe35].Emis->PopOpc = (pops[2] - pops[4]); 00610 TauLines[ipTFe35].Emis->phots = pops[4]*AulEscp[4][2]; 00611 TauLines[ipTFe35].Emis->xIntensity = 00612 TauLines[ipTFe35].Emis->phots*TauLines[ipTFe35].EnergyErg; 00613 00614 TauLines[ipTFe46].Lo->Pop = pops[3]; 00615 TauLines[ipTFe46].Hi->Pop = pops[5]; 00616 if( iteration == 1 ) 00617 { 00618 /* prevent maser due to pumping by Ly-alpha */ 00619 TauLines[ipTFe46].Emis->PopOpc = pops[3]; 00620 } 00621 else 00622 { 00623 TauLines[ipTFe46].Emis->PopOpc = (pops[3] - pops[5]); 00624 } 00625 TauLines[ipTFe46].Emis->PopOpc = (pops[3] - pops[5]); 00626 TauLines[ipTFe46].Emis->phots = pops[5]*AulEscp[5][3]; 00627 TauLines[ipTFe46].Emis->xIntensity = 00628 TauLines[ipTFe46].Emis->phots*TauLines[ipTFe46].EnergyErg; 00629 00630 TauLines[ipTFe56].Lo->Pop = pops[4]; 00631 TauLines[ipTFe56].Hi->Pop = pops[5]; 00632 if( iteration == 1 ) 00633 { 00634 /* prevent maser due to pumping by Ly-alpha */ 00635 TauLines[ipTFe56].Emis->PopOpc = pops[4]; 00636 } 00637 else 00638 { 00639 TauLines[ipTFe56].Emis->PopOpc = (pops[4] - pops[5]); 00640 } 00641 TauLines[ipTFe56].Emis->phots = pops[5]*AulEscp[5][4]; 00642 TauLines[ipTFe56].Emis->xIntensity = 00643 TauLines[ipTFe56].Emis->phots*TauLines[ipTFe56].EnergyErg; 00644 00645 /* Jack's funny FeII lines, data from 00646 * >>refer fe2 energy Johansson, S., Brage, T., Leckrone, D.S., Nave, G. & 00647 * >>refercon Wahlgren, G.M. 1995, ApJ 446, 361 */ 00648 PutCS(10.,&TauLines[ipT191]); 00649 atom_level2(&TauLines[ipT191]); 00650 } 00651 00652 { 00653 /*@-redef@*/ 00654 enum{DEBUG_LOC=false}; 00655 /*@+redef@*/ 00656 if( DEBUG_LOC && iteration>1 && nzone>=4 ) 00657 { 00658 fprintf(ioQQQ,"DEBUG2\t%.2e\t%.2e\t%.2e\n", 00659 phycon.te, 00660 FeII.Fe2_large_cool , 00661 FeII.Fe2_UVsimp_cool ); 00662 } 00663 } 00664 00665 return; 00666 00667 } 00668 00669 /*CoolIron - calculation total cooling due to Fe */ 00670 void CoolIron(void) 00671 { 00672 long int i; 00673 double cs , 00674 cs12, cs13, cs23, 00675 cs2s2p, 00676 cs2s3p; 00677 realnum p2, 00678 rate; 00679 00680 static bool lgFe22First=true; 00681 static long int *ipFe22Pump=NULL, 00682 nFe22Pump=0; 00683 double Fe22_pump_rate; 00684 00685 DEBUG_ENTRY( "CoolIron()" ); 00686 00687 /* cooling by FeI 24m, 34.2 m */ 00688 /* >>chng 03 nov 15, add these lines */ 00692 /*>>refer Fe1 cs Hollenbach & McKee 89 */ 00693 /* the 24.0 micron line */ 00694 rate = (realnum)(1.2e-7 * dense.eden + 00695 /* >>chng 05 jul 05, eden to cdsqte */ 00696 /*8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/ 00697 8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]); 00698 LineConvRate2CS( &TauLines[ipFe1_24m] , rate ); 00699 00700 rate = (realnum)(9.3e-8 * dense.eden + 00701 /* >>chng 05 jul 05, eden to cdsqte */ 00702 /*5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/ 00703 5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]); 00704 LineConvRate2CS( &TauLines[ipFe1_35m] , rate ); 00705 00706 rate = (realnum)(1.2e-7 * dense.eden + 00707 /* >>chng 05 jul 05, eden to cdsqte */ 00708 /*6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/ 00709 6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]); 00710 TauDummy.Hi->g = TauLines[ipFe1_35m].Hi->g; 00711 LineConvRate2CS( &TauDummy , rate ); 00712 /* this says that line is a dummy, not real one */ 00713 TauDummy.Hi->g = 0.; 00714 00715 atom_level3(&TauLines[ipFe1_24m],&TauLines[ipFe1_35m],&TauDummy); 00716 00717 /* series of FeI lines from Dima Verner's list, each 2-lev atom 00718 * 00719 * Fe I 3884 */ 00720 MakeCS(&TauLines[ipFeI3884]); 00721 atom_level2(&TauLines[ipFeI3884]); 00722 00723 /* Fe I 3729 */ 00724 MakeCS(&TauLines[ipFeI3729]); 00725 atom_level2(&TauLines[ipFeI3729]); 00726 00727 /* Fe I 3457 */ 00728 MakeCS(&TauLines[ipFeI3457]); 00729 atom_level2(&TauLines[ipFeI3457]); 00730 00731 /* Fe I 3021 */ 00732 MakeCS(&TauLines[ipFeI3021]); 00733 atom_level2(&TauLines[ipFeI3021]); 00734 00735 /* Fe I 2966 */ 00736 MakeCS(&TauLines[ipFeI2966]); 00737 atom_level2(&TauLines[ipFeI2966]); 00738 00739 /* >>chng 05 dec 03, move Fe2 FeII Fe II cooling into separate routine */ 00740 Fe2_cooling(); 00741 00742 /* lump 3p and 3f together; cs= 00743 * >>refer fe3 as Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384 00744 * A from 00745 * >>refer fe3 as Garstang, R.H., 1957, Vistas in Astronomy, 1, 268 00746 * FE III 5270, is 20.9% of total 00747 * >>chng 05 feb 18, Kevin Blagrave email 00748 * average wavelength is 4823 with statistical weight averaging of upper energy level, 00749 * as per , change 5th number from 2.948 to 2.984, also photon energy 00750 * from 3.78 to 4.12 */ 00751 00752 /* >>chng 05 dec 16, FeIII update by Kevin Blagrave */ 00753 /*CoolHeavy.c5270 = atom_pop2(5.53,25.,30.,0.435,2.984e4,dense.xIonDense[ipIRON][2])* 00754 4.12e-12; 00755 CoolAdd("Fe 3",5270,CoolHeavy.c5270);*/ 00756 /* FeIII 1122 entire multiplet - atomic data=A's Dima, CS = guess */ 00757 PutCS(25.,&TauLines[ipT1122]); 00758 atom_level2(&TauLines[ipT1122]); 00759 00760 /* call 14 level atom */ 00761 Fe3Lev14(); 00762 00763 /* call 12 level atom */ 00764 Fe4Lev12(); 00765 00766 /* FE V 3892 + 3839, data from Shields */ 00767 CoolHeavy.c3892 = atom_pop2(7.4,25.,5.,0.6,3.7e4,dense.xIonDense[ipIRON][4])* 00768 5.11e-12; 00769 CoolAdd("Fe 5",3892,CoolHeavy.c3892); 00770 00771 /* FE VI 5177+5146+4972+4967 00772 * data from 00773 * >>refer fe6 as Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384 */ 00774 CoolHeavy.c5177 = atom_pop2(1.9,28.,18.,0.52,2.78e4,dense.xIonDense[ipIRON][5])* 00775 3.84e-12; 00776 CoolAdd("Fe 6",5177,CoolHeavy.c5177); 00777 00778 /* >>chng 04 nov 04, add multi-level fe7 */ 00779 Fe7Lev8(); 00780 00781 /* Fe IX 245,242 00782 * all atomic data from 00783 * >>refer fe9 all Flower, D.R. 1976, A&A, 56, 451 */ 00784 /* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */ 00785 PutCS(0.123,&TauLines[ipT245]); 00786 /* AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41) 00787 * C245 = AtomSeqBeryllium( .087 ,.038 ,.188 , t245 ,71. ) * 8.14E-11 */ 00788 AtomSeqBeryllium(.087,.038,.188,&TauLines[ipT245],71.); 00789 00790 CoolHeavy.c242 = atoms.PopLevels[3]*8.22e-11*71.; 00791 00792 /* Fe X Fe 10 Fe10 6374, A from 00793 * >>referold fe10 as Mason, H. 1975, MNRAS 170, 651 00794 * >>referold fe10 cs Mohan, M., Hibbert, A., & Kingston, A.E. 1994, ApJ, 434, 389 00795 * >>referold fe10 as Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209 00796 * does not agree with Mohan et al. by factor of 4 00797 * 20x larger than Mason numbers used pre 1994 */ 00798 /*cs = 13.3-2.*MIN2(7.0,phycon.alogte); */ 00799 /*cs = MAX2(0.1,cs); */ 00805 /*cs = 1.;*/ 00806 /* >>chng 00 dec 05, update Fe10 collisions strength to Tayal 2000 */ 00807 /* >>refer fe10 cs Tayal, S.S., 2000, ApJ, 544, 575-580 */ 00808 /* >>chng 04 nov 10, 172.9 was mult rather than div by temp law, 00809 * had no effect due to min on cs that lie below it */ 00810 /*cs = 172.9 /( phycon.te30 * phycon.te05 * phycon.te02 * phycon.te005 );*/ 00811 /* tabulated cs ends at 30,000K, do not allow cs to grow larger than largest 00812 * tabulated value */ 00813 /*cs = MIN2( cs, 3.251 );*/ 00814 /* >>chng 05 dec 19, update As to 00815 * >>refer fe10 As Aggarwal, K.M. & Keenan, F.P. 2004, A&A, 427, 763 00816 * value changed from 54.9 to 68.9 */ 00817 /* >>chng 05 dec 19, update cs to 00818 *>>refer fe10 cs Aggarwal, K.M. & Keenan, F.P. 2005, A&A, 439, 1215 */ 00819 cs = Fe_10_11_13_cs( 00820 /* ion, one of 10, 11, or 13 on physics scale */ 00821 10, 00822 /* initial and final index of levels, with lowest energy 1, highest of 5 00823 * on physics scale */ 00824 1, 00825 2 ); 00826 00827 PutCS(cs,&TauLines[ipFe106375]); 00828 atom_level2(&TauLines[ipFe106375]); 00829 /* c6374 = atom_pop2(cs ,4.,2.,69.5,2.26e4,dense.xIonDense(26,10))*3.12e-12 00830 * dCooldT = dCooldT + c6374 * 2.26e4 * tsq1 00831 * call CoolAdd( 'Fe10' , 6374 , c6374 ) 00832 * 00833 * this is the E1 line that can pump [FeX] */ 00834 cs = 0.85*sexp(0.045*1.259e6/phycon.te); 00835 cs = MAX2(0.05,cs); 00836 PutCS(cs,&TauLines[ipT352]); 00837 atom_level2(&TauLines[ipT352]); 00838 00839 /* Fe XI fe11 fe 11, 7892, 6.08 mic, CS from 00840 * >>refer fe11 as Mason, H. 1975, MNRAS 170, 651 00841 * A from 00842 * >>refer fe11 as Mendoza, C., & Zeippen, C.J. 1983, MNRAS, 202, 981 00843 * following reference has very extensive energy levels and As 00844 * >>refer fe11 as Fritzsche, S., Dong, C.Z., & Traebert, E., 2000, MNRAS, 318, 263 */ 00845 /*cs = 0.27;*/ 00846 /* >>chng 97 jan 31, set cs to unity, above ignored resonances */ 00847 /*cs = 1.;*/ 00848 /* >>chng 00 dec 05, update Fe11 CS to Tayal 2000 */ 00849 /* >>referold fe11 cs Tayal, S.S., 2000, ApJ, 544, 575-580 */ 00850 /* this is the low to mid transition */ 00851 /*cs = 0.0282 * phycon.te30*phycon.te05*phycon.te01;*/ 00852 /* CS is about 2.0 across broad temp range in following reference: 00853 * >>refer Fe11 cs Aggarwal, K.M., & Keenan, F.P. 2003, MNRAS, 338, 412 */ 00854 /*cs = 2.; 00855 PutCS(cs,&TauLines[ipTFe07]);*/ 00856 /* Tayal value is 10x larger than previous values */ 00857 /* Aggarwal & Keenan is about same as Tayal */ 00858 /*cs = 0.48; */ 00859 /*cs = 0.50; 00860 PutCS(cs,&TauLines[ipTFe61]);*/ 00861 /* Tayal value is 3x larger than previous values */ 00862 /*cs = 0.35; tayal number */ 00863 /*cs = 0.55; 00864 PutCS(cs,&TauDummy); 00865 atom_level3(&TauLines[ipTFe07],&TauLines[ipTFe61],&TauDummy);*/ 00866 00867 /* >>refer fe11 cs Kafatos, M., & Lynch, J.P. 1980, ApJS, 42, 611 */ 00868 /*CoolHeavy.c1467 = atom_pop3(9.,5.,1.,.24,.028,.342,100.,1012.,9.3,5.43e4, 00869 6.19e4,&p2,dense.xIonDense[ipIRON][11-1],0.,0.,0.)*1012.*1.36e-11; 00870 CoolHeavy.c2649 = p2*91.0*7.512e-12;*/ 00871 /*CoolAdd("Fe11",1467,CoolHeavy.c1467); 00872 CoolAdd("Fe11",2469,CoolHeavy.c2649);*/ 00873 00874 /* >>chng 05 dec 18, use give level Fe 11 model */ 00875 Fe11Lev5(); 00876 00877 /* A's for 2-1 and 3-1 from 00878 * >>refer fe12 as Tayal, S.S., & Henry, R.J.W. 1986, ApJ, 302, 200 00879 * CS fro 00880 * >>refer fe12 cs Tayal, S.S., Henry, R.J.W., Pradhan, A.K. 1987, ApJ, 319, 951 00881 * a(3-2) scaled from both ref */ 00882 CoolHeavy.c2568 = atom_pop3(4.,10.,6.,0.72,0.69,2.18,8.1e4,1.84e6,1.33e6, 00883 6.37e4,4.91e4,&p2,dense.xIonDense[ipIRON][12-1],0.,0.,0.)*1.33e6*6.79e-12; 00884 CoolAdd("Fe12",2568,CoolHeavy.c2568); 00885 CoolHeavy.c1242 = CoolHeavy.c2568*2.30*1.38; 00886 CoolAdd("Fe12",1242,CoolHeavy.c1242); 00887 CoolHeavy.c2170 = p2*8.09e4*8.82e-12; 00888 CoolAdd("Fe12",2170,CoolHeavy.c2170); 00889 00890 /* [Fe XIII] fe13 fe 13 1.07, 1.08 microns */ 00891 /* >>chng 00 dec 05, update Fe13 CS to Tayal 2000 */ 00892 /* >>refer fe13 cs Tayal, S.S., 2000, ApJ, 544, 575-580 00893 cs = 5.08e-3 * phycon.te30* phycon.te10; 00894 PutCS(cs,&TauLines[ipFe1310]); 00895 PutCS(2.1,&TauLines[ipFe1311]); 00896 PutCS(.46,&TauDummy); 00897 atom_level3(&TauLines[ipFe1310],&TauLines[ipFe1311],&TauDummy); */ 00898 00899 /* Fe13 lines from 1D and 1S - 00900 >>chng 01 aug 07, added these lines */ 00901 /* >>refer fe11 cs Tayal, S.S., 2000, ApJ, 544, 575-580 */ 00902 /* >>refer fe13 as Shirai, T., Sugar, J., Musgrove, A., & Wiese, W.L., 2000., 00903 * >>refercon J Phys Chem Ref Data Monograph 8 */ 00904 /* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) 00905 CoolHeavy.fe13_1216 = atom_pop3( 9.,5.,1., 5.08,0.447,0.678, 103.2 ,1010.,7.99, 00906 48068.,43440.,&p2,dense.xIonDense[ipIRON][13-1],0.,0.,0.)*1010.* 00907 1.64e-11; 00908 CoolHeavy.fe13_2302 = CoolHeavy.fe13_1216*0.528*0.00791; 00909 CoolHeavy.fe13_3000 = p2*103.2*7.72e-12;*/ 00910 /*CoolAdd("Fe13",1216,CoolHeavy.fe13_1216); 00911 CoolAdd("Fe13",3000,CoolHeavy.fe13_3000); 00912 CoolAdd("Fe13",2302,CoolHeavy.fe13_2302);*/ 00913 00914 /* >>chng 05 dec 18, use give level Fe 13 model */ 00915 Fe13Lev5(); 00916 00917 /* Fe XIV 5303, cs from 00918 * >>refer Fe14 cs Storey, P.J., Mason, H.E., Saraph, H.E., 1996, A&A, 309, 677 */ 00919 fe14cs(phycon.alogte,&cs); 00920 /* >>chng 97 jan 31, set cs to unity, above is VERY large, >10 */ 00921 /* >>chng 01 may 30, back to their value, as per discussion at Lexington 2000 */ 00922 /* >>chng 05 aug 04, following line commented out, had set 5303 to constant value */ 00923 /*cs = 3.1;*/ 00924 /* >>chng 05 dec 22, A from 60.3 to 59.7, experimental value 59.7 +/- 0.4 in 00925 * >>refer Fe14 as Trabert, E. 2004, A&A, 415, L39 */ 00926 CoolHeavy.c5303 = atom_pop2(cs,2.,4.,59.7,2.71e4,dense.xIonDense[ipIRON][14-1])* 00927 3.75e-12; 00928 thermal.dCooldT += CoolHeavy.c5303*2.71e4*thermal.tsq1; 00929 CoolAdd("Fe14",5303,CoolHeavy.c5303); 00930 00931 /* Fe 18 974.8A;, cs from 00932 * >>referold fe18 cs Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */ 00933 /* >>refer fe18 cs Berrington,K.A.,Saraph, H.E. & Tully, J.A. 1998, A&AS, 129, 161 */ 00934 /*>>chng 06 jul 19 Changes made-Humeshkar Nemala*/ 00935 /*cs = MIN2(0.143,0.804/(phycon.te10*phycon.te05));*/ 00936 if(phycon.te < 1.29E6) 00937 { 00938 cs = (realnum)(0.3465/((phycon.te10/phycon.te01)*phycon.te001*phycon.te0003)); 00939 } 00940 else if(phycon.te < 5.135E6) 00941 { 00942 cs = (realnum)((1.1062E-02)*phycon.te10*phycon.te05*phycon.te003*phycon.te0005); 00943 } 00944 else 00945 { 00946 cs = (realnum)((60.5728)/(phycon.te40*phycon.te003*phycon.te0001*phycon.te0005)); 00947 } 00948 00949 PutCS(cs,&TauLines[ipFe18975]); 00950 atom_level2(&TauLines[ipFe18975]); 00951 00952 /* O-like Fe19, 3P ground term, 7046.72A vac wl, 1328.90A */ 00953 /* >>refer fe19 cs Butler, K., & Zeippen, C.J., 2001, A&A, 372, 1083 */ 00954 cs12 = 0.0627 / phycon.te03; 00955 cs13 = 0.692 /(phycon.te10*phycon.te01); 00956 cs23 = 0.04; 00957 /*CoolHeavy.c7082 = atom_pop3(5.,1.,3.,0.0132,0.0404,0.0137,0.505,1.46e4,*/ 00958 /* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */ 00959 CoolHeavy.c7082 = atom_pop3(5.,1.,3.,cs12,cs13,cs23,0.505,1.46e4, 00960 41.2,1.083e5,2.03e4,&p2,dense.xIonDense[ipIRON][19-1],0.,0.,0.)*41.2* 00961 2.81e-12; 00962 CoolHeavy.c1118 = CoolHeavy.c7082*354.4*6.335; 00963 CoolHeavy.c1328 = p2*0.505*1.50e-11; 00964 CoolAdd("Fe19",7082,CoolHeavy.c7082); 00965 CoolAdd("Fe19",1118,CoolHeavy.c1118); 00966 CoolAdd("Fe19",1328,CoolHeavy.c1328); 00967 00968 CoolHeavy.c592 = atom_pop2(0.0913,9.,5.,1.64e4,2.428e5,dense.xIonDense[ipIRON][19-1])* 00969 3.36e-11; 00970 CoolAdd("Fe19",592,CoolHeavy.c592); 00971 00972 /* >>chng 01 aug 10, add this line */ 00973 /* Fe20 721.40A, 578*/ 00974 /* >>refer fe20 cs Butler, K., & Zeippen, C.J., 2001, A&A, 372, 1078 */ 00975 /*>>refer fe20 as Merkelis, G., Martinson, I., Kisielius, R., & Vilkas, M.J., 1999, 00976 >>refercon Physica Scripta, 59, 122 */ 00977 cs = 1.17 /(phycon.te20/phycon.te01); 00978 PutCS(cs , &TauLines[ipTFe20_721]); 00979 cs = 0.248 /(phycon.te10/phycon.te01); 00980 PutCS(cs , &TauLines[ipTFe20_578]); 00981 cs = 0.301 /(phycon.te10/phycon.te02); 00982 PutCS(cs , &TauDummy); 00983 atom_level3(&TauLines[ipTFe20_721],&TauLines[ipTFe20_578],&TauDummy); 00984 00985 /* >>chng 00 oct 26, much larger CS for following transition */ 00986 /* Fe 21 1354, 2299 A, cs eval at 1e6K 00987 * >>refer Fe21 cs Aggarwall, K.M., 1991, ApJS 77, 677 */ 00988 PutCS(0.072,&TauLines[ipTFe13]); 00989 PutCS(0.269,&TauLines[ipTFe23]); 00990 PutCS(0.055,&TauDummy); 00991 atom_level3(&TauLines[ipTFe13],&TauLines[ipTFe23],&TauDummy); 00992 00993 /*>>refer Fe22 energy Feldman, U., Curdt, W., Landi, E., Wilhelm, K., 2000, ApJ, 544, 508 */ 00994 /*>>chng 00 oct 26, added these fe 21 lines, removed old CoolHeavy.fs846, the lowest */ 00995 00996 /* one time initialization if first call, and level 2 lines are on */ 00997 if( lgFe22First && nWindLine ) 00998 { 00999 lgFe22First = false; 01000 nFe22Pump = 0; 01001 for( i=0; i<nWindLine; ++i ) 01002 { 01003 /* don't test on nelem==ipIRON since lines on physics, not C, scale */ 01004 if( TauLine2[i].Hi->nelem ==26 && TauLine2[i].Hi->IonStg==22 ) 01005 { 01006 ++nFe22Pump; 01007 } 01008 } 01009 if( nFe22Pump<0 ) 01010 TotalInsanity(); 01011 else if( nFe22Pump > 0 ) 01012 /* create the space - can't malloc 0 bytes */ 01013 ipFe22Pump = (long *)MALLOC((unsigned)(nFe22Pump)*sizeof(long) ); 01014 nFe22Pump = 0; 01015 for( i=0; i<nWindLine; ++i ) 01016 { 01017 /* don't test on nelem==ipIRON since lines on physics, not C, scale */ 01018 if( TauLine2[i].Hi->nelem ==26 && TauLine2[i].Hi->IonStg==22 ) 01019 { 01020 # if 0 01021 DumpLine( &TauLine2[i] ); 01022 # endif 01023 ipFe22Pump[nFe22Pump] = i; 01024 ++nFe22Pump; 01025 } 01026 } 01027 } 01028 else 01029 /* level 2 lines are not enabled */ 01030 nFe22Pump = 0; 01031 01032 /* now sum pump rates */ 01033 Fe22_pump_rate = 0.; 01034 for( i=0; i<nFe22Pump; ++i ) 01035 { 01036 Fe22_pump_rate += TauLine2[ipFe22Pump[i]].Emis->pump; 01037 # if 0 01038 fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n", 01039 i, 01040 TauLine2[ipFe22Pump[i]].WLAng , TauLine2[ipFe22Pump[i]].pump ); 01041 # endif 01042 } 01043 01044 /*AtomSeqBoron compute cooling from 5-level boron sequence model atom */ 01045 /*>>refer fe22 cs Zhang, H.L., Graziani, M., & Pradhan, A.K., 1994, A&A 283, 319*/ 01046 /*>>refer fe22 as Dankwort, W., & Trefftz, E., 1978, A&A 65, 93-98 */ 01047 AtomSeqBoron(&TauLines[ipFe22_846], 01048 &TauLines[ipFe22_247], 01049 &TauLines[ipFe22_217], 01050 &TauLines[ipFe22_348], 01051 &TauLines[ipFe22_292], 01052 &TauLines[ipFe22_253], 01053 /*double cs40, 01054 double cs32, 01055 double cs42, 01056 double cs43, */ 01057 0.00670 , 0.0614 , 0.0438 , 0.122 , Fe22_pump_rate ,"Fe22"); 01058 01059 /* Fe 22 845.6, C.S.=guess, A from 01060 * >>refer fe22 as Froese Fischer, C. 1983, J.Phys. B, 16, 157 01061 CoolHeavy.fs846 = atom_pop2(0.2,2.,4.,1.5e4,1.701e5,dense.xIonDense[ipIRON][22-1])* 01062 2.35e-11; 01063 CoolAdd("Fe22",846,CoolHeavy.fs846);*/ 01064 01066 /* FE 23 262.6, 287A, 1909-LIKE, 01067 * cs from 01068 * >>refer fe23 cs Bhatia, A.K., & Mason, H.E. 1986, A&A, 155, 413 */ 01069 CoolHeavy.c263 = atom_pop2(0.04,1.,9.,1.6e7,5.484e5,dense.xIonDense[ipIRON][23-1])* 01070 7.58e-11; 01071 CoolAdd("Fe23",262,CoolHeavy.c263); 01072 01073 /* FE 24, 255, 192 Li seq doublet 01074 * >>refer fe24 cs Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */ 01075 ligbar(26,&TauLines[ipT192],&TauLines[ipT11],&cs2s2p,&cs2s3p); 01076 01077 PutCS(cs2s2p,&TauLines[ipT192]); 01078 atom_level2(&TauLines[ipT192]); 01079 01080 /* funny factor (should have been 0.5) due to energy change */ 01081 PutCS(cs2s2p*0.376,&TauLines[ipT255]); 01082 atom_level2(&TauLines[ipT255]); 01083 01084 PutCS(cs2s3p,&TauLines[ipT11]); 01085 atom_level2(&TauLines[ipT11]); 01086 01088 TauLines[ipT353].Emis->PopOpc = dense.xIonDense[ipIRON][11-1]; 01089 TauLines[ipT353].Lo->Pop = dense.xIonDense[ipIRON][11-1]; 01090 TauLines[ipT353].Hi->Pop = 0.; 01091 TauLines[ipT347].Emis->PopOpc = dense.xIonDense[ipIRON][14-1]; 01092 TauLines[ipT347].Lo->Pop = dense.xIonDense[ipIRON][14-1]; 01093 TauLines[ipT347].Hi->Pop = 0.; 01094 01095 return; 01096 } 01097 01098 /*fe14cs compute collision strength for forbidden transition 01099 * w/in Fe XIV ground term. From 01100 * >>refer fe14 cs Storey, P.J., Mason, H.E., Saraph, H.E., 1996, A&A, 309, 677 01101 * */ 01102 STATIC void fe14cs(double te1, 01103 double *csfe14) 01104 { 01105 double a, 01106 b, 01107 c, 01108 d, 01109 telog1, 01110 telog2, 01111 telog3; 01112 01113 DEBUG_ENTRY( "fe14cs()" ); 01114 01115 /* limit range in log T: */ 01116 telog1 = te1; 01117 telog1 = MIN2(9.0,telog1); 01118 telog1 = MAX2(4.0,telog1); 01119 01120 /* compute square and cube */ 01121 telog2 = telog1*telog1; 01122 telog3 = telog2*telog1; 01123 01124 /* compute cs: 01125 * */ 01126 if( telog1 <= 5.0 ) 01127 { 01128 a = 557.05536; 01129 b = -324.56109; 01130 c = 63.437974; 01131 d = -4.1365147; 01132 *csfe14 = a + b*telog1 + c*telog2 + d*telog3; 01133 } 01134 else 01135 { 01136 a = 0.19515493; 01137 b = 2.9404407; 01138 c = 4.9578944; 01139 d = 0.79887506; 01140 *csfe14 = a + b*exp(-0.5*((telog1-c)*(telog1-c)/d)); 01141 } 01142 return; 01143 } 01144 01145 /*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */ 01146 STATIC void Fe4Lev12(void) 01147 { 01148 const int NLFE4 = 12; 01149 bool lgZeroPop; 01150 int nNegPop; 01151 long int i, 01152 j; 01153 static bool lgFirst=true; 01154 01155 double dfe4dt; 01156 01157 /*static long int **ipdest; */ 01158 static double 01159 **AulEscp, 01160 **col_str, 01161 **AulDest, 01162 depart[NLFE4], 01163 pops[NLFE4], 01164 destroy[NLFE4], 01165 create[NLFE4], 01166 **CollRate, 01167 **AulPump; 01168 01169 static const double Fe4A[NLFE4][NLFE4] = { 01170 {0.,0.,0.,1.e-5,0.,1.368,.89,0.,1.3e-3,1.8e-4,.056,.028}, 01171 {0.,0.,2.6e-8,0.,0.,0.,0.,0.,1.7e-7,0.,0.,0.}, 01172 {0.,0.,0.,0.,3.5e-7,6.4e-10,0.,0.,6.315e-4,0.,6.7e-7,0.}, 01173 {0.,0.,0.,0.,1.1e-6,6.8e-5,8.6e-6,3.4e-10,7.6e-5,1.e-7,5.8e-4,2.8e-4}, 01174 {0.,0.,0.,0.,0.,1.5e-5,1.3e-9,0.,7.6e-4,0.,1.1e-6,6.0e-7}, 01175 {0.,0.,0.,0.,0.,0.,1.1e-5,1.2e-13,.038,9.9e-7,.022,.018}, 01176 {0.,0.,0.,0.,0.,0.,0.,3.7e-5,2.9e-6,.034,3.5e-3,.039}, 01177 {0.,0.,0.,0.,0.,0.,0.,0.,0.,.058,3.1e-6,1.4e-3}, 01178 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-4,3.1e-14}, 01179 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.9e-19,1.0e-5}, 01180 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-7}, 01181 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.} 01182 }; 01183 static const double Fe4CS[NLFE4][NLFE4] = { 01184 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, 01185 {0.98,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, 01186 {0.8167,3.72,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, 01187 {0.49,0.0475,0.330,0.,0.,0.,0.,0.,0.,0.,0.,0.}, 01188 {0.6533,0.473,2.26,1.64,0.,0.,0.,0.,0.,0.,0.,0.}, 01189 {0.45,0.686,0.446,0.106,0.254,0.,0.,0.,0.,0.,0.,0.}, 01190 {0.30,0.392,0.152,0.269,0.199,0.605,0.,0.,0.,0.,0.,0.}, 01191 {0.15,0.0207,0.190,0.0857,0.166,0.195,0.327,0.,0.,0.,0.,0.}, 01192 {0.512,1.23,0.733,0.174,0.398,0.623,0.335,0.102,0.,0.,0.,0.}, 01193 {0.128,0.0583,0.185,0.200,0.188,0.0835,0.127,0.0498,0.0787,0.,0.,0.}, 01194 {0.384,0.578,0.534,0.363,0.417,0.396,0.210,0.171,0.810,0.101,0.,0.}, 01195 {0.256,0.234,0.306,0.318,0.403,0.209,0.195,0.112,0.195,0.458,0.727,0.} 01196 }; 01197 01198 static const double gfe4[NLFE4]={6.,12.,10.,6.,8.,6.,4.,2.,8.,2.,6.,4.}; 01199 01200 /* excitation energies in Kelvin 01201 static double ex[NLFE4]={0.,46395.8,46464.,46475.6,46483.,50725., 01202 50838.,50945.,55796.,55966.,56021.,56025.};*/ 01203 /*>>refer Fe3 energies version 3 NIST Atomic Spectra Database */ 01204 /*>>chng 05 dec 17, from Kelvin above to excitation energies in wn */ 01205 static const double excit_wn[NLFE4]={0.,32245.5,32292.8,32301.2,32305.7,35253.8, 01206 35333.3,35406.6,38779.4,38896.7,38935.1,38938.2}; 01207 01208 DEBUG_ENTRY( "Fe4Lev12()" ); 01209 01210 if( lgFirst ) 01211 { 01212 /* will never do this again */ 01213 lgFirst = false; 01214 /* allocate the 1D arrays*/ 01215 /* create space for the 2D arrays */ 01216 AulPump = ((double **)MALLOC((NLFE4)*sizeof(double *))); 01217 CollRate = ((double **)MALLOC((NLFE4)*sizeof(double *))); 01218 AulDest = ((double **)MALLOC((NLFE4)*sizeof(double *))); 01219 AulEscp = ((double **)MALLOC((NLFE4)*sizeof(double *))); 01220 col_str = ((double **)MALLOC((NLFE4)*sizeof(double *))); 01221 for( i=0; i<(NLFE4); ++i ) 01222 { 01223 AulPump[i] = ((double *)MALLOC((NLFE4)*sizeof(double ))); 01224 CollRate[i] = ((double *)MALLOC((NLFE4)*sizeof(double ))); 01225 AulDest[i] = ((double *)MALLOC((NLFE4)*sizeof(double ))); 01226 AulEscp[i] = ((double *)MALLOC((NLFE4)*sizeof(double ))); 01227 col_str[i] = ((double *)MALLOC((NLFE4)*sizeof(double ))); 01228 } 01229 } 01230 01231 /* bail if no Fe+3 */ 01232 if( dense.xIonDense[ipIRON][3] <= 0. ) 01233 { 01234 fe.Fe4CoolTot = 0.; 01235 fe.fe40401 = 0.; 01236 fe.fe42836 = 0.; 01237 fe.fe42829 = 0.; 01238 fe.fe42567 = 0.; 01239 fe.fe41207 = 0.; 01240 fe.fe41206 = 0.; 01241 fe.fe41106 = 0.; 01242 fe.fe41007 = 0.; 01243 fe.fe41008 = 0.; 01244 fe.fe40906 = 0.; 01245 CoolAdd("Fe 4",0,0.); 01246 01247 /* level populations */ 01248 /* LIMLEVELN is the dimension of the atoms vectors */ 01249 ASSERT( NLFE4 <= LIMLEVELN); 01250 for( i=0; i < NLFE4; i++ ) 01251 { 01252 atoms.PopLevels[i] = 0.; 01253 atoms.DepLTELevels[i] = 1.; 01254 } 01255 return; 01256 } 01257 /* number of levels in model ion */ 01258 01259 /* these are in wavenumbers 01260 * data excit_wn/ 0., 32245.5, 32293., 32301.2, 01261 * 1 32306., 35254., 35333., 35407., 38779., 38897., 38935., 01262 * 2 38938./ 01263 * excitation energies in Kelvin */ 01264 01265 /* A's are from Garstang, R.H., MNRAS 118, 572 (1958). 01266 * each set is for a lower level indicated by second element in array, 01267 * index runs over upper level 01268 * A's are saved into arrays as data(up,lo) */ 01269 01270 /* collision strengths from Berrington and Pelan Ast Ap S 114, 367. 01271 * order is cs(low,up) */ 01272 01273 /* all elements are used, and must be set to zero if zero */ 01274 for( i=0; i < NLFE4; i++ ) 01275 { 01276 create[i] = 0.; 01277 destroy[i] = 0.; 01278 for( j=0; j < NLFE4; j++ ) 01279 { 01280 /*data[j][i] = 1e33;*/ 01281 col_str[j][i] = 0.; 01282 AulEscp[j][i] = 0.; 01283 AulDest[j][i] = 0.; 01284 AulPump[j][i] = 0.; 01285 } 01286 } 01287 01288 /* fill in einstein As and collision strengths */ 01289 for( i=0; i < NLFE4; i++ ) 01290 { 01291 for( j=i + 1; j < NLFE4; j++ ) 01292 { 01293 /*data[i][j] = Fe4A[i][j];*/ 01294 AulEscp[j][i] = Fe4A[i][j]; 01295 /*data[j][i] = Fe4CS[j][i];*/ 01296 col_str[j][i] = Fe4CS[j][i]; 01297 } 01298 } 01299 01300 /* leveln itself is well-protected against zero abundances, 01301 * low temperatures */ 01302 01303 atom_levelN(NLFE4, 01304 dense.xIonDense[ipIRON][3], 01305 gfe4, 01306 excit_wn, 01307 'w', 01308 pops, 01309 depart, 01310 &AulEscp , 01311 &col_str , 01312 &AulDest, 01313 &AulPump, 01314 &CollRate, 01315 create, 01316 destroy, 01317 /* say atom_levelN should evaluate coll rates from cs */ 01318 false, 01319 &fe.Fe4CoolTot, 01320 &dfe4dt, 01321 "FeIV", 01322 /* nNegPop positive if negative pops occured, negative if too cold */ 01323 &nNegPop, 01324 &lgZeroPop, 01325 false ); 01326 01327 /* LIMLEVELN is the dim of the PopLevels vector */ 01328 ASSERT( NLFE4 <= LIMLEVELN ); 01329 for( i=0; i<NLFE4; ++i) 01330 { 01331 atoms.PopLevels[i] = pops[i]; 01332 atoms.DepLTELevels[i] = depart[i]; 01333 } 01334 01335 if( nNegPop>0 ) 01336 { 01337 fprintf( ioQQQ, " fe4levl2 found negative populations\n" ); 01338 ShowMe(); 01339 cdEXIT(EXIT_FAILURE); 01340 } 01341 01342 CoolAdd("Fe 4",0,fe.Fe4CoolTot); 01343 01344 thermal.dCooldT += dfe4dt; 01345 01346 /* three transitions hst can observe 01347 * 4 -1 3095.9A and 5 -1 3095.9A */ 01348 fe.fe40401 = (pops[3]*Fe4A[0][3]*(excit_wn[3] - excit_wn[0]) + 01349 pops[4]*Fe4A[0][4]*(excit_wn[4] - excit_wn[0]) )*T1CM*BOLTZMANN; 01350 01351 fe.fe42836 = pops[5]*Fe4A[0][5]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN; 01352 01353 fe.fe42829 = pops[6]*Fe4A[0][6]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN; 01354 01355 fe.fe42567 = (pops[10]*Fe4A[0][10]*(excit_wn[10] - excit_wn[0]) + 01356 pops[11]*Fe4A[0][11]*(excit_wn[10] - excit_wn[0]))*T1CM*BOLTZMANN; 01357 01358 fe.fe41207 = pops[11]*Fe4A[6][11]*(excit_wn[11] - excit_wn[6])*T1CM*BOLTZMANN; 01359 fe.fe41206 = pops[11]*Fe4A[5][11]*(excit_wn[11] - excit_wn[5])*T1CM*BOLTZMANN; 01360 fe.fe41106 = pops[10]*Fe4A[5][10]*(excit_wn[10] - excit_wn[5])*T1CM*BOLTZMANN; 01361 fe.fe41007 = pops[9]*Fe4A[6][9]*(excit_wn[9] - excit_wn[6])*T1CM*BOLTZMANN; 01362 fe.fe41008 = pops[9]*Fe4A[7][9]*(excit_wn[9] - excit_wn[7])*T1CM*BOLTZMANN; 01363 fe.fe40906 = pops[8]*Fe4A[5][8]*(excit_wn[8] - excit_wn[5])*T1CM*BOLTZMANN; 01364 return; 01365 } 01366 01367 /*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */ 01368 STATIC void Fe7Lev8(void) 01369 { 01370 bool lgZeroPop; 01371 int nNegPop; 01372 double scale; 01373 long int i, 01374 j; 01375 static bool lgFirst=true; 01376 static long int ipPump=-1; 01377 01378 double dfe7dt, 01379 FUV_pump; 01380 01381 long int ihi , ilo; 01382 static double 01383 *depart, 01384 *pops, 01385 *destroy, 01386 *create , 01387 **AulDest, 01388 **CollRate, 01389 **AulPump, 01390 **AulNet, 01391 **col_str; 01392 /* 01393 following are FeVII lines predicted in limit_laser_2000 01394 Fe 7 5721A -24.399 0.1028 01395 Fe 7 4989A -24.607 0.0637 01396 Fe 7 4894A -24.838 0.0374 01397 Fe 7 4699A -25.693 0.0052 01398 Fe 7 6087A -24.216 0.1566 01399 Fe 7 5159A -24.680 0.0538 01400 Fe 7 4943A -25.048 0.0231 01401 Fe 7 3587A -24.285 0.1336 01402 Fe 7 5277A -25.053 0.0228 01403 Fe 7 3759A -24.139 0.1870 01404 Fe 7 3.384m -25.621 0.0062 01405 Fe 7 2.629m -25.357 0.0113 01406 Fe 7 9.510m -24.467 0.0878 01407 Fe 7 7.810m -24.944 0.0293 01408 */ 01409 01410 /* statistical weights for the n levels */ 01411 static double gfe7[NLFE7]={5.,7.,9.,5.,1.,3.,5.,9.}; 01412 01413 /*refer Fe7 energies Ekbert, J.O. 1981, Phys. Scr 23, 7 */ 01414 /*static double ex[NLFE7]={0.,1047.,2327.,17475.,20037.,20428.,21275.,28915.};*/ 01415 /* excitation energies in wavenumbers 01416 * >>chng 05 dec 15, rest of code had assumed that there were energies in Kelvin 01417 * rather than wavenumbers. Bug caught by Kevin Blagrave 01418 * modified atom_levelN to accept either Kelvin or wavenumbers */ 01419 static double excit_wn[NLFE7]={0. , 1051.5 , 2331.5 , 17475.5 , 20040.3 , 20430.1 , 21278.6 , 28927.3 }; 01420 01421 DEBUG_ENTRY( "Fe7Lev8()" ); 01422 01423 if( lgFirst ) 01424 { 01425 /* will never do this again */ 01426 lgFirst = false; 01427 /* allocate the 1D arrays*/ 01428 /* create space for the 2D arrays */ 01429 depart = ((double *)MALLOC((NLFE7)*sizeof(double))); 01430 pops = ((double *)MALLOC((NLFE7)*sizeof(double))); 01431 destroy = ((double *)MALLOC((NLFE7)*sizeof(double))); 01432 create = ((double *)MALLOC((NLFE7)*sizeof(double))); 01433 /* now the 2-d arrays */ 01434 fe.Fe7_wl = ((double **)MALLOC((NLFE7)*sizeof(double *))); 01435 fe.Fe7_emiss = ((double **)MALLOC((NLFE7)*sizeof(double *))); 01436 AulNet = ((double **)MALLOC((NLFE7)*sizeof(double *))); 01437 col_str = ((double **)MALLOC((NLFE7)*sizeof(double *))); 01438 AulPump = ((double **)MALLOC((NLFE7)*sizeof(double *))); 01439 CollRate = ((double **)MALLOC((NLFE7)*sizeof(double *))); 01440 AulDest = ((double **)MALLOC((NLFE7)*sizeof(double *))); 01441 for( i=0; i<(NLFE7); ++i ) 01442 { 01443 fe.Fe7_wl[i] = ((double *)MALLOC((NLFE7)*sizeof(double ))); 01444 fe.Fe7_emiss[i] = ((double *)MALLOC((NLFE7)*sizeof(double ))); 01445 AulNet[i] = ((double *)MALLOC((NLFE7)*sizeof(double ))); 01446 col_str[i] = ((double *)MALLOC((NLFE7)*sizeof(double ))); 01447 AulPump[i] = ((double *)MALLOC((NLFE7)*sizeof(double ))); 01448 CollRate[i] = ((double *)MALLOC((NLFE7)*sizeof(double ))); 01449 AulDest[i] = ((double *)MALLOC((NLFE7)*sizeof(double ))); 01450 } 01451 01452 /* set some to constant values after zeroing out */ 01453 for( i=0; i<NLFE7; ++i ) 01454 { 01455 create[i] = 0.; 01456 destroy[i] = 0.; 01457 for( j=0; j<NLFE7; ++j ) 01458 { 01459 AulNet[i][j] = 0.; 01460 col_str[i][j] = 0.; 01461 CollRate[i][j] = 0.; 01462 AulDest[i][j] = 0.; 01463 AulPump[i][j] = 0.; 01464 fe.Fe7_wl[i][j] = 0.; 01465 fe.Fe7_emiss[i][j] = 0.; 01466 } 01467 } 01468 set_NaN( fe.Fe7_wl[2][1] ); 01469 set_NaN( fe.Fe7_emiss[2][1] ); 01470 set_NaN( fe.Fe7_wl[1][0] ); 01471 set_NaN( fe.Fe7_emiss[1][0] ); 01472 /* two levels within ground term must never be addressed in this array - 01473 * they are fully transferred */ 01474 for( ilo=0; ilo<NLFE7-1; ++ilo ) 01475 { 01476 /* must not do 1-0 or 2-1, which are transferred lines */ 01477 for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi ) 01478 { 01479 fe.Fe7_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / RefIndex( excit_wn[ihi]-excit_wn[ilo] ); 01480 } 01481 } 01482 01483 /* assume FeVII is optically thin - just use As as net escape */ 01484 /*>>refer Fe7 As Berrington, K.A., Nakazaki, S., & Norrington, P.H. 2000, A&AS, 142, 313 */ 01485 AulNet[1][0] = 0.0325; 01486 AulNet[2][0] = 0.167e-8; 01487 /* following corrected from 3.25 to 0.372 as per Keith Berrington email */ 01488 AulNet[3][0] = 0.372; 01489 AulNet[4][0] = 0.135; 01490 AulNet[5][0] = 0.502e-1; 01491 /* following corrected from 0.174 to 0.0150 as per Keith Berrington email */ 01492 AulNet[6][0] = 0.0150; 01493 AulNet[7][0] = 0.959e-3; 01494 01495 AulNet[2][1] = 0.466e-1; 01496 AulNet[3][1] = 0.603; 01497 AulNet[5][1] = 0.762e-1; 01498 AulNet[6][1] = 0.697e-1; 01499 AulNet[7][1] = 0.343; 01500 01501 AulNet[3][2] = 0.139e-2; 01502 AulNet[6][2] = 0.735e-1; 01503 AulNet[7][2] = 0.503; 01504 01505 AulNet[4][3] = 0.472e-6; 01506 AulNet[5][3] = 0.572e-1; 01507 /* following corrected from 0.191 to 0.182 as per Keith Berrington email */ 01508 AulNet[6][3] = 0.182; 01509 AulNet[7][3] = 0.414e-2; 01510 01511 AulNet[5][4] = 0.115e-2; 01512 AulNet[6][4] = 0.139e-7; 01513 01514 AulNet[6][5] = 0.743e-2; 01515 01516 AulNet[7][6] = 0.454e-4; 01517 01518 /* collision strengths at log T 4.5 */ 01520 /*>>refer Fe7 cs Berrington, K.A., Nakazaki, S., & Norrington, P.H. 2000, A&AS, 142, 313 */ 01521 # if 0 01522 if( fudge(-1) ) 01523 { 01524 fprintf(ioQQQ,"DEBUG fudge call cool_iron\n"); 01525 scale = fudge(0); 01526 } 01527 else 01528 # endif 01529 scale = 1.; 01530 01531 col_str[1][0] = 3.35; 01532 col_str[2][0] = 1.17; 01533 col_str[3][0] = 0.959; 01534 col_str[4][0] = 0.299; 01535 col_str[5][0] = 0.633; 01536 col_str[6][0] = 0.549; 01537 col_str[7][0] = 1.24*scale; 01538 01539 col_str[2][1] = 4.11; 01540 col_str[3][1] = 1.29; 01541 col_str[4][1] = 0.235; 01542 col_str[5][1] = 0.833; 01543 col_str[6][1] = 1.06; 01544 col_str[7][1] = 1.74*scale; 01545 01546 col_str[3][2] = 1.60; 01547 col_str[4][2] = 0.187; 01548 col_str[5][2] = 0.690; 01549 col_str[6][2] = 1.94; 01550 col_str[7][2] = 2.25*scale; 01551 01552 col_str[4][3] = 0.172; 01553 col_str[5][3] = 0.531; 01554 col_str[6][3] = 1.06; 01555 col_str[7][3] = 2.02; 01556 01557 col_str[5][4] = 0.370; 01558 col_str[6][4] = 0.324; 01559 col_str[7][4] = 0.164; 01560 01561 col_str[6][5] = 1.17; 01562 col_str[7][5] = 0.495; 01563 01564 col_str[7][6] = 0.903; 01565 01566 /* check whether level 2 lines are on, and if so, find the driving lines that 01567 * can pump the upper levels of this atom */ 01568 if( nWindLine > 0 ) 01569 { 01570 ipPump = -1; 01571 for( i=0; i<nWindLine; ++i ) 01572 { 01573 /* don't test on nelem==ipIRON since lines on physics, not C, scale */ 01574 if( TauLine2[i].Hi->nelem ==26 && TauLine2[i].Hi->IonStg==7 && 01575 /* >>chng 05 jul 10, move line to wavelength that agrees with nist tables 01576 * here and in level2 data file 01577 * NB must add a few wn to second number to get hit - 01578 * logic is that this is lowest E1 transition */ 01579 (TauLine2[i].EnergyWN-4.27360E+05) < 0. ) 01580 { 01581 ipPump = i; 01582 break; 01583 } 01584 } 01585 if( ipPump<0 ) 01586 { 01587 fprintf(ioQQQ,"PROBLEM Fe7Lev8 cannot identify the FUV driving line.\n"); 01588 TotalInsanity(); 01589 } 01590 } 01591 else 01592 { 01593 ipPump = 0; 01594 } 01595 } 01596 01597 /* bail if no ions */ 01598 if( dense.xIonDense[ipIRON][6] <= 0. ) 01599 { 01600 CoolAdd("Fe 7",0,0.); 01601 01602 fe.Fe7CoolTot = 0.; 01603 for( ilo=0; ilo<NLFE7-1; ++ilo ) 01604 { 01605 /* must not do 1-0 or 2-1, which are transferred lines */ 01606 for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi ) 01607 { 01608 fe.Fe7_emiss[ihi][ilo] = 0.; 01609 } 01610 } 01611 TauLines[ipFe0795].Lo->Pop = 0.; 01612 TauLines[ipFe0795].Emis->PopOpc = 0.; 01613 TauLines[ipFe0795].Hi->Pop = 0.; 01614 TauLines[ipFe0795].Emis->xIntensity = 0.; 01615 TauLines[ipFe0795].Coll.cool = 0.; 01616 TauLines[ipFe0795].Emis->phots = 0.; 01617 TauLines[ipFe0795].Emis->ColOvTot = 0.; 01618 TauLines[ipFe0795].Coll.heat = 0.; 01619 CoolAdd( "Fe 7", TauLines[ipFe0795].WLAng , 0.); 01620 TauLines[ipFe0778].Lo->Pop = 0.; 01621 TauLines[ipFe0778].Emis->PopOpc = 0.; 01622 TauLines[ipFe0778].Hi->Pop = 0.; 01623 TauLines[ipFe0778].Emis->xIntensity = 0.; 01624 TauLines[ipFe0778].Coll.cool = 0.; 01625 TauLines[ipFe0778].Emis->phots = 0.; 01626 TauLines[ipFe0778].Emis->ColOvTot = 0.; 01627 TauLines[ipFe0778].Coll.heat = 0.; 01628 CoolAdd( "Fe 7", TauLines[ipFe0778].WLAng , 0.); 01629 /* level populations */ 01630 /* LIMLEVELN is the dimension of the atoms vectors */ 01631 ASSERT( NLFE7 <= LIMLEVELN); 01632 for( i=0; i < NLFE7; i++ ) 01633 { 01634 atoms.PopLevels[i] = 0.; 01635 atoms.DepLTELevels[i] = 1.; 01636 } 01637 return; 01638 } 01639 01640 /* do pump rate for FUV excitation of 3P (levels 5-7 on physics scale, not C scale) */ 01641 if( ipPump ) 01642 { 01643 FUV_pump = TauLine2[ipPump].Emis->pump * 0.3 /(0.3+TauLine2[ipPump].Emis->Pesc); 01644 } 01645 else 01646 { 01647 FUV_pump = 0.; 01648 } 01649 01650 /* this represents photoexcitation of 3P from ground level 01651 * >>chng 04 nov 22, upper array elements were on physics not c scale, off by one, TE */ 01652 AulPump[0][4] = FUV_pump; 01653 AulPump[1][4] = FUV_pump; 01654 AulPump[2][4] = FUV_pump; 01655 AulPump[0][5] = FUV_pump; 01656 AulPump[1][5] = FUV_pump; 01657 AulPump[2][5] = FUV_pump; 01658 AulPump[0][6] = FUV_pump; 01659 AulPump[1][6] = FUV_pump; 01660 AulPump[2][6] = FUV_pump; 01661 /* these were set in the previous call to atom_levelN as the previous pump times 01662 * ratio of statistical weights, so this is the downward pump rate */ 01663 AulPump[4][0] = 0; 01664 AulPump[4][1] = 0; 01665 AulPump[4][2] = 0; 01666 AulPump[5][0] = 0; 01667 AulPump[5][1] = 0; 01668 AulPump[5][2] = 0; 01669 AulPump[6][0] = 0; 01670 AulPump[6][1] = 0; 01671 AulPump[6][2] = 0; 01672 01673 /* within ground term update to rt results */ 01674 AulNet[1][0] = TauLines[ipFe0795].Emis->Aul*(TauLines[ipFe0795].Emis->Pesc + TauLines[ipFe0795].Emis->Pelec_esc); 01675 AulDest[1][0] = TauLines[ipFe0795].Emis->Aul*TauLines[ipFe0795].Emis->Pdest; 01676 AulPump[0][1] = TauLines[ipFe0795].Emis->pump; 01677 AulPump[1][0] = 0.; 01678 01679 AulNet[2][1] = TauLines[ipFe0778].Emis->Aul*(TauLines[ipFe0778].Emis->Pesc + TauLines[ipFe0778].Emis->Pelec_esc); 01680 AulDest[2][1] = TauLines[ipFe0778].Emis->Aul*TauLines[ipFe0778].Emis->Pdest; 01681 AulPump[1][2] = TauLines[ipFe0778].Emis->pump; 01682 AulPump[2][1] = 0.; 01683 01684 /* nNegPop positive if negative pops occurred, negative if too cold */ 01685 atom_levelN( 01686 /* number of levels */ 01687 NLFE7, 01688 /* the abundance of the ion, cm-3 */ 01689 dense.xIonDense[ipIRON][6], 01690 /* the statistical weights */ 01691 gfe7, 01692 /* the excitation energies in wavenumbers */ 01693 excit_wn, 01694 /* units are wavenumbers */ 01695 'w', 01696 /* the derived populations - cm-3 */ 01697 pops, 01698 /* the derived departure coefficients */ 01699 depart, 01700 /* the net emission rate, Aul * escp prob */ 01701 &AulNet , 01702 /* the collision strengths */ 01703 &col_str , 01704 /* A * destruction prob */ 01705 &AulDest, 01706 /* pumping rate */ 01707 &AulPump, 01708 /* collision rate, s-1, must defined if no collision strengths */ 01709 &CollRate, 01710 /* creation vector */ 01711 create, 01712 /* destruction vector */ 01713 destroy, 01714 /* say atom_levelN should evaluate coll rates from cs */ 01715 false, 01716 &fe.Fe7CoolTot, 01717 &dfe7dt, 01718 "Fe 7", 01719 &nNegPop, 01720 &lgZeroPop, 01721 false ); 01722 01723 /* LIMLEVELN is the dim of the PopLevels vector */ 01724 ASSERT( NLFE7 <= LIMLEVELN ); 01725 for( i=0; i<NLFE7; ++i) 01726 { 01727 atoms.PopLevels[i] = pops[i]; 01728 atoms.DepLTELevels[i] = depart[i]; 01729 } 01730 01731 if( lgZeroPop ) 01732 { 01733 /* this case, too cool to excite upper levels of atom, but will want 01734 * to do IR lines in ground term */ 01735 PutCS(col_str[1][0],&TauLines[ipFe0795]); 01736 PutCS(col_str[2][1],&TauLines[ipFe0778]); 01737 PutCS(col_str[2][0],&TauDummy); 01738 atom_level3(&TauLines[ipFe0795],&TauLines[ipFe0778],&TauDummy); 01739 atoms.PopLevels[0] = TauLines[ipFe0795].Lo->Pop; 01740 atoms.PopLevels[1] = TauLines[ipFe0795].Hi->Pop; 01741 atoms.PopLevels[2] = TauLines[ipFe0778].Hi->Pop; 01742 for( ilo=0; ilo<NLFE7-1; ++ilo ) 01743 { 01744 /* must not do 1-0 or 2-1, which are transferred lines */ 01745 for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi ) 01746 { 01747 fe.Fe7_emiss[ihi][ilo] = 0.; 01748 } 01749 } 01750 } 01751 else 01752 { 01753 /* this case non-zero pops, but must set up vars within transition structure */ 01754 TauLines[ipFe0795].Lo->Pop = atoms.PopLevels[0]; 01755 TauLines[ipFe0795].Hi->Pop = atoms.PopLevels[1]; 01756 TauLines[ipFe0795].Emis->PopOpc = (pops[0] - pops[1]*gfe7[0]/gfe7[1]); 01757 TauLines[ipFe0795].Emis->xIntensity = TauLines[ipFe0795].Emis->phots*TauLines[ipFe0795].EnergyErg; 01758 TauLines[ipFe0795].Emis->phots = TauLines[ipFe0795].Emis->Aul*(TauLines[ipFe0795].Emis->Pesc + TauLines[ipFe0795].Emis->Pelec_esc)*pops[1]; 01759 TauLines[ipFe0795].Emis->ColOvTot = CollRate[0][1]/(CollRate[0][1]+TauLines[ipFe0795].Emis->pump); 01760 TauLines[ipFe0795].Coll.cool = 0.; 01761 TauLines[ipFe0795].Coll.heat = 0.; 01762 01763 TauLines[ipFe0778].Lo->Pop = atoms.PopLevels[1]; 01764 TauLines[ipFe0778].Hi->Pop = atoms.PopLevels[2]; 01765 TauLines[ipFe0778].Emis->PopOpc = (pops[1] - pops[2]*gfe7[1]/gfe7[2]); 01766 TauLines[ipFe0778].Emis->xIntensity = TauLines[ipFe0778].Emis->phots*TauLines[ipFe0778].EnergyErg; 01767 TauLines[ipFe0778].Emis->phots = TauLines[ipFe0778].Emis->Aul*(TauLines[ipFe0778].Emis->Pesc + TauLines[ipFe0778].Emis->Pelec_esc)*pops[2]; 01768 TauLines[ipFe0778].Emis->ColOvTot = CollRate[1][2]/(CollRate[1][2]+TauLines[ipFe0778].Emis->pump); 01769 TauLines[ipFe0778].Coll.heat = 0.; 01770 TauLines[ipFe0778].Coll.cool = 0.; 01771 } 01772 01773 if( nNegPop>0 ) 01774 { 01775 fprintf( ioQQQ, "PROBLEM Fe7Lev8 found negative populations\n" ); 01776 ShowMe(); 01777 cdEXIT(EXIT_FAILURE); 01778 } 01779 01780 /* add cooling then its derivative */ 01781 CoolAdd("Fe 7",0,fe.Fe7CoolTot); 01782 /* derivative of cooling */ 01783 thermal.dCooldT += dfe7dt; 01784 01785 /* find emission in each line */ 01786 for( ilo=0; ilo<NLFE7-1; ++ilo ) 01787 { 01788 /* must not do 1-0 or 2-1, which are transferred lines */ 01789 for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi ) 01790 { 01791 /* emission in these lines */ 01792 fe.Fe7_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*ERG1CM; 01793 } 01794 } 01795 return; 01796 } 01797 01798 /*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion 01799 * >>chng 05 dec 17, code provided by Kevin Blagrave and described in 01800 *>>refer Fe3 model Blagrave, K.P.M., Martin, P.G. & Baldwin, J.A. 01801 *>>refercon 2006, ApJ, in press (astro-ph 0603167) */ 01802 STATIC void Fe3Lev14(void) 01803 { 01804 bool lgZeroPop; 01805 int nNegPop; 01806 long int i, 01807 j; 01808 static bool lgFirst=true; 01809 01810 double dfe3dt; 01811 01812 long int ihi , ilo; 01813 static double 01814 *depart, 01815 *pops, 01816 *destroy, 01817 *create , 01818 **AulDest, 01819 **CollRate, 01820 **AulPump, 01821 **AulNet, 01822 **col_str; 01823 01824 /* statistical weights for the n levels */ 01825 static double gfe3[NLFE3]={9.,7.,5.,3.,1.,5.,13.,11.,9.,3.,1.,9.,7.,5.}; 01826 01827 /*refer Fe3 energies NIST version 3 Atomic Spectra Database */ 01828 /* from smallest to largest energy (not in multiplet groupings) */ 01829 01830 /* energy in wavenumbers */ 01831 static double excit_wn[NLFE3]={ 01832 0.0 , 436.2, 738.9, 932.4, 1027.3, 01833 19404.8, 20051.1, 20300.8, 20481.9, 20688.4, 01834 21208.5, 21462.2, 21699.9, 21857.2 }; 01835 01836 DEBUG_ENTRY( "Fe3Lev14()" ); 01837 01838 if( lgFirst ) 01839 { 01840 /* will never do this again */ 01841 lgFirst = false; 01842 /* allocate the 1D arrays*/ 01843 /* create space for the 2D arrays */ 01844 depart = ((double *)MALLOC((NLFE3)*sizeof(double))); 01845 pops = ((double *)MALLOC((NLFE3)*sizeof(double))); 01846 destroy = ((double *)MALLOC((NLFE3)*sizeof(double))); 01847 create = ((double *)MALLOC((NLFE3)*sizeof(double))); 01848 /* now the 2-d arrays */ 01849 fe.Fe3_wl = ((double **)MALLOC((NLFE3)*sizeof(double *))); 01850 fe.Fe3_emiss = ((double **)MALLOC((NLFE3)*sizeof(double *))); 01851 AulNet = ((double **)MALLOC((NLFE3)*sizeof(double *))); 01852 col_str = ((double **)MALLOC((NLFE3)*sizeof(double *))); 01853 AulPump = ((double **)MALLOC((NLFE3)*sizeof(double *))); 01854 CollRate = ((double **)MALLOC((NLFE3)*sizeof(double *))); 01855 AulDest = ((double **)MALLOC((NLFE3)*sizeof(double *))); 01856 for( i=0; i<(NLFE3); ++i ) 01857 { 01858 fe.Fe3_wl[i] = ((double *)MALLOC((NLFE3)*sizeof(double ))); 01859 fe.Fe3_emiss[i] = ((double *)MALLOC((NLFE3)*sizeof(double ))); 01860 AulNet[i] = ((double *)MALLOC((NLFE3)*sizeof(double ))); 01861 col_str[i] = ((double *)MALLOC((NLFE3)*sizeof(double ))); 01862 AulPump[i] = ((double *)MALLOC((NLFE3)*sizeof(double ))); 01863 CollRate[i] = ((double *)MALLOC((NLFE3)*sizeof(double ))); 01864 AulDest[i] = ((double *)MALLOC((NLFE3)*sizeof(double ))); 01865 } 01866 01867 /* set some to constant values after zeroing out */ 01868 for( i=0; i<NLFE3; ++i ) 01869 { 01870 create[i] = 0.; 01871 destroy[i] = 0.; 01872 for( j=0; j<NLFE3; ++j ) 01873 { 01874 AulNet[i][j] = 0.; 01875 col_str[i][j] = 0.; 01876 CollRate[i][j] = 0.; 01877 AulDest[i][j] = 0.; 01878 AulPump[i][j] = 0.; 01879 fe.Fe3_wl[i][j] = 0.; 01880 fe.Fe3_emiss[i][j] = 0.; 01881 } 01882 } 01883 /* calculates wavelengths of transitions */ 01884 /* dividing by RefIndex converts the vacuum wavelength to air wavelength */ 01885 for( ihi=1; ihi<NLFE3; ++ihi ) 01886 { 01887 for( ilo=0; ilo<ihi; ++ilo ) 01888 { 01889 fe.Fe3_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / 01890 RefIndex( (excit_wn[ihi]-excit_wn[ilo]) ); 01891 } 01892 } 01893 01894 /* assume FeIII is optically thin - just use As as net escape */ 01895 /*>>refer Fe3 as Quinet, P., 1996, A&AS, 116, 573 */ 01896 AulNet[1][0] = 2.8e-3; 01897 AulNet[7][0] = 4.9e-6; 01898 AulNet[8][0] = 5.7e-3; 01899 AulNet[11][0] = 4.5e-1; 01900 AulNet[12][0] = 4.2e-2; 01901 01902 AulNet[2][1] = 1.8e-3; 01903 AulNet[5][1] = 4.2e-1; 01904 AulNet[8][1] = 1.0e-3; 01905 AulNet[11][1] = 8.4e-2; 01906 AulNet[12][1] = 2.5e-1; 01907 AulNet[13][1] = 2.7e-2; 01908 01909 AulNet[3][2] = 7.0e-4; 01910 AulNet[5][2] = 5.1e-5; 01911 AulNet[9][2] = 5.4e-1; 01912 AulNet[12][2] = 8.5e-2; 01913 AulNet[13][2] = 9.8e-2; 01914 01915 AulNet[4][3] = 1.4e-4; 01916 AulNet[5][3] = 3.9e-2; 01917 AulNet[9][3] = 4.1e-5; 01918 AulNet[10][3] = 7.0e-1; 01919 AulNet[13][3] = 4.7e-2; 01920 01921 AulNet[9][4] = 9.3e-2; 01922 01923 AulNet[9][5] = 4.7e-2; 01924 AulNet[12][5] = 2.5e-6; 01925 AulNet[13][5] = 1.7e-5; 01926 01927 AulNet[7][6] = 2.7e-4; 01928 01929 AulNet[8][7] = 1.2e-4; 01930 AulNet[11][7] = 6.6e-4; 01931 01932 AulNet[11][8] = 1.6e-3; 01933 AulNet[12][8] = 7.8e-4; 01934 01935 AulNet[10][9] = 8.4e-3; 01936 AulNet[13][9] = 2.8e-7; 01937 01938 AulNet[12][11] = 3.0e-4; 01939 01940 AulNet[13][12] = 1.4e-4; 01941 01942 /* collision strengths at log T 4 */ 01944 /*>>refer Fe3 cs Zhang, H. 1996, 119, 523 */ 01945 col_str[1][0] = 2.92; 01946 col_str[2][0] = 1.24; 01947 col_str[3][0] = 0.595; 01948 col_str[4][0] = 0.180; 01949 col_str[5][0] = 0.580; 01950 col_str[6][0] = 1.34; 01951 col_str[7][0] = 0.489; 01952 col_str[8][0] = 0.0926; 01953 col_str[9][0] = 0.165; 01954 col_str[10][0] = 0.0213; 01955 col_str[11][0] = 1.07; 01956 col_str[12][0] = 0.435; 01957 col_str[13][0] = 0.157; 01958 01959 col_str[2][1] = 2.06; 01960 col_str[3][1] = 0.799; 01961 col_str[4][1] = 0.225; 01962 col_str[5][1] = 0.335; 01963 col_str[6][1] = 0.555; 01964 col_str[7][1] = 0.609; 01965 col_str[8][1] = 0.367; 01966 col_str[9][1] = 0.195; 01967 col_str[10][1] = 0.0698; 01968 col_str[11][1] = 0.538; 01969 col_str[12][1] = 0.484; 01970 col_str[13][1] = 0.285; 01971 01972 col_str[3][2] = 1.29; 01973 col_str[4][2] = 0.312; 01974 col_str[5][2] = 0.173; 01975 col_str[6][2] = 0.178; 01976 col_str[7][2] = 0.430; 01977 col_str[8][2] = 0.486; 01978 col_str[9][2] = 0.179; 01979 col_str[10][2] = 0.0741; 01980 col_str[11][2] = 0.249; 01981 col_str[12][2] = 0.362; 01982 col_str[13][2] = 0.324; 01983 01984 col_str[4][3] = 0.493; 01985 col_str[5][3] = 0.0767; 01986 col_str[6][3] = 0.0348; 01987 col_str[7][3] = 0.223; 01988 col_str[8][3] = 0.401; 01989 col_str[9][3] = 0.126; 01990 col_str[10][3] = 0.0528; 01991 col_str[11][3] = 0.101; 01992 col_str[12][3] = 0.207; 01993 col_str[13][3] = 0.253; 01994 01995 col_str[5][4] = 0.0211; 01996 col_str[6][4] = 0.00122; 01997 col_str[7][4] = 0.0653; 01998 col_str[8][4] = 0.154; 01999 col_str[9][4] = 0.0453; 02000 col_str[10][4] = 0.0189; 02001 col_str[11][4] = 0.0265; 02002 col_str[12][4] = 0.0654; 02003 col_str[13][4] = 0.0950; 02004 02005 col_str[6][5] = 0.403; 02006 col_str[7][5] = 0.213; 02007 col_str[8][5] = 0.0939; 02008 col_str[9][5] = 1.10; 02009 col_str[10][5] = 0.282; 02010 col_str[11][5] = 0.942; 02011 col_str[12][5] = 0.768; 02012 col_str[13][5] = 0.579; 02013 02014 col_str[7][6] = 2.84; /* 10-9 */ 02015 col_str[8][6] = 0.379; /* 11-9 */ 02016 col_str[9][6] = 0.0876; /* 7-9 */ 02017 col_str[10][6] = 0.00807; /* 8-9 */ 02018 col_str[11][6] = 1.85; /* 12-9 */ 02019 col_str[12][6] = 0.667; /* 13-9 */ 02020 col_str[13][6] = 0.0905; /* 14-9 */ 02021 02022 col_str[8][7] = 3.07; /* 11-10 */ 02023 col_str[9][7] = 0.167; /* 7-10 */ 02024 col_str[10][7] = 0.0526; /* 8-10 */ 02025 col_str[11][7] = 0.814; /* 12-10 */ 02026 col_str[12][7] = 0.837; /* 13-10 */ 02027 col_str[13][7] = 0.626; /* 14-10 */ 02028 02029 col_str[9][8] = 0.181; /* 7-11 */ 02030 col_str[10][8] = 0.0854; /* 8-11 */ 02031 col_str[11][8] = 0.180; /* 12-11 */ 02032 col_str[12][8] = 0.778; /* 13-11 */ 02033 col_str[13][8] = 0.941; /* 14-11 */ 02034 02035 col_str[10][9] = 0.377; /* 8-7 */ 02036 col_str[11][9] = 0.603; /* 12-7 */ 02037 col_str[12][9] = 0.472; /* 13-7 */ 02038 col_str[13][9] = 0.302; /* 14-7 */ 02039 02040 col_str[11][10] = 0.216; /* 12-8 */ 02041 col_str[12][10] = 0.137; /* 13-8 */ 02042 col_str[13][10] = 0.106; /* 14-8 */ 02043 02044 col_str[12][11] = 1.25; 02045 col_str[13][11] = 0.292; 02046 02047 col_str[13][12] = 1.10; 02048 } 02049 02050 /* bail if no ions */ 02051 if( dense.xIonDense[ipIRON][2] <= 0. ) 02052 { 02053 CoolAdd("Fe 3",0,0.); 02054 02055 fe.Fe3CoolTot = 0.; 02056 for( ihi=1; ihi<NLFE3; ++ihi ) 02057 { 02058 for( ilo=0; ilo<ihi; ++ilo ) 02059 { 02060 fe.Fe3_emiss[ihi][ilo] = 0.; 02061 } 02062 } 02063 /* level populations */ 02064 /* LIMLEVELN is the dimension of the atoms vectors */ 02065 ASSERT( NLFE3 <= LIMLEVELN); 02066 for( i=0; i < NLFE3; i++ ) 02067 { 02068 atoms.PopLevels[i] = 0.; 02069 atoms.DepLTELevels[i] = 1.; 02070 } 02071 return; 02072 } 02073 02074 /* nNegPop positive if negative pops occurred, negative if too cold */ 02075 atom_levelN( 02076 /* number of levels */ 02077 NLFE3, 02078 /* the abundance of the ion, cm-3 */ 02079 dense.xIonDense[ipIRON][2], 02080 /* the statistical weights */ 02081 gfe3, 02082 /* the excitation energies */ 02083 excit_wn, 02084 'w', 02085 /* the derived populations - cm-3 */ 02086 pops, 02087 /* the derived departure coefficients */ 02088 depart, 02089 /* the net emission rate, Aul * escape prob */ 02090 &AulNet , 02091 /* the collision strengths */ 02092 &col_str , 02093 /* A * destruction prob */ 02094 &AulDest, 02095 /* pumping rate */ 02096 &AulPump, 02097 /* collision rate, s-1, must defined if no collision strengths */ 02098 &CollRate, 02099 /* creation vector */ 02100 create, 02101 /* destruction vector */ 02102 destroy, 02103 /* say atom_levelN should evaluate coll rates from cs */ 02104 false, 02105 &fe.Fe3CoolTot, 02106 &dfe3dt, 02107 "Fe 3", 02108 &nNegPop, 02109 &lgZeroPop, 02110 false ); 02111 02112 /* LIMLEVELN is the dim of the PopLevels vector */ 02113 ASSERT( NLFE3 <= LIMLEVELN ); 02114 for( i=0; i<NLFE3; ++i) 02115 { 02116 atoms.PopLevels[i] = pops[i]; 02117 atoms.DepLTELevels[i] = depart[i]; 02118 } 02119 02120 if( nNegPop>0 ) 02121 { 02122 fprintf( ioQQQ, " Fe3Lev14 found negative populations\n" ); 02123 ShowMe(); 02124 cdEXIT(EXIT_FAILURE); 02125 } 02126 02127 /* add cooling then its derivative */ 02128 CoolAdd("Fe 3",0,fe.Fe3CoolTot); 02129 /* derivative of cooling */ 02130 thermal.dCooldT += dfe3dt; 02131 02132 /* find emission in each line */ 02133 for( ihi=1; ihi<NLFE3; ++ihi ) 02134 { 02135 for( ilo=0; ilo<ihi; ++ilo ) 02136 { 02137 /* emission in these lines */ 02138 fe.Fe3_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN; 02139 } 02140 } 02141 return; 02142 } 02143 02144 /*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */ 02145 STATIC void Fe11Lev5(void) 02146 { 02147 bool lgZeroPop; 02148 int nNegPop; 02149 long int i, 02150 j; 02151 static bool lgFirst=true; 02152 02153 double dCool_dT; 02154 02155 long int ihi , ilo; 02156 static double 02157 *depart, 02158 *pops, 02159 *destroy, 02160 *create , 02161 **AulDest, 02162 **CollRate, 02163 **AulPump, 02164 **AulNet, 02165 **col_str , 02166 TeUsed=-1.; 02167 02168 /* statistical weights for the n levels */ 02169 static double stat_wght[NLFE11]={5.,3.,1.,5.,1.}; 02170 02171 /*refer Fe11 energies NIST version 3 Atomic Spectra Database */ 02172 /* from smallest to largest energy (not in multiplet groupings) */ 02173 02174 /* energy in wavenumbers */ 02175 static double excit_wn[NLFE11]={ 02176 0.0 , 12667.9 , 14312. , 37743.6 , 80814.7 }; 02177 02178 DEBUG_ENTRY( "Fe11Lev5()" ); 02179 02180 if( lgFirst ) 02181 { 02182 /* will never do this again */ 02183 lgFirst = false; 02184 /* allocate the 1D arrays*/ 02185 /* create space for the 2D arrays */ 02186 depart = ((double *)MALLOC((NLFE11)*sizeof(double))); 02187 pops = ((double *)MALLOC((NLFE11)*sizeof(double))); 02188 destroy = ((double *)MALLOC((NLFE11)*sizeof(double))); 02189 create = ((double *)MALLOC((NLFE11)*sizeof(double))); 02190 /* now the 2-d arrays */ 02191 fe.Fe11_wl = ((double **)MALLOC((NLFE11)*sizeof(double *))); 02192 fe.Fe11_emiss = ((double **)MALLOC((NLFE11)*sizeof(double *))); 02193 AulNet = ((double **)MALLOC((NLFE11)*sizeof(double *))); 02194 col_str = ((double **)MALLOC((NLFE11)*sizeof(double *))); 02195 AulPump = ((double **)MALLOC((NLFE11)*sizeof(double *))); 02196 CollRate = ((double **)MALLOC((NLFE11)*sizeof(double *))); 02197 AulDest = ((double **)MALLOC((NLFE11)*sizeof(double *))); 02198 for( i=0; i<(NLFE11); ++i ) 02199 { 02200 fe.Fe11_wl[i] = ((double *)MALLOC((NLFE11)*sizeof(double ))); 02201 fe.Fe11_emiss[i] = ((double *)MALLOC((NLFE11)*sizeof(double ))); 02202 AulNet[i] = ((double *)MALLOC((NLFE11)*sizeof(double ))); 02203 col_str[i] = ((double *)MALLOC((NLFE11)*sizeof(double ))); 02204 AulPump[i] = ((double *)MALLOC((NLFE11)*sizeof(double ))); 02205 CollRate[i] = ((double *)MALLOC((NLFE11)*sizeof(double ))); 02206 AulDest[i] = ((double *)MALLOC((NLFE11)*sizeof(double ))); 02207 } 02208 02209 /* set some to constant values after zeroing out */ 02210 for( i=0; i<NLFE11; ++i ) 02211 { 02212 create[i] = 0.; 02213 destroy[i] = 0.; 02214 for( j=0; j<NLFE11; ++j ) 02215 { 02216 AulNet[i][j] = 0.; 02217 col_str[i][j] = 0.; 02218 CollRate[i][j] = 0.; 02219 AulDest[i][j] = 0.; 02220 AulPump[i][j] = 0.; 02221 fe.Fe11_wl[i][j] = 0.; 02222 fe.Fe11_emiss[i][j] = 0.; 02223 } 02224 } 02225 /* calculates wavelengths of transitions */ 02226 /* dividing by RefIndex converts the vacuum wavelength to air wavelength */ 02227 for( ihi=1; ihi<NLFE11; ++ihi ) 02228 { 02229 for( ilo=0; ilo<ihi; ++ilo ) 02230 { 02231 fe.Fe11_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / 02232 RefIndex( (excit_wn[ihi]-excit_wn[ilo]) ); 02233 } 02234 } 02235 02236 /*>>refer Fe11 as Mendoza C. & Zeippen, C. J. 1983, MNRAS, 202, 981 */ 02237 AulNet[4][0] = 1.7; 02238 AulNet[4][1] = 990.; 02239 AulNet[4][3] = 8.3; 02240 AulNet[3][0] = 92.3; 02241 AulNet[3][1] = 9.44; 02242 AulNet[3][2] = 1.4e-3; 02243 AulNet[2][0] = 9.9e-3; 02244 AulNet[1][0] = 43.7; 02245 AulNet[2][1] = 0.226; 02246 02247 } 02248 02249 /* bail if no ions */ 02250 if( dense.xIonDense[ipIRON][10] <= 0. ) 02251 { 02252 CoolAdd("Fe11",0,0.); 02253 02254 fe.Fe11CoolTot = 0.; 02255 for( ihi=1; ihi<NLFE11; ++ihi ) 02256 { 02257 for( ilo=0; ilo<ihi; ++ilo ) 02258 { 02259 fe.Fe11_emiss[ihi][ilo] = 0.; 02260 } 02261 } 02262 /* level populations */ 02263 /* LIMLEVELN is the dimension of the atoms vectors */ 02264 ASSERT( NLFE11 <= LIMLEVELN); 02265 for( i=0; i < NLFE11; i++ ) 02266 { 02267 atoms.PopLevels[i] = 0.; 02268 atoms.DepLTELevels[i] = 1.; 02269 } 02270 return; 02271 } 02272 02273 /* evaluate collision strengths if Te changed by too much */ 02274 if( fabs(phycon.te / TeUsed - 1. ) > 0.05 ) 02275 { 02276 TeUsed = phycon.te; 02277 02278 /* collision strengths at current temperature */ 02279 for( ihi=1; ihi<NLFE11; ++ihi ) 02280 { 02281 for( ilo=0; ilo<ihi; ++ilo ) 02282 { 02283 col_str[ihi][ilo] = Fe_10_11_13_cs( 11 , ilo+1 , ihi+1 ); 02284 ASSERT( col_str[ihi][ilo] > 0. ); 02285 } 02286 } 02287 } 02288 02289 /* nNegPop positive if negative pops occurred, negative if too cold */ 02290 atom_levelN( 02291 /* number of levels */ 02292 NLFE11, 02293 /* the abundance of the ion, cm-3 */ 02294 dense.xIonDense[ipIRON][10], 02295 /* the statistical weights */ 02296 stat_wght, 02297 /* the excitation energies */ 02298 excit_wn, 02299 'w', 02300 /* the derived populations - cm-3 */ 02301 pops, 02302 /* the derived departure coefficients */ 02303 depart, 02304 /* the net emission rate, Aul * escp prob */ 02305 &AulNet , 02306 /* the collision strengths */ 02307 &col_str , 02308 /* A * destruction prob */ 02309 &AulDest, 02310 /* pumping rate */ 02311 &AulPump, 02312 /* collision rate, s-1, must defined if no collision strengths */ 02313 &CollRate, 02314 /* creation vector */ 02315 create, 02316 /* destruction vector */ 02317 destroy, 02318 /* say atom_levelN should evaluate coll rates from cs */ 02319 false, 02320 &fe.Fe11CoolTot, 02321 &dCool_dT, 02322 "Fe11", 02323 &nNegPop, 02324 &lgZeroPop, 02325 false ); 02326 02327 /* LIMLEVELN is the dim of the PopLevels vector */ 02328 ASSERT( NLFE11 <= LIMLEVELN ); 02329 for( i=0; i<NLFE11; ++i) 02330 { 02331 atoms.PopLevels[i] = pops[i]; 02332 atoms.DepLTELevels[i] = depart[i]; 02333 } 02334 02335 if( nNegPop>0 ) 02336 { 02337 fprintf( ioQQQ, " Fe11Lev5 found negative populations\n" ); 02338 ShowMe(); 02339 cdEXIT(EXIT_FAILURE); 02340 } 02341 02342 /* add cooling then its derivative */ 02343 CoolAdd("Fe11",0,fe.Fe11CoolTot); 02344 /* derivative of cooling */ 02345 thermal.dCooldT += dCool_dT; 02346 02347 /* find emission in each line */ 02348 for( ihi=1; ihi<NLFE11; ++ihi ) 02349 { 02350 for( ilo=0; ilo<ihi; ++ilo ) 02351 { 02352 /* emission in these lines */ 02353 fe.Fe11_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]* 02354 (excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN; 02355 } 02356 } 02357 return; 02358 } 02359 02360 /*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */ 02361 STATIC void Fe13Lev5(void) 02362 { 02363 bool lgZeroPop; 02364 int nNegPop; 02365 long int i, 02366 j; 02367 static bool lgFirst=true; 02368 02369 double dCool_dT; 02370 02371 long int ihi , ilo; 02372 static double 02373 *depart, 02374 *pops, 02375 *destroy, 02376 *create , 02377 **AulDest, 02378 **CollRate, 02379 **AulPump, 02380 **AulNet, 02381 **col_str , 02382 TeUsed=-1.; 02383 02384 /* statistical weights for the n levels */ 02385 static double stat_wght[NLFE13]={1.,3.,5.,5.,1.}; 02386 02387 /*refer Fe13 energies NIST version 3 Atomic Spectra Database */ 02388 /* energy in wavenumbers */ 02389 static double excit_wn[NLFE13]={ 02390 0.0 , 9302.5 , 18561.0 , 48068. , 91508. }; 02391 02392 DEBUG_ENTRY( "Fe13Lev5()" ); 02393 02394 if( lgFirst ) 02395 { 02396 /* will never do this again */ 02397 lgFirst = false; 02398 /* allocate the 1D arrays*/ 02399 /* create space for the 2D arrays */ 02400 depart = ((double *)MALLOC((NLFE13)*sizeof(double))); 02401 pops = ((double *)MALLOC((NLFE13)*sizeof(double))); 02402 destroy = ((double *)MALLOC((NLFE13)*sizeof(double))); 02403 create = ((double *)MALLOC((NLFE13)*sizeof(double))); 02404 /* now the 2-d arrays */ 02405 fe.Fe13_wl = ((double **)MALLOC((NLFE13)*sizeof(double *))); 02406 fe.Fe13_emiss = ((double **)MALLOC((NLFE13)*sizeof(double *))); 02407 AulNet = ((double **)MALLOC((NLFE13)*sizeof(double *))); 02408 col_str = ((double **)MALLOC((NLFE13)*sizeof(double *))); 02409 AulPump = ((double **)MALLOC((NLFE13)*sizeof(double *))); 02410 CollRate = ((double **)MALLOC((NLFE13)*sizeof(double *))); 02411 AulDest = ((double **)MALLOC((NLFE13)*sizeof(double *))); 02412 for( i=0; i<(NLFE13); ++i ) 02413 { 02414 fe.Fe13_wl[i] = ((double *)MALLOC((NLFE13)*sizeof(double ))); 02415 fe.Fe13_emiss[i] = ((double *)MALLOC((NLFE13)*sizeof(double ))); 02416 AulNet[i] = ((double *)MALLOC((NLFE13)*sizeof(double ))); 02417 col_str[i] = ((double *)MALLOC((NLFE13)*sizeof(double ))); 02418 AulPump[i] = ((double *)MALLOC((NLFE13)*sizeof(double ))); 02419 CollRate[i] = ((double *)MALLOC((NLFE13)*sizeof(double ))); 02420 AulDest[i] = ((double *)MALLOC((NLFE13)*sizeof(double ))); 02421 } 02422 02423 /* set some to constant values after zeroing out */ 02424 for( i=0; i<NLFE13; ++i ) 02425 { 02426 create[i] = 0.; 02427 destroy[i] = 0.; 02428 for( j=0; j<NLFE13; ++j ) 02429 { 02430 AulNet[i][j] = 0.; 02431 col_str[i][j] = 0.; 02432 CollRate[i][j] = 0.; 02433 AulDest[i][j] = 0.; 02434 AulPump[i][j] = 0.; 02435 fe.Fe13_wl[i][j] = 0.; 02436 fe.Fe13_emiss[i][j] = 0.; 02437 } 02438 } 02439 /* calculates wavelengths of transitions */ 02440 /* dividing by RefIndex converts the vacuum wavelength to air wavelength */ 02441 for( ihi=1; ihi<NLFE13; ++ihi ) 02442 { 02443 for( ilo=0; ilo<ihi; ++ilo ) 02444 { 02445 fe.Fe13_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / 02446 RefIndex( (excit_wn[ihi]-excit_wn[ilo]) ); 02447 } 02448 } 02449 02450 /*>>refer Fe13 as Huang, K.-N. 1985, At. Data Nucl. Data Tables 32, 503 */ 02451 AulNet[4][1] = 1010.; 02452 AulNet[4][2] = 3.8; 02453 AulNet[4][3] = 8.1; 02454 AulNet[3][1] = 45.7; 02455 AulNet[3][2] = 57.6; 02456 AulNet[2][0] = 0.0063; 02457 AulNet[1][0] = 14.0; 02458 AulNet[2][1] = 9.87; 02459 02460 } 02461 02462 /* bail if no ions */ 02463 if( dense.xIonDense[ipIRON][12] <= 0. ) 02464 { 02465 CoolAdd("Fe13",0,0.); 02466 02467 fe.Fe13CoolTot = 0.; 02468 for( ihi=1; ihi<NLFE13; ++ihi ) 02469 { 02470 for( ilo=0; ilo<ihi; ++ilo ) 02471 { 02472 fe.Fe13_emiss[ihi][ilo] = 0.; 02473 } 02474 } 02475 /* level populations */ 02476 /* LIMLEVELN is the dimension of the atoms vectors */ 02477 ASSERT( NLFE13 <= LIMLEVELN); 02478 for( i=0; i < NLFE13; i++ ) 02479 { 02480 atoms.PopLevels[i] = 0.; 02481 atoms.DepLTELevels[i] = 1.; 02482 } 02483 return; 02484 } 02485 02486 /* evaluate collision strengths if Te changed by too much */ 02487 if( fabs(phycon.te / TeUsed - 1. ) > 0.05 ) 02488 { 02489 TeUsed = phycon.te; 02490 02491 /* collision strengths at current temperature */ 02492 for( ihi=1; ihi<NLFE13; ++ihi ) 02493 { 02494 for( ilo=0; ilo<ihi; ++ilo ) 02495 { 02496 col_str[ihi][ilo] = Fe_10_11_13_cs( 13 , ilo+1 , ihi+1 ); 02497 ASSERT( col_str[ihi][ilo] > 0. ); 02498 } 02499 } 02500 } 02501 02502 /*fprintf(ioQQQ,"DEBUG Fe13 N %.2e %.2e %.2e %.2e %.2e %.2e\n", 02503 col_str[1][0] , col_str[2][1] , col_str[2][0] , 02504 AulNet[1][0] , AulNet[2][1] , AulNet[2][0]);*/ 02505 02506 /* nNegPop positive if negative pops occurred, negative if too cold */ 02507 atom_levelN( 02508 /* number of levels */ 02509 NLFE13, 02510 /* the abundance of the ion, cm-3 */ 02511 dense.xIonDense[ipIRON][12], 02512 /* the statistical weights */ 02513 stat_wght, 02514 /* the excitation energies */ 02515 excit_wn, 02516 'w', 02517 /* the derived populations - cm-3 */ 02518 pops, 02519 /* the derived departure coefficients */ 02520 depart, 02521 /* the net emission rate, Aul * escape prob */ 02522 &AulNet , 02523 /* the collision strengths */ 02524 &col_str , 02525 /* A * destruction prob */ 02526 &AulDest, 02527 /* pumping rate */ 02528 &AulPump, 02529 /* collision rate, s-1, must defined if no collision strengths */ 02530 &CollRate, 02531 /* creation vector */ 02532 create, 02533 /* destruction vector */ 02534 destroy, 02535 /* say atom_levelN should evaluate coll rates from cs */ 02536 false, 02537 &fe.Fe13CoolTot, 02538 &dCool_dT, 02539 "Fe13", 02540 &nNegPop, 02541 &lgZeroPop, 02542 false ); 02543 02544 /* LIMLEVELN is the dim of the PopLevels vector */ 02545 ASSERT( NLFE13 <= LIMLEVELN ); 02546 for( i=0; i<NLFE13; ++i) 02547 { 02548 atoms.PopLevels[i] = pops[i]; 02549 atoms.DepLTELevels[i] = depart[i]; 02550 } 02551 02552 if( nNegPop>0 ) 02553 { 02554 fprintf( ioQQQ, " Fe13Lev5 found negative populations\n" ); 02555 ShowMe(); 02556 cdEXIT(EXIT_FAILURE); 02557 } 02558 02559 /* add cooling then its derivative */ 02560 CoolAdd("Fe13",0,fe.Fe13CoolTot); 02561 /* derivative of cooling */ 02562 thermal.dCooldT += dCool_dT; 02563 02564 /* find emission in each line */ 02565 for( ihi=1; ihi<NLFE13; ++ihi ) 02566 { 02567 for( ilo=0; ilo<ihi; ++ilo ) 02568 { 02569 /* emission in these lines */ 02570 fe.Fe13_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN; 02571 } 02572 } 02573 return; 02574 }