cloudy  trunk
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 */
30 STATIC void Fe3Lev14(void);
31 
32 /*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
33 STATIC void Fe4Lev12(void);
34 
35 /*Fe2_cooling compute cooling due to FeII emission */
36 STATIC 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 */
85  FeIIIntenZero();
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();
181  FeIILevelPops();
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]);
473  atom_level2(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 */
494 void 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]);
511  LineConvRate2CS( TauLines[ipFe1_24m] , rate );
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]);
517  LineConvRate2CS( TauLines[ipFe1_35m] , rate );
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 
528  atom_level3(TauLines[ipFe1_24m],TauLines[ipFe1_35m],*TauDummy);
529 
530  /* series of FeI lines from Dima Verner's list, each 2-lev atom
531  *
532  * Fe I 3884 */
533  MakeCS(TauLines[ipFeI3884]);
534  atom_level2(TauLines[ipFeI3884]);
535 
536  /* Fe I 3729 */
537  MakeCS(TauLines[ipFeI3729]);
538  atom_level2(TauLines[ipFeI3729]);
539 
540  /* Fe I 3457 */
541  MakeCS(TauLines[ipFeI3457]);
542  atom_level2(TauLines[ipFeI3457]);
543 
544  /* Fe I 3021 */
545  MakeCS(TauLines[ipFeI3021]);
546  atom_level2(TauLines[ipFeI3021]);
547 
548  /* Fe I 2966 */
549  MakeCS(TauLines[ipFeI2966]);
550  atom_level2(TauLines[ipFeI2966]);
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]);
568  atom_level2(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 */
585 STATIC 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();
761  cdEXIT(EXIT_FAILURE);
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) */
793 STATIC 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 */
957  ASSERT( NLFE3 <= LIMLEVELN);
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 */
967  atom_levelN(
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();
1016  cdEXIT(EXIT_FAILURE);
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 
1036 double 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 
1159 double 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 
1187 double 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 }
thermal.h
ipTFe16
long ipTFe16
Definition: atmdat_readin.cpp:87
atom_pop2
double atom_pop2(double omega, double g1, double g2, double a21, double bltz, double abund)
Definition: atom_pop2.cpp:9
t_dense::eden
double eden
Definition: dense.h:190
ipTFe34
long ipTFe34
Definition: atmdat_readin.cpp:87
ipTr48
long ipTr48
Definition: atmdat_readin.cpp:87
ipTuv3
long ipTuv3
Definition: atmdat_readin.cpp:86
dense
t_dense dense
Definition: dense.cpp:24
atoms.h
MakeCS
void MakeCS(const TransitionProxy &t)
Definition: transition.cpp:613
t_FeII::lgPrint
bool lgPrint
Definition: atomfeii.h:204
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
ipTFe46
long ipTFe46
Definition: atmdat_readin.cpp:88
t_hydro::dstfe2lya
realnum dstfe2lya
Definition: hydrogenic.h:60
realnum
float realnum
Definition: cddefines.h:103
conv.h
STATIC
#define STATIC
Definition: cddefines.h:97
t_conv::lgFirstSweepThisZone
bool lgFirstSweepThisZone
Definition: conv.h:155
t_fe::fe42567
double fe42567
Definition: fe.h:20
atom_level3
void atom_level3(const TransitionProxy &t10, const TransitionProxy &t21, const TransitionProxy &t20)
Definition: atom_level3.cpp:15
ipT1122
long ipT1122
Definition: atmdat_readin.cpp:88
ipFeI3884
long ipFeI3884
Definition: atmdat_readin.cpp:85
ipIRON
const int ipIRON
Definition: cddefines.h:330
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
Fe3_cs
double Fe3_cs(long ipLo, long ipHi)
Definition: cool_iron.cpp:1036
ex
static double * ex
Definition: species2.cpp:28
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
T1CM
const UNUSED double T1CM
Definition: physconst.h:167
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
lines_service.h
CoolAdd
void CoolAdd(const char *chLabel, realnum lambda, double cool)
Definition: cool_etc.cpp:13
ipFeI3729
long ipFeI3729
Definition: atmdat_readin.cpp:85
t_FeII::lgFeIILargeOn
bool lgFeIILargeOn
Definition: atomfeii.h:198
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
Fe4_cs
double Fe4_cs(long ipLo, long ipHi)
Definition: cool_iron.cpp:1159
iso.h
atom_levelN
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)
Definition: atom_leveln.cpp:15
t_FeII::ddT_Fe2_UVsimp_cool
double ddT_Fe2_UVsimp_cool
Definition: atomfeii.h:257
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
t_thermal::ctot
double ctot
Definition: thermal.h:112
fe
t_fe fe
Definition: fe.cpp:5
t_FeII::for7
double for7
Definition: atomfeii.h:259
Fe4Lev12
STATIC void Fe4Lev12(void)
Definition: cool_iron.cpp:585
Fe3Lev14
STATIC void Fe3Lev14(void)
Definition: cool_iron.cpp:793
nzone
long int nzone
Definition: cddefines.cpp:14
LIMLEVELN
const long LIMLEVELN
Definition: atoms.h:237
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
AulPump
static double ** AulPump
Definition: species2.cpp:29
col_str
static double ** col_str
Definition: species2.cpp:29
t_FeII::lgSlow
bool lgSlow
Definition: atomfeii.h:201
dense.h
t_fe::fe40401
double fe40401
Definition: fe.h:21
FeIIIntenZero
void FeIIIntenZero(void)
Definition: atom_feii.cpp:1851
coolheavy.h
cooling.h
t_FeII::Fe2_UVsimp_cool
double Fe2_UVsimp_cool
Definition: atomfeii.h:258
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
fe.h
CollRate
static double ** CollRate
Definition: species2.cpp:29
atoms
t_atoms atoms
Definition: atoms.cpp:5
TauDummy
TransitionProxy::iterator TauDummy
Definition: taulines.cpp:60
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_trace::nTrConvg
int nTrConvg
Definition: trace.h:27
t_fe::fe41007
double fe41007
Definition: fe.h:25
CoolIron
void CoolIron(void)
Definition: cool_iron.cpp:494
t_FeII::Fe2_large_heat
double Fe2_large_heat
Definition: atomfeii.h:252
t_thermal::heating
double heating[LIMELM][LIMELM]
Definition: thermal.h:158
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
ipTFe56
long ipTFe56
Definition: atmdat_readin.cpp:88
FeII
t_FeII FeII
Definition: atomfeii.cpp:5
lgFirst
static bool lgFirst
Definition: prt_linesum.cpp:16
ipFeI3457
long ipFeI3457
Definition: atmdat_readin.cpp:85
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:141
PutCS
void PutCS(double cs, const TransitionProxy &t)
Definition: transition.cpp:317
t_fe::Fe3_wl
double ** Fe3_wl
Definition: fe.h:31
MAX2
#define MAX2
Definition: cddefines.h:782
Fe2_cooling
STATIC void Fe2_cooling(void)
Definition: cool_iron.cpp:36
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
ipFe1_24m
long ipFe1_24m
Definition: atmdat_readin.cpp:96
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
iteration
long int iteration
Definition: cddefines.cpp:16
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
t_atoms::DepLTELevels
double DepLTELevels[LIMLEVELN+1]
Definition: atoms.h:271
TauLines
TransitionList TauLines("TauLines", &AnonStates)
ipTFe26
long ipTFe26
Definition: atmdat_readin.cpp:87
hydrogenic.h
t_fe::fe42829
double fe42829
Definition: fe.h:19
t_trace::lgCoolTr
bool lgCoolTr
Definition: trace.h:112
t_fe::Fe4CoolTot
double Fe4CoolTot
Definition: fe.h:17
ipH2p
const int ipH2p
Definition: iso.h:29
t_fe::Fe3CoolTot
double Fe3CoolTot
Definition: fe.h:30
ipTFe35
long ipTFe35
Definition: atmdat_readin.cpp:88
t_fe::fe41008
double fe41008
Definition: fe.h:26
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
AulEscp
static double ** AulEscp
Definition: species2.cpp:29
LineConvRate2CS
void LineConvRate2CS(const TransitionProxy &t, realnum rate)
Definition: transition.cpp:521
physconst.h
FeIILevelPops
void FeIILevelPops(void)
Definition: atom_feii.cpp:718
t_FeII::Fe2_large_cool
double Fe2_large_cool
Definition: atomfeii.h:250
ligbar.h
t_fe::fe41106
double fe41106
Definition: fe.h:24
conv
t_conv conv
Definition: conv.cpp:5
atom_level2
void atom_level2(const TransitionProxy &t)
Definition: atom_level2.cpp:17
t_fe::fe41206
double fe41206
Definition: fe.h:23
t_CoolHeavy::c3892
double c3892
Definition: coolheavy.h:49
t_fe::fe41207
double fe41207
Definition: fe.h:22
t_fe::fe40906
double fe40906
Definition: fe.h:27
ipFeI3021
long ipFeI3021
Definition: atmdat_readin.cpp:86
t_thermal::dCooldT
double dCooldT
Definition: thermal.h:119
ipT191
long ipT191
Definition: atmdat_readin.cpp:89
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
pops
static double * pops
Definition: species2.cpp:28
CoolHeavy
t_CoolHeavy CoolHeavy
Definition: coolheavy.cpp:5
phycon.h
ShowMe
void ShowMe(void)
Definition: service.cpp:181
depart
static double * depart
Definition: species2.cpp:28
Fe5_cs
double Fe5_cs(long ipLo, long ipHi)
Definition: cool_iron.cpp:1187
ipH1s
const int ipH1s
Definition: iso.h:27
t_conv::lgLastSweepThisZone
bool lgLastSweepThisZone
Definition: conv.h:157
t_atoms::PopLevels
double PopLevels[LIMLEVELN+1]
Definition: atoms.h:270
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
NLFE3
#define NLFE3
Definition: fe.h:10
t_conv::lgSearch
bool lgSearch
Definition: conv.h:175
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
BOLTZMANN
const UNUSED double BOLTZMANN
Definition: physconst.h:97
t_phycon::te
double te
Definition: phycon.h:11
ipFeI2966
long ipFeI2966
Definition: atmdat_readin.cpp:86
FeII_RT_Make
void FeII_RT_Make(void)
Definition: atom_feii.cpp:1538
AulDest
static double ** AulDest
Definition: species2.cpp:29
t_FeII::ddT_Fe2_large_cool
double ddT_Fe2_large_cool
Definition: atomfeii.h:251
ipFe1_35m
long ipFe1_35m
Definition: atmdat_readin.cpp:96
t_fe::Fe3_emiss
double ** Fe3_emiss
Definition: fe.h:31
atomfeii.h
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
t_fe::fe42836
double fe42836
Definition: fe.h:18
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17