cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cool_iron.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 /*Fe_10_11_13_cs evaluate collision strength for Fe */
5 /*fe14cs compute collision strengths for forbidden transitions */
6 /*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */
7 /*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
8 /*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
9 /*Fe2_cooling compute cooling due to FeII emission */
10 /*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */
11 /*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */
12 #include "cddefines.h"
13 #include "physconst.h"
14 #include "dense.h"
15 #include "coolheavy.h"
16 #include "taulines.h"
17 #include "phycon.h"
18 #include "iso.h"
19 #include "conv.h"
20 #include "trace.h"
21 #include "hydrogenic.h"
22 #include "ligbar.h"
23 #include "cooling.h"
24 #include "thermal.h"
25 #include "lines_service.h"
26 #include "atoms.h"
27 #include "atomfeii.h"
28 #include "fe.h"
29 #define NLFE2 6
30 
31 /*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */
32 STATIC void Fe11Lev5(void);
33 
34 /*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */
35 STATIC void Fe13Lev5(void);
36 
37 /*fe14cs compute collision strengths for forbidden transitions */
38 STATIC void fe14cs(double te1,
39  double *csfe14);
40 
41 /*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
42 STATIC void Fe7Lev8(void);
43 
44 /*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */
45 STATIC void Fe3Lev14(void);
46 
47 /*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
48 STATIC void Fe4Lev12(void);
49 
50 /*Fe_10_11_13_cs evaluate collision strength for Fe */
52  /* ion, one of 10, 11, or 13 on physics scale */
53  int ion,
54  /* initial and final index of levels, with lowest energy 1, highest of 5
55  * on physics scale */
56  int init,
57  int final )
58 {
59 # define N 10
60  static double Fe10cs[6][6][2];
61  static double Fe11cs[6][6][2];
62  static double Fe13cs[6][6][2];
63  int i, j;
64  double cs;
65  int index = 0;
66  double temp_max, temp_min = 4;
67  double temp_log = phycon.alogte;
68  static bool lgFirstTime = true;
69 
70  DEBUG_ENTRY( "Fe_10_11_13_cs()" );
71 
72  if( lgFirstTime )
73  {
74  /* fit coefficients for collision strengths */
75  double aFe10[N] = {10.859,-1.1541,11.593,22.333,-0.4283,7.5663,3.087,1.0937,0.8261,59.678};
76  double bFe10[N] = {-1.4804,0.4956,-2.1096,-4.1288,0.1929,-1.3525,-0.5531,-0.1748,-0.1286,-11.081};
77  double aFe11[N] = {5.7269,1.2885,4.0877,0.4571,1.2911,2.2339,0.3621,0.7972,0.2225,1.1021};
78  double bFe11[N] = {-0.7559,-0.1671,-0.5678,-0.0653,-0.1589,-0.2924,-0.0506,-0.1038,-0.0302,-0.1062};
79  double aFe13[N] = {2.9102,1.8682,-0.353,0.0622,14.229,-4.3845,0.0375,-6.9222,0.688,-0.0609};
80  double bFe13[N] = {-0.4158,-0.242,0.1417,0.0023,-2.0643,1.2573,0.0286,2.0919,-0.083,0.1487};
81  /* do one-time initialization
82  * assigning array initially to zero */
83  for(i=0; i<6; i++)
84  {
85  for(j=0; j<6; j++)
86  {
87  set_NaN( Fe10cs[i][j], 2L );
88  set_NaN( Fe11cs[i][j], 2L );
89  set_NaN( Fe13cs[i][j], 2L );
90  }
91  }
92 
93  /* reading coefficients into 3D array */
94  for(i=1; i<6; i++)
95  {
96  for(j=i+1; j<6; j++)
97  {
98  Fe10cs[i][j][0] = aFe10[index];
99  Fe10cs[i][j][1] = bFe10[index];
100  Fe11cs[i][j][0] = aFe11[index];
101  Fe11cs[i][j][1] = bFe11[index];
102  Fe13cs[i][j][0] = aFe13[index];
103  Fe13cs[i][j][1] = bFe13[index];
104  index++;
105  }
106  }
107  lgFirstTime = false;
108  }
109 
110  /* Invalid entries returns '-1':the initial indices are smaller than
111  * the final indices */
112  if(init >= final)
113  {
114  cs = -1;
115  }
116  /* Invalid returns '-1': the indices are greater than 5 or smaller than 0 */
117  else if(init < 1 || init > 4 || final < 2 || final > 5)
118  {
119  cs = -1;
120  }
121  else
122  {
123  /* cs = a + b*log(T)
124  * if temp is out of range, return boundary values */
125  if(ion == 10)
126  {
127  temp_max = 5;
128  temp_log = MAX2(temp_log, temp_min);
129  temp_log = MIN2(temp_log, temp_max);
130  cs = Fe10cs[init][final][0] + Fe10cs[init][final][1]*temp_log;
131  }
132  else if(ion == 11)
133  {
134  temp_max = 6.7;
135  temp_log = MAX2(temp_log, temp_min);
136  temp_log = MIN2(temp_log, temp_max);
137  cs = Fe11cs[init][final][0] + Fe11cs[init][final][1]*temp_log;
138  }
139  else if(ion ==13)
140  {
141  temp_max = 5;
142  temp_log = MAX2(temp_log, temp_min);
143  temp_log = MIN2(temp_log, temp_max);
144  cs = Fe13cs[init][final][0] + Fe13cs[init][final][1]*temp_log;
145  }
146  else
147  /* this can't happen */
148  TotalInsanity();
149  }
150 
151  return cs;
152 
153 # undef N
154 }
155 
156 /*Fe2_cooling compute cooling due to FeII emission */
157 STATIC void Fe2_cooling( void )
158 {
159  long int i , j;
160  int nNegPop;
161 
162  static double **AulPump,
163  **CollRate,
164  **AulEscp,
165  **col_str ,
166  **AulDest,
167  *depart,
168  *pops,
169  *destroy,
170  *create;
171 
172  static bool lgFirst=true;
173  bool lgZeroPop;
174 
175  /* stat weights for Fred's 6 level model FeII atom */
176  static double gFe2[NLFE2]={1.,1.,1.,1.,1.,1.};
177  /* excitation energies (Kelvin) for Fred's 6 level model FeII atom */
178  static double ex[NLFE2]={0.,3.32e4,5.68e4,6.95e4,1.15e5,1.31e5};
179 
180  /* these are used to only evaluated FeII one time per temperature, zone,
181  * and abundance */
182  static double TUsed = 0.;
183  static realnum AbunUsed = 0.;
184  /* remember which zone last evaluated with */
185  static long int nZUsed=-1,
186  /* make sure at least two calls per zone */
187  nCall=0;
188 
189  DEBUG_ENTRY( "Fe2_cooling()" );
190 
191  /* return if nothing to do */
192  if( dense.xIonDense[ipIRON][1] == 0. )
193  {
194  /* zero abundance, do nothing */
195  /* heating cooling and derivative from large atom */
196  FeII.Fe2_large_cool = 0.;
197  FeII.Fe2_large_heat = 0.;
199 
200  /* cooling and derivative from simple UV atom */
201  FeII.Fe2_UVsimp_cool = 0.;
203 
204  /* now zero out intensities of all FeII lines and level populations */
205  FeIIIntenZero();
206  return;
207  }
208 
209  /* this can be the large 371 level FeII atom
210  * >>chng 05 dec 04, always call this but with default of 16 levels only
211  * >>chng 00 jan 06, total rewrite mostly done
212  * >>chng 97 jan 17, evaluate large FeII atom cooling every time temp changes
213  * >>chng 97 jan 31, added check on zone since must reevaluate even const temp
214  * >>chng 99 may 21, reevaluate when zone or temperature changes, but not when
215  * abundance changes, since we can simply rescale cooling */
216 
217  /* totally reevaluate large atom if new zone, or cooling is significant
218  * and temperature changed, we are in search phase,
219  * lgSlow option set true with atom FeII slow, forces constant
220  * evaluation of atom */
221  if( FeII.lgSlow ||
222  nzone < 1 ||
223  nzone != nZUsed ||
224  /* on new call, nCall is now set at previous zone's number of calls.
225  * it is set to zero below, so on second call, nCall is 0. On
226  * third call nCall is 1. Check for <1 is to insure at least two calls */
227  nCall < 1 ||
228  /* check whether things have changed on later calls */
229  ( ! fp_equal( phycon.te, TUsed ) && fabs(FeII.Fe2_large_cool/thermal.ctot)> 0.002 &&
230  fabs(dense.xIonDense[ipIRON][1]-AbunUsed)/SDIV(AbunUsed)> 0.002 ) ||
231  ( ! fp_equal( phycon.te, TUsed ) && fabs(FeII.Fe2_large_cool/thermal.ctot)> 0.01) )
232  {
233 
234  if( nZUsed == nzone )
235  {
236  /* not first call, increment, check above to make sure at least
237  * two evaluations */
238  ++nCall;
239  }
240  else
241  {
242  /* first call this zone set nCall to zero*/
243  nCall = 0;
244  }
245 
246  /* option to trace convergence and FeII calls */
247  if( trace.nTrConvg>=5 )
248  {
249  fprintf( ioQQQ, " CoolIron5 calling FeIILevelPops since ");
250  if( ! fp_equal( phycon.te, TUsed ) )
251  {
252  fprintf( ioQQQ,
253  "temperature changed, old new are %g %g, nCall %li ",
254  TUsed, phycon.te , nCall);
255  }
256  else if( nzone != nZUsed )
257  {
258  fprintf( ioQQQ,
259  "new zone, nCall %li ", nCall );
260  }
261  else if( FeII.lgSlow )
262  {
263  fprintf( ioQQQ,
264  "FeII.lgSlow set %li", nCall );
265  }
266  else if( conv.lgSearch )
267  {
268  fprintf( ioQQQ,
269  " in search phase %li", nCall );
270  }
271  else if( nCall < 2 )
272  {
273  fprintf( ioQQQ,
274  "not second nCall %li " , nCall );
275  }
276  else if( ! fp_equal( phycon.te, TUsed ) && FeII.Fe2_large_cool/thermal.ctot > 0.001 )
277  {
278  fprintf( ioQQQ,
279  "temp or cooling changed, new are %g %g nCall %li ",
280  phycon.te, FeII.Fe2_large_cool, nCall );
281  }
282  else
283  {
284  fprintf(ioQQQ, "????");
285  }
286  fprintf(ioQQQ, "\n");
287  }
288 
289  /* remember parameters for current conditions */
290  TUsed = phycon.te;
291  AbunUsed = dense.xIonDense[ipIRON][1];
292  nZUsed = nzone;
293 
294  /* this print turned on with atom FeII print command */
295  if( FeII.lgPrint )
296  {
297  fprintf(ioQQQ,
298  " FeIILevelPops called zone %4li te %5f abun %10e c(fe/tot):%6f nCall %li\n",
299  nzone,phycon.te,AbunUsed,FeII.Fe2_large_cool/thermal.ctot,nCall);
300  }
301 
302  /* this solves the multi-level problem,
303  * sets FeII.Fe2_large_cool,
304  FeII.Fe2_large_heat, &
305  FeII.ddT_Fe2_large_cool
306  but does nothing with them */
307  FeIILevelPops();
308  {
309  /*@-redef@*/
310  enum{DEBUG_LOC=false};
311  /*@+redef@*/
312  if( DEBUG_LOC && iteration>1 && nzone>=4 )
313  {
314  fprintf(ioQQQ,"DEBUG1\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
315  nzone,
316  phycon.te,
318  dense.eden,
322  thermal.ctot );
323  }
324  }
325 
326  if( trace.nTrConvg>=5 || FeII.lgPrint)
327  {
328  /* spacing needed to get proper trace convergence output */
329  fprintf( ioQQQ, " FeIILevelPops5 returned cool=%.2e heat=%.2e derivative=%.2e\n ",
331  }
332 
333  }
334  else if( ! fp_equal( dense.xIonDense[ipIRON][1], AbunUsed ) )
335  {
336  realnum ratio;
337  /* this branch, same zone and temperature, but small change in abundance, so just
338  * rescale cooling and derivative by this change. assumption is that very small changes
339  * in abundance occurs as ots rates damp out */
340  if( trace.nTrConvg>=5 )
341  {
342  fprintf( ioQQQ,
343  " CoolIron rescaling FeIILevelPops since small change, CFe2=%.2e CTOT=%.2e\n",
345  }
346  ratio = dense.xIonDense[ipIRON][1]/AbunUsed;
347  FeII.Fe2_large_cool *= ratio;
348  FeII.ddT_Fe2_large_cool *= ratio;
349  FeII.Fe2_large_heat *= ratio;
350  AbunUsed = dense.xIonDense[ipIRON][1];
351  }
352  else
353  {
354  /* this is case where temp is unchanged, so heating and cooling same too */
355  if( trace.nTrConvg>=5 )
356  {
357  fprintf( ioQQQ, " CoolIron NOT calling FeIILevelPops\n");
358  }
359  }
360 
361  /* evaluate some strong lines that would have been evaluated by the 16 level atom */
362  FeIIFillLow16();
363 
364  /* now update total cooling and its derivative
365  * all paths flow through here */
366  /* FeII.Fecool = FeII.Fe2_large_cool; */
367  CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_large_cool));
368 
369  /* add negative cooling to heating stack */
370  thermal.heating[25][27] = MAX2(0.,FeII.Fe2_large_heat);
371 
372  /* counts as heating derivative if negative cooling */
373  if( FeII.Fe2_large_cool > 0. )
374  {
375  /* >>chng 01 mar 16, add factor of 3 due to conv problems after changing damper */
377  }
378 
379  if( trace.lgTrace && trace.lgCoolTr )
380  {
381  fprintf( ioQQQ, " Large FeII returns te, cooling, dc=%11.3e%11.3e%11.3e\n",
383  }
384 
385  /* >>chng 05 nov 29, still do simple UV atom if only ground term is done */
386  if( !FeII.lgFeIILargeOn )
387  {
388 
389  /* following treatment of Fe II follows
390  * >>refer fe2 model Wills, B.J., Wills, D., Netzer, H. 1985, ApJ, 288, 143
391  * all elements are used, and must be set to zero if zero */
392 
393  /* set up space for simple model of UV FeII emission */
394  if( lgFirst )
395  {
396  /* will never do this again in this core load */
397  lgFirst = false;
398  /* allocate the 1D arrays*/
399  pops = (double *)MALLOC( sizeof(double)*(NLFE2) );
400  create = (double *)MALLOC( sizeof(double)*(NLFE2) );
401  destroy = (double *)MALLOC( sizeof(double)*(NLFE2) );
402  depart = (double *)MALLOC( sizeof(double)*(NLFE2) );
403  /* create space for the 2D arrays */
404  AulPump = ((double **)MALLOC((NLFE2)*sizeof(double *)));
405  CollRate = ((double **)MALLOC((NLFE2)*sizeof(double *)));
406  AulDest = ((double **)MALLOC((NLFE2)*sizeof(double *)));
407  AulEscp = ((double **)MALLOC((NLFE2)*sizeof(double *)));
408  col_str = ((double **)MALLOC((NLFE2)*sizeof(double *)));
409 
410  for( i=0; i<(NLFE2); ++i )
411  {
412  AulPump[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
413  CollRate[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
414  AulDest[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
415  AulEscp[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
416  col_str[i] = ((double *)MALLOC((NLFE2)*sizeof(double )));
417  }
418  }
419 
420  /*zero out all arrays, then check that upper diagonal remains zero below */
421  for( i=0; i < NLFE2; i++ )
422  {
423  create[i] = 0.;
424  destroy[i] = 0.;
425  for( j=0; j < NLFE2; j++ )
426  {
427  /*data[j][i] = 0.;*/
428  col_str[j][i] = 0.;
429  AulEscp[j][i] = 0.;
430  AulDest[j][i] = 0.;
431  AulPump[j][i] = 0.;
432  }
433  }
434 
435  /* now put in real data for lines */
436  AulEscp[1][0] = 1.;
437  AulEscp[2][0] = ( TauLines[ipTuv3].Emis->Pesc + TauLines[ipTuv3].Emis->Pelec_esc)*TauLines[ipTuv3].Emis->Aul;
438  AulDest[2][0] = TauLines[ipTuv3].Emis->Pdest*TauLines[ipTuv3].Emis->Aul;
439  AulPump[0][2] = TauLines[ipTuv3].Emis->pump;
440 
442  AulDest[5][0] = TauLines[ipTFe16].Emis->Pdest*TauLines[ipTFe16].Emis->Aul;
443  /* continuum pumping of n=6 */
444  AulPump[0][5] = TauLines[ipTFe16].Emis->pump;
445  /* Ly-alpha pumping */
446 
447  double PumpLyaFeII = StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*dense.xIonDense[ipHYDROGEN][1]*
448  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul*
450 
451  if( iteration == 1 )
452  PumpLyaFeII = 0.;
453 
454  AulPump[0][5] += PumpLyaFeII;
455 
456  AulEscp[2][1] = (TauLines[ipTr48].Emis->Pesc + TauLines[ipTr48].Emis->Pelec_esc)*TauLines[ipTr48].Emis->Aul;
457  AulDest[2][1] = TauLines[ipTr48].Emis->Pdest*TauLines[ipTr48].Emis->Aul;
458  AulPump[1][2] = TauLines[ipTr48].Emis->pump;
459 
461  AulDest[5][1] = TauLines[ipTFe26].Emis->Pdest*TauLines[ipTFe26].Emis->Aul;
462  AulPump[1][5] = TauLines[ipTFe26].Emis->pump;
463 
464  AulEscp[3][2] = (TauLines[ipTFe34].Emis->Pesc + TauLines[ipTFe34].Emis->Pelec_esc)*TauLines[ipTFe34].Emis->Aul;
465  AulDest[3][2] = TauLines[ipTFe34].Emis->Pdest*TauLines[ipTFe34].Emis->Aul;
466  AulPump[2][3] = TauLines[ipTFe34].Emis->pump;
467 
469  AulDest[4][2] = TauLines[ipTFe35].Emis->Pdest*TauLines[ipTFe35].Emis->Aul;
470  AulPump[2][4] = TauLines[ipTFe35].Emis->pump;
471 
472  AulEscp[5][3] = (TauLines[ipTFe46].Emis->Pesc + TauLines[ipTFe46].Emis->Pelec_esc)*TauLines[ipTFe46].Emis->Aul;
473  AulDest[5][3] = TauLines[ipTFe46].Emis->Pdest*TauLines[ipTFe46].Emis->Aul;
474  AulPump[3][5] = TauLines[ipTFe46].Emis->pump;
475 
477  AulDest[5][4] = TauLines[ipTFe56].Emis->Pdest*TauLines[ipTFe56].Emis->Aul;
478  AulPump[4][5] = TauLines[ipTFe56].Emis->pump;
479 
480  /* these are collision strengths */
481  col_str[1][0] = 1.;
482  col_str[2][0] = 12.;
483  col_str[3][0] = 1.;
484  col_str[4][0] = 1.;
485  col_str[5][0] = 12.;
486  col_str[2][1] = 6.;
487  col_str[3][1] = 1.;
488  col_str[4][1] = 1.;
489  col_str[5][1] = 12.;
490  col_str[3][2] = 6.;
491  col_str[4][2] = 12.;
492  col_str[5][2] = 1.;
493  col_str[4][3] = 1.;
494  col_str[5][3] = 12.;
495  col_str[5][4] = 6.;
496 
497  /*void atom_levelN(long,long,realnum,double[],double[],double[],double*,
498  double*,double*,long*,realnum*,realnum*,STRING,int*);*/
499  atom_levelN(NLFE2,
500  dense.xIonDense[ipIRON][1],
501  gFe2,
502  ex,
503  'K',
504  pops,
505  depart,
506  &AulEscp ,
507  &col_str,
508  &AulDest,
509  &AulPump,
510  &CollRate,
511  create,
512  destroy,
513  false,/* say atom_levelN should evaluate coll rates from cs */
514  /*&&ipdest,*/
517  "FeII",
518  &nNegPop,
519  &lgZeroPop,
520  false );
521 
522  /* nNegPop positive if negative pops occurred, negative if too cold */
523  if( nNegPop>0 /*negative if too cold - that is not negative and is OK */ )
524  {
525  fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population for simple UV FeII.\n");
526  }
527 
528  /* add heating - cooling by this process */;
529  CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_UVsimp_cool));
530  thermal.heating[25][27] = MAX2(0.,-FeII.Fe2_UVsimp_cool);
532 
533  /* LIMLEVELN is the dim of the PopLevels vector */
534  ASSERT( NLFE2 <= LIMLEVELN );
535  for( i=0; i<NLFE2; ++i)
536  {
537  atoms.PopLevels[i] = pops[i];
538  atoms.DepLTELevels[i] = depart[i];
539  }
540 
541  TauLines[ipTuv3].Lo->Pop = pops[0];
542  TauLines[ipTuv3].Hi->Pop = pops[2];
543  TauLines[ipTuv3].Emis->PopOpc = (pops[0] - pops[2]);
544  TauLines[ipTuv3].Emis->phots = pops[2]*AulEscp[2][0];
547 
548  TauLines[ipTr48].Lo->Pop = pops[1];
549  TauLines[ipTr48].Hi->Pop = pops[2];
550  TauLines[ipTr48].Emis->PopOpc = (pops[1] - pops[2]);
551  TauLines[ipTr48].Emis->phots = pops[2]*AulEscp[2][1];
554 
555  FeII.for7 = pops[1]*AulEscp[1][0]*4.65e-12;
556 
557  TauLines[ipTFe16].Lo->Pop = pops[0];
558  TauLines[ipTFe16].Hi->Pop = pops[5];
559  /* Lyman alpha optical depths are not known on first iteration,
560  * inward optical depths used, so line trapping overestimated,
561  * this can cause artificial maser in FeII - prevent by not
562  * including stimulated emission correction on first iteration */
563  if( iteration == 1 )
564  {
565  /* prevent maser due to pumping by Ly-alpha */
566  TauLines[ipTFe16].Emis->PopOpc = pops[0];
567  }
568  else
569  {
570  TauLines[ipTFe16].Emis->PopOpc = (pops[0] - pops[5]);
571  }
572  TauLines[ipTFe16].Emis->phots = pops[5]*AulEscp[5][0];
575 
576  TauLines[ipTFe26].Lo->Pop = pops[1];
577  TauLines[ipTFe26].Hi->Pop = pops[5];
578  if( iteration == 1 )
579  {
580  /* prevent maser due to pumping by Ly-alpha */
581  TauLines[ipTFe26].Emis->PopOpc = pops[1];
582  }
583  else
584  {
585  TauLines[ipTFe26].Emis->PopOpc = (pops[1] - pops[5]);
586  }
587  TauLines[ipTFe26].Emis->phots = pops[5]*AulEscp[5][1];
590 
591  TauLines[ipTFe34].Lo->Pop = pops[2];
592  TauLines[ipTFe34].Hi->Pop = pops[3];
593  if( iteration == 1 )
594  {
595  /* prevent maser due to pumping by Ly-alpha */
596  /* this line is indirectly pumped after decays from 6 to 4 */
597  TauLines[ipTFe34].Emis->PopOpc = pops[2];
598  }
599  else
600  {
601  TauLines[ipTFe34].Emis->PopOpc = (pops[2] - pops[3]);
602  }
603  TauLines[ipTFe34].Emis->phots = pops[3]*AulEscp[3][2];
606 
607  TauLines[ipTFe35].Lo->Pop = pops[2];
608  TauLines[ipTFe35].Hi->Pop = pops[4];
609  TauLines[ipTFe35].Emis->PopOpc = (pops[2] - pops[4]);
610  TauLines[ipTFe35].Emis->phots = pops[4]*AulEscp[4][2];
613 
614  TauLines[ipTFe46].Lo->Pop = pops[3];
615  TauLines[ipTFe46].Hi->Pop = pops[5];
616  if( iteration == 1 )
617  {
618  /* prevent maser due to pumping by Ly-alpha */
619  TauLines[ipTFe46].Emis->PopOpc = pops[3];
620  }
621  else
622  {
623  TauLines[ipTFe46].Emis->PopOpc = (pops[3] - pops[5]);
624  }
625  TauLines[ipTFe46].Emis->PopOpc = (pops[3] - pops[5]);
626  TauLines[ipTFe46].Emis->phots = pops[5]*AulEscp[5][3];
629 
630  TauLines[ipTFe56].Lo->Pop = pops[4];
631  TauLines[ipTFe56].Hi->Pop = pops[5];
632  if( iteration == 1 )
633  {
634  /* prevent maser due to pumping by Ly-alpha */
635  TauLines[ipTFe56].Emis->PopOpc = pops[4];
636  }
637  else
638  {
639  TauLines[ipTFe56].Emis->PopOpc = (pops[4] - pops[5]);
640  }
641  TauLines[ipTFe56].Emis->phots = pops[5]*AulEscp[5][4];
644 
645  /* Jack's funny FeII lines, data from
646  * >>refer fe2 energy Johansson, S., Brage, T., Leckrone, D.S., Nave, G. &
647  * >>refercon Wahlgren, G.M. 1995, ApJ 446, 361 */
648  PutCS(10.,&TauLines[ipT191]);
649  atom_level2(&TauLines[ipT191]);
650  }
651 
652  {
653  /*@-redef@*/
654  enum{DEBUG_LOC=false};
655  /*@+redef@*/
656  if( DEBUG_LOC && iteration>1 && nzone>=4 )
657  {
658  fprintf(ioQQQ,"DEBUG2\t%.2e\t%.2e\t%.2e\n",
659  phycon.te,
662  }
663  }
664 
665  return;
666 
667 }
668 
669 /*CoolIron - calculation total cooling due to Fe */
670 void CoolIron(void)
671 {
672  long int i;
673  double cs ,
674  cs12, cs13, cs23,
675  cs2s2p,
676  cs2s3p;
677  realnum p2,
678  rate;
679 
680  static bool lgFe22First=true;
681  static long int *ipFe22Pump=NULL,
682  nFe22Pump=0;
683  double Fe22_pump_rate;
684 
685  DEBUG_ENTRY( "CoolIron()" );
686 
687  /* cooling by FeI 24m, 34.2 m */
688  /* >>chng 03 nov 15, add these lines */
692  /*>>refer Fe1 cs Hollenbach & McKee 89 */
693  /* the 24.0 micron line */
694  rate = (realnum)(1.2e-7 * dense.eden +
695  /* >>chng 05 jul 05, eden to cdsqte */
696  /*8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
697  8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
698  LineConvRate2CS( &TauLines[ipFe1_24m] , rate );
699 
700  rate = (realnum)(9.3e-8 * dense.eden +
701  /* >>chng 05 jul 05, eden to cdsqte */
702  /*5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
703  5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
704  LineConvRate2CS( &TauLines[ipFe1_35m] , rate );
705 
706  rate = (realnum)(1.2e-7 * dense.eden +
707  /* >>chng 05 jul 05, eden to cdsqte */
708  /*6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
709  6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
711  LineConvRate2CS( &TauDummy , rate );
712  /* this says that line is a dummy, not real one */
713  TauDummy.Hi->g = 0.;
714 
715  atom_level3(&TauLines[ipFe1_24m],&TauLines[ipFe1_35m],&TauDummy);
716 
717  /* series of FeI lines from Dima Verner's list, each 2-lev atom
718  *
719  * Fe I 3884 */
721  atom_level2(&TauLines[ipFeI3884]);
722 
723  /* Fe I 3729 */
725  atom_level2(&TauLines[ipFeI3729]);
726 
727  /* Fe I 3457 */
729  atom_level2(&TauLines[ipFeI3457]);
730 
731  /* Fe I 3021 */
733  atom_level2(&TauLines[ipFeI3021]);
734 
735  /* Fe I 2966 */
737  atom_level2(&TauLines[ipFeI2966]);
738 
739  /* >>chng 05 dec 03, move Fe2 FeII Fe II cooling into separate routine */
740  Fe2_cooling();
741 
742  /* lump 3p and 3f together; cs=
743  * >>refer fe3 as Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384
744  * A from
745  * >>refer fe3 as Garstang, R.H., 1957, Vistas in Astronomy, 1, 268
746  * FE III 5270, is 20.9% of total
747  * >>chng 05 feb 18, Kevin Blagrave email
748  * average wavelength is 4823 with statistical weight averaging of upper energy level,
749  * as per , change 5th number from 2.948 to 2.984, also photon energy
750  * from 3.78 to 4.12 */
751 
752  /* >>chng 05 dec 16, FeIII update by Kevin Blagrave */
753  /*CoolHeavy.c5270 = atom_pop2(5.53,25.,30.,0.435,2.984e4,dense.xIonDense[ipIRON][2])*
754  4.12e-12;
755  CoolAdd("Fe 3",5270,CoolHeavy.c5270);*/
756  /* FeIII 1122 entire multiplet - atomic data=A's Dima, CS = guess */
757  PutCS(25.,&TauLines[ipT1122]);
758  atom_level2(&TauLines[ipT1122]);
759 
760  /* call 14 level atom */
761  Fe3Lev14();
762 
763  /* call 12 level atom */
764  Fe4Lev12();
765 
766  /* FE V 3892 + 3839, data from Shields */
767  CoolHeavy.c3892 = atom_pop2(7.4,25.,5.,0.6,3.7e4,dense.xIonDense[ipIRON][4])*
768  5.11e-12;
769  CoolAdd("Fe 5",3892,CoolHeavy.c3892);
770 
771  /* FE VI 5177+5146+4972+4967
772  * data from
773  * >>refer fe6 as Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384 */
774  CoolHeavy.c5177 = atom_pop2(1.9,28.,18.,0.52,2.78e4,dense.xIonDense[ipIRON][5])*
775  3.84e-12;
776  CoolAdd("Fe 6",5177,CoolHeavy.c5177);
777 
778  /* >>chng 04 nov 04, add multi-level fe7 */
779  Fe7Lev8();
780 
781  /* Fe IX 245,242
782  * all atomic data from
783  * >>refer fe9 all Flower, D.R. 1976, A&A, 56, 451 */
784  /* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */
785  PutCS(0.123,&TauLines[ipT245]);
786  /* AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41)
787  * C245 = AtomSeqBeryllium( .087 ,.038 ,.188 , t245 ,71. ) * 8.14E-11 */
788  AtomSeqBeryllium(.087,.038,.188,&TauLines[ipT245],71.);
789 
790  CoolHeavy.c242 = atoms.PopLevels[3]*8.22e-11*71.;
791 
792  /* Fe X Fe 10 Fe10 6374, A from
793  * >>referold fe10 as Mason, H. 1975, MNRAS 170, 651
794  * >>referold fe10 cs Mohan, M., Hibbert, A., & Kingston, A.E. 1994, ApJ, 434, 389
795  * >>referold fe10 as Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209
796  * does not agree with Mohan et al. by factor of 4
797  * 20x larger than Mason numbers used pre 1994 */
798  /*cs = 13.3-2.*MIN2(7.0,phycon.alogte); */
799  /*cs = MAX2(0.1,cs); */
805  /*cs = 1.;*/
806  /* >>chng 00 dec 05, update Fe10 collisions strength to Tayal 2000 */
807  /* >>refer fe10 cs Tayal, S.S., 2000, ApJ, 544, 575-580 */
808  /* >>chng 04 nov 10, 172.9 was mult rather than div by temp law,
809  * had no effect due to min on cs that lie below it */
810  /*cs = 172.9 /( phycon.te30 * phycon.te05 * phycon.te02 * phycon.te005 );*/
811  /* tabulated cs ends at 30,000K, do not allow cs to grow larger than largest
812  * tabulated value */
813  /*cs = MIN2( cs, 3.251 );*/
814  /* >>chng 05 dec 19, update As to
815  * >>refer fe10 As Aggarwal, K.M. & Keenan, F.P. 2004, A&A, 427, 763
816  * value changed from 54.9 to 68.9 */
817  /* >>chng 05 dec 19, update cs to
818  *>>refer fe10 cs Aggarwal, K.M. & Keenan, F.P. 2005, A&A, 439, 1215 */
819  cs = Fe_10_11_13_cs(
820  /* ion, one of 10, 11, or 13 on physics scale */
821  10,
822  /* initial and final index of levels, with lowest energy 1, highest of 5
823  * on physics scale */
824  1,
825  2 );
826 
827  PutCS(cs,&TauLines[ipFe106375]);
828  atom_level2(&TauLines[ipFe106375]);
829  /* c6374 = atom_pop2(cs ,4.,2.,69.5,2.26e4,dense.xIonDense(26,10))*3.12e-12
830  * dCooldT = dCooldT + c6374 * 2.26e4 * tsq1
831  * call CoolAdd( 'Fe10' , 6374 , c6374 )
832  *
833  * this is the E1 line that can pump [FeX] */
834  cs = 0.85*sexp(0.045*1.259e6/phycon.te);
835  cs = MAX2(0.05,cs);
836  PutCS(cs,&TauLines[ipT352]);
837  atom_level2(&TauLines[ipT352]);
838 
839  /* Fe XI fe11 fe 11, 7892, 6.08 mic, CS from
840  * >>refer fe11 as Mason, H. 1975, MNRAS 170, 651
841  * A from
842  * >>refer fe11 as Mendoza, C., & Zeippen, C.J. 1983, MNRAS, 202, 981
843  * following reference has very extensive energy levels and As
844  * >>refer fe11 as Fritzsche, S., Dong, C.Z., & Traebert, E., 2000, MNRAS, 318, 263 */
845  /*cs = 0.27;*/
846  /* >>chng 97 jan 31, set cs to unity, above ignored resonances */
847  /*cs = 1.;*/
848  /* >>chng 00 dec 05, update Fe11 CS to Tayal 2000 */
849  /* >>referold fe11 cs Tayal, S.S., 2000, ApJ, 544, 575-580 */
850  /* this is the low to mid transition */
851  /*cs = 0.0282 * phycon.te30*phycon.te05*phycon.te01;*/
852  /* CS is about 2.0 across broad temp range in following reference:
853  * >>refer Fe11 cs Aggarwal, K.M., & Keenan, F.P. 2003, MNRAS, 338, 412 */
854  /*cs = 2.;
855  PutCS(cs,&TauLines[ipTFe07]);*/
856  /* Tayal value is 10x larger than previous values */
857  /* Aggarwal & Keenan is about same as Tayal */
858  /*cs = 0.48; */
859  /*cs = 0.50;
860  PutCS(cs,&TauLines[ipTFe61]);*/
861  /* Tayal value is 3x larger than previous values */
862  /*cs = 0.35; tayal number */
863  /*cs = 0.55;
864  PutCS(cs,&TauDummy);
865  atom_level3(&TauLines[ipTFe07],&TauLines[ipTFe61],&TauDummy);*/
866 
867  /* >>refer fe11 cs Kafatos, M., & Lynch, J.P. 1980, ApJS, 42, 611 */
868  /*CoolHeavy.c1467 = atom_pop3(9.,5.,1.,.24,.028,.342,100.,1012.,9.3,5.43e4,
869  6.19e4,&p2,dense.xIonDense[ipIRON][11-1],0.,0.,0.)*1012.*1.36e-11;
870  CoolHeavy.c2649 = p2*91.0*7.512e-12;*/
871  /*CoolAdd("Fe11",1467,CoolHeavy.c1467);
872  CoolAdd("Fe11",2469,CoolHeavy.c2649);*/
873 
874  /* >>chng 05 dec 18, use give level Fe 11 model */
875  Fe11Lev5();
876 
877  /* A's for 2-1 and 3-1 from
878  * >>refer fe12 as Tayal, S.S., & Henry, R.J.W. 1986, ApJ, 302, 200
879  * CS fro
880  * >>refer fe12 cs Tayal, S.S., Henry, R.J.W., Pradhan, A.K. 1987, ApJ, 319, 951
881  * a(3-2) scaled from both ref */
882  CoolHeavy.c2568 = atom_pop3(4.,10.,6.,0.72,0.69,2.18,8.1e4,1.84e6,1.33e6,
883  6.37e4,4.91e4,&p2,dense.xIonDense[ipIRON][12-1],0.,0.,0.)*1.33e6*6.79e-12;
884  CoolAdd("Fe12",2568,CoolHeavy.c2568);
885  CoolHeavy.c1242 = CoolHeavy.c2568*2.30*1.38;
886  CoolAdd("Fe12",1242,CoolHeavy.c1242);
887  CoolHeavy.c2170 = p2*8.09e4*8.82e-12;
888  CoolAdd("Fe12",2170,CoolHeavy.c2170);
889 
890  /* [Fe XIII] fe13 fe 13 1.07, 1.08 microns */
891  /* >>chng 00 dec 05, update Fe13 CS to Tayal 2000 */
892  /* >>refer fe13 cs Tayal, S.S., 2000, ApJ, 544, 575-580
893  cs = 5.08e-3 * phycon.te30* phycon.te10;
894  PutCS(cs,&TauLines[ipFe1310]);
895  PutCS(2.1,&TauLines[ipFe1311]);
896  PutCS(.46,&TauDummy);
897  atom_level3(&TauLines[ipFe1310],&TauLines[ipFe1311],&TauDummy); */
898 
899  /* Fe13 lines from 1D and 1S -
900  >>chng 01 aug 07, added these lines */
901  /* >>refer fe11 cs Tayal, S.S., 2000, ApJ, 544, 575-580 */
902  /* >>refer fe13 as Shirai, T., Sugar, J., Musgrove, A., & Wiese, W.L., 2000.,
903  * >>refercon J Phys Chem Ref Data Monograph 8 */
904  /* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2)
905  CoolHeavy.fe13_1216 = atom_pop3( 9.,5.,1., 5.08,0.447,0.678, 103.2 ,1010.,7.99,
906  48068.,43440.,&p2,dense.xIonDense[ipIRON][13-1],0.,0.,0.)*1010.*
907  1.64e-11;
908  CoolHeavy.fe13_2302 = CoolHeavy.fe13_1216*0.528*0.00791;
909  CoolHeavy.fe13_3000 = p2*103.2*7.72e-12;*/
910  /*CoolAdd("Fe13",1216,CoolHeavy.fe13_1216);
911  CoolAdd("Fe13",3000,CoolHeavy.fe13_3000);
912  CoolAdd("Fe13",2302,CoolHeavy.fe13_2302);*/
913 
914  /* >>chng 05 dec 18, use give level Fe 13 model */
915  Fe13Lev5();
916 
917  /* Fe XIV 5303, cs from
918  * >>refer Fe14 cs Storey, P.J., Mason, H.E., Saraph, H.E., 1996, A&A, 309, 677 */
919  fe14cs(phycon.alogte,&cs);
920  /* >>chng 97 jan 31, set cs to unity, above is VERY large, >10 */
921  /* >>chng 01 may 30, back to their value, as per discussion at Lexington 2000 */
922  /* >>chng 05 aug 04, following line commented out, had set 5303 to constant value */
923  /*cs = 3.1;*/
924  /* >>chng 05 dec 22, A from 60.3 to 59.7, experimental value 59.7 +/- 0.4 in
925  * >>refer Fe14 as Trabert, E. 2004, A&A, 415, L39 */
926  CoolHeavy.c5303 = atom_pop2(cs,2.,4.,59.7,2.71e4,dense.xIonDense[ipIRON][14-1])*
927  3.75e-12;
929  CoolAdd("Fe14",5303,CoolHeavy.c5303);
930 
931  /* Fe 18 974.8A;, cs from
932  * >>referold fe18 cs Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */
933  /* >>refer fe18 cs Berrington,K.A.,Saraph, H.E. & Tully, J.A. 1998, A&AS, 129, 161 */
934  /*>>chng 06 jul 19 Changes made-Humeshkar Nemala*/
935  /*cs = MIN2(0.143,0.804/(phycon.te10*phycon.te05));*/
936  if(phycon.te < 1.29E6)
937  {
938  cs = (realnum)(0.3465/((phycon.te10/phycon.te01)*phycon.te001*phycon.te0003));
939  }
940  else if(phycon.te < 5.135E6)
941  {
942  cs = (realnum)((1.1062E-02)*phycon.te10*phycon.te05*phycon.te003*phycon.te0005);
943  }
944  else
945  {
946  cs = (realnum)((60.5728)/(phycon.te40*phycon.te003*phycon.te0001*phycon.te0005));
947  }
948 
949  PutCS(cs,&TauLines[ipFe18975]);
950  atom_level2(&TauLines[ipFe18975]);
951 
952  /* O-like Fe19, 3P ground term, 7046.72A vac wl, 1328.90A */
953  /* >>refer fe19 cs Butler, K., & Zeippen, C.J., 2001, A&A, 372, 1083 */
954  cs12 = 0.0627 / phycon.te03;
955  cs13 = 0.692 /(phycon.te10*phycon.te01);
956  cs23 = 0.04;
957  /*CoolHeavy.c7082 = atom_pop3(5.,1.,3.,0.0132,0.0404,0.0137,0.505,1.46e4,*/
958  /* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
959  CoolHeavy.c7082 = atom_pop3(5.,1.,3.,cs12,cs13,cs23,0.505,1.46e4,
960  41.2,1.083e5,2.03e4,&p2,dense.xIonDense[ipIRON][19-1],0.,0.,0.)*41.2*
961  2.81e-12;
962  CoolHeavy.c1118 = CoolHeavy.c7082*354.4*6.335;
963  CoolHeavy.c1328 = p2*0.505*1.50e-11;
964  CoolAdd("Fe19",7082,CoolHeavy.c7082);
965  CoolAdd("Fe19",1118,CoolHeavy.c1118);
966  CoolAdd("Fe19",1328,CoolHeavy.c1328);
967 
968  CoolHeavy.c592 = atom_pop2(0.0913,9.,5.,1.64e4,2.428e5,dense.xIonDense[ipIRON][19-1])*
969  3.36e-11;
970  CoolAdd("Fe19",592,CoolHeavy.c592);
971 
972  /* >>chng 01 aug 10, add this line */
973  /* Fe20 721.40A, 578*/
974  /* >>refer fe20 cs Butler, K., & Zeippen, C.J., 2001, A&A, 372, 1078 */
975  /*>>refer fe20 as Merkelis, G., Martinson, I., Kisielius, R., & Vilkas, M.J., 1999,
976  >>refercon Physica Scripta, 59, 122 */
977  cs = 1.17 /(phycon.te20/phycon.te01);
978  PutCS(cs , &TauLines[ipTFe20_721]);
979  cs = 0.248 /(phycon.te10/phycon.te01);
980  PutCS(cs , &TauLines[ipTFe20_578]);
981  cs = 0.301 /(phycon.te10/phycon.te02);
982  PutCS(cs , &TauDummy);
983  atom_level3(&TauLines[ipTFe20_721],&TauLines[ipTFe20_578],&TauDummy);
984 
985  /* >>chng 00 oct 26, much larger CS for following transition */
986  /* Fe 21 1354, 2299 A, cs eval at 1e6K
987  * >>refer Fe21 cs Aggarwall, K.M., 1991, ApJS 77, 677 */
988  PutCS(0.072,&TauLines[ipTFe13]);
989  PutCS(0.269,&TauLines[ipTFe23]);
990  PutCS(0.055,&TauDummy);
991  atom_level3(&TauLines[ipTFe13],&TauLines[ipTFe23],&TauDummy);
992 
993  /*>>refer Fe22 energy Feldman, U., Curdt, W., Landi, E., Wilhelm, K., 2000, ApJ, 544, 508 */
994  /*>>chng 00 oct 26, added these fe 21 lines, removed old CoolHeavy.fs846, the lowest */
995 
996  /* one time initialization if first call, and level 2 lines are on */
997  if( lgFe22First && nWindLine )
998  {
999  lgFe22First = false;
1000  nFe22Pump = 0;
1001  for( i=0; i<nWindLine; ++i )
1002  {
1003  /* don't test on nelem==ipIRON since lines on physics, not C, scale */
1004  if( TauLine2[i].Hi->nelem ==26 && TauLine2[i].Hi->IonStg==22 )
1005  {
1006  ++nFe22Pump;
1007  }
1008  }
1009  if( nFe22Pump<0 )
1010  TotalInsanity();
1011  else if( nFe22Pump > 0 )
1012  /* create the space - can't malloc 0 bytes */
1013  ipFe22Pump = (long *)MALLOC((unsigned)(nFe22Pump)*sizeof(long) );
1014  nFe22Pump = 0;
1015  for( i=0; i<nWindLine; ++i )
1016  {
1017  /* don't test on nelem==ipIRON since lines on physics, not C, scale */
1018  if( TauLine2[i].Hi->nelem ==26 && TauLine2[i].Hi->IonStg==22 )
1019  {
1020 # if 0
1021  DumpLine( &TauLine2[i] );
1022 # endif
1023  ipFe22Pump[nFe22Pump] = i;
1024  ++nFe22Pump;
1025  }
1026  }
1027  }
1028  else
1029  /* level 2 lines are not enabled */
1030  nFe22Pump = 0;
1031 
1032  /* now sum pump rates */
1033  Fe22_pump_rate = 0.;
1034  for( i=0; i<nFe22Pump; ++i )
1035  {
1036  Fe22_pump_rate += TauLine2[ipFe22Pump[i]].Emis->pump;
1037 # if 0
1038  fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n",
1039  i,
1040  TauLine2[ipFe22Pump[i]].WLAng , TauLine2[ipFe22Pump[i]].pump );
1041 # endif
1042  }
1043 
1044  /*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
1045  /*>>refer fe22 cs Zhang, H.L., Graziani, M., & Pradhan, A.K., 1994, A&A 283, 319*/
1046  /*>>refer fe22 as Dankwort, W., & Trefftz, E., 1978, A&A 65, 93-98 */
1048  &TauLines[ipFe22_247],
1049  &TauLines[ipFe22_217],
1050  &TauLines[ipFe22_348],
1051  &TauLines[ipFe22_292],
1052  &TauLines[ipFe22_253],
1053  /*double cs40,
1054  double cs32,
1055  double cs42,
1056  double cs43, */
1057  0.00670 , 0.0614 , 0.0438 , 0.122 , Fe22_pump_rate ,"Fe22");
1058 
1059  /* Fe 22 845.6, C.S.=guess, A from
1060  * >>refer fe22 as Froese Fischer, C. 1983, J.Phys. B, 16, 157
1061  CoolHeavy.fs846 = atom_pop2(0.2,2.,4.,1.5e4,1.701e5,dense.xIonDense[ipIRON][22-1])*
1062  2.35e-11;
1063  CoolAdd("Fe22",846,CoolHeavy.fs846);*/
1064 
1066  /* FE 23 262.6, 287A, 1909-LIKE,
1067  * cs from
1068  * >>refer fe23 cs Bhatia, A.K., & Mason, H.E. 1986, A&A, 155, 413 */
1069  CoolHeavy.c263 = atom_pop2(0.04,1.,9.,1.6e7,5.484e5,dense.xIonDense[ipIRON][23-1])*
1070  7.58e-11;
1071  CoolAdd("Fe23",262,CoolHeavy.c263);
1072 
1073  /* FE 24, 255, 192 Li seq doublet
1074  * >>refer fe24 cs Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
1075  ligbar(26,&TauLines[ipT192],&TauLines[ipT11],&cs2s2p,&cs2s3p);
1076 
1077  PutCS(cs2s2p,&TauLines[ipT192]);
1078  atom_level2(&TauLines[ipT192]);
1079 
1080  /* funny factor (should have been 0.5) due to energy change */
1081  PutCS(cs2s2p*0.376,&TauLines[ipT255]);
1082  atom_level2(&TauLines[ipT255]);
1083 
1084  PutCS(cs2s3p,&TauLines[ipT11]);
1085  atom_level2(&TauLines[ipT11]);
1086 
1089  TauLines[ipT353].Lo->Pop = dense.xIonDense[ipIRON][11-1];
1090  TauLines[ipT353].Hi->Pop = 0.;
1092  TauLines[ipT347].Lo->Pop = dense.xIonDense[ipIRON][14-1];
1093  TauLines[ipT347].Hi->Pop = 0.;
1094 
1095  return;
1096 }
1097 
1098 /*fe14cs compute collision strength for forbidden transition
1099  * w/in Fe XIV ground term. From
1100  * >>refer fe14 cs Storey, P.J., Mason, H.E., Saraph, H.E., 1996, A&A, 309, 677
1101  * */
1102 STATIC void fe14cs(double te1,
1103  double *csfe14)
1104 {
1105  double a,
1106  b,
1107  c,
1108  d,
1109  telog1,
1110  telog2,
1111  telog3;
1112 
1113  DEBUG_ENTRY( "fe14cs()" );
1114 
1115  /* limit range in log T: */
1116  telog1 = te1;
1117  telog1 = MIN2(9.0,telog1);
1118  telog1 = MAX2(4.0,telog1);
1119 
1120  /* compute square and cube */
1121  telog2 = telog1*telog1;
1122  telog3 = telog2*telog1;
1123 
1124  /* compute cs:
1125  * */
1126  if( telog1 <= 5.0 )
1127  {
1128  a = 557.05536;
1129  b = -324.56109;
1130  c = 63.437974;
1131  d = -4.1365147;
1132  *csfe14 = a + b*telog1 + c*telog2 + d*telog3;
1133  }
1134  else
1135  {
1136  a = 0.19515493;
1137  b = 2.9404407;
1138  c = 4.9578944;
1139  d = 0.79887506;
1140  *csfe14 = a + b*exp(-0.5*((telog1-c)*(telog1-c)/d));
1141  }
1142  return;
1143 }
1144 
1145 /*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
1146 STATIC void Fe4Lev12(void)
1147 {
1148  const int NLFE4 = 12;
1149  bool lgZeroPop;
1150  int nNegPop;
1151  long int i,
1152  j;
1153  static bool lgFirst=true;
1154 
1155  double dfe4dt;
1156 
1157  /*static long int **ipdest; */
1158  static double
1159  **AulEscp,
1160  **col_str,
1161  **AulDest,
1162  depart[NLFE4],
1163  pops[NLFE4],
1164  destroy[NLFE4],
1165  create[NLFE4],
1166  **CollRate,
1167  **AulPump;
1168 
1169  static const double Fe4A[NLFE4][NLFE4] = {
1170  {0.,0.,0.,1.e-5,0.,1.368,.89,0.,1.3e-3,1.8e-4,.056,.028},
1171  {0.,0.,2.6e-8,0.,0.,0.,0.,0.,1.7e-7,0.,0.,0.},
1172  {0.,0.,0.,0.,3.5e-7,6.4e-10,0.,0.,6.315e-4,0.,6.7e-7,0.},
1173  {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},
1174  {0.,0.,0.,0.,0.,1.5e-5,1.3e-9,0.,7.6e-4,0.,1.1e-6,6.0e-7},
1175  {0.,0.,0.,0.,0.,0.,1.1e-5,1.2e-13,.038,9.9e-7,.022,.018},
1176  {0.,0.,0.,0.,0.,0.,0.,3.7e-5,2.9e-6,.034,3.5e-3,.039},
1177  {0.,0.,0.,0.,0.,0.,0.,0.,0.,.058,3.1e-6,1.4e-3},
1178  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-4,3.1e-14},
1179  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.9e-19,1.0e-5},
1180  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-7},
1181  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}
1182  };
1183  static const double Fe4CS[NLFE4][NLFE4] = {
1184  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1185  {0.98,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1186  {0.8167,3.72,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1187  {0.49,0.0475,0.330,0.,0.,0.,0.,0.,0.,0.,0.,0.},
1188  {0.6533,0.473,2.26,1.64,0.,0.,0.,0.,0.,0.,0.,0.},
1189  {0.45,0.686,0.446,0.106,0.254,0.,0.,0.,0.,0.,0.,0.},
1190  {0.30,0.392,0.152,0.269,0.199,0.605,0.,0.,0.,0.,0.,0.},
1191  {0.15,0.0207,0.190,0.0857,0.166,0.195,0.327,0.,0.,0.,0.,0.},
1192  {0.512,1.23,0.733,0.174,0.398,0.623,0.335,0.102,0.,0.,0.,0.},
1193  {0.128,0.0583,0.185,0.200,0.188,0.0835,0.127,0.0498,0.0787,0.,0.,0.},
1194  {0.384,0.578,0.534,0.363,0.417,0.396,0.210,0.171,0.810,0.101,0.,0.},
1195  {0.256,0.234,0.306,0.318,0.403,0.209,0.195,0.112,0.195,0.458,0.727,0.}
1196  };
1197 
1198  static const double gfe4[NLFE4]={6.,12.,10.,6.,8.,6.,4.,2.,8.,2.,6.,4.};
1199 
1200  /* excitation energies in Kelvin
1201  static double ex[NLFE4]={0.,46395.8,46464.,46475.6,46483.,50725.,
1202  50838.,50945.,55796.,55966.,56021.,56025.};*/
1203  /*>>refer Fe3 energies version 3 NIST Atomic Spectra Database */
1204  /*>>chng 05 dec 17, from Kelvin above to excitation energies in wn */
1205  static const double excit_wn[NLFE4]={0.,32245.5,32292.8,32301.2,32305.7,35253.8,
1206  35333.3,35406.6,38779.4,38896.7,38935.1,38938.2};
1207 
1208  DEBUG_ENTRY( "Fe4Lev12()" );
1209 
1210  if( lgFirst )
1211  {
1212  /* will never do this again */
1213  lgFirst = false;
1214  /* allocate the 1D arrays*/
1215  /* create space for the 2D arrays */
1216  AulPump = ((double **)MALLOC((NLFE4)*sizeof(double *)));
1217  CollRate = ((double **)MALLOC((NLFE4)*sizeof(double *)));
1218  AulDest = ((double **)MALLOC((NLFE4)*sizeof(double *)));
1219  AulEscp = ((double **)MALLOC((NLFE4)*sizeof(double *)));
1220  col_str = ((double **)MALLOC((NLFE4)*sizeof(double *)));
1221  for( i=0; i<(NLFE4); ++i )
1222  {
1223  AulPump[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
1224  CollRate[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
1225  AulDest[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
1226  AulEscp[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
1227  col_str[i] = ((double *)MALLOC((NLFE4)*sizeof(double )));
1228  }
1229  }
1230 
1231  /* bail if no Fe+3 */
1232  if( dense.xIonDense[ipIRON][3] <= 0. )
1233  {
1234  fe.Fe4CoolTot = 0.;
1235  fe.fe40401 = 0.;
1236  fe.fe42836 = 0.;
1237  fe.fe42829 = 0.;
1238  fe.fe42567 = 0.;
1239  fe.fe41207 = 0.;
1240  fe.fe41206 = 0.;
1241  fe.fe41106 = 0.;
1242  fe.fe41007 = 0.;
1243  fe.fe41008 = 0.;
1244  fe.fe40906 = 0.;
1245  CoolAdd("Fe 4",0,0.);
1246 
1247  /* level populations */
1248  /* LIMLEVELN is the dimension of the atoms vectors */
1249  ASSERT( NLFE4 <= LIMLEVELN);
1250  for( i=0; i < NLFE4; i++ )
1251  {
1252  atoms.PopLevels[i] = 0.;
1253  atoms.DepLTELevels[i] = 1.;
1254  }
1255  return;
1256  }
1257  /* number of levels in model ion */
1258 
1259  /* these are in wavenumbers
1260  * data excit_wn/ 0., 32245.5, 32293., 32301.2,
1261  * 1 32306., 35254., 35333., 35407., 38779., 38897., 38935.,
1262  * 2 38938./
1263  * excitation energies in Kelvin */
1264 
1265  /* A's are from Garstang, R.H., MNRAS 118, 572 (1958).
1266  * each set is for a lower level indicated by second element in array,
1267  * index runs over upper level
1268  * A's are saved into arrays as data(up,lo) */
1269 
1270  /* collision strengths from Berrington and Pelan Ast Ap S 114, 367.
1271  * order is cs(low,up) */
1272 
1273  /* all elements are used, and must be set to zero if zero */
1274  for( i=0; i < NLFE4; i++ )
1275  {
1276  create[i] = 0.;
1277  destroy[i] = 0.;
1278  for( j=0; j < NLFE4; j++ )
1279  {
1280  /*data[j][i] = 1e33;*/
1281  col_str[j][i] = 0.;
1282  AulEscp[j][i] = 0.;
1283  AulDest[j][i] = 0.;
1284  AulPump[j][i] = 0.;
1285  }
1286  }
1287 
1288  /* fill in einstein As and collision strengths */
1289  for( i=0; i < NLFE4; i++ )
1290  {
1291  for( j=i + 1; j < NLFE4; j++ )
1292  {
1293  /*data[i][j] = Fe4A[i][j];*/
1294  AulEscp[j][i] = Fe4A[i][j];
1295  /*data[j][i] = Fe4CS[j][i];*/
1296  col_str[j][i] = Fe4CS[j][i];
1297  }
1298  }
1299 
1300  /* leveln itself is well-protected against zero abundances,
1301  * low temperatures */
1302 
1303  atom_levelN(NLFE4,
1304  dense.xIonDense[ipIRON][3],
1305  gfe4,
1306  excit_wn,
1307  'w',
1308  pops,
1309  depart,
1310  &AulEscp ,
1311  &col_str ,
1312  &AulDest,
1313  &AulPump,
1314  &CollRate,
1315  create,
1316  destroy,
1317  /* say atom_levelN should evaluate coll rates from cs */
1318  false,
1319  &fe.Fe4CoolTot,
1320  &dfe4dt,
1321  "FeIV",
1322  /* nNegPop positive if negative pops occured, negative if too cold */
1323  &nNegPop,
1324  &lgZeroPop,
1325  false );
1326 
1327  /* LIMLEVELN is the dim of the PopLevels vector */
1328  ASSERT( NLFE4 <= LIMLEVELN );
1329  for( i=0; i<NLFE4; ++i)
1330  {
1331  atoms.PopLevels[i] = pops[i];
1332  atoms.DepLTELevels[i] = depart[i];
1333  }
1334 
1335  if( nNegPop>0 )
1336  {
1337  fprintf( ioQQQ, " fe4levl2 found negative populations\n" );
1338  ShowMe();
1339  cdEXIT(EXIT_FAILURE);
1340  }
1341 
1342  CoolAdd("Fe 4",0,fe.Fe4CoolTot);
1343 
1344  thermal.dCooldT += dfe4dt;
1345 
1346  /* three transitions hst can observe
1347  * 4 -1 3095.9A and 5 -1 3095.9A */
1348  fe.fe40401 = (pops[3]*Fe4A[0][3]*(excit_wn[3] - excit_wn[0]) +
1349  pops[4]*Fe4A[0][4]*(excit_wn[4] - excit_wn[0]) )*T1CM*BOLTZMANN;
1350 
1351  fe.fe42836 = pops[5]*Fe4A[0][5]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN;
1352 
1353  fe.fe42829 = pops[6]*Fe4A[0][6]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN;
1354 
1355  fe.fe42567 = (pops[10]*Fe4A[0][10]*(excit_wn[10] - excit_wn[0]) +
1356  pops[11]*Fe4A[0][11]*(excit_wn[10] - excit_wn[0]))*T1CM*BOLTZMANN;
1357 
1358  fe.fe41207 = pops[11]*Fe4A[6][11]*(excit_wn[11] - excit_wn[6])*T1CM*BOLTZMANN;
1359  fe.fe41206 = pops[11]*Fe4A[5][11]*(excit_wn[11] - excit_wn[5])*T1CM*BOLTZMANN;
1360  fe.fe41106 = pops[10]*Fe4A[5][10]*(excit_wn[10] - excit_wn[5])*T1CM*BOLTZMANN;
1361  fe.fe41007 = pops[9]*Fe4A[6][9]*(excit_wn[9] - excit_wn[6])*T1CM*BOLTZMANN;
1362  fe.fe41008 = pops[9]*Fe4A[7][9]*(excit_wn[9] - excit_wn[7])*T1CM*BOLTZMANN;
1363  fe.fe40906 = pops[8]*Fe4A[5][8]*(excit_wn[8] - excit_wn[5])*T1CM*BOLTZMANN;
1364  return;
1365 }
1366 
1367 /*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
1368 STATIC void Fe7Lev8(void)
1369 {
1370  bool lgZeroPop;
1371  int nNegPop;
1372  double scale;
1373  long int i,
1374  j;
1375  static bool lgFirst=true;
1376  static long int ipPump=-1;
1377 
1378  double dfe7dt,
1379  FUV_pump;
1380 
1381  long int ihi , ilo;
1382  static double
1383  *depart,
1384  *pops,
1385  *destroy,
1386  *create ,
1387  **AulDest,
1388  **CollRate,
1389  **AulPump,
1390  **AulNet,
1391  **col_str;
1392  /*
1393  following are FeVII lines predicted in limit_laser_2000
1394  Fe 7 5721A -24.399 0.1028
1395  Fe 7 4989A -24.607 0.0637
1396  Fe 7 4894A -24.838 0.0374
1397  Fe 7 4699A -25.693 0.0052
1398  Fe 7 6087A -24.216 0.1566
1399  Fe 7 5159A -24.680 0.0538
1400  Fe 7 4943A -25.048 0.0231
1401  Fe 7 3587A -24.285 0.1336
1402  Fe 7 5277A -25.053 0.0228
1403  Fe 7 3759A -24.139 0.1870
1404  Fe 7 3.384m -25.621 0.0062
1405  Fe 7 2.629m -25.357 0.0113
1406  Fe 7 9.510m -24.467 0.0878
1407  Fe 7 7.810m -24.944 0.0293
1408  */
1409 
1410  /* statistical weights for the n levels */
1411  static double gfe7[NLFE7]={5.,7.,9.,5.,1.,3.,5.,9.};
1412 
1413  /*refer Fe7 energies Ekbert, J.O. 1981, Phys. Scr 23, 7 */
1414  /*static double ex[NLFE7]={0.,1047.,2327.,17475.,20037.,20428.,21275.,28915.};*/
1415  /* excitation energies in wavenumbers
1416  * >>chng 05 dec 15, rest of code had assumed that there were energies in Kelvin
1417  * rather than wavenumbers. Bug caught by Kevin Blagrave
1418  * modified atom_levelN to accept either Kelvin or wavenumbers */
1419  static double excit_wn[NLFE7]={0. , 1051.5 , 2331.5 , 17475.5 , 20040.3 , 20430.1 , 21278.6 , 28927.3 };
1420 
1421  DEBUG_ENTRY( "Fe7Lev8()" );
1422 
1423  if( lgFirst )
1424  {
1425  /* will never do this again */
1426  lgFirst = false;
1427  /* allocate the 1D arrays*/
1428  /* create space for the 2D arrays */
1429  depart = ((double *)MALLOC((NLFE7)*sizeof(double)));
1430  pops = ((double *)MALLOC((NLFE7)*sizeof(double)));
1431  destroy = ((double *)MALLOC((NLFE7)*sizeof(double)));
1432  create = ((double *)MALLOC((NLFE7)*sizeof(double)));
1433  /* now the 2-d arrays */
1434  fe.Fe7_wl = ((double **)MALLOC((NLFE7)*sizeof(double *)));
1435  fe.Fe7_emiss = ((double **)MALLOC((NLFE7)*sizeof(double *)));
1436  AulNet = ((double **)MALLOC((NLFE7)*sizeof(double *)));
1437  col_str = ((double **)MALLOC((NLFE7)*sizeof(double *)));
1438  AulPump = ((double **)MALLOC((NLFE7)*sizeof(double *)));
1439  CollRate = ((double **)MALLOC((NLFE7)*sizeof(double *)));
1440  AulDest = ((double **)MALLOC((NLFE7)*sizeof(double *)));
1441  for( i=0; i<(NLFE7); ++i )
1442  {
1443  fe.Fe7_wl[i] = ((double *)MALLOC((NLFE7)*sizeof(double )));
1444  fe.Fe7_emiss[i] = ((double *)MALLOC((NLFE7)*sizeof(double )));
1445  AulNet[i] = ((double *)MALLOC((NLFE7)*sizeof(double )));
1446  col_str[i] = ((double *)MALLOC((NLFE7)*sizeof(double )));
1447  AulPump[i] = ((double *)MALLOC((NLFE7)*sizeof(double )));
1448  CollRate[i] = ((double *)MALLOC((NLFE7)*sizeof(double )));
1449  AulDest[i] = ((double *)MALLOC((NLFE7)*sizeof(double )));
1450  }
1451 
1452  /* set some to constant values after zeroing out */
1453  for( i=0; i<NLFE7; ++i )
1454  {
1455  create[i] = 0.;
1456  destroy[i] = 0.;
1457  for( j=0; j<NLFE7; ++j )
1458  {
1459  AulNet[i][j] = 0.;
1460  col_str[i][j] = 0.;
1461  CollRate[i][j] = 0.;
1462  AulDest[i][j] = 0.;
1463  AulPump[i][j] = 0.;
1464  fe.Fe7_wl[i][j] = 0.;
1465  fe.Fe7_emiss[i][j] = 0.;
1466  }
1467  }
1468  set_NaN( fe.Fe7_wl[2][1] );
1469  set_NaN( fe.Fe7_emiss[2][1] );
1470  set_NaN( fe.Fe7_wl[1][0] );
1471  set_NaN( fe.Fe7_emiss[1][0] );
1472  /* two levels within ground term must never be addressed in this array -
1473  * they are fully transferred */
1474  for( ilo=0; ilo<NLFE7-1; ++ilo )
1475  {
1476  /* must not do 1-0 or 2-1, which are transferred lines */
1477  for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
1478  {
1479  fe.Fe7_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / RefIndex( excit_wn[ihi]-excit_wn[ilo] );
1480  }
1481  }
1482 
1483  /* assume FeVII is optically thin - just use As as net escape */
1484  /*>>refer Fe7 As Berrington, K.A., Nakazaki, S., & Norrington, P.H. 2000, A&AS, 142, 313 */
1485  AulNet[1][0] = 0.0325;
1486  AulNet[2][0] = 0.167e-8;
1487  /* following corrected from 3.25 to 0.372 as per Keith Berrington email */
1488  AulNet[3][0] = 0.372;
1489  AulNet[4][0] = 0.135;
1490  AulNet[5][0] = 0.502e-1;
1491  /* following corrected from 0.174 to 0.0150 as per Keith Berrington email */
1492  AulNet[6][0] = 0.0150;
1493  AulNet[7][0] = 0.959e-3;
1494 
1495  AulNet[2][1] = 0.466e-1;
1496  AulNet[3][1] = 0.603;
1497  AulNet[5][1] = 0.762e-1;
1498  AulNet[6][1] = 0.697e-1;
1499  AulNet[7][1] = 0.343;
1500 
1501  AulNet[3][2] = 0.139e-2;
1502  AulNet[6][2] = 0.735e-1;
1503  AulNet[7][2] = 0.503;
1504 
1505  AulNet[4][3] = 0.472e-6;
1506  AulNet[5][3] = 0.572e-1;
1507  /* following corrected from 0.191 to 0.182 as per Keith Berrington email */
1508  AulNet[6][3] = 0.182;
1509  AulNet[7][3] = 0.414e-2;
1510 
1511  AulNet[5][4] = 0.115e-2;
1512  AulNet[6][4] = 0.139e-7;
1513 
1514  AulNet[6][5] = 0.743e-2;
1515 
1516  AulNet[7][6] = 0.454e-4;
1517 
1518  /* collision strengths at log T 4.5 */
1520  /*>>refer Fe7 cs Berrington, K.A., Nakazaki, S., & Norrington, P.H. 2000, A&AS, 142, 313 */
1521 # if 0
1522  if( fudge(-1) )
1523  {
1524  fprintf(ioQQQ,"DEBUG fudge call cool_iron\n");
1525  scale = fudge(0);
1526  }
1527  else
1528 # endif
1529  scale = 1.;
1530 
1531  col_str[1][0] = 3.35;
1532  col_str[2][0] = 1.17;
1533  col_str[3][0] = 0.959;
1534  col_str[4][0] = 0.299;
1535  col_str[5][0] = 0.633;
1536  col_str[6][0] = 0.549;
1537  col_str[7][0] = 1.24*scale;
1538 
1539  col_str[2][1] = 4.11;
1540  col_str[3][1] = 1.29;
1541  col_str[4][1] = 0.235;
1542  col_str[5][1] = 0.833;
1543  col_str[6][1] = 1.06;
1544  col_str[7][1] = 1.74*scale;
1545 
1546  col_str[3][2] = 1.60;
1547  col_str[4][2] = 0.187;
1548  col_str[5][2] = 0.690;
1549  col_str[6][2] = 1.94;
1550  col_str[7][2] = 2.25*scale;
1551 
1552  col_str[4][3] = 0.172;
1553  col_str[5][3] = 0.531;
1554  col_str[6][3] = 1.06;
1555  col_str[7][3] = 2.02;
1556 
1557  col_str[5][4] = 0.370;
1558  col_str[6][4] = 0.324;
1559  col_str[7][4] = 0.164;
1560 
1561  col_str[6][5] = 1.17;
1562  col_str[7][5] = 0.495;
1563 
1564  col_str[7][6] = 0.903;
1565 
1566  /* check whether level 2 lines are on, and if so, find the driving lines that
1567  * can pump the upper levels of this atom */
1568  if( nWindLine > 0 )
1569  {
1570  ipPump = -1;
1571  for( i=0; i<nWindLine; ++i )
1572  {
1573  /* don't test on nelem==ipIRON since lines on physics, not C, scale */
1574  if( TauLine2[i].Hi->nelem ==26 && TauLine2[i].Hi->IonStg==7 &&
1575  /* >>chng 05 jul 10, move line to wavelength that agrees with nist tables
1576  * here and in level2 data file
1577  * NB must add a few wn to second number to get hit -
1578  * logic is that this is lowest E1 transition */
1579  (TauLine2[i].EnergyWN-4.27360E+05) < 0. )
1580  {
1581  ipPump = i;
1582  break;
1583  }
1584  }
1585  if( ipPump<0 )
1586  {
1587  fprintf(ioQQQ,"PROBLEM Fe7Lev8 cannot identify the FUV driving line.\n");
1588  TotalInsanity();
1589  }
1590  }
1591  else
1592  {
1593  ipPump = 0;
1594  }
1595  }
1596 
1597  /* bail if no ions */
1598  if( dense.xIonDense[ipIRON][6] <= 0. )
1599  {
1600  CoolAdd("Fe 7",0,0.);
1601 
1602  fe.Fe7CoolTot = 0.;
1603  for( ilo=0; ilo<NLFE7-1; ++ilo )
1604  {
1605  /* must not do 1-0 or 2-1, which are transferred lines */
1606  for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
1607  {
1608  fe.Fe7_emiss[ihi][ilo] = 0.;
1609  }
1610  }
1611  TauLines[ipFe0795].Lo->Pop = 0.;
1612  TauLines[ipFe0795].Emis->PopOpc = 0.;
1613  TauLines[ipFe0795].Hi->Pop = 0.;
1615  TauLines[ipFe0795].Coll.cool = 0.;
1616  TauLines[ipFe0795].Emis->phots = 0.;
1617  TauLines[ipFe0795].Emis->ColOvTot = 0.;
1618  TauLines[ipFe0795].Coll.heat = 0.;
1619  CoolAdd( "Fe 7", TauLines[ipFe0795].WLAng , 0.);
1620  TauLines[ipFe0778].Lo->Pop = 0.;
1621  TauLines[ipFe0778].Emis->PopOpc = 0.;
1622  TauLines[ipFe0778].Hi->Pop = 0.;
1624  TauLines[ipFe0778].Coll.cool = 0.;
1625  TauLines[ipFe0778].Emis->phots = 0.;
1626  TauLines[ipFe0778].Emis->ColOvTot = 0.;
1627  TauLines[ipFe0778].Coll.heat = 0.;
1628  CoolAdd( "Fe 7", TauLines[ipFe0778].WLAng , 0.);
1629  /* level populations */
1630  /* LIMLEVELN is the dimension of the atoms vectors */
1631  ASSERT( NLFE7 <= LIMLEVELN);
1632  for( i=0; i < NLFE7; i++ )
1633  {
1634  atoms.PopLevels[i] = 0.;
1635  atoms.DepLTELevels[i] = 1.;
1636  }
1637  return;
1638  }
1639 
1640  /* do pump rate for FUV excitation of 3P (levels 5-7 on physics scale, not C scale) */
1641  if( ipPump )
1642  {
1643  FUV_pump = TauLine2[ipPump].Emis->pump * 0.3 /(0.3+TauLine2[ipPump].Emis->Pesc);
1644  }
1645  else
1646  {
1647  FUV_pump = 0.;
1648  }
1649 
1650  /* this represents photoexcitation of 3P from ground level
1651  * >>chng 04 nov 22, upper array elements were on physics not c scale, off by one, TE */
1652  AulPump[0][4] = FUV_pump;
1653  AulPump[1][4] = FUV_pump;
1654  AulPump[2][4] = FUV_pump;
1655  AulPump[0][5] = FUV_pump;
1656  AulPump[1][5] = FUV_pump;
1657  AulPump[2][5] = FUV_pump;
1658  AulPump[0][6] = FUV_pump;
1659  AulPump[1][6] = FUV_pump;
1660  AulPump[2][6] = FUV_pump;
1661  /* these were set in the previous call to atom_levelN as the previous pump times
1662  * ratio of statistical weights, so this is the downward pump rate */
1663  AulPump[4][0] = 0;
1664  AulPump[4][1] = 0;
1665  AulPump[4][2] = 0;
1666  AulPump[5][0] = 0;
1667  AulPump[5][1] = 0;
1668  AulPump[5][2] = 0;
1669  AulPump[6][0] = 0;
1670  AulPump[6][1] = 0;
1671  AulPump[6][2] = 0;
1672 
1673  /* within ground term update to rt results */
1675  AulDest[1][0] = TauLines[ipFe0795].Emis->Aul*TauLines[ipFe0795].Emis->Pdest;
1676  AulPump[0][1] = TauLines[ipFe0795].Emis->pump;
1677  AulPump[1][0] = 0.;
1678 
1680  AulDest[2][1] = TauLines[ipFe0778].Emis->Aul*TauLines[ipFe0778].Emis->Pdest;
1681  AulPump[1][2] = TauLines[ipFe0778].Emis->pump;
1682  AulPump[2][1] = 0.;
1683 
1684  /* nNegPop positive if negative pops occurred, negative if too cold */
1685  atom_levelN(
1686  /* number of levels */
1687  NLFE7,
1688  /* the abundance of the ion, cm-3 */
1689  dense.xIonDense[ipIRON][6],
1690  /* the statistical weights */
1691  gfe7,
1692  /* the excitation energies in wavenumbers */
1693  excit_wn,
1694  /* units are wavenumbers */
1695  'w',
1696  /* the derived populations - cm-3 */
1697  pops,
1698  /* the derived departure coefficients */
1699  depart,
1700  /* the net emission rate, Aul * escp prob */
1701  &AulNet ,
1702  /* the collision strengths */
1703  &col_str ,
1704  /* A * destruction prob */
1705  &AulDest,
1706  /* pumping rate */
1707  &AulPump,
1708  /* collision rate, s-1, must defined if no collision strengths */
1709  &CollRate,
1710  /* creation vector */
1711  create,
1712  /* destruction vector */
1713  destroy,
1714  /* say atom_levelN should evaluate coll rates from cs */
1715  false,
1716  &fe.Fe7CoolTot,
1717  &dfe7dt,
1718  "Fe 7",
1719  &nNegPop,
1720  &lgZeroPop,
1721  false );
1722 
1723  /* LIMLEVELN is the dim of the PopLevels vector */
1724  ASSERT( NLFE7 <= LIMLEVELN );
1725  for( i=0; i<NLFE7; ++i)
1726  {
1727  atoms.PopLevels[i] = pops[i];
1728  atoms.DepLTELevels[i] = depart[i];
1729  }
1730 
1731  if( lgZeroPop )
1732  {
1733  /* this case, too cool to excite upper levels of atom, but will want
1734  * to do IR lines in ground term */
1735  PutCS(col_str[1][0],&TauLines[ipFe0795]);
1736  PutCS(col_str[2][1],&TauLines[ipFe0778]);
1737  PutCS(col_str[2][0],&TauDummy);
1738  atom_level3(&TauLines[ipFe0795],&TauLines[ipFe0778],&TauDummy);
1742  for( ilo=0; ilo<NLFE7-1; ++ilo )
1743  {
1744  /* must not do 1-0 or 2-1, which are transferred lines */
1745  for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
1746  {
1747  fe.Fe7_emiss[ihi][ilo] = 0.;
1748  }
1749  }
1750  }
1751  else
1752  {
1753  /* this case non-zero pops, but must set up vars within transition structure */
1756  TauLines[ipFe0795].Emis->PopOpc = (pops[0] - pops[1]*gfe7[0]/gfe7[1]);
1759  TauLines[ipFe0795].Emis->ColOvTot = CollRate[0][1]/(CollRate[0][1]+TauLines[ipFe0795].Emis->pump);
1760  TauLines[ipFe0795].Coll.cool = 0.;
1761  TauLines[ipFe0795].Coll.heat = 0.;
1762 
1765  TauLines[ipFe0778].Emis->PopOpc = (pops[1] - pops[2]*gfe7[1]/gfe7[2]);
1768  TauLines[ipFe0778].Emis->ColOvTot = CollRate[1][2]/(CollRate[1][2]+TauLines[ipFe0778].Emis->pump);
1769  TauLines[ipFe0778].Coll.heat = 0.;
1770  TauLines[ipFe0778].Coll.cool = 0.;
1771  }
1772 
1773  if( nNegPop>0 )
1774  {
1775  fprintf( ioQQQ, "PROBLEM Fe7Lev8 found negative populations\n" );
1776  ShowMe();
1777  cdEXIT(EXIT_FAILURE);
1778  }
1779 
1780  /* add cooling then its derivative */
1781  CoolAdd("Fe 7",0,fe.Fe7CoolTot);
1782  /* derivative of cooling */
1783  thermal.dCooldT += dfe7dt;
1784 
1785  /* find emission in each line */
1786  for( ilo=0; ilo<NLFE7-1; ++ilo )
1787  {
1788  /* must not do 1-0 or 2-1, which are transferred lines */
1789  for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
1790  {
1791  /* emission in these lines */
1792  fe.Fe7_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*ERG1CM;
1793  }
1794  }
1795  return;
1796 }
1797 
1798 /*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion
1799  * >>chng 05 dec 17, code provided by Kevin Blagrave and described in
1800  *>>refer Fe3 model Blagrave, K.P.M., Martin, P.G. & Baldwin, J.A.
1801  *>>refercon 2006, ApJ, in press (astro-ph 0603167) */
1802 STATIC void Fe3Lev14(void)
1803 {
1804  bool lgZeroPop;
1805  int nNegPop;
1806  long int i,
1807  j;
1808  static bool lgFirst=true;
1809 
1810  double dfe3dt;
1811 
1812  long int ihi , ilo;
1813  static double
1814  *depart,
1815  *pops,
1816  *destroy,
1817  *create ,
1818  **AulDest,
1819  **CollRate,
1820  **AulPump,
1821  **AulNet,
1822  **col_str;
1823 
1824  /* statistical weights for the n levels */
1825  static double gfe3[NLFE3]={9.,7.,5.,3.,1.,5.,13.,11.,9.,3.,1.,9.,7.,5.};
1826 
1827  /*refer Fe3 energies NIST version 3 Atomic Spectra Database */
1828  /* from smallest to largest energy (not in multiplet groupings) */
1829 
1830  /* energy in wavenumbers */
1831  static double excit_wn[NLFE3]={
1832  0.0 , 436.2, 738.9, 932.4, 1027.3,
1833  19404.8, 20051.1, 20300.8, 20481.9, 20688.4,
1834  21208.5, 21462.2, 21699.9, 21857.2 };
1835 
1836  DEBUG_ENTRY( "Fe3Lev14()" );
1837 
1838  if( lgFirst )
1839  {
1840  /* will never do this again */
1841  lgFirst = false;
1842  /* allocate the 1D arrays*/
1843  /* create space for the 2D arrays */
1844  depart = ((double *)MALLOC((NLFE3)*sizeof(double)));
1845  pops = ((double *)MALLOC((NLFE3)*sizeof(double)));
1846  destroy = ((double *)MALLOC((NLFE3)*sizeof(double)));
1847  create = ((double *)MALLOC((NLFE3)*sizeof(double)));
1848  /* now the 2-d arrays */
1849  fe.Fe3_wl = ((double **)MALLOC((NLFE3)*sizeof(double *)));
1850  fe.Fe3_emiss = ((double **)MALLOC((NLFE3)*sizeof(double *)));
1851  AulNet = ((double **)MALLOC((NLFE3)*sizeof(double *)));
1852  col_str = ((double **)MALLOC((NLFE3)*sizeof(double *)));
1853  AulPump = ((double **)MALLOC((NLFE3)*sizeof(double *)));
1854  CollRate = ((double **)MALLOC((NLFE3)*sizeof(double *)));
1855  AulDest = ((double **)MALLOC((NLFE3)*sizeof(double *)));
1856  for( i=0; i<(NLFE3); ++i )
1857  {
1858  fe.Fe3_wl[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
1859  fe.Fe3_emiss[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
1860  AulNet[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
1861  col_str[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
1862  AulPump[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
1863  CollRate[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
1864  AulDest[i] = ((double *)MALLOC((NLFE3)*sizeof(double )));
1865  }
1866 
1867  /* set some to constant values after zeroing out */
1868  for( i=0; i<NLFE3; ++i )
1869  {
1870  create[i] = 0.;
1871  destroy[i] = 0.;
1872  for( j=0; j<NLFE3; ++j )
1873  {
1874  AulNet[i][j] = 0.;
1875  col_str[i][j] = 0.;
1876  CollRate[i][j] = 0.;
1877  AulDest[i][j] = 0.;
1878  AulPump[i][j] = 0.;
1879  fe.Fe3_wl[i][j] = 0.;
1880  fe.Fe3_emiss[i][j] = 0.;
1881  }
1882  }
1883  /* calculates wavelengths of transitions */
1884  /* dividing by RefIndex converts the vacuum wavelength to air wavelength */
1885  for( ihi=1; ihi<NLFE3; ++ihi )
1886  {
1887  for( ilo=0; ilo<ihi; ++ilo )
1888  {
1889  fe.Fe3_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) /
1890  RefIndex( (excit_wn[ihi]-excit_wn[ilo]) );
1891  }
1892  }
1893 
1894  /* assume FeIII is optically thin - just use As as net escape */
1895  /*>>refer Fe3 as Quinet, P., 1996, A&AS, 116, 573 */
1896  AulNet[1][0] = 2.8e-3;
1897  AulNet[7][0] = 4.9e-6;
1898  AulNet[8][0] = 5.7e-3;
1899  AulNet[11][0] = 4.5e-1;
1900  AulNet[12][0] = 4.2e-2;
1901 
1902  AulNet[2][1] = 1.8e-3;
1903  AulNet[5][1] = 4.2e-1;
1904  AulNet[8][1] = 1.0e-3;
1905  AulNet[11][1] = 8.4e-2;
1906  AulNet[12][1] = 2.5e-1;
1907  AulNet[13][1] = 2.7e-2;
1908 
1909  AulNet[3][2] = 7.0e-4;
1910  AulNet[5][2] = 5.1e-5;
1911  AulNet[9][2] = 5.4e-1;
1912  AulNet[12][2] = 8.5e-2;
1913  AulNet[13][2] = 9.8e-2;
1914 
1915  AulNet[4][3] = 1.4e-4;
1916  AulNet[5][3] = 3.9e-2;
1917  AulNet[9][3] = 4.1e-5;
1918  AulNet[10][3] = 7.0e-1;
1919  AulNet[13][3] = 4.7e-2;
1920 
1921  AulNet[9][4] = 9.3e-2;
1922 
1923  AulNet[9][5] = 4.7e-2;
1924  AulNet[12][5] = 2.5e-6;
1925  AulNet[13][5] = 1.7e-5;
1926 
1927  AulNet[7][6] = 2.7e-4;
1928 
1929  AulNet[8][7] = 1.2e-4;
1930  AulNet[11][7] = 6.6e-4;
1931 
1932  AulNet[11][8] = 1.6e-3;
1933  AulNet[12][8] = 7.8e-4;
1934 
1935  AulNet[10][9] = 8.4e-3;
1936  AulNet[13][9] = 2.8e-7;
1937 
1938  AulNet[12][11] = 3.0e-4;
1939 
1940  AulNet[13][12] = 1.4e-4;
1941 
1942  /* collision strengths at log T 4 */
1944  /*>>refer Fe3 cs Zhang, H. 1996, 119, 523 */
1945  col_str[1][0] = 2.92;
1946  col_str[2][0] = 1.24;
1947  col_str[3][0] = 0.595;
1948  col_str[4][0] = 0.180;
1949  col_str[5][0] = 0.580;
1950  col_str[6][0] = 1.34;
1951  col_str[7][0] = 0.489;
1952  col_str[8][0] = 0.0926;
1953  col_str[9][0] = 0.165;
1954  col_str[10][0] = 0.0213;
1955  col_str[11][0] = 1.07;
1956  col_str[12][0] = 0.435;
1957  col_str[13][0] = 0.157;
1958 
1959  col_str[2][1] = 2.06;
1960  col_str[3][1] = 0.799;
1961  col_str[4][1] = 0.225;
1962  col_str[5][1] = 0.335;
1963  col_str[6][1] = 0.555;
1964  col_str[7][1] = 0.609;
1965  col_str[8][1] = 0.367;
1966  col_str[9][1] = 0.195;
1967  col_str[10][1] = 0.0698;
1968  col_str[11][1] = 0.538;
1969  col_str[12][1] = 0.484;
1970  col_str[13][1] = 0.285;
1971 
1972  col_str[3][2] = 1.29;
1973  col_str[4][2] = 0.312;
1974  col_str[5][2] = 0.173;
1975  col_str[6][2] = 0.178;
1976  col_str[7][2] = 0.430;
1977  col_str[8][2] = 0.486;
1978  col_str[9][2] = 0.179;
1979  col_str[10][2] = 0.0741;
1980  col_str[11][2] = 0.249;
1981  col_str[12][2] = 0.362;
1982  col_str[13][2] = 0.324;
1983 
1984  col_str[4][3] = 0.493;
1985  col_str[5][3] = 0.0767;
1986  col_str[6][3] = 0.0348;
1987  col_str[7][3] = 0.223;
1988  col_str[8][3] = 0.401;
1989  col_str[9][3] = 0.126;
1990  col_str[10][3] = 0.0528;
1991  col_str[11][3] = 0.101;
1992  col_str[12][3] = 0.207;
1993  col_str[13][3] = 0.253;
1994 
1995  col_str[5][4] = 0.0211;
1996  col_str[6][4] = 0.00122;
1997  col_str[7][4] = 0.0653;
1998  col_str[8][4] = 0.154;
1999  col_str[9][4] = 0.0453;
2000  col_str[10][4] = 0.0189;
2001  col_str[11][4] = 0.0265;
2002  col_str[12][4] = 0.0654;
2003  col_str[13][4] = 0.0950;
2004 
2005  col_str[6][5] = 0.403;
2006  col_str[7][5] = 0.213;
2007  col_str[8][5] = 0.0939;
2008  col_str[9][5] = 1.10;
2009  col_str[10][5] = 0.282;
2010  col_str[11][5] = 0.942;
2011  col_str[12][5] = 0.768;
2012  col_str[13][5] = 0.579;
2013 
2014  col_str[7][6] = 2.84; /* 10-9 */
2015  col_str[8][6] = 0.379; /* 11-9 */
2016  col_str[9][6] = 0.0876; /* 7-9 */
2017  col_str[10][6] = 0.00807; /* 8-9 */
2018  col_str[11][6] = 1.85; /* 12-9 */
2019  col_str[12][6] = 0.667; /* 13-9 */
2020  col_str[13][6] = 0.0905; /* 14-9 */
2021 
2022  col_str[8][7] = 3.07; /* 11-10 */
2023  col_str[9][7] = 0.167; /* 7-10 */
2024  col_str[10][7] = 0.0526; /* 8-10 */
2025  col_str[11][7] = 0.814; /* 12-10 */
2026  col_str[12][7] = 0.837; /* 13-10 */
2027  col_str[13][7] = 0.626; /* 14-10 */
2028 
2029  col_str[9][8] = 0.181; /* 7-11 */
2030  col_str[10][8] = 0.0854; /* 8-11 */
2031  col_str[11][8] = 0.180; /* 12-11 */
2032  col_str[12][8] = 0.778; /* 13-11 */
2033  col_str[13][8] = 0.941; /* 14-11 */
2034 
2035  col_str[10][9] = 0.377; /* 8-7 */
2036  col_str[11][9] = 0.603; /* 12-7 */
2037  col_str[12][9] = 0.472; /* 13-7 */
2038  col_str[13][9] = 0.302; /* 14-7 */
2039 
2040  col_str[11][10] = 0.216; /* 12-8 */
2041  col_str[12][10] = 0.137; /* 13-8 */
2042  col_str[13][10] = 0.106; /* 14-8 */
2043 
2044  col_str[12][11] = 1.25;
2045  col_str[13][11] = 0.292;
2046 
2047  col_str[13][12] = 1.10;
2048  }
2049 
2050  /* bail if no ions */
2051  if( dense.xIonDense[ipIRON][2] <= 0. )
2052  {
2053  CoolAdd("Fe 3",0,0.);
2054 
2055  fe.Fe3CoolTot = 0.;
2056  for( ihi=1; ihi<NLFE3; ++ihi )
2057  {
2058  for( ilo=0; ilo<ihi; ++ilo )
2059  {
2060  fe.Fe3_emiss[ihi][ilo] = 0.;
2061  }
2062  }
2063  /* level populations */
2064  /* LIMLEVELN is the dimension of the atoms vectors */
2065  ASSERT( NLFE3 <= LIMLEVELN);
2066  for( i=0; i < NLFE3; i++ )
2067  {
2068  atoms.PopLevels[i] = 0.;
2069  atoms.DepLTELevels[i] = 1.;
2070  }
2071  return;
2072  }
2073 
2074  /* nNegPop positive if negative pops occurred, negative if too cold */
2075  atom_levelN(
2076  /* number of levels */
2077  NLFE3,
2078  /* the abundance of the ion, cm-3 */
2079  dense.xIonDense[ipIRON][2],
2080  /* the statistical weights */
2081  gfe3,
2082  /* the excitation energies */
2083  excit_wn,
2084  'w',
2085  /* the derived populations - cm-3 */
2086  pops,
2087  /* the derived departure coefficients */
2088  depart,
2089  /* the net emission rate, Aul * escape prob */
2090  &AulNet ,
2091  /* the collision strengths */
2092  &col_str ,
2093  /* A * destruction prob */
2094  &AulDest,
2095  /* pumping rate */
2096  &AulPump,
2097  /* collision rate, s-1, must defined if no collision strengths */
2098  &CollRate,
2099  /* creation vector */
2100  create,
2101  /* destruction vector */
2102  destroy,
2103  /* say atom_levelN should evaluate coll rates from cs */
2104  false,
2105  &fe.Fe3CoolTot,
2106  &dfe3dt,
2107  "Fe 3",
2108  &nNegPop,
2109  &lgZeroPop,
2110  false );
2111 
2112  /* LIMLEVELN is the dim of the PopLevels vector */
2113  ASSERT( NLFE3 <= LIMLEVELN );
2114  for( i=0; i<NLFE3; ++i)
2115  {
2116  atoms.PopLevels[i] = pops[i];
2117  atoms.DepLTELevels[i] = depart[i];
2118  }
2119 
2120  if( nNegPop>0 )
2121  {
2122  fprintf( ioQQQ, " Fe3Lev14 found negative populations\n" );
2123  ShowMe();
2124  cdEXIT(EXIT_FAILURE);
2125  }
2126 
2127  /* add cooling then its derivative */
2128  CoolAdd("Fe 3",0,fe.Fe3CoolTot);
2129  /* derivative of cooling */
2130  thermal.dCooldT += dfe3dt;
2131 
2132  /* find emission in each line */
2133  for( ihi=1; ihi<NLFE3; ++ihi )
2134  {
2135  for( ilo=0; ilo<ihi; ++ilo )
2136  {
2137  /* emission in these lines */
2138  fe.Fe3_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN;
2139  }
2140  }
2141  return;
2142 }
2143 
2144 /*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */
2145 STATIC void Fe11Lev5(void)
2146 {
2147  bool lgZeroPop;
2148  int nNegPop;
2149  long int i,
2150  j;
2151  static bool lgFirst=true;
2152 
2153  double dCool_dT;
2154 
2155  long int ihi , ilo;
2156  static double
2157  *depart,
2158  *pops,
2159  *destroy,
2160  *create ,
2161  **AulDest,
2162  **CollRate,
2163  **AulPump,
2164  **AulNet,
2165  **col_str ,
2166  TeUsed=-1.;
2167 
2168  /* statistical weights for the n levels */
2169  static double stat_wght[NLFE11]={5.,3.,1.,5.,1.};
2170 
2171  /*refer Fe11 energies NIST version 3 Atomic Spectra Database */
2172  /* from smallest to largest energy (not in multiplet groupings) */
2173 
2174  /* energy in wavenumbers */
2175  static double excit_wn[NLFE11]={
2176  0.0 , 12667.9 , 14312. , 37743.6 , 80814.7 };
2177 
2178  DEBUG_ENTRY( "Fe11Lev5()" );
2179 
2180  if( lgFirst )
2181  {
2182  /* will never do this again */
2183  lgFirst = false;
2184  /* allocate the 1D arrays*/
2185  /* create space for the 2D arrays */
2186  depart = ((double *)MALLOC((NLFE11)*sizeof(double)));
2187  pops = ((double *)MALLOC((NLFE11)*sizeof(double)));
2188  destroy = ((double *)MALLOC((NLFE11)*sizeof(double)));
2189  create = ((double *)MALLOC((NLFE11)*sizeof(double)));
2190  /* now the 2-d arrays */
2191  fe.Fe11_wl = ((double **)MALLOC((NLFE11)*sizeof(double *)));
2192  fe.Fe11_emiss = ((double **)MALLOC((NLFE11)*sizeof(double *)));
2193  AulNet = ((double **)MALLOC((NLFE11)*sizeof(double *)));
2194  col_str = ((double **)MALLOC((NLFE11)*sizeof(double *)));
2195  AulPump = ((double **)MALLOC((NLFE11)*sizeof(double *)));
2196  CollRate = ((double **)MALLOC((NLFE11)*sizeof(double *)));
2197  AulDest = ((double **)MALLOC((NLFE11)*sizeof(double *)));
2198  for( i=0; i<(NLFE11); ++i )
2199  {
2200  fe.Fe11_wl[i] = ((double *)MALLOC((NLFE11)*sizeof(double )));
2201  fe.Fe11_emiss[i] = ((double *)MALLOC((NLFE11)*sizeof(double )));
2202  AulNet[i] = ((double *)MALLOC((NLFE11)*sizeof(double )));
2203  col_str[i] = ((double *)MALLOC((NLFE11)*sizeof(double )));
2204  AulPump[i] = ((double *)MALLOC((NLFE11)*sizeof(double )));
2205  CollRate[i] = ((double *)MALLOC((NLFE11)*sizeof(double )));
2206  AulDest[i] = ((double *)MALLOC((NLFE11)*sizeof(double )));
2207  }
2208 
2209  /* set some to constant values after zeroing out */
2210  for( i=0; i<NLFE11; ++i )
2211  {
2212  create[i] = 0.;
2213  destroy[i] = 0.;
2214  for( j=0; j<NLFE11; ++j )
2215  {
2216  AulNet[i][j] = 0.;
2217  col_str[i][j] = 0.;
2218  CollRate[i][j] = 0.;
2219  AulDest[i][j] = 0.;
2220  AulPump[i][j] = 0.;
2221  fe.Fe11_wl[i][j] = 0.;
2222  fe.Fe11_emiss[i][j] = 0.;
2223  }
2224  }
2225  /* calculates wavelengths of transitions */
2226  /* dividing by RefIndex converts the vacuum wavelength to air wavelength */
2227  for( ihi=1; ihi<NLFE11; ++ihi )
2228  {
2229  for( ilo=0; ilo<ihi; ++ilo )
2230  {
2231  fe.Fe11_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) /
2232  RefIndex( (excit_wn[ihi]-excit_wn[ilo]) );
2233  }
2234  }
2235 
2236  /*>>refer Fe11 as Mendoza C. & Zeippen, C. J. 1983, MNRAS, 202, 981 */
2237  AulNet[4][0] = 1.7;
2238  AulNet[4][1] = 990.;
2239  AulNet[4][3] = 8.3;
2240  AulNet[3][0] = 92.3;
2241  AulNet[3][1] = 9.44;
2242  AulNet[3][2] = 1.4e-3;
2243  AulNet[2][0] = 9.9e-3;
2244  AulNet[1][0] = 43.7;
2245  AulNet[2][1] = 0.226;
2246 
2247  }
2248 
2249  /* bail if no ions */
2250  if( dense.xIonDense[ipIRON][10] <= 0. )
2251  {
2252  CoolAdd("Fe11",0,0.);
2253 
2254  fe.Fe11CoolTot = 0.;
2255  for( ihi=1; ihi<NLFE11; ++ihi )
2256  {
2257  for( ilo=0; ilo<ihi; ++ilo )
2258  {
2259  fe.Fe11_emiss[ihi][ilo] = 0.;
2260  }
2261  }
2262  /* level populations */
2263  /* LIMLEVELN is the dimension of the atoms vectors */
2264  ASSERT( NLFE11 <= LIMLEVELN);
2265  for( i=0; i < NLFE11; i++ )
2266  {
2267  atoms.PopLevels[i] = 0.;
2268  atoms.DepLTELevels[i] = 1.;
2269  }
2270  return;
2271  }
2272 
2273  /* evaluate collision strengths if Te changed by too much */
2274  if( fabs(phycon.te / TeUsed - 1. ) > 0.05 )
2275  {
2276  TeUsed = phycon.te;
2277 
2278  /* collision strengths at current temperature */
2279  for( ihi=1; ihi<NLFE11; ++ihi )
2280  {
2281  for( ilo=0; ilo<ihi; ++ilo )
2282  {
2283  col_str[ihi][ilo] = Fe_10_11_13_cs( 11 , ilo+1 , ihi+1 );
2284  ASSERT( col_str[ihi][ilo] > 0. );
2285  }
2286  }
2287  }
2288 
2289  /* nNegPop positive if negative pops occurred, negative if too cold */
2290  atom_levelN(
2291  /* number of levels */
2292  NLFE11,
2293  /* the abundance of the ion, cm-3 */
2294  dense.xIonDense[ipIRON][10],
2295  /* the statistical weights */
2296  stat_wght,
2297  /* the excitation energies */
2298  excit_wn,
2299  'w',
2300  /* the derived populations - cm-3 */
2301  pops,
2302  /* the derived departure coefficients */
2303  depart,
2304  /* the net emission rate, Aul * escp prob */
2305  &AulNet ,
2306  /* the collision strengths */
2307  &col_str ,
2308  /* A * destruction prob */
2309  &AulDest,
2310  /* pumping rate */
2311  &AulPump,
2312  /* collision rate, s-1, must defined if no collision strengths */
2313  &CollRate,
2314  /* creation vector */
2315  create,
2316  /* destruction vector */
2317  destroy,
2318  /* say atom_levelN should evaluate coll rates from cs */
2319  false,
2320  &fe.Fe11CoolTot,
2321  &dCool_dT,
2322  "Fe11",
2323  &nNegPop,
2324  &lgZeroPop,
2325  false );
2326 
2327  /* LIMLEVELN is the dim of the PopLevels vector */
2328  ASSERT( NLFE11 <= LIMLEVELN );
2329  for( i=0; i<NLFE11; ++i)
2330  {
2331  atoms.PopLevels[i] = pops[i];
2332  atoms.DepLTELevels[i] = depart[i];
2333  }
2334 
2335  if( nNegPop>0 )
2336  {
2337  fprintf( ioQQQ, " Fe11Lev5 found negative populations\n" );
2338  ShowMe();
2339  cdEXIT(EXIT_FAILURE);
2340  }
2341 
2342  /* add cooling then its derivative */
2343  CoolAdd("Fe11",0,fe.Fe11CoolTot);
2344  /* derivative of cooling */
2345  thermal.dCooldT += dCool_dT;
2346 
2347  /* find emission in each line */
2348  for( ihi=1; ihi<NLFE11; ++ihi )
2349  {
2350  for( ilo=0; ilo<ihi; ++ilo )
2351  {
2352  /* emission in these lines */
2353  fe.Fe11_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*
2354  (excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN;
2355  }
2356  }
2357  return;
2358 }
2359 
2360 /*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */
2361 STATIC void Fe13Lev5(void)
2362 {
2363  bool lgZeroPop;
2364  int nNegPop;
2365  long int i,
2366  j;
2367  static bool lgFirst=true;
2368 
2369  double dCool_dT;
2370 
2371  long int ihi , ilo;
2372  static double
2373  *depart,
2374  *pops,
2375  *destroy,
2376  *create ,
2377  **AulDest,
2378  **CollRate,
2379  **AulPump,
2380  **AulNet,
2381  **col_str ,
2382  TeUsed=-1.;
2383 
2384  /* statistical weights for the n levels */
2385  static double stat_wght[NLFE13]={1.,3.,5.,5.,1.};
2386 
2387  /*refer Fe13 energies NIST version 3 Atomic Spectra Database */
2388  /* energy in wavenumbers */
2389  static double excit_wn[NLFE13]={
2390  0.0 , 9302.5 , 18561.0 , 48068. , 91508. };
2391 
2392  DEBUG_ENTRY( "Fe13Lev5()" );
2393 
2394  if( lgFirst )
2395  {
2396  /* will never do this again */
2397  lgFirst = false;
2398  /* allocate the 1D arrays*/
2399  /* create space for the 2D arrays */
2400  depart = ((double *)MALLOC((NLFE13)*sizeof(double)));
2401  pops = ((double *)MALLOC((NLFE13)*sizeof(double)));
2402  destroy = ((double *)MALLOC((NLFE13)*sizeof(double)));
2403  create = ((double *)MALLOC((NLFE13)*sizeof(double)));
2404  /* now the 2-d arrays */
2405  fe.Fe13_wl = ((double **)MALLOC((NLFE13)*sizeof(double *)));
2406  fe.Fe13_emiss = ((double **)MALLOC((NLFE13)*sizeof(double *)));
2407  AulNet = ((double **)MALLOC((NLFE13)*sizeof(double *)));
2408  col_str = ((double **)MALLOC((NLFE13)*sizeof(double *)));
2409  AulPump = ((double **)MALLOC((NLFE13)*sizeof(double *)));
2410  CollRate = ((double **)MALLOC((NLFE13)*sizeof(double *)));
2411  AulDest = ((double **)MALLOC((NLFE13)*sizeof(double *)));
2412  for( i=0; i<(NLFE13); ++i )
2413  {
2414  fe.Fe13_wl[i] = ((double *)MALLOC((NLFE13)*sizeof(double )));
2415  fe.Fe13_emiss[i] = ((double *)MALLOC((NLFE13)*sizeof(double )));
2416  AulNet[i] = ((double *)MALLOC((NLFE13)*sizeof(double )));
2417  col_str[i] = ((double *)MALLOC((NLFE13)*sizeof(double )));
2418  AulPump[i] = ((double *)MALLOC((NLFE13)*sizeof(double )));
2419  CollRate[i] = ((double *)MALLOC((NLFE13)*sizeof(double )));
2420  AulDest[i] = ((double *)MALLOC((NLFE13)*sizeof(double )));
2421  }
2422 
2423  /* set some to constant values after zeroing out */
2424  for( i=0; i<NLFE13; ++i )
2425  {
2426  create[i] = 0.;
2427  destroy[i] = 0.;
2428  for( j=0; j<NLFE13; ++j )
2429  {
2430  AulNet[i][j] = 0.;
2431  col_str[i][j] = 0.;
2432  CollRate[i][j] = 0.;
2433  AulDest[i][j] = 0.;
2434  AulPump[i][j] = 0.;
2435  fe.Fe13_wl[i][j] = 0.;
2436  fe.Fe13_emiss[i][j] = 0.;
2437  }
2438  }
2439  /* calculates wavelengths of transitions */
2440  /* dividing by RefIndex converts the vacuum wavelength to air wavelength */
2441  for( ihi=1; ihi<NLFE13; ++ihi )
2442  {
2443  for( ilo=0; ilo<ihi; ++ilo )
2444  {
2445  fe.Fe13_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) /
2446  RefIndex( (excit_wn[ihi]-excit_wn[ilo]) );
2447  }
2448  }
2449 
2450  /*>>refer Fe13 as Huang, K.-N. 1985, At. Data Nucl. Data Tables 32, 503 */
2451  AulNet[4][1] = 1010.;
2452  AulNet[4][2] = 3.8;
2453  AulNet[4][3] = 8.1;
2454  AulNet[3][1] = 45.7;
2455  AulNet[3][2] = 57.6;
2456  AulNet[2][0] = 0.0063;
2457  AulNet[1][0] = 14.0;
2458  AulNet[2][1] = 9.87;
2459 
2460  }
2461 
2462  /* bail if no ions */
2463  if( dense.xIonDense[ipIRON][12] <= 0. )
2464  {
2465  CoolAdd("Fe13",0,0.);
2466 
2467  fe.Fe13CoolTot = 0.;
2468  for( ihi=1; ihi<NLFE13; ++ihi )
2469  {
2470  for( ilo=0; ilo<ihi; ++ilo )
2471  {
2472  fe.Fe13_emiss[ihi][ilo] = 0.;
2473  }
2474  }
2475  /* level populations */
2476  /* LIMLEVELN is the dimension of the atoms vectors */
2477  ASSERT( NLFE13 <= LIMLEVELN);
2478  for( i=0; i < NLFE13; i++ )
2479  {
2480  atoms.PopLevels[i] = 0.;
2481  atoms.DepLTELevels[i] = 1.;
2482  }
2483  return;
2484  }
2485 
2486  /* evaluate collision strengths if Te changed by too much */
2487  if( fabs(phycon.te / TeUsed - 1. ) > 0.05 )
2488  {
2489  TeUsed = phycon.te;
2490 
2491  /* collision strengths at current temperature */
2492  for( ihi=1; ihi<NLFE13; ++ihi )
2493  {
2494  for( ilo=0; ilo<ihi; ++ilo )
2495  {
2496  col_str[ihi][ilo] = Fe_10_11_13_cs( 13 , ilo+1 , ihi+1 );
2497  ASSERT( col_str[ihi][ilo] > 0. );
2498  }
2499  }
2500  }
2501 
2502  /*fprintf(ioQQQ,"DEBUG Fe13 N %.2e %.2e %.2e %.2e %.2e %.2e\n",
2503  col_str[1][0] , col_str[2][1] , col_str[2][0] ,
2504  AulNet[1][0] , AulNet[2][1] , AulNet[2][0]);*/
2505 
2506  /* nNegPop positive if negative pops occurred, negative if too cold */
2507  atom_levelN(
2508  /* number of levels */
2509  NLFE13,
2510  /* the abundance of the ion, cm-3 */
2511  dense.xIonDense[ipIRON][12],
2512  /* the statistical weights */
2513  stat_wght,
2514  /* the excitation energies */
2515  excit_wn,
2516  'w',
2517  /* the derived populations - cm-3 */
2518  pops,
2519  /* the derived departure coefficients */
2520  depart,
2521  /* the net emission rate, Aul * escape prob */
2522  &AulNet ,
2523  /* the collision strengths */
2524  &col_str ,
2525  /* A * destruction prob */
2526  &AulDest,
2527  /* pumping rate */
2528  &AulPump,
2529  /* collision rate, s-1, must defined if no collision strengths */
2530  &CollRate,
2531  /* creation vector */
2532  create,
2533  /* destruction vector */
2534  destroy,
2535  /* say atom_levelN should evaluate coll rates from cs */
2536  false,
2537  &fe.Fe13CoolTot,
2538  &dCool_dT,
2539  "Fe13",
2540  &nNegPop,
2541  &lgZeroPop,
2542  false );
2543 
2544  /* LIMLEVELN is the dim of the PopLevels vector */
2545  ASSERT( NLFE13 <= LIMLEVELN );
2546  for( i=0; i<NLFE13; ++i)
2547  {
2548  atoms.PopLevels[i] = pops[i];
2549  atoms.DepLTELevels[i] = depart[i];
2550  }
2551 
2552  if( nNegPop>0 )
2553  {
2554  fprintf( ioQQQ, " Fe13Lev5 found negative populations\n" );
2555  ShowMe();
2556  cdEXIT(EXIT_FAILURE);
2557  }
2558 
2559  /* add cooling then its derivative */
2560  CoolAdd("Fe13",0,fe.Fe13CoolTot);
2561  /* derivative of cooling */
2562  thermal.dCooldT += dCool_dT;
2563 
2564  /* find emission in each line */
2565  for( ihi=1; ihi<NLFE13; ++ihi )
2566  {
2567  for( ilo=0; ilo<ihi; ++ilo )
2568  {
2569  /* emission in these lines */
2570  fe.Fe13_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN;
2571  }
2572  }
2573  return;
2574 }

Generated for cloudy by doxygen 1.8.3.1