cloudy trunk
Loading...
Searching...
No Matches
cool_iron.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*CoolIron compute iron cooling */
4/*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */
5/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
6/*Fe2_cooling compute cooling due to FeII emission */
7/*Fe3_cs return collision strengths for low level Fe3 lines */
8/*Fe4_cs return collision strengths for low level Fe4 lines */
9/*Fe5_cs return collision strengths for low level Fe5 lines */
10#include "cddefines.h"
11#include "physconst.h"
12#include "dense.h"
13#include "coolheavy.h"
14#include "taulines.h"
15#include "phycon.h"
16#include "iso.h"
17#include "conv.h"
18#include "cooling.h"
19#include "trace.h"
20#include "hydrogenic.h"
21#include "ligbar.h"
22#include "cooling.h"
23#include "thermal.h"
24#include "lines_service.h"
25#include "atoms.h"
26#include "atomfeii.h"
27#include "fe.h"
28
29/*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */
30STATIC void Fe3Lev14(void);
31
32/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
33STATIC void Fe4Lev12(void);
34
35/*Fe2_cooling compute cooling due to FeII emission */
36STATIC void Fe2_cooling( void )
37{
38 const int NLFE2 = 6;
39 long int i, j;
40 int nNegPop;
41
42 static double **AulPump,
43 **CollRate,
44 **AulEscp,
45 **col_str ,
46 **AulDest,
47 *depart,
48 *pops,
49 *destroy,
50 *create;
51
52 static bool lgFirst=true;
53 bool lgZeroPop;
54
55 /* stat weights for Fred's 6 level model FeII atom */
56 static double gFe2[NLFE2]={1.,1.,1.,1.,1.,1.};
57 /* excitation energies (Kelvin) for Fred's 6 level model FeII atom */
58 static double ex[NLFE2]={0.,3.32e4,5.68e4,6.95e4,1.15e5,1.31e5};
59
60 /* these are used to only evaluated FeII one time per temperature, zone,
61 * and abundance */
62 static double TUsed = 0.;
63 static double AbunUsed = 0.;
64 /* remember which zone last evaluated with */
65 static long int nZUsed=-1,
66 /* make sure at least two calls per zone */
67 nCall=0;
68
69 DEBUG_ENTRY( "Fe2_cooling()" );
70
71 /* return if nothing to do */
72 if( dense.xIonDense[ipIRON][1] == 0. )
73 {
74 /* zero abundance, do nothing */
75 /* heating cooling and derivative from large atom */
76 FeII.Fe2_large_cool = 0.;
77 FeII.Fe2_large_heat = 0.;
78 FeII.ddT_Fe2_large_cool = 0.;
79
80 /* cooling and derivative from simple UV atom */
81 FeII.Fe2_UVsimp_cool = 0.;
82 FeII.ddT_Fe2_UVsimp_cool = 0.;
83
84 /* now zero out intensities of all FeII lines and level populations */
86 return;
87 }
88
89 /* evaluate FeII model atom, by default only the lowest 16 levels
90 * but number of levels can be increased with atom feii command */
91
92 /* totally reevaluate FeII atom if new zone, or cooling is significant
93 * and temperature changed, we are in search phase,
94 * lgSlow option set true with atom FeII slow, forces constant
95 * evaluation of atom */
96 if( FeII.lgSlow ||
97 conv.lgFirstSweepThisZone || conv.lgLastSweepThisZone ||
98 /* check whether things have changed on later calls */
99 ( !fp_equal( phycon.te, TUsed ) && fabs(FeII.Fe2_large_cool/thermal.ctot) > 0.002 &&
100 fabs(dense.xIonDense[ipIRON][1]-AbunUsed)/SDIV(AbunUsed) > 0.002 ) ||
101 ( !fp_equal( phycon.te, TUsed ) && fabs(FeII.Fe2_large_cool/thermal.ctot) > 0.01) )
102 {
103
104 if( conv.lgFirstSweepThisZone )
105 /* first call this zone set nCall to zero*/
106 nCall = 0;
107 else
108 /* not first call, increment, check above to make sure at least
109 * two evaluations */
110 ++nCall;
111
112 /* option to trace convergence and FeII calls */
113 if( trace.nTrConvg >= 5 )
114 {
115 fprintf( ioQQQ, " CoolIron5 calling FeIILevelPops since ");
116 if( conv.lgFirstSweepThisZone )
117 {
118 fprintf( ioQQQ,
119 "first sweep this zone." );
120 }
121 else if( ! fp_equal( phycon.te, TUsed ) )
122 {
123 fprintf( ioQQQ,
124 "temperature changed, old new are %g %g, nCall %li ",
125 TUsed, phycon.te , nCall);
126 }
127 else if( nzone != nZUsed )
128 {
129 fprintf( ioQQQ,
130 "new zone, nCall %li ", nCall );
131 }
132 else if( FeII.lgSlow )
133 {
134 fprintf( ioQQQ,
135 "FeII.lgSlow set %li", nCall );
136 }
137 else if( conv.lgSearch )
138 {
139 fprintf( ioQQQ,
140 " in search phase %li", nCall );
141 }
142 else if( nCall < 2 )
143 {
144 fprintf( ioQQQ,
145 "not second nCall %li " , nCall );
146 }
147 else if( ! fp_equal( phycon.te, TUsed ) && FeII.Fe2_large_cool/thermal.ctot > 0.001 )
148 {
149 fprintf( ioQQQ,
150 "temp or cooling changed, new are %g %g nCall %li ",
151 phycon.te, FeII.Fe2_large_cool, nCall );
152 }
153 else
154 {
155 fprintf(ioQQQ, "????");
156 }
157 fprintf(ioQQQ, "\n");
158 }
159
160 /* remember parameters for current conditions */
161 TUsed = phycon.te;
162 AbunUsed = dense.xIonDense[ipIRON][1];
163 nZUsed = nzone;
164
165 /* this print turned on with atom FeII print command */
166 if( FeII.lgPrint )
167 {
168 fprintf(ioQQQ,
169 " FeIILevelPops called zone %4li te %5f abun %10e c(fe/tot):%6f nCall %li\n",
170 nzone,phycon.te,AbunUsed,FeII.Fe2_large_cool/thermal.ctot,nCall);
171 }
172
173 /* this solves the multi-level problem,
174 * sets FeII.Fe2_large_cool,
175 FeII.Fe2_large_heat, &
176 FeII.ddT_Fe2_large_cool
177 but does nothing with them */
178
179 // line RT update, followed by level population solver
180 FeII_RT_Make();
182 {
183 enum{DEBUG_LOC=false};
184 if( DEBUG_LOC && iteration > 1 && nzone >=4 )
185 {
186 fprintf(ioQQQ,"DEBUG1\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
187 nzone,
188 phycon.te,
189 dense.gas_phase[ipHYDROGEN],
190 dense.eden,
191 FeII.Fe2_large_cool ,
192 FeII.ddT_Fe2_large_cool ,
193 FeII.Fe2_large_cool/dense.eden/dense.gas_phase[ipHYDROGEN] ,
194 thermal.ctot );
195 }
196 }
197
198 if( trace.nTrConvg >= 5 || FeII.lgPrint)
199 {
200 /* spacing needed to get proper trace convergence output */
201 fprintf( ioQQQ, " FeIILevelPops5 returned cool=%.2e heat=%.2e derivative=%.2e\n",
202 FeII.Fe2_large_cool,FeII.Fe2_large_heat ,FeII.ddT_Fe2_large_cool);
203 }
204
205 }
206 else if( ! fp_equal( dense.xIonDense[ipIRON][1], AbunUsed ) )
207 {
208 realnum ratio;
209 /* this branch, same zone and temperature, but small change in abundance, so just
210 * rescale cooling and derivative by this change. assumption is that very small changes
211 * in abundance occurs as ots rates damp out */
212 if( trace.nTrConvg >= 5 )
213 {
214 fprintf( ioQQQ,
215 " CoolIron rescaling FeIILevelPops since small change, CFe2=%.2e CTOT=%.2e\n",
216 FeII.Fe2_large_cool,thermal.ctot);
217 }
218 ratio = dense.xIonDense[ipIRON][1]/AbunUsed;
219 FeII.Fe2_large_cool *= ratio;
220 FeII.ddT_Fe2_large_cool *= ratio;
221 FeII.Fe2_large_heat *= ratio;
222 AbunUsed = dense.xIonDense[ipIRON][1];
223 }
224 else
225 {
226 /* this is case where temp is unchanged, so heating and cooling same too */
227 if( trace.nTrConvg >= 5 )
228 {
229 fprintf( ioQQQ, " CoolIron NOT calling FeIILevelPops\n");
230 }
231 }
232
233 /* now update total cooling and its derivative
234 * all paths flow through here */
235 /* FeII.Fecool = FeII.Fe2_large_cool; */
236 CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_large_cool));
237
238 /* add negative cooling to heating stack */
239 thermal.heating[25][27] = MAX2(0.,FeII.Fe2_large_heat);
240
241 /* counts as heating derivative if negative cooling */
242 if( FeII.Fe2_large_cool > 0. )
243 {
244 /* >>chng 01 mar 16, add factor of 3 due to conv problems after changing damper */
245 thermal.dCooldT += 3.*FeII.ddT_Fe2_large_cool;
246 }
247
248 if( trace.lgTrace && trace.lgCoolTr )
249 {
250 fprintf( ioQQQ, " Large FeII returns te, cooling, dc=%11.3e%11.3e%11.3e\n",
251 phycon.te, FeII.Fe2_large_cool, FeII.ddT_Fe2_large_cool );
252 }
253
254 /* >>chng 05 nov 29, still do simple UV atom if only ground term is done */
255 if( !FeII.lgFeIILargeOn )
256 {
257
258 /* following treatment of Fe II follows
259 * >>refer fe2 model Wills, B.J., Wills, D., Netzer, H. 1985, ApJ, 288, 143
260 * all elements are used, and must be set to zero if zero */
261
262 /* set up space for simple model of UV FeII emission */
263 if( lgFirst )
264 {
265 /* will never do this again in this core load */
266 lgFirst = false;
267 /* allocate the 1D arrays*/
268 pops = (double *)MALLOC( sizeof(double)*(NLFE2) );
269 create = (double *)MALLOC( sizeof(double)*(NLFE2) );
270 destroy = (double *)MALLOC( sizeof(double)*(NLFE2) );
271 depart = (double *)MALLOC( sizeof(double)*(NLFE2) );
272 /* create space for the 2D arrays */
273 AulPump = ((double **)MALLOC((NLFE2)*sizeof(double *)));
274 CollRate = ((double **)MALLOC((NLFE2)*sizeof(double *)));
275 AulDest = ((double **)MALLOC((NLFE2)*sizeof(double *)));
276 AulEscp = ((double **)MALLOC((NLFE2)*sizeof(double *)));
277 col_str = ((double **)MALLOC((NLFE2)*sizeof(double *)));
278
279 for( i=0; i < NLFE2; ++i )
280 {
281 AulPump[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
282 CollRate[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
283 AulDest[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
284 AulEscp[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
285 col_str[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
286 }
287 }
288
289 /*zero out all arrays, then check that upper diagonal remains zero below */
290 for( i=0; i < NLFE2; i++ )
291 {
292 create[i] = 0.;
293 destroy[i] = 0.;
294 for( j=0; j < NLFE2; j++ )
295 {
296 /*data[j][i] = 0.;*/
297 col_str[j][i] = 0.;
298 AulEscp[j][i] = 0.;
299 AulDest[j][i] = 0.;
300 AulPump[j][i] = 0.;
301 }
302 }
303
304 /* now put in real data for lines */
305 AulEscp[1][0] = 1.;
306 AulEscp[2][0] = ( TauLines[ipTuv3].Emis().Pesc() + TauLines[ipTuv3].Emis().Pelec_esc())*TauLines[ipTuv3].Emis().Aul();
307 AulDest[2][0] = TauLines[ipTuv3].Emis().Pdest()*TauLines[ipTuv3].Emis().Aul();
308 AulPump[0][2] = TauLines[ipTuv3].Emis().pump();
309
310 AulEscp[5][0] = (TauLines[ipTFe16].Emis().Pesc() + TauLines[ipTFe16].Emis().Pelec_esc())*TauLines[ipTFe16].Emis().Aul();
311 AulDest[5][0] = TauLines[ipTFe16].Emis().Pdest()*TauLines[ipTFe16].Emis().Aul();
312 /* continuum pumping of n=6 */
313 AulPump[0][5] = TauLines[ipTFe16].Emis().pump();
314 /* Ly-alpha pumping */
315
316 double PumpLyaFeII = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()*
317 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul()*
318 hydro.dstfe2lya/SDIV(dense.xIonDense[ipIRON][1]);
319
320 AulPump[0][5] += PumpLyaFeII;
321
322 AulEscp[2][1] = (TauLines[ipTr48].Emis().Pesc() + TauLines[ipTr48].Emis().Pelec_esc())*TauLines[ipTr48].Emis().Aul();
323 AulDest[2][1] = TauLines[ipTr48].Emis().Pdest()*TauLines[ipTr48].Emis().Aul();
324 AulPump[1][2] = TauLines[ipTr48].Emis().pump();
325
326 AulEscp[5][1] = (TauLines[ipTFe26].Emis().Pesc() + TauLines[ipTFe26].Emis().Pelec_esc())*TauLines[ipTFe26].Emis().Aul();
327 AulDest[5][1] = TauLines[ipTFe26].Emis().Pdest()*TauLines[ipTFe26].Emis().Aul();
328 AulPump[1][5] = TauLines[ipTFe26].Emis().pump();
329
330 AulEscp[3][2] = (TauLines[ipTFe34].Emis().Pesc() + TauLines[ipTFe34].Emis().Pelec_esc())*TauLines[ipTFe34].Emis().Aul();
331 AulDest[3][2] = TauLines[ipTFe34].Emis().Pdest()*TauLines[ipTFe34].Emis().Aul();
332 AulPump[2][3] = TauLines[ipTFe34].Emis().pump();
333
334 AulEscp[4][2] = (TauLines[ipTFe35].Emis().Pesc() + TauLines[ipTFe35].Emis().Pelec_esc())*TauLines[ipTFe35].Emis().Aul();
335 AulDest[4][2] = TauLines[ipTFe35].Emis().Pdest()*TauLines[ipTFe35].Emis().Aul();
336 AulPump[2][4] = TauLines[ipTFe35].Emis().pump();
337
338 AulEscp[5][3] = (TauLines[ipTFe46].Emis().Pesc() + TauLines[ipTFe46].Emis().Pelec_esc())*TauLines[ipTFe46].Emis().Aul();
339 AulDest[5][3] = TauLines[ipTFe46].Emis().Pdest()*TauLines[ipTFe46].Emis().Aul();
340 AulPump[3][5] = TauLines[ipTFe46].Emis().pump();
341
342 AulEscp[5][4] = (TauLines[ipTFe56].Emis().Pesc() + TauLines[ipTFe56].Emis().Pelec_esc())*TauLines[ipTFe56].Emis().Aul();
343 AulDest[5][4] = TauLines[ipTFe56].Emis().Pdest()*TauLines[ipTFe56].Emis().Aul();
344 AulPump[4][5] = TauLines[ipTFe56].Emis().pump();
345
346 /* these are collision strengths */
347 col_str[1][0] = 1.;
348 col_str[2][0] = 12.;
349 col_str[3][0] = 1.;
350 col_str[4][0] = 1.;
351 col_str[5][0] = 12.;
352 col_str[2][1] = 6.;
353 col_str[3][1] = 1.;
354 col_str[4][1] = 1.;
355 col_str[5][1] = 12.;
356 col_str[3][2] = 6.;
357 col_str[4][2] = 12.;
358 col_str[5][2] = 1.;
359 col_str[4][3] = 1.;
360 col_str[5][3] = 12.;
361 col_str[5][4] = 6.;
362
363 /*void atom_levelN(long,long,realnum,double[],double[],double[],double*,
364 double*,double*,long*,realnum*,realnum*,STRING,int*);*/
365 atom_levelN(NLFE2,
366 dense.xIonDense[ipIRON][1],
367 gFe2,
368 ex,
369 'K',
370 pops,
371 depart,
372 &AulEscp ,
373 &col_str,
374 &AulDest,
375 &AulPump,
376 &CollRate,
377 create,
378 destroy,
379 false,/* say atom_levelN should evaluate coll rates from cs */
380 /*&&ipdest,*/
381 &FeII.Fe2_UVsimp_cool,
382 &FeII.ddT_Fe2_UVsimp_cool,
383 "FeII",
384 &nNegPop,
385 &lgZeroPop,
386 false );
387
388 /* nNegPop positive if negative pops occurred, negative if too cold */
389 if( nNegPop > 0 /*negative if too cold - that is not negative and is OK */ )
390 {
391 fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population for simple UV FeII.\n");
392 }
393
394 /* add heating - cooling by this process */;
395 CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_UVsimp_cool));
396 thermal.heating[25][27] = MAX2(0.,-FeII.Fe2_UVsimp_cool);
397 thermal.dCooldT += FeII.ddT_Fe2_UVsimp_cool;
398
399 /* LIMLEVELN is the dim of the PopLevels vector */
400 ASSERT( NLFE2 <= LIMLEVELN );
401 for( i=0; i < NLFE2; ++i )
402 {
403 atoms.PopLevels[i] = pops[i];
404 atoms.DepLTELevels[i] = depart[i];
405 }
406
407 (*TauLines[ipTuv3].Lo()).Pop() = pops[0];
408 (*TauLines[ipTuv3].Hi()).Pop() = pops[2];
409 TauLines[ipTuv3].Emis().PopOpc() = (pops[0] - pops[2]);
410 TauLines[ipTuv3].Emis().phots() = pops[2]*AulEscp[2][0];
411 TauLines[ipTuv3].Emis().xIntensity() =
412 TauLines[ipTuv3].Emis().phots()*TauLines[ipTuv3].EnergyErg();
413
414 (*TauLines[ipTr48].Lo()).Pop() = pops[1];
415 (*TauLines[ipTr48].Hi()).Pop() = pops[2];
416 TauLines[ipTr48].Emis().PopOpc() = (pops[1] - pops[2]);
417 TauLines[ipTr48].Emis().phots() = pops[2]*AulEscp[2][1];
418 TauLines[ipTr48].Emis().xIntensity() =
419 TauLines[ipTr48].Emis().phots()*TauLines[ipTr48].EnergyErg();
420
421 FeII.for7 = pops[1]*AulEscp[1][0]*4.65e-12;
422
423 (*TauLines[ipTFe16].Lo()).Pop() = pops[0];
424 (*TauLines[ipTFe16].Hi()).Pop() = pops[5];
425 /* Lyman alpha optical depths are not known on first iteration,
426 * inward optical depths used, so line trapping overestimated,
427 * this can cause artificial maser in FeII - prevent by not
428 * including stimulated emission correction on first iteration */
429 TauLines[ipTFe16].Emis().PopOpc() = (pops[0] - pops[5]);
430 TauLines[ipTFe16].Emis().phots() = pops[5]*AulEscp[5][0];
431 TauLines[ipTFe16].Emis().xIntensity() =
432 TauLines[ipTFe16].Emis().phots()*TauLines[ipTFe16].EnergyErg();
433
434 (*TauLines[ipTFe26].Lo()).Pop() = pops[1];
435 (*TauLines[ipTFe26].Hi()).Pop() = pops[5];
436 TauLines[ipTFe26].Emis().PopOpc() = (pops[1] - pops[5]);
437 TauLines[ipTFe26].Emis().phots() = pops[5]*AulEscp[5][1];
438 TauLines[ipTFe26].Emis().xIntensity() =
439 TauLines[ipTFe26].Emis().phots()*TauLines[ipTFe26].EnergyErg();
440
441 (*TauLines[ipTFe34].Lo()).Pop() = pops[2];
442 (*TauLines[ipTFe34].Hi()).Pop() = pops[3];
443 TauLines[ipTFe34].Emis().PopOpc() = (pops[2] - pops[3]);
444 TauLines[ipTFe34].Emis().phots() = pops[3]*AulEscp[3][2];
445 TauLines[ipTFe34].Emis().xIntensity() =
446 TauLines[ipTFe34].Emis().phots()*TauLines[ipTFe34].EnergyErg();
447
448 (*TauLines[ipTFe35].Lo()).Pop() = pops[2];
449 (*TauLines[ipTFe35].Hi()).Pop() = pops[4];
450 TauLines[ipTFe35].Emis().PopOpc() = (pops[2] - pops[4]);
451 TauLines[ipTFe35].Emis().phots() = pops[4]*AulEscp[4][2];
452 TauLines[ipTFe35].Emis().xIntensity() =
453 TauLines[ipTFe35].Emis().phots()*TauLines[ipTFe35].EnergyErg();
454
455 (*TauLines[ipTFe46].Lo()).Pop() = pops[3];
456 (*TauLines[ipTFe46].Hi()).Pop() = pops[5];
457 TauLines[ipTFe46].Emis().PopOpc() = (pops[3] - pops[5]);
458 TauLines[ipTFe46].Emis().phots() = pops[5]*AulEscp[5][3];
459 TauLines[ipTFe46].Emis().xIntensity() =
460 TauLines[ipTFe46].Emis().phots()*TauLines[ipTFe46].EnergyErg();
461
462 (*TauLines[ipTFe56].Lo()).Pop() = pops[4];
463 (*TauLines[ipTFe56].Hi()).Pop() = pops[5];
464 TauLines[ipTFe56].Emis().PopOpc() = (pops[4] - pops[5]);
465 TauLines[ipTFe56].Emis().phots() = pops[5]*AulEscp[5][4];
466 TauLines[ipTFe56].Emis().xIntensity() =
467 TauLines[ipTFe56].Emis().phots()*TauLines[ipTFe56].EnergyErg();
468
469 /* Jack's funny FeII lines, data from
470 * >>refer fe2 energy Johansson, S., Brage, T., Leckrone, D.S., Nave, G. &
471 * >>refercon Wahlgren, G.M. 1995, ApJ 446, 361 */
472 PutCS(10.,TauLines[ipT191]);
474 }
475
476 {
477 /*@-redef@*/
478 enum{DEBUG_LOC=false};
479 /*@+redef@*/
480 if( DEBUG_LOC && iteration > 1 && nzone >=4 )
481 {
482 fprintf(ioQQQ,"DEBUG2\t%.2e\t%.2e\t%.2e\n",
483 phycon.te,
484 FeII.Fe2_large_cool ,
485 FeII.Fe2_UVsimp_cool );
486 }
487 }
488
489 return;
490
491}
492
493/*CoolIron - calculation total cooling due to Fe */
494void CoolIron(void)
495{
496 realnum rate;
497
498 DEBUG_ENTRY( "CoolIron()" );
499
500 /* cooling by FeI 24m, 34.2 m */
501 /* >>chng 03 nov 15, add these lines */
505 /*>>refer Fe1 cs Hollenbach & McKee 89 */
506 /* the 24.0 micron line */
507 rate = (realnum)(1.2e-7 * dense.eden +
508 /* >>chng 05 jul 05, eden to cdsqte */
509 /*8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
510 8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
512
513 rate = (realnum)(9.3e-8 * dense.eden +
514 /* >>chng 05 jul 05, eden to cdsqte */
515 /*5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
516 5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
518
519 rate = (realnum)(1.2e-7 * dense.eden +
520 /* >>chng 05 jul 05, eden to cdsqte */
521 /*6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
522 6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
523 (*(*TauDummy).Hi()).g() = (*TauLines[ipFe1_35m].Hi()).g();
524 LineConvRate2CS( *TauDummy , rate );
525 /* this says that line is a dummy, not real one */
526 (*(*TauDummy).Hi()).g() = 0.;
527
529
530 /* series of FeI lines from Dima Verner's list, each 2-lev atom
531 *
532 * Fe I 3884 */
535
536 /* Fe I 3729 */
539
540 /* Fe I 3457 */
543
544 /* Fe I 3021 */
547
548 /* Fe I 2966 */
551
552 /* >>chng 05 dec 03, move Fe2 FeII Fe II cooling into separate routine */
553 Fe2_cooling();
554
555 /* lump 3p and 3f together; cs=
556 * >>refer fe3 as Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384
557 * A from
558 * >>refer fe3 as Garstang, R.H., 1957, Vistas in Astronomy, 1, 268
559 * FE III 5270, is 20.9% of total
560 * >>chng 05 feb 18, Kevin Blagrave email
561 * average wavelength is 4823 with statistical weight averaging of upper energy level,
562 * as per , change 5th number from 2.948 to 2.984, also photon energy
563 * from 3.78 to 4.12 */
564
565 /* >>chng 05 dec 16, FeIII update by Kevin Blagrave */
566 /* FeIII 1122 entire multiplet - atomic data=A's Dima, CS = guess */
567 PutCS(25.,TauLines[ipT1122]);
569
570 /* call 14 level atom */
571 Fe3Lev14();
572
573 /* call 12 level atom */
574 Fe4Lev12();
575
576 /* FE V 3892 + 3839, data from Shields */
577 CoolHeavy.c3892 = atom_pop2(7.4,25.,5.,0.6,3.7e4,dense.xIonDense[ipIRON][4])*
578 5.11e-12;
579 CoolAdd("Fe 5",3892,CoolHeavy.c3892);
580
581 return;
582}
583
584/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
585STATIC void Fe4Lev12(void)
586{
587 const int NLFE4 = 12;
588 bool lgZeroPop;
589 int nNegPop;
590 long int i,
591 j;
592 static bool lgFirst=true;
593
594 double dfe4dt;
595
596 /*static long int **ipdest; */
597 static double
598 **AulEscp,
599 **col_str,
600 **AulDest,
601 depart[NLFE4],
602 pops[NLFE4],
603 destroy[NLFE4],
604 create[NLFE4],
605 **CollRate,
606 **AulPump;
607
608 static const double Fe4A[NLFE4][NLFE4] = {
609 {0.,0.,0.,1.e-5,0.,1.368,.89,0.,1.3e-3,1.8e-4,.056,.028},
610 {0.,0.,2.6e-8,0.,0.,0.,0.,0.,1.7e-7,0.,0.,0.},
611 {0.,0.,0.,0.,3.5e-7,6.4e-10,0.,0.,6.315e-4,0.,6.7e-7,0.},
612 {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},
613 {0.,0.,0.,0.,0.,1.5e-5,1.3e-9,0.,7.6e-4,0.,1.1e-6,6.0e-7},
614 {0.,0.,0.,0.,0.,0.,1.1e-5,1.2e-13,.038,9.9e-7,.022,.018},
615 {0.,0.,0.,0.,0.,0.,0.,3.7e-5,2.9e-6,.034,3.5e-3,.039},
616 {0.,0.,0.,0.,0.,0.,0.,0.,0.,.058,3.1e-6,1.4e-3},
617 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-4,3.1e-14},
618 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.9e-19,1.0e-5},
619 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-7},
620 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}
621 };
622
623 static const double gfe4[NLFE4]={6.,12.,10.,6.,8.,6.,4.,2.,8.,2.,6.,4.};
624
625 /* excitation energies in Kelvin
626 static double ex[NLFE4]={0.,46395.8,46464.,46475.6,46483.,50725.,
627 50838.,50945.,55796.,55966.,56021.,56025.};*/
628 /*>>refer Fe3 energies version 3 NIST Atomic Spectra Database */
629 /*>>chng 05 dec 17, from Kelvin above to excitation energies in wn */
630 static const double excit_wn[NLFE4]={0.,32245.5,32292.8,32301.2,32305.7,35253.8,
631 35333.3,35406.6,38779.4,38896.7,38935.1,38938.2};
632
633 DEBUG_ENTRY( "Fe4Lev12()" );
634
635 if( lgFirst )
636 {
637 /* will never do this again */
638 lgFirst = false;
639 /* allocate the 1D arrays*/
640 /* create space for the 2D arrays */
641 AulPump = ((double **)MALLOC((NLFE4)*sizeof(double *)));
642 CollRate = ((double **)MALLOC((NLFE4)*sizeof(double *)));
643 AulDest = ((double **)MALLOC((NLFE4)*sizeof(double *)));
644 AulEscp = ((double **)MALLOC((NLFE4)*sizeof(double *)));
645 col_str = ((double **)MALLOC((NLFE4)*sizeof(double *)));
646 for( i=0; i < NLFE4; ++i )
647 {
648 AulPump[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
649 CollRate[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
650 AulDest[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
651 AulEscp[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
652 col_str[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
653 }
654 }
655
656 /* bail if no Fe+3 */
657 if( dense.xIonDense[ipIRON][3] <= 0. )
658 {
659 fe.Fe4CoolTot = 0.;
660 fe.fe40401 = 0.;
661 fe.fe42836 = 0.;
662 fe.fe42829 = 0.;
663 fe.fe42567 = 0.;
664 fe.fe41207 = 0.;
665 fe.fe41206 = 0.;
666 fe.fe41106 = 0.;
667 fe.fe41007 = 0.;
668 fe.fe41008 = 0.;
669 fe.fe40906 = 0.;
670 CoolAdd("Fe 4",0,0.);
671
672 /* level populations */
673 /* LIMLEVELN is the dimension of the atoms vectors */
674 ASSERT( NLFE4 <= LIMLEVELN);
675 for( i=0; i < NLFE4; i++ )
676 {
677 atoms.PopLevels[i] = 0.;
678 atoms.DepLTELevels[i] = 1.;
679 }
680 return;
681 }
682 /* number of levels in model ion */
683
684 /* these are in wavenumbers
685 * data excit_wn/ 0., 32245.5, 32293., 32301.2,
686 * 1 32306., 35254., 35333., 35407., 38779., 38897., 38935.,
687 * 2 38938./
688 * excitation energies in Kelvin */
689
690 /* A's are from Garstang, R.H., MNRAS 118, 572 (1958).
691 * each set is for a lower level indicated by second element in array,
692 * index runs over upper level
693 * A's are saved into arrays as data(up,lo) */
694
695 /* collision strengths from Berrington and Pelan Ast Ap S 114, 367.
696 * order is cs(low,up) */
697
698 /* all elements are used, and must be set to zero if zero */
699 for( i=0; i < NLFE4; i++ )
700 {
701 create[i] = 0.;
702 destroy[i] = 0.;
703 for( j=0; j < NLFE4; j++ )
704 {
705 col_str[j][i] = 0.;
706 AulEscp[j][i] = 0.;
707 AulDest[j][i] = 0.;
708 AulPump[j][i] = 0.;
709 }
710 }
711
712 /* fill in Einstein As and collision strengths */
713 for( long ipHi=1; ipHi < NLFE4; ipHi++ )
714 {
715 for( long ipLo=0; ipLo < ipHi; ipLo++ )
716 {
717 AulEscp[ipHi][ipLo] = Fe4A[ipLo][ipHi];
718 col_str[ipHi][ipLo] = Fe4_cs(ipLo, ipHi);
719 }
720 }
721
722 /* leveln itself is well-protected against zero abundances,
723 * low temperatures */
724
725 atom_levelN(NLFE4,
726 dense.xIonDense[ipIRON][3],
727 gfe4,
728 excit_wn,
729 'w',
730 pops,
731 depart,
732 &AulEscp ,
733 &col_str ,
734 &AulDest,
735 &AulPump,
736 &CollRate,
737 create,
738 destroy,
739 /* say atom_levelN should evaluate coll rates from cs */
740 false,
741 &fe.Fe4CoolTot,
742 &dfe4dt,
743 "FeIV",
744 /* nNegPop positive if negative pops occured, negative if too cold */
745 &nNegPop,
746 &lgZeroPop,
747 false );
748
749 /* LIMLEVELN is the dim of the PopLevels vector */
750 ASSERT( NLFE4 <= LIMLEVELN );
751 for( i=0; i < NLFE4; ++i )
752 {
753 atoms.PopLevels[i] = pops[i];
754 atoms.DepLTELevels[i] = depart[i];
755 }
756
757 if( nNegPop > 0 )
758 {
759 fprintf( ioQQQ, " fe4levl2 found negative populations\n" );
760 ShowMe();
762 }
763
764 CoolAdd("Fe 4",0,fe.Fe4CoolTot);
765
766 thermal.dCooldT += dfe4dt;
767
768 /* three transitions hst can observe
769 * 4 -1 3095.9A and 5 -1 3095.9A */
770 fe.fe40401 = (pops[3]*Fe4A[0][3]*(excit_wn[3] - excit_wn[0]) +
771 pops[4]*Fe4A[0][4]*(excit_wn[4] - excit_wn[0]) )*T1CM*BOLTZMANN;
772
773 fe.fe42836 = pops[5]*Fe4A[0][5]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN;
774
775 fe.fe42829 = pops[6]*Fe4A[0][6]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN;
776
777 fe.fe42567 = (pops[10]*Fe4A[0][10]*(excit_wn[10] - excit_wn[0]) +
778 pops[11]*Fe4A[0][11]*(excit_wn[10] - excit_wn[0]))*T1CM*BOLTZMANN;
779
780 fe.fe41207 = pops[11]*Fe4A[6][11]*(excit_wn[11] - excit_wn[6])*T1CM*BOLTZMANN;
781 fe.fe41206 = pops[11]*Fe4A[5][11]*(excit_wn[11] - excit_wn[5])*T1CM*BOLTZMANN;
782 fe.fe41106 = pops[10]*Fe4A[5][10]*(excit_wn[10] - excit_wn[5])*T1CM*BOLTZMANN;
783 fe.fe41007 = pops[9]*Fe4A[6][9]*(excit_wn[9] - excit_wn[6])*T1CM*BOLTZMANN;
784 fe.fe41008 = pops[9]*Fe4A[7][9]*(excit_wn[9] - excit_wn[7])*T1CM*BOLTZMANN;
785 fe.fe40906 = pops[8]*Fe4A[5][8]*(excit_wn[8] - excit_wn[5])*T1CM*BOLTZMANN;
786 return;
787}
788
789/*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion
790 * >>chng 05 dec 17, code provided by Kevin Blagrave and described in
791 *>>refer Fe3 model Blagrave, K.P.M., Martin, P.G. & Baldwin, J.A.
792 *>>refercon 2006, ApJ, 644, 1006B (astro-ph 0603167) */
793STATIC void Fe3Lev14(void)
794{
795 bool lgZeroPop;
796 int nNegPop;
797 long int i,
798 j;
799 static bool lgFirst=true;
800
801 double dfe3dt;
802
803 long int ihi , ilo;
804 static double
805 *depart,
806 *pops,
807 *destroy,
808 *create ,
809 **AulDest,
810 **CollRate,
811 **AulPump,
812 **AulNet,
813 **col_str;
814
815 /* statistical weights for the n levels */
816 static double gfe3[NLFE3]={9.,7.,5.,3.,1.,5.,13.,11.,9.,3.,1.,9.,7.,5.};
817
818 /*refer Fe3 energies NIST version 3 Atomic Spectra Database */
819 /* from smallest to largest energy (not in multiplet groupings) */
820
821 /* energy in wavenumbers */
822 static double excit_wn[NLFE3]={
823 0.0 , 436.2, 738.9, 932.4, 1027.3,
824 19404.8, 20051.1, 20300.8, 20481.9, 20688.4,
825 21208.5, 21462.2, 21699.9, 21857.2 };
826
827 DEBUG_ENTRY( "Fe3Lev14()" );
828
829 if( lgFirst )
830 {
831 /* will never do this again */
832 lgFirst = false;
833 /* allocate the 1D arrays*/
834 /* create space for the 2D arrays */
835 depart = ((double *)MALLOC((NLFE3)*sizeof(double)));
836 pops = ((double *)MALLOC((NLFE3)*sizeof(double)));
837 destroy = ((double *)MALLOC((NLFE3)*sizeof(double)));
838 create = ((double *)MALLOC((NLFE3)*sizeof(double)));
839 /* now the 2-d arrays */
840 fe.Fe3_wl = ((double **)MALLOC((NLFE3)*sizeof(double *)));
841 fe.Fe3_emiss = ((double **)MALLOC((NLFE3)*sizeof(double *)));
842 AulNet = ((double **)MALLOC((NLFE3)*sizeof(double *)));
843 col_str = ((double **)MALLOC((NLFE3)*sizeof(double *)));
844 AulPump = ((double **)MALLOC((NLFE3)*sizeof(double *)));
845 CollRate = ((double **)MALLOC((NLFE3)*sizeof(double *)));
846 AulDest = ((double **)MALLOC((NLFE3)*sizeof(double *)));
847 for( i=0; i < NLFE3; ++i )
848 {
849 fe.Fe3_wl[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
850 fe.Fe3_emiss[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
851 AulNet[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
852 col_str[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
853 AulPump[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
854 CollRate[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
855 AulDest[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
856 }
857
858 /* set some to constant values after zeroing out */
859 for( i=0; i < NLFE3; ++i )
860 {
861 create[i] = 0.;
862 destroy[i] = 0.;
863 for( j=0; j < NLFE3; ++j )
864 {
865 AulNet[i][j] = 0.;
866 col_str[i][j] = 0.;
867 CollRate[i][j] = 0.;
868 AulDest[i][j] = 0.;
869 AulPump[i][j] = 0.;
870 fe.Fe3_wl[i][j] = 0.;
871 fe.Fe3_emiss[i][j] = 0.;
872 }
873 }
874 /* calculates wavelengths of transitions */
875 /* dividing by RefIndex converts the vacuum wavelength to air wavelength */
876 for( ihi=1; ihi < NLFE3; ++ihi )
877 {
878 for( ilo=0; ilo < ihi; ++ilo )
879 {
880 fe.Fe3_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) /
881 RefIndex( (excit_wn[ihi]-excit_wn[ilo]) );
882 }
883 }
884
885 /* assume FeIII is optically thin - just use As as net escape */
886 /*>>refer Fe3 as Quinet, P., 1996, A&AS, 116, 573 */
887 AulNet[1][0] = 2.8e-3;
888 AulNet[7][0] = 4.9e-6;
889 AulNet[8][0] = 5.7e-3;
890 AulNet[11][0] = 4.5e-1;
891 AulNet[12][0] = 4.2e-2;
892
893 AulNet[2][1] = 1.8e-3;
894 AulNet[5][1] = 4.2e-1;
895 AulNet[8][1] = 1.0e-3;
896 AulNet[11][1] = 8.4e-2;
897 AulNet[12][1] = 2.5e-1;
898 AulNet[13][1] = 2.7e-2;
899
900 AulNet[3][2] = 7.0e-4;
901 AulNet[5][2] = 5.1e-5;
902 AulNet[9][2] = 5.4e-1;
903 AulNet[12][2] = 8.5e-2;
904 AulNet[13][2] = 9.8e-2;
905
906 AulNet[4][3] = 1.4e-4;
907 AulNet[5][3] = 3.9e-2;
908 AulNet[9][3] = 4.1e-5;
909 AulNet[10][3] = 7.0e-1;
910 AulNet[13][3] = 4.7e-2;
911
912 AulNet[9][4] = 9.3e-2;
913
914 AulNet[9][5] = 4.7e-2;
915 AulNet[12][5] = 2.5e-6;
916 AulNet[13][5] = 1.7e-5;
917
918 AulNet[7][6] = 2.7e-4;
919
920 AulNet[8][7] = 1.2e-4;
921 AulNet[11][7] = 6.6e-4;
922
923 AulNet[11][8] = 1.6e-3;
924 AulNet[12][8] = 7.8e-4;
925
926 AulNet[10][9] = 8.4e-3;
927 AulNet[13][9] = 2.8e-7;
928
929 AulNet[12][11] = 3.0e-4;
930
931 AulNet[13][12] = 1.4e-4;
932
933 for( int ipHi = 1; ipHi < NLFE3; ipHi++)
934 {
935 for( int ipLo = 0; ipLo < ipHi; ipLo++)
936 {
937 col_str[ipHi][ipLo] = Fe3_cs(ipLo,ipHi);
938 }
939 }
940 }
941
942 /* bail if no ions */
943 if( dense.xIonDense[ipIRON][2] <= 0. )
944 {
945 CoolAdd("Fe 3",0,0.);
946
947 fe.Fe3CoolTot = 0.;
948 for( ihi=1; ihi < NLFE3; ++ihi )
949 {
950 for( ilo=0; ilo < ihi; ++ilo )
951 {
952 fe.Fe3_emiss[ihi][ilo] = 0.;
953 }
954 }
955 /* level populations */
956 /* LIMLEVELN is the dimension of the atoms vectors */
958 for( i=0; i < NLFE3; i++ )
959 {
960 atoms.PopLevels[i] = 0.;
961 atoms.DepLTELevels[i] = 1.;
962 }
963 return;
964 }
965
966 /* nNegPop positive if negative pops occurred, negative if too cold */
968 /* number of levels */
969 NLFE3,
970 /* the abundance of the ion, cm-3 */
971 dense.xIonDense[ipIRON][2],
972 /* the statistical weights */
973 gfe3,
974 /* the excitation energies */
975 excit_wn,
976 'w',
977 /* the derived populations - cm-3 */
978 pops,
979 /* the derived departure coefficients */
980 depart,
981 /* the net emission rate, Aul * escape prob */
982 &AulNet ,
983 /* the collision strengths */
984 &col_str ,
985 /* A * destruction prob */
986 &AulDest,
987 /* pumping rate */
988 &AulPump,
989 /* collision rate, s-1, must defined if no collision strengths */
990 &CollRate,
991 /* creation vector */
992 create,
993 /* destruction vector */
994 destroy,
995 /* say atom_levelN should evaluate coll rates from cs */
996 false,
997 &fe.Fe3CoolTot,
998 &dfe3dt,
999 "Fe 3",
1000 &nNegPop,
1001 &lgZeroPop,
1002 false );
1003
1004 /* LIMLEVELN is the dim of the PopLevels vector */
1005 ASSERT( NLFE3 <= LIMLEVELN );
1006 for( i=0; i < NLFE3; ++i )
1007 {
1008 atoms.PopLevels[i] = pops[i];
1009 atoms.DepLTELevels[i] = depart[i];
1010 }
1011
1012 if( nNegPop > 0 )
1013 {
1014 fprintf( ioQQQ, " Fe3Lev14 found negative populations\n" );
1015 ShowMe();
1017 }
1018
1019 /* add cooling then its derivative */
1020 CoolAdd("Fe 3",0,fe.Fe3CoolTot);
1021 /* derivative of cooling */
1022 thermal.dCooldT += dfe3dt;
1023
1024 /* find emission in each line */
1025 for( ihi=1; ihi < NLFE3; ++ihi )
1026 {
1027 for( ilo=0; ilo < ihi; ++ilo )
1028 {
1029 /* emission in these lines */
1030 fe.Fe3_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN;
1031 }
1032 }
1033 return;
1034}
1035
1036double Fe3_cs(long ipLo, long ipHi)
1037{
1038 static double col_str[14][14];
1039
1040 static double lgOneTimeMustInit=true;
1041 if( lgOneTimeMustInit )
1042 {
1043 lgOneTimeMustInit = false;
1044 /* collision strengths at log T 4 */
1046 /*>>refer Fe3 cs Zhang, H. 1996, 119, 523 */
1047 col_str[1][0] = 2.92;
1048 col_str[2][0] = 1.24;
1049 col_str[3][0] = 0.595;
1050 col_str[4][0] = 0.180;
1051 col_str[5][0] = 0.580;
1052 col_str[6][0] = 1.34;
1053 col_str[7][0] = 0.489;
1054 col_str[8][0] = 0.0926;
1055 col_str[9][0] = 0.165;
1056 col_str[10][0] = 0.0213;
1057 col_str[11][0] = 1.07;
1058 col_str[12][0] = 0.435;
1059 col_str[13][0] = 0.157;
1060
1061 col_str[2][1] = 2.06;
1062 col_str[3][1] = 0.799;
1063 col_str[4][1] = 0.225;
1064 col_str[5][1] = 0.335;
1065 col_str[6][1] = 0.555;
1066 col_str[7][1] = 0.609;
1067 col_str[8][1] = 0.367;
1068 col_str[9][1] = 0.195;
1069 col_str[10][1] = 0.0698;
1070 col_str[11][1] = 0.538;
1071 col_str[12][1] = 0.484;
1072 col_str[13][1] = 0.285;
1073
1074 col_str[3][2] = 1.29;
1075 col_str[4][2] = 0.312;
1076 col_str[5][2] = 0.173;
1077 col_str[6][2] = 0.178;
1078 col_str[7][2] = 0.430;
1079 col_str[8][2] = 0.486;
1080 col_str[9][2] = 0.179;
1081 col_str[10][2] = 0.0741;
1082 col_str[11][2] = 0.249;
1083 col_str[12][2] = 0.362;
1084 col_str[13][2] = 0.324;
1085
1086 col_str[4][3] = 0.493;
1087 col_str[5][3] = 0.0767;
1088 col_str[6][3] = 0.0348;
1089 col_str[7][3] = 0.223;
1090 col_str[8][3] = 0.401;
1091 col_str[9][3] = 0.126;
1092 col_str[10][3] = 0.0528;
1093 col_str[11][3] = 0.101;
1094 col_str[12][3] = 0.207;
1095 col_str[13][3] = 0.253;
1096
1097 col_str[5][4] = 0.0211;
1098 col_str[6][4] = 0.00122;
1099 col_str[7][4] = 0.0653;
1100 col_str[8][4] = 0.154;
1101 col_str[9][4] = 0.0453;
1102 col_str[10][4] = 0.0189;
1103 col_str[11][4] = 0.0265;
1104 col_str[12][4] = 0.0654;
1105 col_str[13][4] = 0.0950;
1106
1107 col_str[6][5] = 0.403;
1108 col_str[7][5] = 0.213;
1109 col_str[8][5] = 0.0939;
1110 col_str[9][5] = 1.10;
1111 col_str[10][5] = 0.282;
1112 col_str[11][5] = 0.942;
1113 col_str[12][5] = 0.768;
1114 col_str[13][5] = 0.579;
1115
1116 col_str[7][6] = 2.84; /* 10-9 */
1117 col_str[8][6] = 0.379; /* 11-9 */
1118 col_str[9][6] = 0.0876; /* 7-9 */
1119 col_str[10][6] = 0.00807; /* 8-9 */
1120 col_str[11][6] = 1.85; /* 12-9 */
1121 col_str[12][6] = 0.667; /* 13-9 */
1122 col_str[13][6] = 0.0905; /* 14-9 */
1123
1124 col_str[8][7] = 3.07; /* 11-10 */
1125 col_str[9][7] = 0.167; /* 7-10 */
1126 col_str[10][7] = 0.0526; /* 8-10 */
1127 col_str[11][7] = 0.814; /* 12-10 */
1128 col_str[12][7] = 0.837; /* 13-10 */
1129 col_str[13][7] = 0.626; /* 14-10 */
1130
1131 col_str[9][8] = 0.181; /* 7-11 */
1132 col_str[10][8] = 0.0854; /* 8-11 */
1133 col_str[11][8] = 0.180; /* 12-11 */
1134 col_str[12][8] = 0.778; /* 13-11 */
1135 col_str[13][8] = 0.941; /* 14-11 */
1136
1137 col_str[10][9] = 0.377; /* 8-7 */
1138 col_str[11][9] = 0.603; /* 12-7 */
1139 col_str[12][9] = 0.472; /* 13-7 */
1140 col_str[13][9] = 0.302; /* 14-7 */
1141
1142 col_str[11][10] = 0.216; /* 12-8 */
1143 col_str[12][10] = 0.137; /* 13-8 */
1144 col_str[13][10] = 0.106; /* 14-8 */
1145
1146 col_str[12][11] = 1.25;
1147 col_str[13][11] = 0.292;
1148
1149 col_str[13][12] = 1.10;
1150 }
1151
1152 ASSERT( ipHi > ipLo );
1153 double CollisionStrength = col_str[ipHi][ipLo];
1154 ASSERT( CollisionStrength >0. );
1155
1156 return( CollisionStrength );
1157}
1158
1159double Fe4_cs(long ipLo, long ipHi)
1160{
1161 const int NLFE4=12;
1162
1163 /* collision strengths from Berrington and Pelan Ast Ap S 114, 367.
1164 * order is cs(low,up) */
1165 static const double Fe4CS[NLFE4][NLFE4] = {
1166 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1167 {0.98,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1168 {0.8167,3.72,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1169 {0.49,0.0475,0.330,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1170 {0.6533,0.473,2.26,1.64,0.,0.,0.,0.,0.,0.,0.,0.},
1171 {0.45,0.686,0.446,0.106,0.254,0.,0.,0.,0.,0.,0.,0.},
1172 {0.30,0.392,0.152,0.269,0.199,0.605,0.,0.,0.,0.,0.,0.},
1173 {0.15,0.0207,0.190,0.0857,0.166,0.195,0.327,0.,0.,0.,0.,0.},
1174 {0.512,1.23,0.733,0.174,0.398,0.623,0.335,0.102,0.,0.,0.,0.},
1175 {0.128,0.0583,0.185,0.200,0.188,0.0835,0.127,0.0498,0.0787,0.,0.,0.},
1176 {0.384,0.578,0.534,0.363,0.417,0.396,0.210,0.171,0.810,0.101,0.,0.},
1177 {0.256,0.234,0.306,0.318,0.403,0.209,0.195,0.112,0.195,0.458,0.727,0.}
1178 };
1179
1180 ASSERT( ipHi > ipLo );
1181 double CollisionStrength = Fe4CS[ipHi][ipLo];
1182 ASSERT( CollisionStrength >0. );
1183
1184 return( CollisionStrength );
1185}
1186
1187double Fe5_cs(long ipLo, long ipHi)
1188{
1189 const int NLFE5 = 14;
1190 static double col_str[NLFE5][NLFE5];
1191 /*>>refer Fe5 cs Shields ApJ 219, 559. */
1192
1193
1194 static double lgOneTimeMustInit=true;
1195 if( lgOneTimeMustInit )
1196 {
1197 lgOneTimeMustInit = false;
1198 for( int i = 0;i < NLFE5;i++)
1199 {
1200 for( int j = 0;j < NLFE5;j++)
1201 {
1202 col_str[i][j] = 1.;
1203 }
1204 }
1205
1206 col_str[10][3] = 1.4; //3896A
1207 col_str[7][2] = 1.1; //4072A
1208 col_str[13][4] = 3.7; //3892A
1209 col_str[12][3] = 3.7; //3839A
1210 col_str[11][2] = 2.0; //3795A
1211 }
1212
1213 ASSERT( ipHi > ipLo );
1214 double CollisionStrength = col_str[ipHi][ipLo];
1215 ASSERT( CollisionStrength >0. );
1216
1217 return( CollisionStrength );
1218}
long ipTFe34
long ipTFe26
long ipFe1_24m
long ipFe1_35m
long ipT1122
long ipT191
long ipTuv3
long ipTFe46
long ipTFe56
long ipFeI2966
long ipFeI3021
long ipFeI3729
long ipTFe35
long ipFeI3884
long ipFeI3457
long ipTr48
long ipTFe16
void FeIIIntenZero(void)
void FeII_RT_Make(void)
void FeIILevelPops(void)
void atom_level2(const TransitionProxy &t)
void atom_level3(const TransitionProxy &t10, const TransitionProxy &t21, const TransitionProxy &t20)
void atom_levelN(long int nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double ***AulEscp, double ***col_str, double ***AulDest, double ***AulPump, double ***CollRate, const double source[], const double sink[], bool lgCollRateDone, double *cooltl, double *coolder, const char *chLabel, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE, multi_arr< double, 2 > *Cool, multi_arr< double, 2 > *dCooldT)
double atom_pop2(double omega, double g1, double g2, double a21, double bltz, double abund)
Definition atom_pop2.cpp:9
t_FeII FeII
Definition atomfeii.cpp:5
t_atoms atoms
Definition atoms.cpp:5
const long LIMLEVELN
Definition atoms.h:237
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int ipIRON
Definition cddefines.h:330
#define MALLOC(exp)
Definition cddefines.h:501
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
t_conv conv
Definition conv.cpp:5
void CoolAdd(const char *chLabel, realnum lambda, double cool)
Definition cool_etc.cpp:13
STATIC void Fe3Lev14(void)
double Fe3_cs(long ipLo, long ipHi)
double Fe5_cs(long ipLo, long ipHi)
void CoolIron(void)
double Fe4_cs(long ipLo, long ipHi)
STATIC void Fe2_cooling(void)
Definition cool_iron.cpp:36
STATIC void Fe4Lev12(void)
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
t_dense dense
Definition dense.cpp:24
t_fe fe
Definition fe.cpp:5
#define NLFE3
Definition fe.h:10
t_hydro hydro
Definition hydrogenic.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
double RefIndex(double EnergyWN)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double T1CM
Definition physconst.h:167
static bool lgFirst
static double ** AulEscp
Definition species2.cpp:29
static double * ex
Definition species2.cpp:28
static double * pops
Definition species2.cpp:28
static double ** col_str
Definition species2.cpp:29
static double * depart
Definition species2.cpp:28
static double ** AulPump
Definition species2.cpp:29
static double ** AulDest
Definition species2.cpp:29
static double ** CollRate
Definition species2.cpp:29
TransitionProxy::iterator TauDummy
Definition taulines.cpp:60
TransitionList TauLines("TauLines", &AnonStates)
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
void LineConvRate2CS(const TransitionProxy &t, realnum rate)
void MakeCS(const TransitionProxy &t)
void PutCS(double cs, const TransitionProxy &t)