cloudy  trunk
atom_feii.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 /*FeII_Colden maintain H2 column densities within X */
4 /*FeIILevelPops main FeII routine, called by CoolIron to evaluate iron cooling */
5 /*FeIICreate read in needed data from file
6  * convert form of FeII data, as read in from file within routine FeIICreate
7  * into physical form. called by atmdat_readin */
8 /*FeIIPunPop - save level populations */
9 /*AssertFeIIDep called by assert FeII depart coef command */
10 /*FeIIPrint print FeII information */
11 /*FeIICollRatesBoltzmann evaluate collision strenths,
12  * both interpolating on r-mat and creating g-bar
13  * find Boltzmann factors, evaluate collisional rate coefficients */
14 /*FeIIPrint print output from large FeII atom, called by prtzone */
15 /*FeIISumBand sum up large FeII emission over certain bands, called in lineset4 */
16 /*FeII_RT_TauInc called once per zone in RT_tau_inc to increment large FeII atom line optical depths */
17 /*FeII_RT_tau_reset reset optical depths for large FeII atom, called by update after each iteration */
18 /*FeIIPoint called by ContCreatePointers to create pointers for lines in large FeII atom */
19 /*FeIIAccel called by rt_line_driving to compute radiative acceleration due to FeII lines */
20 /*FeIIAddLines save accumulated FeII intensities called by lineset4 */
21 /*FeIISaveLines save FeII lines at end of calculation, if save verner set, called by punch_do*/
22 /*FeIIPunchOpticalDepth save FeII line optical depth at end of calculation, called by punch_do*/
23 /*FeIIPunchLevels save FeII levels and energies */
24 /*FeII_LineZero zero out storage for large FeII atom, called by tauout */
25 /*FeIIIntenZero zero out intensity of FeII atom */
26 /*FeIIZero initialize some variables, called by zero one time before commands parsed */
27 /*FeIIReset reset some variables, called by zero */
28 /*FeIIPunData save line data */
29 /*FeIIPunDepart save some departure coef for large atom,
30  * set with save FeII departure command*/
31 /*FeIIPun1Depart send the departure coef for physical level nPUN to unit ioPUN */
32 /*FeIIContCreate create FeII continuum bins to add lines into ncell cells between wavelengths lambda low and high,
33  * returns number of cells used */
34 /*FeIIBandsCreate returns number of FeII bands */
35 /*FeII_RT_Out do outward rates for FeII, called by RT_diffuse */
36 /*FeII_OTS do ots rates for FeII, called by RT_OTS */
37 /*FeII_RT_Make called by RT_line_all, does large FeII atom radiative transfer */
38 /*FeIILyaPump find rate of Lya excitation of the FeII atom */
39 /*ParseAtomFeII parse the atom FeII command */
40 /*FeIIPunchLineStuff include FeII lines in punched optical depths, etc, called from SaveLineStuff */
41 #include "cddefines.h"
42 #include "cddrive.h"
43 #include "thermal.h"
44 #include "physconst.h"
45 #include "doppvel.h"
46 #include "taulines.h"
47 #include "dense.h"
48 #include "rfield.h"
49 #include "radius.h"
50 #include "lines_service.h"
51 #include "ipoint.h"
52 #include "thirdparty.h"
53 #include "hydrogenic.h"
54 #include "lines.h"
55 #include "rt.h"
56 #include "trace.h"
57 #include "save.h"
58 #include "phycon.h"
59 #include "atomfeii.h"
60 #include "iso.h"
61 #include "pressure.h"
62 #include "parser.h"
63 
64 /* add FeII lines into ncell cells between wavelengths lambda low and high */
65 STATIC void FeIIContCreate(double xLamLow , double xLamHigh , long int ncell );
66 
67 /* read in the FeII Bands file, and sets nFeIIBands, the number of bands,
68  * if argument is "" then use default file with bands,
69  * if filename within "" then use it instead,
70  * return value is 0 if success, 1 if failed */
71 STATIC int FeIIBandsCreate( const char chFile[] );
72 
73 /* this will be the smallest collision strength we will permit with the gbar
74  * collision strengths, or for the data that have zeroes */
75 /* >>chng 99 jul 24, this was 1e-9 in the old fortran version and in baldwin et al. 96,
76  * and has been changed several times, and affects results. this is the smallest
77  * non-zero collision strength in the r-matrix calculations */
78 realnum CS2SMALL = (realnum)1e-5;
79 /* routines used only within this file */
80 
81 /*FeIICollRatesBoltzmann evaluate collision strenths,
82  * both interpolating on r-mat and creating g-bar
83  * find Boltzmann factors, evaluate collisional rate coefficients */
84 STATIC void FeIICollRatesBoltzmann(void);
85 /* find rate of Lya excitation of the FeII atom */
86 STATIC void FeIILyaPump(void);
87 
88 /*extern realnum Fe2LevN[NFE2LEVN][NFE2LEVN][NTA];*/
89 /*extern realnum Fe2LevN[ipHi][ipLo].t[NTA];*/
90 /*realnum ***Fe2LevN;*/
91 /* >>chng 06 mar 01, boost to global namespace */
92 /*transition **Fe2LevN;*/
93 /* flag for the collision strength */
94 int **ncs1;
95 
96 /* all following variables have file scope
97 #define NFEIILINES 68635 */
98 
99 /* this is size of nnPradDat array */
100 #define NPRADDAT 159
101 
102 /* band wavelength, lower and upper bounds, in vacuum Angstroms */
103 /* FeII_Bands[n][3], where n is the number of bands in fe2Bands.dat
104  * these bands are defined in fe2Bands.dat and read in at startup
105  * of calculation */
106 realnum **FeII_Bands;
107 
108 // FeII_Cont[n][0] gives lower and upper bounds continuum wavelengths
109 // in vacuum Angstroms,
110 // [1] and [2] are inward and outward integrated intensity
111 // n is the number of bands, one less that number of cells needed
112 // these bands are defined in FeIIContCreate
113 realnum **FeII_Cont;
114 
115 /* this is the number of bands read in from FeII_bands.ini */
116 long int nFeIIBands;
117 
118 /* number of bands in continuum array */
119 long int nFeIIConBins;
120 
121 /* the dim of this vector this needs one extra since there are 20 numbers per line,
122  * highest not ever used for anything */
123 /*long int nnPradDat[NPRADDAT+1];*/
124 static long int *nnPradDat;
125 
126 /* malloced in feiidata */
127 /* realnum sPradDat[NPRADDAT][NPRADDAT][8];*/
128 /* realnum sPradDat[ipHi][ipLo].t[8];*/
129 static realnum ***sPradDat;
130 
131 /* array used to integrate FeII line intensities over model,
132  * Fe2SavN[upper][lower],
133  *static realnum Fe2SavN[NFE2LEVN][NFE2LEVN];*/
134 static double **Fe2SavN;
135 
136 /* save effective transition rates */
137 static double **Fe2A;
138 
139 /* induced transition rates */
140 static double **Fe2LPump, **Fe2CPump;
141 
142 /* energies read in from fe2energies.dat data file */
143 realnum *Fe2Energies;
144 
145 /* collision rates */
146 static realnum **Fe2Coll;
147 
148 /* Fe2DepCoef[NFE2LEVN];
149 realnum cli[NFEIILINES],
150  cfe[NFEIILINES];*/
151 static double
152  /* departure coefficients */
153  *Fe2DepCoef ,
154  /* level populations */
155  *Fe2LevelPop ,
156  /* column densities */
157  *Fe2ColDen ,
158  /* this will become array of Boltzmann factors */
159  *FeIIBoltzmann;
160  /*FeIIBoltzmann[NFE2LEVN] ,*/
161 
162 static double EnerLyaProf1,
163  EnerLyaProf4,
164  PhotOccNumLyaCenter;
165 static double
166  /* the yVector - will become level populations after matrix inversion */
167  *yVector,
168  /* this is used to call matrix routines */
169  /*xMatrix[NFE2LEVN][NFE2LEVN] ,*/
170  **xMatrix ,
171  /* this will become the very large 1-D array that
172  * is passed to the matrix inversion routine*/
173  *amat;
174 
175 
176 /*FeII_Colden maintain H2 column densities within X */
177 void FeII_Colden( const char *chLabel )
178 {
179  long int n;
180 
181  DEBUG_ENTRY( "FeII_Colden()" );
182 
183  if( strcmp(chLabel,"ZERO") == 0 )
184  {
185  /* zero out column density */
186  for( n=0; n < FeII.nFeIILevel_malloc; ++n )
187  {
188  /* space for the rotation quantum number */
189  Fe2ColDen[n] = 0.;
190  }
191  }
192 
193  else if( strcmp(chLabel,"ADD ") == 0 )
194  {
195  /* add together column densities */
196  for( n=0; n < FeII.nFeIILevel_local; ++n )
197  {
198  /* state-specific FeII column density */
199  Fe2ColDen[n] += Fe2LevelPop[n]*radius.drad_x_fillfac;
200  }
201  }
202 
203  /* check for the print option, which we can't do, else we have a problem */
204  else if( strcmp(chLabel,"PRIN") != 0 )
205  {
206  fprintf( ioQQQ, " FeII_Colden does not understand the label %s\n",
207  chLabel );
208  cdEXIT(EXIT_FAILURE);
209  }
210 
211  return;
212 }
213 
214 /*
215  *=====================================================================
216  */
217 /* FeIICreate read in FeII data from files on disk. called by atmdat_readin
218  * but only if FeII. lgFeIION is true, set with atom FeII verner command */
219 void FeIICreate(void)
220 {
221  FILE *ioDATA;
222 
223  char chLine[FILENAME_PATH_LENGTH_2];
224 
225  long int i,
226  ipHi ,
227  ipLo,
228  lo,
229  ihi,
230  k,
231  m1,
232  m2,
233  m3;
234 
235  DEBUG_ENTRY( "FeIICreate()" );
236 
237  if( lgFeIIMalloc )
238  {
239  /* we have already been called one time, just bail out */
240 
241  return;
242  }
243 
244  /* now set flag so never done again - this is set false when initi
245  * when this is true it is no longer possible to change the number of levels
246  * in the model atom with the atom FeII levels command */
247  lgFeIIMalloc = true;
248 
249  /* remember how many levels this was, so that in future calculations
250  * we can reset the atom to full value */
251  FeII.nFeIILevel_malloc = FeII.nFeIILevel_local;
252 
253  /* set up array to save FeII transition probabilities */
254  Fe2A = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
255 
256  /* second dimension, lower level, for line save array */
257  for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
258  {
259  Fe2A[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel_malloc);
260  }
261 
262  /* set up array to save FeII pumping rates */
263  Fe2CPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
264 
265  /* set up array to save FeII pumping rates */
266  Fe2LPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
267 
268  /* second dimension, lower level, for line save array */
269  for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
270  {
271  Fe2CPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel_malloc);
272 
273  Fe2LPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel_malloc);
274  }
275 
276  /* set up array to save FeII collision rates */
277  Fe2Energies = (realnum *)MALLOC(sizeof(realnum)*(unsigned long)FeII.nFeIILevel_malloc );
278 
279  /* set up array to save FeII collision rates */
280  Fe2Coll = (realnum **)MALLOC(sizeof(realnum *)*(unsigned long)FeII.nFeIILevel_malloc );
281 
282  /* second dimension, lower level, for line save array */
283  for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
284  {
285  Fe2Coll[ipHi]=(realnum*)MALLOC(sizeof(realnum )*(unsigned long)FeII.nFeIILevel_malloc);
286  }
287 
288  /* MALLOC space for the 2D matrix array */
289  xMatrix = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
290 
291  /* now do the second dimension */
292  for( i=0; i<FeII.nFeIILevel_malloc; ++i )
293  {
294  xMatrix[i] = (double *)MALLOC(sizeof(double)*(unsigned long)FeII.nFeIILevel_malloc );
295  }
296  /* MALLOC space for the 1-yVector array */
297  amat=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevel_malloc*FeII.nFeIILevel_malloc) ));
298 
299  /* MALLOC space for the 1-yVector array */
300  yVector=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevel_malloc) ));
301 
302  /* set up array to save FeII line intensities */
303  Fe2SavN = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel_malloc );
304 
305  /* second dimension, lower level, for line save array */
306  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
307  {
308  Fe2SavN[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)ipHi);
309  }
310 
311  /* now MALLOC space for energy level table*/
312  nnPradDat = (long*)MALLOC( (NPRADDAT+1)*sizeof(long) );
313 
314  /*Fe2DepCoef[NFE2LEVN];*/
315  Fe2DepCoef = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
316 
317  /*Fe2LevelPop[NFE2LEVN];*/
318  Fe2LevelPop = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
319 
320  /*Fe2ColDen[NFE2LEVN];*/
321  Fe2ColDen = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
322 
323  /*FeIIBoltzmann[NFE2LEVN];*/
324  FeIIBoltzmann = (double*)MALLOC( (unsigned long)FeII.nFeIILevel_malloc*sizeof(double) );
325 
326 
327  /* MALLOC the realnum sPradDat[NPRADDAT][NPRADDAT][8] array */
328  /* MALLOC the realnum sPradDat[ipHi][ipLo].t[8] array */
329  sPradDat = ((realnum ***)MALLOC(NPRADDAT*sizeof(realnum **)));
330 
331  /* >>chng 00 dec 06, changed lower limit of loop to 1, Tru64 alpha's will not allocate 0 bytes!, PvH */
332  sPradDat[0] = NULL;
333  for(ipHi=1; ipHi < NPRADDAT; ipHi++)
334  {
335  /* >>chng 00 dec 06, changed sizeof(realnum) to sizeof(realnum*), PvH */
336  sPradDat[ipHi] = (realnum **)MALLOC((unsigned long)ipHi*sizeof(realnum *));
337 
338  /* now make space for the second dimension */
339  for( ipLo=0; ipLo< ipHi; ipLo++ )
340  {
341  sPradDat[ipHi][ipLo] = (realnum *)MALLOC(8*sizeof(realnum ));
342  }
343  }
344 
345  /* now set junk to make sure reset before used */
346  for(ipHi=0; ipHi < NPRADDAT; ipHi++)
347  {
348  for( ipLo=0; ipLo< ipHi; ipLo++ )
349  {
350  for( k=0; k<8; ++k )
351  {
352  sPradDat[ipHi][ipLo][k] = -FLT_MAX;
353  }
354  }
355  }
356 
357  /* Set arrays to -2 to know which are used */
358  for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
359  {
360  for( ipLo=0; ipLo < FeII.nFeIILevel_malloc; ipLo++ )
361  {
362  FeII.FeIIAul[ipHi][ipLo] = -2.;
363  for( int k=0; k<8; k++ )
364  {
365  FeII.FeIIColl[ipHi][ipLo][k] = -2.;
366  }
367  }
368  FeII.FeIINRGs[ipHi] = -2.;
369  FeII.FeIISTWT[ipHi] = -2.;
370  }
371 
372  /* now create the main emission line array and a helper array for the cs flag */
373  ipFe2LevN.reserve(FeII.nFeIILevel_malloc);
374  ncs1=(int**)MALLOC(sizeof(int *)*(unsigned long)FeII.nFeIILevel_malloc);
375 
376  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
377  {
378  ipFe2LevN.reserve(ipHi,ipHi);
379  ncs1[ipHi]=(int*)MALLOC(sizeof(int)*(unsigned long)ipHi);
380  }
381 
382  ipFe2LevN.alloc();
383  Fe2LevN.resize(ipFe2LevN.size());
384  AllTransitions.push_back(Fe2LevN);
385 
386  int nLine=0;
387  /* now that its created, set to junk */
388  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
389  {
390  for( ipLo=0; ipLo < ipHi; ipLo++ )
391  {
392  ipFe2LevN[ipHi][ipLo] = nLine;
393  Fe2LevN[nLine].Junk();
394  ++nLine;
395  }
396  }
397 
398  /* now assign state and Emis pointers */
399  Fe2LevN.states()->resize(FeII.nFeIILevel_malloc);
400  /* now that its created, set to junk */
401  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
402  {
403  /* add this upper level */
404  for( ipLo=0; ipLo < ipHi; ipLo++ )
405  {
406  TransitionProxy tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
407  tr.setLo(ipLo);
408  tr.setHi(ipHi);
409  tr.AddLine2Stack();
410  tr.Zero();
411  }
412  }
413 
414  if( trace.lgTrace )
415  {
416  fprintf( ioQQQ," FeIICreate opening fe2nn.dat:");
417  }
418 
419  ioDATA = open_data( "fe2nn.dat", "r" );
420 
421  ASSERT( ioDATA !=NULL );
422  /* read in the fe2nn.dat file - this gives the Zheng and Pradhan number of level
423  * for every cloudy level. So nnPradDat[1] is the index in the cloudy level
424  * counting for level 1 of Zheng & Pradan
425  * note that the order of some levels is different, the nnPradDat file goes
426  * 21 23 22 - also that many levels are missing, the file goes 95 99 94 93 116
427  */
428  for( i=0; i < 8; i++ )
429  {
430  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
431  {
432  fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
433  cdEXIT(EXIT_FAILURE);
434  }
435 
436  /* get these integers from fe2nn.dat */
437  k = 20*i;
438  /* NPRADDAT is size of nnPradDat array, 159+1, make sure we do not exceed it */
439  ASSERT( k+19 < NPRADDAT+1 );
440  sscanf( chLine ,
441  "%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld",
442  &nnPradDat[k+0], &nnPradDat[k+1], &nnPradDat[k+2], &nnPradDat[k+3], &nnPradDat[k+4],
443  &nnPradDat[k+5], &nnPradDat[k+6], &nnPradDat[k+7], &nnPradDat[k+8], &nnPradDat[k+9],
444  &nnPradDat[k+10],&nnPradDat[k+11], &nnPradDat[k+12],&nnPradDat[k+13],&nnPradDat[k+14],
445  &nnPradDat[k+15],&nnPradDat[k+16], &nnPradDat[k+17],&nnPradDat[k+18],&nnPradDat[k+19]
446  );
447 # if !defined(NDEBUG)
448  for( m1=0; m1<20; ++m1 )
449  {
450  ASSERT( nnPradDat[k+m1] >= 0 && nnPradDat[k+m1] <= NFE2LEVN );
451  }
452 # endif
453  }
454  fclose(ioDATA);
455 
456  /* now get energies */
457  if( trace.lgTrace )
458  {
459  fprintf( ioQQQ," FeIICreate opening fe2energies.dat:");
460  }
461 
462  ioDATA = open_data( "fe2energies.dat", "r" );
463 
464  /* file now open, read the data */
465  for( ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
466  {
467  /* keep reading until have non-comment line, one that does not
468  * start with # */
469  chLine[0] = '#';
470  while( chLine[0] == '#' )
471  {
472  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
473  {
474  fprintf( ioQQQ, " fe2energies.dat error reading data\n" );
475  cdEXIT(EXIT_FAILURE);
476  }
477  }
478 
479  /* first and last number on this line */
480  double help;
481  sscanf( chLine, "%lf", &help );
482  Fe2Energies[ipHi] = (realnum)help;
483  FeII.FeIINRGs[ipHi] = Fe2Energies[ipHi];
484  }
485  fclose(ioDATA);
486 
487  /* transition probabilities */
488 
489  if( trace.lgTrace )
490  fprintf( ioQQQ," FeIICreate opening fe2rad.dat:");
491 
492  ioDATA = open_data( "fe2rad.dat", "r" );
493 
494  /* get the first line, this is a version number */
495  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
496  {
497  fprintf( ioQQQ, " fe2rad.dat error reading data\n" );
498  cdEXIT(EXIT_FAILURE);
499  }
500  /* scan off three integers */
501  sscanf( chLine ,"%ld%ld%ld",&lo, &ihi, &m1);
502  const int nYrFe2Rad=8, nMoFe2Rad=8, nDyFe2Rad=24;
503  if( lo!=nYrFe2Rad || ihi!=nMoFe2Rad || m1!=nDyFe2Rad )
504  {
505  fprintf( ioQQQ, "DISASTER fe2rad.dat has the wrong magic numbers, expected "
506  "%2i %2i %2i and got %2li %2li %2li\n" ,
507  nYrFe2Rad, nMoFe2Rad, nDyFe2Rad,
508  lo, ihi, m1);
509  cdEXIT(EXIT_FAILURE);
510  }
511 
512  /* file now open, read the data */
513  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
514  {
515  for( ipLo=0; ipLo < ipHi; ipLo++ )
516  {
517  /* following double since used in sscanf */
518  double save[2];
519  /* keep reading until have non-comment line, one that does not
520  * start with # */
521  chLine[0] = '#';
522  while( chLine[0] == '#' )
523  {
524  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
525  {
526  fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
527  cdEXIT(EXIT_FAILURE);
528  }
529  }
530 
531  /* first and last number on this line */
532  sscanf( chLine ,
533  "%ld%ld%ld%ld%lf%lf%ld",
534  &lo, &ihi, &m1, &m2 ,
535  &save[0], &save[1] , &m3);
536  /* the indices ihi and lo within this array were
537  * read in from the line above */
538  const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ihi-1][lo-1]];
539  (*tr.Lo()).g() = (realnum)m1;
540  FeII.FeIISTWT[lo-1] = (*tr.Lo()).g();
541  (*tr.Hi()).g() = (realnum)m2;
542  FeII.FeIISTWT[ihi-1] = (*tr.Hi()).g();
543  tr.Emis().Aul() = (realnum)save[0];
544  FeII.FeIIAul[ihi-1][lo-1] = tr.Emis().Aul();
545 
546  /*>>chng 06 apr 10, option to use table of energies */
547 # define USE_OLD true
548  if( USE_OLD )
549  tr.EnergyWN() = (realnum)save[1];
550  else
551  tr.EnergyWN() = Fe2Energies[ihi-1]-Fe2Energies[lo-1];
552 
553  /* Aul == -1 indicates intercombination line with no real Aul */
554  if( fp_equal( tr.Emis().Aul() , (realnum)-1.0f ) )
555  {
556  /* these are made-up intercombination lines, set gf to 1e-5 */
557  tr.Emis().gf() = 1e-5f;
558 
559  /* get corresponding A */
560  tr.Emis().Aul() = tr.Emis().gf()*(realnum)TRANS_PROB_CONST*
561  POW2(tr.EnergyWN())*(*tr.Lo()).g()/(*tr.Hi()).g();
562  }
563 
564  /* the last column of fe2rad.dat, and is 1, 2, or 3.
565  * 1 signifies that transition is permitted,
566  * 2 is semi-forbidden
567  * 3 forbidden, within lowest 63 levels are forbidden, first permitted
568  * transition is from 64 */
569  ncs1[ihi-1][lo-1] = (int)m3;
570  /* Use above to determine E1,M1,E2 ...etc */
571  }
572  }
573  fclose(ioDATA);
574 
575  /* now read collision data in fe2col.dat
576  * These are from the following sources
577  >>refer Fe2 CS Zhang, H. L., & Pradhan, A. K. 1995, A&A 293, 953
578  >>refer Fe2 CS Bautista, M., (private communication),
579  >>refer Fe2 CS Mewe, R. 1972, A&A, 20, 215
580  */
581 
582  if( trace.lgTrace )
583  {
584  fprintf( ioQQQ," FeIICreate opening fe2col.dat:");
585  }
586 
587  ioDATA = open_data( "fe2col.dat", "r" );
588 
589  ASSERT( ioDATA !=NULL);
590  for( ipHi=1; ipHi<NPRADDAT; ++ipHi )
591  {
592  for( ipLo=0; ipLo<ipHi; ++ipLo )
593  {
594  /* double since used in sscanf */
595  double save[8];
596  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
597  {
598  fprintf( ioQQQ, " fe2col.dat error reading data\n" );
599  cdEXIT(EXIT_FAILURE);
600  }
601  sscanf( chLine ,
602  "%ld%ld%lf%lf%lf%lf%lf%lf%lf%lf",
603  &m1, &m2,
604  &save[0], &save[1] , &save[2],&save[3], &save[4] , &save[5],
605  &save[6], &save[7]
606  );
607  for( k=0; k<8; ++k )
608  {
609  /* the max is here because there are some zeroes in the data file.
610  * this is unphysical but is part of their distribution. As a result
611  * don't let the cs fall below 0.01 */
612  sPradDat[m2-1][m1-1][k] = max(CS2SMALL , (realnum)save[k] );
613 
614  long ipHiFe2 = MAX2( nnPradDat[m2-1] , nnPradDat[m1-1] );
615  long ipLoFe2 = MIN2( nnPradDat[m2-1] , nnPradDat[m1-1] );
616  FeII.FeIIColl[ipHiFe2-1][ipLoFe2-1][k] = save[k];
617  }
618  }
619  }
620 
621  fclose( ioDATA );
622 
623  /*generate needed opacity data for the large FeII atom */
624 
625  /* this routine is called only one time when cloudy is init
626  * for the very first time. It converts the FeII data stored
627  * in the large FeII arrays into the array storage needed by cloudy
628  * for its line optical depth arrays
629  */
630 
631  /* convert large FeII line arrays into standard heavy el ar */
632  for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
633  {
634  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
635  {
636  /* pull information out of existing data arrays */
637  const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
638  ASSERT( tr.EnergyWN() > 0. );
639  ASSERT( tr.Emis().Aul()> 0. );
640 
641  /* now put into standard internal line format
642  Fe2LevN[ipHi][ipLo].WLAng = (realnum)(1.e8/Fe2LevN[ipHi][ipLo].EnergyWN); */
643  /* >>chng 03 oct 28, above neglected index of refraction of air -
644  * convert to below */
645  tr.WLAng() =
646  (realnum)(1.0e8/
647  tr.EnergyWN()/
648  RefIndex( tr.EnergyWN() ));
649 
650  /* generate gf from A */
651  tr.Emis().gf() =
652  (realnum)(tr.Emis().Aul()*(*tr.Hi()).g()/
653  TRANS_PROB_CONST/POW2(tr.EnergyWN()));
654 
655  (*tr.Hi()).IonStg() = 2;
656  (*tr.Hi()).nelem() = 26;
657  /* which redistribution function??
658  * For resonance line use K2 (-1),
659  * for subordinate lines use CRD with wings */
660  /* >>chng 01 mar 09, all had been 1, inc redis with wings */
661  /* to reset this, so that the code works as it did pre 01 mar 29,
662  * use command
663  * atom FeII redistribution resonance pdr
664  * atom FeII redistribution subordinate pdr */
665  if( ipLo == 0 )
666  {
667  tr.Emis().iRedisFun() = FeII.ipRedisFcnResonance;
668  }
669  else
670  {
671  /* >>chng 01 feb 27, had been -1, crd with core only,
672  * change to crd with wings as per discussion with Ivan Hubeny */
673  tr.Emis().iRedisFun() = FeII.ipRedisFcnSubordinate;
674  }
675  tr.Emis().phots() = 0.;
676  tr.Emis().ots() = 0.;
677  tr.Emis().FracInwd() = 1.;
678  }
679  }
680 
681  /* finally get FeII bands, this sets */
682  if( FeIIBandsCreate("") )
683  {
684  fprintf( ioQQQ," failed to open FeII bands file \n");
685  cdEXIT(EXIT_FAILURE);
686  }
687 
688  /* now establish the FeII continuum, these are set to
689  * 1000, 7000, and 1000 in FeIIZero in this file, and
690  * are reset with the atom FeII continuum command */
691  FeIIContCreate( FeII.feconwlLo , FeII.feconwlHi , FeII.nfe2con );
692 
693  /* this must be true */
694  ASSERT( FeII.nFeIILevel_malloc == FeII.nFeIILevel_local );
695 
696  return;
697 }
698 
699 /*
700  *=====================================================================
701  */
702 /***********************************************************************
703  *** version of FeIILevelPops.f with overlap in procces 05.28.97 ooooooo
704  ******change in common block *te* sqrte 05.28.97
705  *******texc is fixed 03.28.97
706  *********this version has work on overlap
707  *********this version has # of zones (ZoneCnt) 07.03.97
708  *********taux - optical depth depends on iter correction 03.03.97
709  *********calculate cooling (Fe2_large_cool) and output fecool from Cloudy 01.29.97god
710  *********and cooling derivative (ddT_Fe2_large_cool)
711  ************ this program for 371 level iron model 12/14/1995
712  ************ this program for 371 level iron model 1/11/1996
713  ************ this program for 371 level iron model 3/21/1996
714  ************ this program without La 3/27/1996
715  ************ this program for 371 level iron model 4/9/1996
716  ************ includes:FeIICollRatesBoltzmann;cooling;overlapping in lines */
717 STATIC bool lgFeIIEverCalled=false;
718 void FeIILevelPops( void )
719 {
720  long int i,
721  ipHi ,
722  ipLo ,
723  n;
724  /* used in test for non-positive level populations */
725  bool lgPopNeg;
726 
727  double
728  EnergyWN,
729  pop ,
730  sum;
731 
732  int32 info;
733  int32 ipiv[NFE2LEVN];
734 
735  DEBUG_ENTRY( "FeIILevelPops()" );
736 
737  if( trace.lgTrace )
738  {
739  fprintf( ioQQQ," FeIILevelPops fe2 pops called\n");
740  }
741 
742  /* FeII.lgSimulate was set true with simulate flag on atom FeII command,
743  * for bebugging without actually calling the routine */
744  if( FeII.lgSimulate )
745  return;
746 
747  lgFeIIEverCalled = true;
748  /* zero out some arrays */
749  for( n=0; n<FeII.nFeIILevel_local; ++n)
750  {
751  for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo )
752  {
753  Fe2CPump[ipLo][n] = 0.;
754  Fe2LPump[ipLo][n] = 0.;
755  Fe2A[ipLo][n] = 0.;
756  xMatrix[ipLo][n] = 0.;
757  }
758  }
759 
760  /* make the g-bar collision strengths and do linear interpolation on r-matrix data.
761  * this also sets Boltzmann factors for all levels,
762  * sets values of FeColl used below, but only if temp has changed */
763  FeIICollRatesBoltzmann();
764 
765  /* pumping and spontantous decays */
766  for( n=0; n<FeII.nFeIILevel_local; ++n)
767  {
768  for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
769  {
770  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][n]];
771  /* continuum pumping rate from n to upper ipHi */
772  Fe2CPump[n][ipHi] = tr.Emis().pump();
773 
774  /* continuum pumping rate from ipHi to lower n */
775  Fe2CPump[ipHi][n] = tr.Emis().pump()*
776  (*tr.Hi()).g()/(*tr.Lo()).g();
777 
778  /* spontaneous decays */
779  Fe2A[ipHi][n] = tr.Emis().Aul()*(tr.Emis().Pesc() + tr.Emis().Pelec_esc() +
780  tr.Emis().Pdest());
781  }
782  }
783 
784  /* now do pumping of atom by Lya */
785  FeIILyaPump();
786 
787  /* **************************************************************************
788  *
789  * final rates into matrix
790  *
791  ***************************************************************************/
792 
793  /* fill in xMatrix with matrix elements */
794  for( n=0; n<FeII.nFeIILevel_local; ++n)
795  {
796  /* all processes leaving level n going down*/
797  for( ipLo=0; ipLo<n; ++ipLo )
798  {
799  xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipLo] + Fe2LPump[n][ipLo]+ Fe2A[n][ipLo] +
800  Fe2Coll[n][ipLo]*dense.eden;
801  }
802  /* all processes leaving level n going up*/
803  for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
804  {
805  xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipHi] + Fe2LPump[n][ipHi] + Fe2Coll[n][ipHi]*dense.eden;
806  }
807  /* all processes entering level n from below*/
808  for( ipLo=0; ipLo<n; ++ipLo )
809  {
810  xMatrix[ipLo][n] = xMatrix[ipLo][n] - Fe2CPump[ipLo][n] - Fe2LPump[ipLo][n] - Fe2Coll[ipLo][n]*dense.eden;
811  }
812  /* all processes entering level n from above*/
813  for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
814  {
815  xMatrix[ipHi][n] = xMatrix[ipHi][n] - Fe2CPump[ipHi][n] - Fe2LPump[ipHi][n] - Fe2Coll[ipHi][n]*dense.eden -
816  Fe2A[ipHi][n];
817  }
818 
819  /* the top row of the matrix is the sum of level populations */
820  xMatrix[n][0] = 1.0;
821  }
822 
823  {
824  /* option to print out entire matrix */
825  enum {DEBUG_LOC=false};
826  if( DEBUG_LOC )
827  {
828  /* print the matrices */
829  for( n=0; n<FeII.nFeIILevel_local; ++n)
830  {
831  /*fprintf(ioQQQ,"\n");*/
832  /* now print the matrix*/
833  for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo )
834  {
835  fprintf(ioQQQ," %.1e", xMatrix[n][ipLo]);
836  }
837  fprintf(ioQQQ,"\n");
838  }
839  }
840  }
841 
842  /* define the Y Vector. The oth element is the sum of all level populations
843  * adding up to the total population. The remaining elements are the level
844  * balance equations adding up to zero */
845  yVector[0] = 1.0;
846  for( n=1; n < FeII.nFeIILevel_local; n++ )
847  {
848  yVector[n] = 0.0;
849  }
850 
851  /* create the 1-yVector array that will save vector,
852  * this is the macro trick */
853 # ifdef AMAT
854 # undef AMAT
855 # endif
856 # define AMAT(I_,J_) (*(amat+(I_)*FeII.nFeIILevel_local+(J_)))
857 
858  /* copy current contents of xMatrix array over to special amat array*/
859  for( ipHi=0; ipHi < FeII.nFeIILevel_local; ipHi++ )
860  {
861  for( i=0; i < FeII.nFeIILevel_local; i++ )
862  {
863  AMAT(i,ipHi) = xMatrix[i][ipHi];
864  }
865  }
866 
867  info = 0;
868 
869  /* do the linear algebra to find the level populations */
870  getrf_wrapper(FeII.nFeIILevel_local, FeII.nFeIILevel_local, amat, FeII.nFeIILevel_local, ipiv, &info);
871  getrs_wrapper('N', FeII.nFeIILevel_local, 1, amat, FeII.nFeIILevel_local, ipiv, yVector, FeII.nFeIILevel_local, &info);
872 
873  if( info != 0 )
874  {
875  fprintf( ioQQQ, "DISASTER FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" );
876  cdEXIT(EXIT_FAILURE);
877  }
878 
879  /* yVector now contains the level populations */
880 
881  /* this better be false after this loop - if not then non-positive level pops */
882  lgPopNeg = false;
883  /* copy all level pops over to Fe2LevelPop */
884  for( ipLo=0; ipLo < FeII.nFeIILevel_local; ipLo++ )
885  {
886  if(yVector[ipLo] < 0. )
887  {
888  lgPopNeg = true;
889  fprintf(ioQQQ,"PROBLEM FeIILevelPops finds non-positive level population, level is %ld pop is %g\n",
890  ipLo , yVector[ipLo] );
891  }
892  /* this is now correct level population, cm^-3 */
893  Fe2LevelPop[ipLo] = yVector[ipLo] * dense.xIonDense[ipIRON][1];
894  }
895  if( lgPopNeg )
896  {
897  /* this is here to use the lgPopNeg value for something, if not here then
898  * lint and some compilers will note that var is never used */
899  fprintf(ioQQQ , "PROBLEM FeIILevelPops exits with negative level populations.\n");
900  }
901 
902  /* >>chng 06 jun 24, make sure remainder of populations up through max
903  * limit are zero - this makes safe the case where the number
904  * of levels actually computed has been trimmed down from largest
905  * possible number of levels, for instance, in cool gas */
906  for( ipLo=FeII.nFeIILevel_local; ipLo < FeII.nFeIILevel_malloc; ++ipLo )
907  {
908  Fe2LevelPop[ipLo] = 0.;
909  }
910 
911  /* now set line opacities, intensities, and level populations
912  * >>chng 06 jun ipLo should go up to FeII.nFeIILevel_local-1 since this
913  * is the largest lower level with non-zero population */
914  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
915  {
916  /* >>chng 06 jun 24, upper level should go to limit
917  * of all that were allocated */
918  for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
919  {
920  /* >>chng 06 jun 24, in all of these the product
921  * yVector[ipHi]*dense.xIonDense[ipIRON][1] has been replaced
922  * with Fe2LevelPop[ipLo] - should always have been this way,
923  * and saves a mult */
924  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
925  tr.Emis().PopOpc() = (Fe2LevelPop[ipLo] -
926  Fe2LevelPop[ipHi]*(*tr.Lo()).g()/(*tr.Hi()).g());
927 
928  (*tr.Lo()).Pop() = Fe2LevelPop[ipLo];
929 
930  (*tr.Hi()).Pop() = Fe2LevelPop[ipHi];
931 
932  tr.Emis().phots() = Fe2LevelPop[ipHi]*
933  tr.Emis().Aul()*(tr.Emis().Pesc() + tr.Emis().Pelec_esc() );
934 
935  tr.Emis().xIntensity() = tr.Emis().phots()*
936  tr.EnergyErg();
937 
938  /* ratio of collisional (new) to pumped excitations */
939  /* >>chng 02 mar 04, add test MAX2 to prevent div by zero */
940  tr.Emis().ColOvTot() = (realnum)(Fe2Coll[ipLo][ipHi]*dense.eden /
941  SDIV( Fe2Coll[ipLo][ipHi]*dense.eden + Fe2CPump[ipLo][ipHi] + Fe2LPump[ipLo][ipHi] ) );
942  }
943  }
944 
945  /* only do this if large atom is on and more than ground terms computed -
946  * do not if only lowest levels are computed */
947  if( FeII.lgFeIILargeOn )
948  {
949  /* the hydrogen Lya destruction rate, then probability */
950  hydro.dstfe2lya = 0.;
951  EnergyWN = 0.;
952  /* count how many photons were removed-added */
953  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
954  {
955  for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_local; ipHi++ )
956  {
957  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
958  EnergyWN += Fe2LPump[ipLo][ipHi] + Fe2LPump[ipHi][ipLo];
959  hydro.dstfe2lya += (realnum)(
960  (*tr.Lo()).Pop()*Fe2LPump[ipLo][ipHi] -
961  (*tr.Hi()).Pop()*Fe2LPump[ipHi][ipLo]);
962  }
963  }
964  /* the destruction prob comes from
965  * dest rate = n(2p) * A(21) * PDest */
966  pop = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop();
967  if( pop > SMALLFLOAT )
968  {
969  hydro.dstfe2lya /= (realnum)(pop * iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul());
970  }
971  else
972  {
973  hydro.dstfe2lya = 0.;
974  }
975  /* NB - do not update hydro.dstfe2lya here if large FeII not on since
976  * done in FeII overlap */
977  }
978 
979  {
980  enum {DEBUG_LOC=false};
981  if( DEBUG_LOC)
982  {
983  /*lint -e644 EnergyWN not init */
984  fprintf(ioQQQ," sum all %.1e dest rate%.1e escR= %.1e\n",
985  EnergyWN,hydro.dstfe2lya,
986  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc());
987  /*lint +e644 EnergyWN not init */
988  }
989  }
990 
991  /* next two blocks determine departure coefficients for the atom */
992 
993  /* first sum up partition function for the model atom */
994  Fe2DepCoef[0] = 1.0;
995  sum = 1.0;
996  for( i=1; i < FeII.nFeIILevel_local; i++ )
997  {
998  /* Boltzmann factor relative to ground times ratio of statistical weights */
999  Fe2DepCoef[i] = Fe2DepCoef[0]*FeIIBoltzmann[i]*
1000  (*Fe2LevN[ipFe2LevN[i][0]].Hi()).g()/ (*Fe2LevN[ipFe2LevN[1][0]].Lo()).g();
1001  /* this sum is the partition function for the atom */
1002  sum += Fe2DepCoef[i];
1003  }
1004 
1005  /* now renormalize departure coefficients with partition function */
1006  for( i=0; i < FeII.nFeIILevel_local; i++ )
1007  {
1008  /* divide by total partition function, Fe2DepCoef is now the fraction of the
1009  * population that is in this level in TE */
1010  Fe2DepCoef[i] /= sum;
1011 
1012  /* convert to true departure coefficient */
1013  if( Fe2DepCoef[i]>SMALLFLOAT )
1014  {
1015  Fe2DepCoef[i] = yVector[i]/Fe2DepCoef[i];
1016  }
1017  else
1018  {
1019  Fe2DepCoef[i] = 0.;
1020  }
1021  }
1022 
1023  /*calculate cooling, heating, and cooling derivative */
1024  /* this is net cooling, cooling minus heating */
1025  FeII.Fe2_large_cool = 0.0f;
1026  /* this is be heating, not heating minus cooling */
1027  FeII.Fe2_large_heat = 0.f;
1028  /* derivative of cooling */
1029  FeII.ddT_Fe2_large_cool = 0.0f;
1030 
1031  /* compute heating and cooling due to model atom */
1032  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
1033  {
1034  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
1035  {
1036  double OneLine;
1037  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1038 
1039  /* net cooling due to single line */
1040  OneLine = (Fe2Coll[ipLo][ipHi]*Fe2LevelPop[ipLo] - Fe2Coll[ipHi][ipLo]*Fe2LevelPop[ipHi])*
1041  dense.eden*tr.EnergyErg();
1042 
1043  /* net cooling due to this line */
1044  FeII.Fe2_large_cool += (realnum)MAX2(0., OneLine);
1045 
1046  /* net heating due to line */
1047  FeII.Fe2_large_heat += (realnum)MAX2(0., -OneLine);
1048 
1049  /* derivative of FeII cooling */
1050  if( OneLine >= 0. )
1051  {
1052  /* net coolant, exponential dominates deriv */
1053  FeII.ddT_Fe2_large_cool += (realnum)OneLine*
1054  (Fe2LevN[ipFe2LevN[ipHi][0]].EnergyK()*thermal.tsq1 - thermal.halfte);
1055  }
1056  else
1057  {
1058  /* net heating, te^-1/2 dominates */
1059  FeII.ddT_Fe2_large_cool -= (realnum)OneLine*thermal.halfte;
1060  }
1061  }
1062  }
1063 
1064  return;
1065 }
1066 
1067 /*FeIICollRatesBoltzmann evaluate collision strenths,
1068  * both interpolating on r-mat and creating g-bar
1069  * find Boltzmann factors, evaluate collisional rate coefficients */
1070 /* >>05 dec 03, verified that this rountine only changes temperature dependent
1071  * quantities - nothing with temperature dependence is done */
1072 STATIC void FeIICollRatesBoltzmann(void)
1073 {
1074  /* this will be used to only reevaluate cs when temp changes */
1075  static double OldTemp = -1.;
1076  long int i,
1077  ipLo ,
1078  ipHi,
1079  ipTrim;
1080  realnum ag,
1081  cg,
1082  dg,
1083  gb,
1084  y;
1085  realnum FracLowTe , FracHighTe;
1086  static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
1087  long int ipTemp,
1088  ipTempp1 ,
1089  ipLoFe2,
1090  ipHiFe2;
1091  static long int nFeII_old = -1;
1092 
1093  DEBUG_ENTRY( "FeIICollRatesBoltzmann()" );
1094 
1095  /* return if temperature has not changed */
1096  /* >>chng 06 feb 14, add test on number of levels changing */
1097  if( fp_equal( phycon.te, OldTemp ) && FeII.nFeIILevel_local == nFeII_old )
1098  {
1099 
1100  return;
1101  }
1102  OldTemp = phycon.te;
1103  nFeII_old = FeII.nFeIILevel_local;
1104 
1105  /* find g-bar collision strength for all levels
1106  * expression for g-bar taken from
1107  *>>refer Fe2 gbar Mewe, R. 1972, A&A, 20, 215 */
1108  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
1109  {
1110  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
1111  {
1112  /* these coefficients are from page 221 of Mewe 1972 */
1113  if( ncs1[ipHi][ipLo] == 3 )
1114  {
1115  /************************forbidden tr**************************/
1116  ag = 0.15f;
1117  cg = 0.f;
1118  dg = 0.f;
1119  }
1120  /************************allowed tr*******************************/
1121  else if( ncs1[ipHi][ipLo] == 1 )
1122  {
1123  ag = 0.6f;
1124  cg = 0.f;
1125  dg = 0.28f;
1126  }
1127  /************************semiforbidden****************************/
1128  else if( ncs1[ipHi][ipLo] == 2 )
1129  {
1130  ag = 0.0f;
1131  cg = 0.1f;
1132  dg = 0.0f;
1133  }
1134  else
1135  {
1136  /* this branch is impossible, since cs must be 1, 2, or 3 */
1137  ag = 0.0f;
1138  cg = 0.1f;
1139  dg = 0.0f;
1140  fprintf(ioQQQ,">>>PROBLEM impossible ncs1 in FeIILevelPops\n");
1141  }
1142 
1143  /*y = 1.438826f*Fe2LevN[ipHi][ipLo].EnergyWN/ phycon.te;*/
1144  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1145  y = tr.EnergyWN()/ (realnum)phycon.te_wn;
1146 
1147  gb = (realnum)(ag + (-cg*POW2(y) + dg)*(log(1.0+1.0/y) - 0.4/
1148  POW2(y + 1.0)) + cg*y);
1149 
1150  tr.Coll().col_str() = 23.861f*1e5f*gb*
1151  tr.Emis().Aul()*(*tr.Hi()).g()/
1152  POW3(tr.EnergyWN());
1153 
1154  /* confirm that collision strength is positive */
1155  ASSERT( tr.Coll().col_str() > 0.);
1156 
1157  /* g-bar cs becomes unphysically small for forbidden transitions -
1158  * this sets a lower limit to the final cs - CS2SMALL is defined above */
1159  tr.Coll().col_str() = MAX2( CS2SMALL, tr.Coll().col_str());
1160  /* this was the logic used in the old fortran version,
1161  * and reproduces the results in Baldwin et al '96
1162  if( Fe2LevN[ipHi][ipLo].cs < 1e-10 )
1163  {
1164  Fe2LevN[ipHi][ipLo].cs = 0.01f;
1165  }
1166  */
1167  }
1168  }
1169  /* end loop setting Mewe 1972 g-bar approximation */
1170 
1171  /* we will interpolate on the set of listed collision strengths -
1172  * where in this set are we? */
1173  if( phycon.te <= tt[0] )
1174  {
1175  /* temperature is lower than lowest tabulated, use the
1176  * lowest tabulated point */
1177  /* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher,
1178  * here both are lowest */
1179  ipTemp = 0;
1180  ipTempp1 = 0;
1181  FracHighTe = 0.;
1182  }
1183  else if( phycon.te > tt[7] )
1184  {
1185  /* temperature is higher than highest tabulated, use the
1186  * highest tabulated point */
1187  /* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher,
1188  * here both are highest */
1189  ipTemp = 7;
1190  ipTempp1 = 7;
1191  FracHighTe = 0.;
1192  }
1193  else
1194  {
1195  /* where in the array is the temperature we need? */
1196  ipTemp = -1;
1197  for( i=0; i < 8-1; i++ )
1198  {
1199  if( phycon.te <= tt[i+1] )
1200  {
1201  ipTemp = i;
1202  break;
1203  }
1204 
1205  }
1206  /* this cannot possibly happen */
1207  if( ipTemp < 0 )
1208  {
1209  fprintf( ioQQQ, " Insanity while looking for temperature in coll str array, te=%g.\n",
1210  phycon.te );
1211  cdEXIT(EXIT_FAILURE);
1212  }
1213  /* ipTemp points to the cell cooler than te, ipTemp+1 to one higher,
1214  * do linear interpolation between */
1215  ipTemp = i;
1216  ipTempp1 = i+1;
1217  FracHighTe = ((realnum)phycon.te - tt[ipTemp])/(tt[ipTempp1] - tt[ipTemp]);
1218  }
1219 
1220  /* this is fraction of final linear interpolated collision strength that
1221  * is weighted by the lower bound cs */
1222  FracLowTe = 1.f - FracHighTe;
1223 
1224  for( ipHi=1; ipHi < NPRADDAT; ipHi++ )
1225  {
1226  for( ipLo=0; ipLo < ipHi; ipLo++ )
1227  {
1228  /* ipHiFe2 should point to upper level of this pair, and
1229  * ipLoFe2 should point to lower level */
1230  ipHiFe2 = MAX2( nnPradDat[ipHi] , nnPradDat[ipLo] );
1231  ipLoFe2 = MIN2( nnPradDat[ipHi] , nnPradDat[ipLo] );
1232  ASSERT( ipHiFe2-1 < NFE2LEVN );
1233  ASSERT( ipHiFe2-1 >= 0 );
1234  ASSERT( ipLoFe2-1 < NFE2LEVN );
1235  ASSERT( ipLoFe2-1 >= 0 );
1236 
1237  /* do linear interpolation for CS, these are dimensioned NPRADDAT = 159 */
1238  if( ipHiFe2-1 < FeII.nFeIILevel_local )
1239  {
1240  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHiFe2-1][ipLoFe2-1]];
1241  /*fprintf(ioQQQ,"DEBUG\t%.2e", Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs);*/
1242  /* do linear interpolation */
1243  tr.Coll().col_str() =
1244  sPradDat[ipHi][ipLo][ipTemp] * FracLowTe +
1245  sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe;
1246  /*fprintf(ioQQQ,"\t%.2e\n", Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs);*/
1247 
1248  /* confirm that this is positive */
1249  ASSERT( tr.Coll().col_str() > 0. );
1250  }
1251  }
1252  }
1253 
1254  /* create Boltzmann factors for all levels */
1255  FeIIBoltzmann[0] = 1.0;
1256  for( ipHi=1; ipHi < FeII.nFeIILevel_local; ipHi++ )
1257  {
1258  /*FeIIBoltzmann[ipHi] = sexp( 1.438826f*Fe2LevN[ipHi][0].EnergyWN/phycon.te );*/
1259  /* >>chng 99 may 21, from above to following slight different number from phyconst.h
1260  * >>chng 05 dec 01, from sexp to dsexp to avoid all 0 Boltzmann factors in low temps
1261  * now that FeII is on all the time */
1262  FeIIBoltzmann[ipHi] = dsexp( Fe2LevN[ipFe2LevN[ipHi][0]].EnergyWN()/phycon.te_wn );
1263  }
1264 
1265  /* now possibly trim down atom if Boltzmann factors for upper levels are zero */
1266  ipTrim = FeII.nFeIILevel_local - 1;
1267  while( FeIIBoltzmann[ipTrim] == 0. && ipTrim > 0 )
1268  {
1269  --ipTrim;
1270  }
1271  /* ipTrim now is the highest level with finite Boltzmann factor -
1272  * use this as the number of levels
1273  *>>chng 05 dec 01, from <=1 to <=0 - func_map does 10K, and large FeII had never
1274  * been tested in that limit. 1 is ok. */
1275  ASSERT( ipTrim > 0 );
1276  /* trim FeII.nFeIILevel_local so that FeIIBoltzmann is positive
1277  * in nearly all cases this does nothing since ipTrim and FeII.nFeIILevel_local
1278  * are equal . . . */
1279  if( ipTrim+1 < FeII.nFeIILevel_local )
1280  {
1281  /* >>chng 06 jun 27, zero out collision data */
1282  for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo )
1283  {
1284  for( ipHi=ipTrim; ipHi<FeII.nFeIILevel_local; ++ipHi )
1285  {
1286  Fe2Coll[ipLo][ipHi] = 0.;
1287  Fe2Coll[ipHi][ipLo] = 0.;
1288  }
1289  }
1290  }
1291  FeII.nFeIILevel_local = ipTrim+1;
1292  /*fprintf(ioQQQ," levels reset to %li\n",FeII.nFeIILevel_local);*/
1293 
1294  /* debug code to print out the collision strengths for some levels */
1295  {
1296  enum {DEBUG_LOC=false};
1297  if( DEBUG_LOC)
1298  {
1299  for( ipLo=0; ipLo < 52; ipLo++ )
1300  {
1301  fprintf(ioQQQ,"%e %e\n",
1302  Fe2LevN[ipFe2LevN[51][ipLo]].Coll().col_str(),Fe2LevN[ipFe2LevN[52][ipLo]].Coll().col_str());
1303  }
1304  cdEXIT(EXIT_FAILURE);
1305  }
1306  }
1307 
1308  /* collisional excitation and deexcitation */
1309  for( ipLo=0; ipLo<FeII.nFeIILevel_local; ++ipLo)
1310  {
1311  for( ipHi=ipLo+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
1312  {
1313  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1314  /* collisional deexcitation rate coefficient from ipHi to lower ipLo
1315  * note that it needs eden to become rate
1316  * units cm3 s-1 */
1317  Fe2Coll[ipHi][ipLo] = (realnum)(COLL_CONST/phycon.sqrte*tr.Coll().col_str()/
1318  (*tr.Hi()).g());
1319  /*fprintf(ioQQQ,"DEBUG fe2coll %li %li %.2e", ipHi , ipLo , Fe2Coll[ipHi][ipLo] );*/
1320 
1321  /* collisional excitation rate coefficient from ipLo to upper ipHi
1322  * units cm3 s-1
1323  * FeIIBoltzmann guaranteed to be > 0 by nFeIILevel_local trimming */
1324  Fe2Coll[ipLo][ipHi] = (realnum)(Fe2Coll[ipHi][ipLo]*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]*
1325  (*tr.Hi()).g()/(*tr.Lo()).g());
1326  /*fprintf(ioQQQ,"DEBUG fe2coll %.2e\n", Fe2Coll[ipLo][ipHi] );*/
1327  }
1328  }
1329 
1330  return;
1331 }
1332 
1333 /*
1334  *=====================================================================
1335  */
1336 /*subroutine FeIIPrint PhotOccNum_at_nu raspechatki naselennostej v cloudy.out ili v
1337  * otdel'nom file unit=33
1338  *!!nado takzhe vklyuchit raspechatku iz perekrytiya linii */
1339 /*FeIIPrint - print output from large FeII atom, called by prtzone */
1340 void FeIIPrint(void)
1341 {
1342 
1343  DEBUG_ENTRY( "FeIIPrint()" );
1344 
1345  return;
1346 }
1347 
1348 /*
1349  *=====================================================================
1350  */
1351 /*FeIISumBand, sum up large FeII emission over certain bands
1352  * units are erg s-1 cm-3, same units as xIntensity in line structure */
1353 /* >>chng 06 apr 11, from using erg as energy unit to vacuum angstroms,
1354  * this fixes bug in physical constant and also uses air wavelengths for wl > 2000A */
1355 double FeIISumBand(
1356  /* lower and upper range to wavelength in Angstroms */
1357  realnum wl1,
1358  realnum wl2,
1359  double *SumBandInward )
1360 {
1361  long int ipHi,
1362  ipLo;
1363  double SumBandFe2_v;
1364 
1365  DEBUG_ENTRY( "FeIISumBand()" );
1366  /*sum up large FeII emission over certain bands */
1367 
1368  SumBandFe2_v = 0.;
1369  *SumBandInward = 0.;
1370  if( dense.xIonDense[ipIRON][1] > SMALLFLOAT )
1371  {
1372  /* line energy in wavenumber */
1373  ASSERT( wl2 > wl1 );
1374  for( ipHi=1; ipHi < FeII.nFeIILevel_local; ++ipHi )
1375  {
1376  for( ipLo=0; ipLo < ipHi; ++ipLo )
1377  {
1378  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1379  if( tr.WLAng() >= wl1 &&
1380  tr.WLAng() < wl2 )
1381  {
1382  SumBandFe2_v += tr.Emis().xIntensity();
1383  *SumBandInward += tr.Emis().xIntensity() *
1384  tr.Emis().FracInwd();
1385  }
1386  }
1387  }
1388  }
1389  return( SumBandFe2_v );
1390 }
1391 
1392 /*
1393  *=====================================================================
1394  */
1395 /*FeII_RT_TauInc called once per zone in RT_tau_inc to increment large FeII atom line optical depths */
1396 void FeII_RT_TauInc(void)
1397 {
1398  long int ipHi,
1399  ipLo;
1400 
1401  DEBUG_ENTRY( "FeII_RT_TauInc()" );
1402 
1403  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
1404  {
1405  /* >>chng 06 jun 24, for upper level go all the way to the
1406  * largest possible number of levels so that optical depths
1407  * of UV transitions are correct even for very cold gas where
1408  * the high level populations are not computed */
1409  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1410  {
1411  /* >>chng 04 aug 28, add test on bogus line */
1412  const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1413  if( tr.ipCont() > 0 )
1414  RT_line_one_tauinc( tr, -8, -8, ipHi, ipLo, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
1415  }
1416  }
1417 
1418  return;
1419 }
1420 
1421 /*
1422  *=====================================================================
1423  */
1424 /*FeII_RT_tau_reset reset optical depths for large FeII atom, called by update after each iteration */
1425 void FeII_RT_tau_reset(void)
1426 {
1427  long int ipHi,
1428  ipLo;
1429 
1430  DEBUG_ENTRY( "FeII_RT_tau_reset()" );
1431 
1432  /*fprintf(ioQQQ,"DEBUG FeIITauAver1 %li %.2e %.2e nFeIILevel_local %li \n",
1433  iteration,
1434  Fe2LevN[9][0].Emis().TauIn() ,
1435  Fe2LevN[9][0].Emis().TauTot(),
1436  FeII.nFeIILevel_local);*/
1437 
1438  /* called at end of iteration */
1439  /* >>chng 05 dec 07, had been FeII.nFeIILevel_local but this may have been trimmed down
1440  * if previous iteration went very deep - not reset until FeIIReset called */
1441  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1442  {
1443  for( ipLo=0; ipLo < ipHi; ipLo++ )
1444  {
1445  RT_line_one_tau_reset( Fe2LevN[ipFe2LevN[ipHi][ipLo]] );
1446  }
1447  }
1448 
1449  return;
1450 }
1451 
1452 /*
1453  *=====================================================================
1454  */
1455 /*FeIIPoint called by ContCreatePointers to create pointers for lines in large FeII atom */
1456 void FeIIPoint(void)
1457 {
1458  long int ipHi,
1459  ip,
1460  ipLo;
1461 
1462  DEBUG_ENTRY( "FeIIPoint()" );
1463 
1464  /* routine called when cloudy sets continuum array indices for Fe2 lines,
1465  * once per coreload */
1466  for( ipLo=0; ipLo < FeII.nFeIILevel_malloc-1; ipLo++ )
1467  {
1468  for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1469  {
1470 
1471  /* >>chng 02 feb 11, set continuum index to negative value for fake transition */
1472  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1473  if( fabs(tr.Emis().Aul()- 1e-5 ) > 1e-8 )
1474  {
1475  ip = ipoint(tr.EnergyRyd());
1476  tr.ipCont() = ip;
1477 
1478  /* do not over write other pointers with FeII since FeII is everywhere */
1479  if( strcmp( rfield.chLineLabel[ip-1], " " ) == 0 )
1480  strcpy( rfield.chLineLabel[ip-1], "FeII" );
1481 
1482  tr.Emis().ipFine() = ipFineCont( tr.EnergyRyd());
1483  }
1484  else
1485  {
1486  tr.ipCont() = -1;
1487  tr.Emis().ipFine() = -1;
1488  }
1489 
1490  tr.Emis().dampXvel() =
1491  (realnum)(tr.Emis().Aul()/
1492  tr.EnergyWN()/PI4);
1493 
1494  /* derive the abs coef, call to function is gf, wl (A), g_low */
1495  tr.Emis().opacity() =
1496  (realnum)abscf(tr.Emis().gf(),
1497  tr.EnergyWN(),
1498  (*tr.Lo()).g());
1499  }
1500  }
1501 
1502  return;
1503 }
1504 
1505 /*
1506  *=====================================================================
1507  */
1508 /*FeIIAccel called by rt_line_driving to compute radiative acceleration due to FeII lines */
1509 void FeIIAccel(double *fe2drive)
1510 {
1511  long int ipHi,
1512  ipLo;
1513 
1514  DEBUG_ENTRY( "FeIIAccel()" );
1515  /*compute acceleration due to large FeII atom */
1516 
1517  /* this routine computes the line driven radiative acceleration
1518  * due to large FeII atom*/
1519 
1520  *fe2drive = 0.;
1521  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
1522  {
1523  for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel_local; ipHi++ )
1524  {
1525  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1526  *fe2drive += tr.Emis().pump()*
1527  tr.EnergyErg()*tr.Emis().PopOpc();
1528  }
1529  }
1530 
1531  return;
1532 }
1533 
1534 /*
1535  *=====================================================================
1536  */
1537 /*FeII_RT_Make called by RT_line_all, does large FeII atom radiative transfer */
1538 void FeII_RT_Make( void )
1539 {
1540  long int ipHi,
1541  ipLo;
1542 
1543  DEBUG_ENTRY( "FeII_RT_Make()" );
1544 
1545  if( trace.lgTrace )
1546  {
1547  fprintf( ioQQQ," FeII_RT_Make called\n");
1548  }
1549 
1550  /* this routine drives calls to make RT relations with large FeII atom */
1551  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
1552  {
1553  /* >>chng 06 jun 24, go up to number allocated to keep UV resonance lines */
1554  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1555  {
1556  /* only evaluate real transitions */
1557  const TransitionProxy &tr=Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1558  if( tr.ipCont() > 0 )
1559  {
1560  /*RT_line_one do rt for emission line structure -
1561  * calls RT_line_escape or RT_line_wind */
1562  RT_line_one( tr, true, 0.f, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
1563  }
1564  }
1565  }
1566 
1567  return;
1568 }
1569 
1570 /*
1571  *=====================================================================
1572  */
1573 /* called in LineSet4 to add FeII lines to save array */
1574 void FeIIAddLines( void )
1575 {
1576 
1577  DEBUG_ENTRY( "FeIIAddLines()" );
1578 
1579  /* this routine puts the emission from the large FeII atom
1580  * into an array to save and integrate them*/
1581 
1582  /* add lines option called from lines, add intensities into storage array */
1583 
1584  /* routine is called three different ways, ipass < 0 before space allocated,
1585  * =0 when time to generate labels (and we zero out array here), and ipass>0
1586  * when time to save intensities */
1587  if( LineSave.ipass == 0 )
1588  {
1589  /* we were called by lines, and we want to zero out Fe2SavN */
1590  for( long ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
1591  {
1592  for( long ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1593  {
1594  Fe2SavN[ipHi][ipLo] = 0.;
1595  }
1596  }
1597  }
1598 
1599  /* this call during calculations, save intensities */
1600  else if( LineSave.ipass == 1 )
1601  {
1602  /* total emission from vol */
1603  for( long ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
1604  {
1605  for( long ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
1606  {
1607  Fe2SavN[ipHi][ipLo] +=
1608  radius.dVeffAper*Fe2LevN[ipFe2LevN[ipHi][ipLo]].Emis().xIntensity();
1609  }
1610  }
1611  }
1612 
1613  {
1614  enum {DEBUG_LOC=false};
1615  if( DEBUG_LOC /*&& (iteration==2)*/ )
1616  {
1617  fprintf(ioQQQ," 69-1\t%li\t%e\n", nzone , Fe2SavN[68][0] );
1618  }
1619  }
1620 
1621  return;
1622 }
1623 
1624 /*FeIIPunchLevels save level energies and stat weights */
1625 void FeIIPunchLevels(
1626  /* file we will save to */
1627  FILE *ioPUN )
1628 {
1629 
1630  long int ipHi;
1631 
1632  DEBUG_ENTRY( "FeIIPunchLevels()" );
1633 
1634  fprintf(ioPUN , "%.2f\t%li\n", 0., (long)(*Fe2LevN[ipFe2LevN[1][0]].Lo()).g() );
1635  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ++ipHi )
1636  {
1637  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][0]];
1638  fprintf(ioPUN , "%.2f\t%li\n",
1639  tr.EnergyWN() ,
1640  (long)(*tr.Hi()).g());
1641  }
1642 
1643  return;
1644 }
1645 
1646 /*FeIIPunchColden save level energies, stat weights, column density */
1647 void FeIIPunchColden(
1648  /* file we will save to */
1649  FILE *ioPUN )
1650 {
1651 
1652  long int n;
1653 
1654  DEBUG_ENTRY( "FeIIPunchColden()" );
1655 
1656  fprintf(ioPUN , "%.2f\t%li\t%.3e\n", 0., (long)(*Fe2LevN[ipFe2LevN[1][0]].Lo()).g() , Fe2ColDen[0]);
1657  for( n=1; n < FeII.nFeIILevel_malloc; ++n )
1658  {
1659  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[n][0]];
1660  fprintf(ioPUN , "%.2f\t%li\t%.3e\n",
1661  tr.EnergyWN() ,
1662  (long)(*tr.Hi()).g(),
1663  Fe2ColDen[n]);
1664  }
1665 
1666  return;
1667 }
1668 
1669 
1670 /*FeIIPunchOpticalDepth save FeII optical depths,
1671  * called by punch_do to save them out,
1672  * save turned on with save lines command */
1673 void FeIIPunchOpticalDepth(
1674  /* file we will save to */
1675  FILE *ioPUN )
1676 {
1677  long int
1678  ipHi,
1679  ipLo;
1680 
1681  DEBUG_ENTRY( "FeIIPunchOpticalDepth()" );
1682 
1683  /* this routine punches the optical depths predicted by the large FeII atom */
1684  /*>>chng 06 may 19, must use total number of levels allocated not current
1685  * number since this decreases as gas grows colder with depth nFeIILevel_local->nFeIILevel_malloc*/
1686  for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
1687  {
1688  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1689  {
1690  /* fe2ener(1) and (2) are lower and upper limits to range of
1691  * wavelengths of interest. entered in ryd with
1692  * save FeII verner command, where they are converted
1693  * to wavenumbers, */
1694  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1695  fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.2e\n",
1696  ipHi+1,
1697  ipLo+1,
1698  tr.WLAng() ,
1699  tr.Emis().TauIn() );
1700  }
1701  }
1702 
1703  return;
1704 }
1705 
1706 /*FeIISaveLines save accumulated FeII intensities, at end of calculation,
1707  * called by punch_do to save them out,
1708  * save turned on with save lines command */
1709 void FeIISaveLines(
1710  /* file we will save to */
1711  FILE *ioPUN )
1712 {
1713  long int MaseHi,
1714  MaseLow,
1715  ipHi,
1716  ipLo;
1717  double hbeta, absint , renorm;
1718  /* >>chng 00 jun 02, demoted next two to realnum, PvH */
1719  realnum TauMase,
1720  thresh;
1721 
1722  DEBUG_ENTRY( "FeIISaveLines()" );
1723 
1724  /* this routine puts the emission from the large FeII atom
1725  * into a line array, and eventually will save it out */
1726 
1727  /* get the normalization line */
1728  if( LineSv[LineSave.ipNormWavL].SumLine[0] > 0. )
1729  renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].SumLine[0];
1730  else
1731  renorm = 1.;
1732 
1733  fprintf( ioPUN, " up low log I, I/I(LineSave), Tau\n" );
1734 
1735  /* first look for any masing lines */
1736  MaseLow = -1;
1737  MaseHi = -1;
1738  TauMase = 0.f;
1739  for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
1740  {
1741  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1742  {
1743  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1744  if( tr.Emis().TauIn() < TauMase )
1745  {
1746  TauMase = tr.Emis().TauIn();
1747  MaseLow = ipLo;
1748  MaseHi = ipHi;
1749  }
1750  }
1751  }
1752 
1753  if( TauMase < 0.f )
1754  {
1755  fprintf( ioPUN, " Most negative optical depth was %4ld%4ld%10.2e\n",
1756  MaseHi, MaseLow, TauMase );
1757  }
1758 
1759  /* now print actual line intensities, Hbeta first */
1760  if( cdLine("TOTL", 4861.36f , &hbeta , &absint)<=0 )
1761  {
1762  fprintf( ioQQQ, " FeIILevelPops could not find Hbeta with cdLine.\n" );
1763  cdEXIT(EXIT_FAILURE);
1764  }
1765 
1766  fprintf( ioPUN, "Hbeta=%7.3f %10.2e\n",
1767  absint ,
1768  hbeta );
1769 
1770  if( renorm > SMALLFLOAT )
1771  /* this is threshold for faintest line, normally 0, set with
1772  * number on save feii lines command */
1773  thresh = FeII.fe2thresh/(realnum)renorm;
1774  else
1775  thresh = 0.f;
1776 
1777  /* must use total number of levels allocated not current
1778  * number since this decreases as gas grows colder with depth nFeIILevel_malloc*/
1779  for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
1780  {
1781  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1782  {
1783  /* fe2ener(1) and (2) are lower and upper limits to range of
1784  * wavelengths of interest. entered in ryd with
1785  * save FeII verner command, where they are converted
1786  * to wavenumbers, */
1787  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1788  if( (Fe2SavN[ipHi][ipLo] > thresh &&
1789  tr.EnergyWN() > FeII.fe2ener[0]) &&
1790  tr.EnergyWN() < FeII.fe2ener[1] )
1791  {
1792  if( FeII.lgShortFe2 )
1793  {
1794  /* short option on save FeII line command
1795  * does not include rel intensity or optical dep */
1796  fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.3f\n",
1797  ipHi+1,
1798  ipLo+1,
1799  tr.WLAng() ,
1800  log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten );
1801  }
1802  else
1803  {
1804  /* long printout does */
1805  fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.3f\t%.2e\t%.2e\n",
1806  ipHi+1,
1807  ipLo+1,
1808  tr.WLAng() ,
1809  log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten,
1810  Fe2SavN[ipHi][ipLo]*renorm ,
1811  tr.Emis().TauIn() );
1812  }
1813  }
1814  }
1815  }
1816 
1817  return;
1818 }
1819 
1820 
1821 /*
1822  *=====================================================================
1823  */
1824 /*FeII_LineZero zero out storage for large FeII atom, called in tauout */
1825 void FeII_LineZero(void)
1826 {
1827  long int ipHi,
1828  ipLo;
1829 
1830  DEBUG_ENTRY( "FeII_LineZero()" );
1831 
1832  /* this routine is called in routine zero and it
1833  * zero's out various elements of the FeII arrays
1834  * it is called on every iteration
1835  * */
1836  for( ipLo=0; ipLo < (FeII.nFeIILevel_malloc - 1); ipLo++ )
1837  {
1838  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1839  {
1840  /* >>chng 03 feb 14, use EmLineZero rather than explicit sets */
1841  Fe2LevN[ipFe2LevN[ipHi][ipLo]].Zero();
1842  }
1843  }
1844 
1845  return;
1846 }
1847 /*
1848  *=====================================================================
1849  */
1850 /*FeIIIntenZero zero out storage for large FeII atom, called in tauout */
1851 void FeIIIntenZero(void)
1852 {
1853  long int ipHi,
1854  ipLo;
1855 
1856  DEBUG_ENTRY( "FeIIIntenZero()" );
1857 
1858  /* this routine is called by routine cool_iron and it
1859  * zero's out various elements of the FeII arrays
1860  * */
1861  Fe2LevelPop[0] = 0.;
1862  for( ipHi=1; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
1863  {
1864  /* >>chng 05 dec 14, zero Fe2LevelPop */
1865  Fe2LevelPop[ipHi] = 0.;
1866  for( ipLo=0; ipLo < ipHi; ipLo++ )
1867  {
1868  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
1869 
1870  /* population of lower level with correction for stim emission */
1871  tr.Emis().PopOpc() = 0.;
1872 
1873  /* population of lower level */
1874  (*tr.Lo()).Pop() = 0.;
1875 
1876  /* population of upper level */
1877  (*tr.Hi()).Pop() = 0.;
1878 
1879  /* following two heat exchange excitation, deexcitation */
1880  tr.Coll().cool() = 0.;
1881  tr.Coll().heat() = 0.;
1882 
1883  /* intensity of line */
1884  tr.Emis().xIntensity() = 0.;
1885 
1886  tr.Emis().phots() = 0.;
1887  tr.Emis().ots() = 0.;
1888  tr.Emis().ColOvTot() = 0.;
1889  }
1890  }
1891 
1892  (*TauLines[ipTuv3].Lo()).Pop() = 0;
1893  (*TauLines[ipTuv3].Hi()).Pop() = 0;
1894  TauLines[ipTuv3].Emis().PopOpc() = 0;
1895  TauLines[ipTuv3].Emis().phots() = 0;
1896  TauLines[ipTuv3].Emis().xIntensity() = 0;
1897 
1898  (*TauLines[ipTr48].Lo()).Pop() = 0;
1899  (*TauLines[ipTr48].Hi()).Pop() = 0;
1900  TauLines[ipTr48].Emis().PopOpc() = 0;
1901  TauLines[ipTr48].Emis().phots() = 0;
1902  TauLines[ipTr48].Emis().xIntensity() = 0;
1903 
1904  FeII.for7 = 0;
1905 
1906  (*TauLines[ipTFe16].Lo()).Pop() = 0;
1907  (*TauLines[ipTFe16].Hi()).Pop() = 0;
1908  TauLines[ipTFe16].Emis().PopOpc() = 0;
1909  TauLines[ipTFe16].Emis().phots() = 0;
1910  TauLines[ipTFe16].Emis().xIntensity() = 0;
1911 
1912  (*TauLines[ipTFe26].Lo()).Pop() = 0;
1913  (*TauLines[ipTFe26].Hi()).Pop() = 0;
1914  TauLines[ipTFe26].Emis().PopOpc() = 0;
1915  TauLines[ipTFe26].Emis().phots() = 0;
1916  TauLines[ipTFe26].Emis().xIntensity() = 0;
1917 
1918  (*TauLines[ipTFe34].Lo()).Pop() = 0;
1919  (*TauLines[ipTFe34].Hi()).Pop() = 0;
1920  TauLines[ipTFe34].Emis().PopOpc() = 0;
1921  TauLines[ipTFe34].Emis().phots() = 0;
1922  TauLines[ipTFe34].Emis().xIntensity() = 0;
1923 
1924  (*TauLines[ipTFe35].Lo()).Pop() = 0;
1925  (*TauLines[ipTFe35].Hi()).Pop() = 0;
1926  TauLines[ipTFe35].Emis().PopOpc() = 0;
1927  TauLines[ipTFe35].Emis().phots() = 0;
1928  TauLines[ipTFe35].Emis().xIntensity() = 0;
1929 
1930  (*TauLines[ipTFe46].Lo()).Pop() = 0;
1931  (*TauLines[ipTFe46].Hi()).Pop() = 0;
1932  TauLines[ipTFe46].Emis().PopOpc() = 0;
1933  TauLines[ipTFe46].Emis().phots() = 0;
1934  TauLines[ipTFe46].Emis().xIntensity() = 0;
1935 
1936  (*TauLines[ipTFe56].Lo()).Pop() = 0;
1937  (*TauLines[ipTFe56].Hi()).Pop() = 0;
1938  TauLines[ipTFe56].Emis().PopOpc() = 0;
1939  TauLines[ipTFe56].Emis().phots() = 0;
1940  TauLines[ipTFe56].Emis().xIntensity() = 0;
1941 
1942  (*TauLines[ipT191].Lo()).Pop() = 0;
1943  (*TauLines[ipT191].Hi()).Pop() = 0;
1944  TauLines[ipT191].Emis().PopOpc() = 0;
1945  TauLines[ipT191].Emis().phots() = 0;
1946  TauLines[ipT191].Emis().xIntensity() = 0;
1947 
1948  return;
1949 }
1950 
1951 /*
1952  *=====================================================================
1953  * save line data for FeII atom */
1954 void FeIIPunData(
1955  /* io unit for save */
1956  FILE* ioPUN ,
1957  /* save all levels if true, only subset if false */
1958  bool lgDoAll )
1959 {
1960  long int ipLo , ipHi;
1961 
1962  DEBUG_ENTRY( "FeIIPunData()" );
1963 
1964  if( lgDoAll )
1965  {
1966  fprintf( ioQQQ,
1967  " FeIIPunData ALL option not implemented yet 1\n" );
1968  cdEXIT(EXIT_FAILURE);
1969  }
1970  else if( lgFeIIEverCalled )
1971  {
1972  long int nSkip=0, limit=MIN2(64, FeII.nFeIILevel_local);
1973 
1974  /* false, only save subset in above init */
1975  /* first 64 do all lines */
1976  //fprintf( ioPUN, "#Lo\tHi\tIon\tWL\tgl\tgu\tgf\tA\tCS\tn(crt)\tdamp\n" );
1977  bool lgPrint = true;
1978  for( ipHi=1; ipHi<limit; ++ipHi )
1979  {
1980  for( ipLo=0; ipLo<ipHi; ++ipLo )
1981  {
1982  //fprintf(ioPUN,"%4li\t%4li\t",ipLo,ipHi );
1983  Save1LineData( Fe2LevN[ipFe2LevN[ipHi][ipLo]] , ioPUN, false , lgPrint);
1984  }
1985  }
1986  fprintf( ioPUN , "\n");
1987 
1988  if( limit==64 )
1989  {
1990  /* higher than 64 only do real transitions - the majority of those above
1991  * n = 64 have no radiative transitions but still exist to hold collision,
1992  * energy, and other line data - they will have Aul == 1e-5 */
1993  for( ipHi=limit; ipHi<FeII.nFeIILevel_local; ++ipHi )
1994  {
1995  /*>>chng 06 jun 23, ipLo had ranged from limit up to ipHi,
1996  * so never printed lines from ipHi to ipLo<64 */
1997  for( ipLo=0; ipLo<ipHi; ++ipLo )
1998  {
1999  /*fprintf(ioPUN , "%li %li\n", ipHi , ipLo );
2000  if( ipHi==70 && ipLo==0 )
2001  Save1LineData( Fe2LevN.begin()+ipHi][ipLo] , ioPUN , false); */
2002  /* ncs1 of 3 and Aul = 1e-5 indicate transition that is not optically
2003  * allowed with fake cs */
2004  const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
2005  if( ncs1[ipHi][ipLo] == 3 &&
2006  (fabs(tr.Emis().Aul()-1e-5) < 1e-8 ) )
2007  {
2008  ++nSkip;
2009  }
2010  else
2011  {
2012  /* add one to ipLo and ipHi so that printed number is on atomic,
2013  * not C, scale */
2014  //fprintf(ioPUN,"%4li\t%4li\t",ipLo+1,ipHi+1 );
2015  Save1LineData( tr , ioPUN , false , lgPrint);
2016  }
2017  }
2018  }
2019  fprintf( ioPUN , " %li lines skipped\n",nSkip);
2020  }
2021  }
2022 
2023  return;
2024 }
2025 
2026 /*
2027  *=====================================================================
2028  */
2029 void FeIIPunDepart(
2030  /* io unit for save */
2031  FILE* ioPUN ,
2032  /* save all levels if true, only subset if false */
2033  bool lgDoAll )
2034 {
2035  /* numer of levels with dep coef that we will save out */
2036 # define NLEVDEP 11
2037  /* these are the levels on the physical, not c, scale (count from 1) */
2038  const int LevDep[NLEVDEP]={1,10,25,45,64,124,206,249,295,347,371};
2039  long int i;
2040  static bool lgFIRST=true;
2041 
2042  DEBUG_ENTRY( "FeIIPunDepart()" );
2043 
2044  /* on first call only, print levels that we will save later */
2045  if( lgFIRST && !lgDoAll )
2046  {
2047  /* but all the rest do */
2048  for( i=0; i<NLEVDEP; ++i )
2049  {
2050  fprintf( ioPUN , "%i\t", LevDep[i] );
2051  }
2052  fprintf( ioPUN , "\n");
2053  lgFIRST = false;
2054  }
2055 
2056  if( lgDoAll )
2057  {
2058  /* true, save all levels, one per line */
2059  for( i=1; i<=FeII.nFeIILevel_local; ++i )
2060  {
2061  FeIIPun1Depart( ioPUN , i );
2062  fprintf( ioPUN , "\n");
2063  }
2064  }
2065  else
2066  {
2067  /* false, only save subset in above init */
2068  for( i=0; i<NLEVDEP; ++i )
2069  {
2070  FeIIPun1Depart( ioPUN , LevDep[i] );
2071  fprintf( ioPUN , "\t");
2072  }
2073  fprintf( ioPUN , "\n");
2074  }
2075 
2076  return;
2077 }
2078 
2079 /*
2080  *=====================================================================
2081  */
2082 void FeIIPun1Depart(
2083  /* the io unit where the print should be directed */
2084  FILE * ioPUN ,
2085  /* the physical (not c) number of the level */
2086  long int nPUN )
2087 {
2088 
2089  DEBUG_ENTRY( "FeIIPun1Depart()" );
2090 
2091  ASSERT( nPUN > 0 );
2092  /* need real assert to keep lint happy */
2093  assert( ioPUN != NULL );
2094 
2095  /* print the level departure coef on ioPUN if the level was computed,
2096  * print a zero if it was not */
2097  if( nPUN <= FeII.nFeIILevel_local )
2098  {
2099  fprintf( ioPUN , "%e ",Fe2DepCoef[nPUN-1] );
2100  }
2101  else
2102  {
2103  fprintf( ioPUN , "%e ",0. );
2104  }
2105 
2106  return;
2107 }
2108 
2109 /*
2110  *=====================================================================
2111  */
2112 void FeIIReset(void)
2113 {
2114 
2115  DEBUG_ENTRY( "FeIIReset()" );
2116 
2117  /* space has been allocated, so reset number of levels to whatever it was */
2118  FeII.nFeIILevel_local = FeII.nFeIILevel_malloc;
2119 
2120  return;
2121 }
2122 
2123 /*
2124  *=====================================================================
2125  */
2126 /*FeIIZero initialize some variables, called by zero one time before commands parsed*/
2127 void FeIIZero(void)
2128 {
2129 
2130  DEBUG_ENTRY( "FeIIZero()" );
2131 
2132  /* heating, cooling, and deriv wrt temperature */
2133  FeII.Fe2_large_cool = 0.;
2134  FeII.ddT_Fe2_large_cool = 0.;
2135  FeII.Fe2_large_heat = 0.;
2136 
2137  /* flag saying that lya pumping of FeII in large atom is turned on */
2138  FeII.lgLyaPumpOn = true;
2139 
2140  /*FeII. lgFeIION = false;*/
2141  /* >>chng 05 nov 29, test effects of always including FeII ground term with large atom
2142  * but if ground term only is done, still also do simple UV approximation */
2143  /*FeII. lgFeIION = true;*/
2144  /* will not compute large FeII atom */
2145  FeII.lgFeIILargeOn = false;
2146 
2147  /* energy range of large FeII atom is zero to infinity */
2148  FeII.fe2ener[0] = 0.;
2149  FeII.fe2ener[1] = 1e8;
2150 
2151  /* pre mar 01, these had both been 1, ipPRD */
2152  /* resonance lines, ipCRD is -1 */
2153  FeII.ipRedisFcnResonance = ipCRD;
2154  /* subordinate lines, ipCRDW is 2 */
2155  FeII.ipRedisFcnSubordinate = ipCRDW;
2156 
2157  /* set zero for the threshold of weakest large FeII atom line to save */
2158  FeII.fe2thresh = 0.;
2159 
2160  /* normally do not constantly reevaluate the atom, set true with
2161  * SLOW key on atom FeII command */
2162  FeII.lgSlow = false;
2163 
2164  /* option to print each call to FeIILevelPops, set with print option on atom FeII */
2165  FeII.lgPrint = false;
2166 
2167  /* option to only simulate calls to FeIILevelPops */
2168  FeII.lgSimulate = false;
2169 
2170  /* set number of levels for the atom
2171  * changed with the atom FeII levels command */
2172  if( lgFeIIMalloc )
2173  {
2174  /* space has been allocated, so reset number of levels to whatever it was */
2175  FeII.nFeIILevel_local = FeII.nFeIILevel_malloc;
2176  }
2177  else
2178  {
2179  /* always include FeII ground term with large atom
2180  * but if ground term only is done, still also do simple UV approximation
2181  * set this to only ground term, - will reset to NFE2LEVN when atom FeII parsed if levels not set */
2182  FeII.nFeIILevel_local = 16;
2183  }
2184 
2185  return;
2186 }
2187 
2188 /*FeIIPunPop - save level populations */
2189 void FeIIPunPop(
2190  /* io unit for save */
2191  FILE* ioPUN ,
2192  /* save range of levels if true, only selected subset if false */
2193  bool lgPunchRange ,
2194  /* if ipPunchRange is true, this is lower bound of range on C scale */
2195  long int ipRangeLo ,
2196  /* if ipPunchRange is true, this is upper bound of range on C scale */
2197  long int ipRangeHi ,
2198  /* flag saying whether to save density (cm-3, true) or relative population (flase) */
2199  bool lgPunchDensity )
2200 {
2201  /* numer of levels with dep coef that we will save out */
2202 # define NLEVPOP 11
2203  /* these are the levels on the physical, not c, scale (count from 1) */
2204  const int LevPop[NLEVPOP]={1,10,25,45,64,124,206,249,295,347,371};
2205  long int i;
2206  static bool lgFIRST=true;
2207  realnum denominator = 1.f;
2208 
2209  DEBUG_ENTRY( "FeIIPunPop()" );
2210 
2211  /* this implements the relative option on save FeII populations command */
2212  if( !lgPunchDensity )
2213  denominator = SDIV( dense.xIonDense[ipIRON][1] );
2214 
2215  /* on first call only, print levels that we will save later,
2216  * but only if we will only save selected levels*/
2217  if( lgFIRST && !lgPunchRange )
2218  {
2219  /* but all the rest do */
2220  for( i=0; i<NLEVPOP; ++i )
2221  {
2222  /* indices for particular levels */
2223  fprintf( ioPUN , "%i\t", LevPop[i] );
2224  }
2225  fprintf( ioPUN , "\n");
2226  lgFIRST = false;
2227  }
2228 
2229  if( lgPunchRange )
2230  {
2231  /* confirm sane level indices */
2232  ASSERT( ipRangeLo>=0 && ipRangeLo<ipRangeHi );
2233 
2234  /* true, save all levels across line,
2235  * both call with physical level so that list is physical */
2236  long nHigh = MIN2(ipRangeHi, FeII.nFeIILevel_local );
2237  for( i=ipRangeLo; i<nHigh; ++i )
2238  {
2239  /* routine takes levels on physics scale */
2240  fprintf( ioPUN , "%.3e\t", Fe2LevelPop[i]/denominator );
2241  }
2242  fprintf( ioPUN , "\n");
2243  }
2244  else
2245  {
2246  /* false, only save subset in above init,
2247  * both call with physical level so that list is physical */
2248  for( i=0; i<NLEVPOP; ++i )
2249  {
2250  fprintf( ioPUN , "%.3e\t", Fe2LevelPop[LevPop[i]-1]/denominator );
2251  }
2252  fprintf( ioPUN , "\n");
2253  }
2254 
2255  return;
2256 }
2257 
2258 /*
2259  *=====================================================================
2260  */
2261 #if 0
2262 void FeIIPun1Pop(
2263  /* the io unit where the print should be directed */
2264  FILE * ioPUN ,
2265  /* the physical (not c) number of the level */
2266  long int nPUN )
2267 {
2268  DEBUG_ENTRY( "FeIIPun1Pop()" );
2269 
2270  ASSERT( nPUN > 0 );
2271  /* need real assert to keep lint happy */
2272  assert( ioPUN != NULL );
2273 
2274  /* print the level population on ioPUN if the level was computed,
2275  * print a zero if it was not,
2276  * note that nPUN is on physical scale, so test is <=*/
2277  if( nPUN <= FeII.nFeIILevel_local )
2278  {
2279  fprintf( ioPUN , "%e ",Fe2LevelPop[nPUN-1] );
2280  }
2281  else
2282  {
2283  fprintf( ioPUN , "%e ",0. );
2284  }
2285 
2286  return;
2287 }
2288 #endif
2289 
2290 /*
2291  *=====================================================================
2292  */
2293 STATIC int FeIIBandsCreate(
2294  /* chFile is optional filename, if void then use default bands,
2295  * if not void then use file specified,
2296  * return value is 0 for success, 1 for failure */
2297  const char chFile[] )
2298 {
2299 
2300  char chLine[FILENAME_PATH_LENGTH_2];
2301  const char* chFile1;
2302  FILE *ioDATA;
2303 
2304  bool lgEOL;
2305  long int i,k;
2306 
2307  /* keep track of whether we have been called - want to be
2308  * called a total of one time */
2309  static bool lgCalled=false;
2310 
2311  DEBUG_ENTRY( "FeIIBandsCreate()" );
2312 
2313  /* return previous number of bands if this is second or later call*/
2314  if( lgCalled )
2315  {
2316  /* success */
2317  return 0;
2318  }
2319  lgCalled = true;
2320 
2321  /* use default filename if void string, else use file specified */
2322  chFile1 = ( strlen(chFile) == 0 ) ? "FeII_bands.ini" : chFile;
2323 
2324  /* get FeII band data */
2325  if( trace.lgTrace )
2326  {
2327  fprintf( ioQQQ, " FeIICreate opening %s:", chFile1 );
2328  }
2329 
2330  ioDATA = open_data( chFile1, "r" );
2331 
2332  /* now count how many bands are in the file */
2333  nFeIIBands = 0;
2334 
2335  /* first line is a version number and does not count */
2336  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
2337  {
2338  fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 );
2339  return 1;
2340  }
2341  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
2342  {
2343  /* we want to count the lines that do not start with #
2344  * since these contain data */
2345  if( chLine[0] != '#')
2346  ++nFeIIBands;
2347  }
2348 
2349  /* now rewind the file so we can read it a second time*/
2350  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
2351  {
2352  fprintf( ioQQQ, " FeIICreate could not rewind %s.\n", chFile1 );
2353  return 1;
2354  }
2355 
2356  FeII_Bands = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(nFeIIBands) );
2357 
2358  /* now make second dim, id wavelength, and lower and upper bounds */
2359  for( i=0; i<nFeIIBands; ++i )
2360  {
2361  FeII_Bands[i] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(3) );
2362  }
2363 
2364  /* first line is a version number - now confirm that it is valid */
2365  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
2366  {
2367  fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 );
2368  return 1;
2369  }
2370  {
2371  i = 1;
2372  const long int iyr = 9, imo=6 , idy = 11;
2373  long iyrread, imoread , idyread;
2374  iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
2375  imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
2376  idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
2377 
2378  if(( iyrread != iyr ) ||
2379  ( imoread != imo ) ||
2380  ( idyread != idy ) )
2381  {
2382  fprintf( ioQQQ,
2383  " PROBLEM FeIIBandsCreate: the version of %s is not the current version.\n", chFile1 );
2384  fprintf( ioQQQ,
2385  " FeIIBandsCreate: I expected the magic numbers %li %li %li but found %li %li %li.\n",
2386  iyr, imo , idy ,iyrread, imoread , idyread );
2387  return 1;
2388  }
2389  }
2390 
2391  /* now read in data again, but save it this time */
2392  k = 0;
2393  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
2394  {
2395  /* we want to count the lines that do not start with #
2396  * since these contain data */
2397  if( chLine[0] != '#')
2398  {
2399  i = 1;
2400  FeII_Bands[k][0] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
2401  if( lgEOL )
2402  {
2403  fprintf( ioQQQ, " There should have been a number on this band line 1. Sorry.\n" );
2404  fprintf( ioQQQ, "string==%s==\n" ,chLine );
2405  return 1;
2406  }
2407  FeII_Bands[k][1] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
2408  if( lgEOL )
2409  {
2410  fprintf( ioQQQ, " There should have been a number on this band line 2. Sorry.\n" );
2411  fprintf( ioQQQ, "string==%s==\n" ,chLine );
2412  return 1;
2413  }
2414  FeII_Bands[k][2] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
2415  if( lgEOL )
2416  {
2417  fprintf( ioQQQ, " There should have been a number on this band line 3. Sorry.\n" );
2418  fprintf( ioQQQ, "string==%s==\n" ,chLine );
2419  return 1;
2420  }
2421  /*fprintf(ioQQQ,
2422  " band data %f %f %f \n", FeII_Bands[k][0],FeII_Bands[k][1],FeII_Bands[k][2]);*/
2423  ++k;
2424  }
2425  }
2426  /* now validate this incoming data */
2427  for( i=0; i<nFeIIBands; ++i )
2428  {
2429  /* make sure all are positive */
2430  if( FeII_Bands[i][0] <=0. || FeII_Bands[i][1] <=0. || FeII_Bands[i][2] <=0. )
2431  {
2432  fprintf( ioQQQ, " FeII band %li has none positive entry.\n",i );
2433  return 1;
2434  }
2435  /* make sure bands bounds are in correct order, shorter - longer wavelength*/
2436  if( FeII_Bands[i][1] >= FeII_Bands[i][2] )
2437  {
2438  fprintf( ioQQQ, " FeII band %li has improper bounds.\n" ,i);
2439  return 1;
2440  }
2441 
2442  }
2443 
2444  fclose(ioDATA);
2445 
2446  /* return success */
2447  return 0;
2448 }
2449 /*
2450  *=====================================================================
2451  */
2452 void AssertFeIIDep( double *pred , double *BigError , double *StdDev )
2453 {
2454  long int n;
2455  double arg ,
2456  error ,
2457  sum2;
2458 
2459  DEBUG_ENTRY( "AssertFeIIDep()" );
2460 
2461  if( FeII.lgSimulate || !lgFeIIEverCalled )
2462  {
2463 
2464  *pred = 0.;
2465  *BigError = 0.;
2466  *StdDev = 0.;
2467 
2468  return;
2469  }
2470 
2471  /* sanity check */
2472  ASSERT( FeII.nFeIILevel_local > 0 );
2473 
2474  /* find sum of deviation of departure coef from unity */
2475  *BigError = 0.;
2476  *pred = 0.;
2477  sum2 = 0;
2478  for( n=0; n<FeII.nFeIILevel_local; ++n )
2479  {
2480  /* get mean */
2481  *pred += Fe2DepCoef[n];
2482 
2483  /* error is the largest deviation from unity for any single level*/
2484  error = fabs( Fe2DepCoef[n] -1. );
2485  /* remember biggest deviation */
2486  *BigError = MAX2( *BigError , error );
2487 
2488  /* get standard deviation */
2489  sum2 += POW2( Fe2DepCoef[n] );
2490  }
2491 
2492  /* now get standard deviation */
2493  arg = sum2 - POW2( *pred ) / (double)FeII.nFeIILevel_local;
2494  ASSERT( (arg >= 0.) );
2495  *StdDev = sqrt( arg / (double)(FeII.nFeIILevel_local - 1.) );
2496 
2497  /* this is average value, should be unity */
2498  *pred /= (double)(FeII.nFeIILevel_local);
2499 
2500  return;
2501 }
2502 
2503 /* do ots rates for FeII, called by RT_OTS */
2504 void FeII_OTS( void )
2505 {
2506  long int ipLo ,
2507  ipHi;
2508 
2509  DEBUG_ENTRY( "FeII_OTS()" );
2510 
2511  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
2512  {
2513  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
2514  {
2515  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
2516  /* >>chng 02 feb 11, skip bogus transitions */
2517  if( tr.ipCont() < 1)
2518  continue;
2519 
2520  /* ots rates, the destp prob was set in hydropesc */
2521  tr.Emis().ots() =
2522  tr.Emis().Aul()*
2523  (*tr.Hi()).Pop()*
2524  tr.Emis().Pdest();
2525 
2526  ASSERT( tr.Emis().ots() >= 0. );
2527 
2528  /* finally dump the ots rate into the stack */
2529  RT_OTS_AddLine(tr.Emis().ots(),
2530  tr.ipCont() );
2531  }
2532  }
2533 
2534  return;
2535 }
2536 
2537 /*
2538  *=====================================================================
2539  */
2540 /*FeII_RT_Out - do outward rates for FeII,
2541  * called by RT_diffuse, which is called by cloudy */
2542 void FeII_RT_Out(void)
2543 {
2544  long int ipLo , ipHi;
2545 
2546  DEBUG_ENTRY( "FeIIRTOut()" );
2547 
2548  /* only do this if Fe+ exists */
2549  if( dense.xIonDense[ipIRON][1] > 0. )
2550  {
2551  /* outward line photons */
2552  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
2553  {
2554  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
2555  {
2556  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
2557  /* >>chng 02 feb 11, skip bogus transitions */
2558  if( tr.ipCont() < 1)
2559  continue;
2560  /* >>chng 03 sep 09, change to outline, ouutlin is
2561  * not exactly the same in the two - tmn missing in outline */
2562  tr.outline_resonance();
2563 
2564  }
2565  }
2566  }
2567 
2568  return;
2569 }
2570 
2571 /*
2572  *=====================================================================
2573  */
2574 STATIC void FeIIContCreate(
2575  /* wavelength of low-lambda end */
2576  double xLamLow ,
2577  /* wavelength of high end */
2578  double xLamHigh ,
2579  /* number of cells to break this into */
2580  long int ncell )
2581 {
2582 
2583  double step , wl1;
2584  long int i;
2585 
2586  /* keep track of whether we have been called - want to be
2587  * called a total of one time */
2588  static bool lgCalled=false;
2589 
2590  DEBUG_ENTRY( "FeIIContCreate()" );
2591 
2592  /* return previous number of bands if this is second or later call*/
2593  if( lgCalled )
2594  {
2595  /* return value of number of bands, may be used by calling program*/
2596  return;
2597  }
2598  lgCalled = true;
2599 
2600  /* how many cells will be needed to go from xLamLow to xLamHigh in ncell steps */
2601  nFeIIConBins = ncell;
2602 
2603  FeII_Cont = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(nFeIIConBins+1) );
2604 
2605  /* now make second dim, id wavelength, and lower and upper bounds */
2606  for( i=0; i<nFeIIConBins+1; ++i )
2607  FeII_Cont[i] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(3) );
2608 
2609  step = log10( xLamHigh/xLamLow)/ncell;
2610  wl1 = log10( xLamLow);
2611  for( i=0; i<nFeIIConBins+1; ++i)
2612  // [n][0] are bounds of cells in Angstroms
2613  FeII_Cont[i][0] = (realnum)pow(10. , ( wl1 + i*step ) );
2614 
2615  return;
2616 }
2617 
2618 /*ParseAtomFeII parse the atom FeII command */
2619 void ParseAtomFeII( Parser &p )
2620 {
2621  DEBUG_ENTRY( "ParseAtomFeII()" );
2622 
2623  /* turn on the large verner atom */
2624  FeII.lgFeIILargeOn = true;
2625  /* >>chng 05 nov 29, reset number of levels when called, so increased above default of 16 */
2626  if( lgFeIIMalloc )
2627  {
2628  /* space has been allocated, so reset number of levels to whatever it was */
2629  FeII.nFeIILevel_local = FeII.nFeIILevel_malloc;
2630  }
2631  else
2632  {
2633  /* space not allocated yet, set to largest possible number of levels */
2634  FeII.nFeIILevel_local = NFE2LEVN;
2635  }
2636 
2637  /* levels keyword is to adjust number of levels. But this only has effect
2638  * BEFORE space is allocated for the FeII arrays */
2639  if( p.nMatch("LEVE") )
2640  {
2641  /* do option only if space not yet allocated */
2642  if( !lgFeIIMalloc )
2643  {
2644  /* number of levels for hydrogen at, 2s is this plus one */
2645  FeII.nFeIILevel_local = (long int)p.getNumberCheck("LEVEL");
2646  if( FeII.nFeIILevel_local <16 )
2647  {
2648  fprintf( ioQQQ, " This would be too few levels, must have at least 16.\n" );
2649  cdEXIT(EXIT_FAILURE);
2650  }
2651  else if( FeII.nFeIILevel_local > NFE2LEVN )
2652  {
2653  fprintf( ioQQQ, " This would be too many levels.\n" );
2654  cdEXIT(EXIT_FAILURE);
2655  }
2656  }
2657  }
2658 
2659  /* slow keyword means do not try to avoid evaluating atom */
2660  else if( p.nMatch("SLOW") )
2661  {
2662  FeII.lgSlow = true;
2663  }
2664 
2665  /* redistribution keyword changes form of redistribution function */
2666  else if( p.nMatch("REDI") )
2667  {
2668  int ipRedis=0;
2669  /* there are three functions, PRD_, CRD_, and CRDW,
2670  * representing partial redistribution,
2671  * complete redistribution (doppler core only, no wings)
2672  * and complete with wings */
2673  /* partial redistribution */
2674  if( p.nMatch(" PRD") )
2675  {
2676  ipRedis = ipPRD;
2677  }
2678  /* complete redistribution */
2679  else if( p.nMatch(" CRD") )
2680  {
2681  ipRedis = ipCRD;
2682  }
2683  /* complete redistribution with wings */
2684  else if( p.nMatch("CRDW") )
2685  {
2686  ipRedis = ipCRDW;
2687  }
2688 
2689  /* if not SHOW option (handled below) then we have a problem */
2690  else if( !p.nMatch("SHOW") )
2691  {
2692  fprintf(ioQQQ," There should have been a second keyword on this command.\n");
2693  fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n");
2694  cdEXIT(EXIT_FAILURE);
2695  }
2696 
2697  /* resonance lines */
2698  if( p.nMatch("RESO") )
2699  {
2700  FeII.ipRedisFcnResonance = ipRedis;
2701  }
2702  /* subordinate lines */
2703  else if( p.nMatch("SUBO") )
2704  {
2705  FeII.ipRedisFcnSubordinate = ipRedis;
2706  }
2707  /* the show option, say what we are assuming */
2708  else if( p.nMatch("SHOW") )
2709  {
2710  fprintf(ioQQQ," FeII resonance lines are ");
2711  if( FeII.ipRedisFcnResonance ==ipCRDW )
2712  {
2713  fprintf(ioQQQ,"complete redistribution with wings\n");
2714  }
2715  else if( FeII.ipRedisFcnResonance ==ipCRD )
2716  {
2717  fprintf(ioQQQ,"complete redistribution with core only.\n");
2718  }
2719  else if( FeII.ipRedisFcnResonance ==ipPRD )
2720  {
2721  fprintf(ioQQQ,"partial redistribution.\n");
2722  }
2723  else
2724  {
2725  fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnResonance.\n");
2726  TotalInsanity();
2727  }
2728 
2729  fprintf(ioQQQ," FeII subordinate lines are ");
2730  if( FeII.ipRedisFcnSubordinate ==ipCRDW )
2731  {
2732  fprintf(ioQQQ,"complete redistribution with wings\n");
2733  }
2734  else if( FeII.ipRedisFcnSubordinate ==ipCRD )
2735  {
2736  fprintf(ioQQQ,"complete redistribution with core only.\n");
2737  }
2738  else if( FeII.ipRedisFcnSubordinate ==ipPRD )
2739  {
2740  fprintf(ioQQQ,"partial redistribution.\n");
2741  }
2742  else
2743  {
2744  fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnSubordinate.\n");
2745  TotalInsanity();
2746  }
2747  }
2748  else
2749  {
2750  fprintf(ioQQQ," here should have been a second keyword on this command.\n");
2751  fprintf(ioQQQ," Options are RESONANCE, SUBORDINATE. Sorry.\n");
2752  cdEXIT(EXIT_FAILURE);
2753  }
2754  }
2755 
2756  /* trace keyword - print info for each call to FeIILevelPops */
2757  else if( p.nMatch("TRAC") )
2758  {
2759  FeII.lgPrint = true;
2760  }
2761 
2762  /* only simulate the FeII atom, do not actually call it */
2763  else if( p.nMatch("SIMU") )
2764  {
2765  /* option to only simulate calls to FeIILevelPops */
2766  FeII.lgSimulate = true;
2767  }
2768 
2769  /* only simulate the FeII atom, do not actually call it */
2770  else if( p.nMatch("CONT") )
2771  {
2772  /* the continuum output with the save FeII continuum command */
2773 
2774  /* number of levels for hydrogen at, 2s is this plus one */
2775  FeII.feconwlLo = (realnum)p.getNumberCheck("low wavelength");
2776  FeII.feconwlHi = (realnum)p.getNumberCheck("high wavelength");
2777  FeII.nfe2con = (long int)p.getNumberCheck("number of intervals");
2778  /* check that all are positive */
2779  if( FeII.feconwlLo<=0. || FeII.feconwlHi<=0. || FeII.nfe2con<= 0 )
2780  {
2781  fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
2782  fprintf(ioQQQ," all three must be greater than zero, sorry.\n");
2783  cdEXIT(EXIT_FAILURE);
2784  }
2785  /* make sure that wl1 is less than wl2 */
2786  if( FeII.feconwlLo>= FeII.feconwlHi )
2787  {
2788  fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
2789  fprintf(ioQQQ," the second wavelength must be greater than the first, sorry.\n");
2790  cdEXIT(EXIT_FAILURE);
2791  }
2792  }
2793  /* note no fall-through error since routine can be called with no options,
2794  * to turn on the large atom */
2795 
2796  return;
2797 }
2798 
2799 void PunFeII( FILE * io )
2800 {
2801  long int n, ipHi;
2802 
2803  DEBUG_ENTRY( "PunFeII()" );
2804 
2805  for( n=0; n<FeII.nFeIILevel_local-1; ++n)
2806  {
2807  for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
2808  {
2809  const TransitionProxy&tr = Fe2LevN[ipFe2LevN[ipHi][n]];
2810  if( tr.ipCont() > 0 )
2811  fprintf(io,"%li\t%li\t%.2e\n", n , ipHi ,
2812  tr.Emis().TauIn() );
2813  }
2814  }
2815 
2816  return;
2817 }
2818 
2819 /* include FeII lines in punched optical depths, etc, called from SaveLineStuff */
2820 void FeIIPunchLineStuff( FILE * io , realnum xLimit , long index)
2821 {
2822  long int n, ipHi;
2823 
2824  DEBUG_ENTRY( "FeIIPunchLineStuff()" );
2825 
2826  for( n=0; n<FeII.nFeIILevel_local-1; ++n)
2827  {
2828  for( ipHi=n+1; ipHi<FeII.nFeIILevel_local; ++ipHi )
2829  {
2830  Save1Line( Fe2LevN[ipFe2LevN[ipHi][n]] , io , xLimit , index, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
2831  }
2832  }
2833 
2834  return;
2835 }
2836 
2837 /* rad pre due to FeII lines called in PresTotCurrent*/
2838 double FeIIRadPress(void)
2839 {
2840 
2841  /* will be used to check on size of opacity, was capped at this value */
2842  realnum smallfloat=SMALLFLOAT*10.f;
2843  double press,
2844  RadPres1;
2845 # undef DEBUGFE
2846 # ifdef DEBUGFE
2847  double RadMax;
2848  long ipLoMax , ipHiMax;
2849 # endif
2850  long int ipLo, ipHi;
2851 
2852  DEBUG_ENTRY( "FeIIRadPress()" );
2853 
2854  press = 0.;
2855  if( !lgFeIIEverCalled )
2856  return 0.;
2857 
2858 # ifdef DEBUGFE
2859  RadMax = 0;
2860  ipLoMax = -1;
2861  ipHiMax = -1;
2862 # endif
2863  /* loop over all levels to find radiation pressure */
2864  for( ipHi=1; ipHi<FeII.nFeIILevel_local; ++ipHi )
2865  {
2866  for( ipLo=0; ipLo<ipHi; ++ipLo)
2867  {
2868  const TransitionProxy &tr = Fe2LevN[ipFe2LevN[ipHi][ipLo]];
2869  /* >>chng 04 aug 27, add test on ipCont() for bogus line */
2870  if( tr.ipCont() > 0 &&
2871  (*tr.Hi()).Pop() > 1e-30 )
2872  {
2873  if( (*tr.Hi()).Pop() > smallfloat &&
2874  tr.Emis().PopOpc() > smallfloat )
2875  {
2876  RadPres1 = PressureRadiationLine( tr, GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
2877 
2878 # ifdef DEBUGFE
2879  if( RadPres1 > RadMax )
2880  {
2881  RadMax = RadPres1;
2882  ipLoMax = ipLo;
2883  ipHiMax = ipHi;
2884  }
2885 # endif
2886  press += RadPres1;
2887  }
2888  }
2889  }
2890  }
2891 
2892 # ifdef DEBUGFE
2893  /* option to print radiation pressure */
2894  if( iteration > 1 || nzone > 1558 )
2895  {
2896  fprintf(ioQQQ,"DEBUG FeIIRadPress finds P(FeII) = %.2e %.2e %li %li width %.2e\n",
2897  press ,
2898  RadMax ,
2899  ipLoMax ,
2900  ipHiMax ,
2901  RT_LineWidth(Fe2LevN[ipFe2LevN[9][0]],GetDopplerWidth(dense.AtomicWeight[ipIRON])) );
2902  DumpLine( Fe2LevN.begin()+ipFe2LevN[9][0] );
2903  }
2904 # endif
2905 # undef DEBUGFE
2906 
2907  return press;
2908 }
2909 
2910 #if 0
2911 /* internal energy of FeII */
2912 double FeII_InterEnergy(void)
2913 {
2914  double energy;
2915 
2916  DEBUG_ENTRY( "FeII_InterEnergy()" );
2917 
2918  /* There is no stack of levels, so we access all levels
2919  * uniquely by via the line stack with Fe2LevN[ipHi][0].Hi */
2920 
2921  energy = 0.;
2922  for( long ipHi=1; ipHi<FeII.nFeIILevel_local; ++ipHi )
2923  {
2924  const TransitionList::iterator&tr = Fe2LevN.begin()+ipFe2LevN[ipHi][0];
2925  if( (*(*tr).Hi()).Pop() > 1e-30 )
2926  {
2927  energy +=
2928  (*(*tr).Hi()).Pop() * (*tr).EnergyErg;
2929  }
2930  }
2931 
2932  return energy;
2933 }
2934 #endif
2935 
2936 /* HP cc cannot compile this routine with any optimization */
2937 #if defined(__HP_aCC)
2938 # pragma OPT_LEVEL 1
2939 #endif
2940 /*FeIILyaPump find rate of Lya excitation of the FeII atom */
2941 STATIC void FeIILyaPump(void)
2942 {
2943 
2944  long int ipLo ,
2945  ipHi;
2946  double EnerLyaProf2,
2947  EnerLyaProf3,
2948  EnergyWN,
2949  Gup_ov_Glo,
2950  PhotOccNum_at_nu,
2951  PumpRate,
2952  de,
2953  FeIILineWidthHz,
2954  taux;
2955 
2956  DEBUG_ENTRY( "FeIILyaPump()" );
2957 
2958  /* lgLyaPumpOn is false if no Lya pumping, with no FeII command */
2959  /* get rates FeII atom is pumped */
2960 
2961  /* elsewhere in this file the dest prob hydro.dstfe2lya is defined from
2962  * quantites derived here, and the resulting populations */
2963  if( FeII.lgLyaPumpOn )
2964  {
2965 
2966  /*************trapeze form La profile:de,EnerLyaProf1,EnerLyaProf2,EnerLyaProf3,EnerLyaProf4*************************
2967  * */
2968  /* width of Lya in cm^-1 */
2969  /* HLineWidth has units of cm/s, as was evaluated in PresTotCurrent */
2970  /* the factor is 1/2 of E(Lya, cm^-1_/c */
2971  de = 1.37194e-06*hydro.HLineWidth*2.0/3.0;
2972  /* 82259 is energy of Lya in wavenumbers, so these are the form of the trapezoid */
2973  EnerLyaProf1 = 82259.0 - de*2.0;
2974  EnerLyaProf2 = 82259.0 - de;
2975  EnerLyaProf3 = 82259.0 + de;
2976  EnerLyaProf4 = 82259.0 + de*2.0;
2977 
2978  /* find Lya photon occupation number */
2979  if( iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop() > SMALLFLOAT )
2980  {
2981  /* This is the photon occupation number at the Lya line center */
2982  PhotOccNumLyaCenter =
2983  MAX2(0.,1.0- iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pelec_esc() -
2984  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc())/
2985  (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()/iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()*3. - 1.0);
2986  }
2987  else
2988  {
2989  /* lya excitation temperature not available */
2990  PhotOccNumLyaCenter = 0.;
2991  }
2992  }
2993  else
2994  {
2995  PhotOccNumLyaCenter = 0.;
2996  de = 0.;
2997  EnerLyaProf1 = 0.;
2998  EnerLyaProf2 = 0.;
2999  EnerLyaProf3 = 0.;
3000  EnerLyaProf4 = 0.;
3001  }
3002 
3003  /* is Lya pumping enabled? */
3004  if( FeII.lgLyaPumpOn )
3005  {
3006  for( ipLo=0; ipLo < (FeII.nFeIILevel_local - 1); ipLo++ )
3007  {
3008  for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel_local; ipHi++ )
3009  {
3010  const TransitionProxy&tr=Fe2LevN[ipFe2LevN[ipHi][ipLo]];
3011  /* on first iteration optical depth in line is inward only, on later
3012  * iterations is total optical depth */
3013  if( iteration == 1 )
3014  {
3015  taux = tr.Emis().TauIn();
3016  }
3017  else
3018  {
3019  taux = tr.Emis().TauTot();
3020  }
3021 
3022  /* Gup_ov_Glo is ratio of g values */
3023  Gup_ov_Glo = (*tr.Hi()).g()/(*tr.Lo()).g();
3024 
3025  /* the energy of the FeII line */
3026  EnergyWN = tr.EnergyWN();
3027 
3028  if( EnergyWN >= EnerLyaProf1 && EnergyWN <= EnerLyaProf4 && taux > 1 )
3029  {
3030  /* this branch, line is within the Lya profile */
3031 
3032  /*
3033  * Lya source function, at peak is PhotOccNumLyaCenter,
3034  *
3035  * Prof2 Prof3
3036  * ----------
3037  * / \
3038  * / \
3039  * / \
3040  * ======================
3041  * Prof1 Prof4
3042  *
3043  */
3044 
3045  if( EnergyWN < EnerLyaProf2 )
3046  {
3047  /* linear interpolation on edge of trapazoid */
3048  PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnergyWN - EnerLyaProf1)/ de;
3049  }
3050  else if( EnergyWN < EnerLyaProf3 )
3051  {
3052  /* this is the central plateau */
3053  PhotOccNum_at_nu = PhotOccNumLyaCenter;
3054  }
3055  else
3056  {
3057  /* linear interpolation on edge of trapazoid */
3058  PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnerLyaProf4 - EnergyWN)/de;
3059  }
3060 
3061  /* at this point Lya source function at FeII line energy is defined, but
3062  * we need to multiply by line width in Hz,
3063  * >>refer Fe2 pump Netzer, H., Elitzur, M., & Ferland, G. J. 1985, ApJ, 299, 752-762*/
3064 
3066  /* width of FeII line in Hz */
3067  FeIILineWidthHz = 1.e8/(EnergyWN*299792.5)*sqrt(log(taux))*GetDopplerWidth(dense.AtomicWeight[ipIRON]);
3068 
3069  /* final Lya pumping rate, s^-1*/
3070  PumpRate = FeIILineWidthHz * PhotOccNum_at_nu * tr.Emis().Aul()*
3071  powi(82259.0f/EnergyWN,3);
3072  /* above must be bogus, use just occ num times A */
3073  PumpRate = tr.Emis().Aul()* PhotOccNum_at_nu;
3074 
3075  /* Lya pumping rate from ipHi to lower n */
3076  Fe2LPump[ipHi][ipLo] += PumpRate;
3077 
3078  /* Lya pumping rate from n to upper ipHi */
3079  Fe2LPump[ipLo][ipHi] += PumpRate*Gup_ov_Glo;
3080  }
3081  }
3082  }
3083  }
3084 
3085  return;
3086 }
3087 
3088 /* end work around bugs in HP compiler */
3089 #if defined(__HP_aCC)
3090 #pragma OPTIMIZE OFF
3091 #pragma OPTIMIZE ON
3092 #endif
ipFineCont
long ipFineCont(double energy_ryd)
Definition: cont_ipoint.cpp:226
thermal.h
TransitionProxy::EnergyErg
realnum EnergyErg() const
Definition: transition.h:78
FeIIBoltzmann
static double * FeIIBoltzmann
Definition: atom_feii.cpp:159
FeIILyaPump
STATIC void FeIILyaPump(void)
Definition: atom_feii.cpp:2941
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
ipTFe16
long ipTFe16
Definition: atmdat_readin.cpp:87
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
lines.h
CollisionProxy::cool
double & cool() const
Definition: collision.h:190
Save1LineData
void Save1LineData(const TransitionProxy &t, FILE *io, bool lgCS_2, bool &lgPrint)
Definition: save_linedata.cpp:278
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
DumpLine
void DumpLine(const TransitionProxy &t)
Definition: transition.cpp:100
ipTuv3
long ipTuv3
Definition: atmdat_readin.cpp:86
t_FeII::fe2ener
realnum fe2ener[2]
Definition: atomfeii.h:233
dense
t_dense dense
Definition: dense.cpp:24
Fe2A
static double ** Fe2A
Definition: atom_feii.cpp:137
t_FeII::lgShortFe2
bool lgShortFe2
Definition: atomfeii.h:225
NPRADDAT
#define NPRADDAT
Definition: atom_feii.cpp:100
Fe2LPump
static double ** Fe2LPump
Definition: atom_feii.cpp:140
FeIISaveLines
void FeIISaveLines(FILE *ioPUN)
Definition: atom_feii.cpp:1709
rfield
t_rfield rfield
Definition: rfield.cpp:8
t_FeII::lgPrint
bool lgPrint
Definition: atomfeii.h:204
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
FeII_Colden
void FeII_Colden(const char *chLabel)
Definition: atom_feii.cpp:177
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:8
EmissionProxy::xIntensity
double & xIntensity() const
Definition: emission.h:483
FeIICollRatesBoltzmann
STATIC void FeIICollRatesBoltzmann(void)
Definition: atom_feii.cpp:1072
EmissionProxy::ipFine
long int & ipFine() const
Definition: emission.h:413
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
ipTFe46
long ipTFe46
Definition: atmdat_readin.cpp:88
t_radius::Conv2PrtInten
double Conv2PrtInten
Definition: radius.h:147
t_hydro::dstfe2lya
realnum dstfe2lya
Definition: hydrogenic.h:60
PunFeII
void PunFeII(FILE *io)
Definition: atom_feii.cpp:2799
realnum
float realnum
Definition: cddefines.h:103
t_radius::dVeffAper
double dVeffAper
Definition: radius.h:87
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
rfield.h
STATIC
#define STATIC
Definition: cddefines.h:97
t_FeII::FeIIColl
double FeIIColl[NFE2LEVN][NFE2LEVN][8]
Definition: atomfeii.h:216
multi_arr::reserve
void reserve(size_type i1)
Definition: container_classes.h:1080
PressureRadiationLine
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
Definition: pressure.h:18
RT_LineWidth
double RT_LineWidth(const TransitionProxy &t, realnum DopplerWidth)
Definition: rt_escprob.cpp:918
EmissionProxy::PopOpc
double & PopOpc() const
Definition: emission.h:603
TransitionList::resize
void resize(size_t newsize)
Definition: transition.h:285
EmissionProxy::iRedisFun
int & iRedisFun() const
Definition: emission.h:403
RT_line_one_tauinc
void RT_line_one_tauinc(const TransitionProxy &t, long int mas_species, long int mas_ion, long int mas_hi, long int mas_lo, realnum DopplerWidth)
Definition: rt_line_one_tauinc.cpp:16
t_FeII::lgLyaPumpOn
bool lgLyaPumpOn
Definition: atomfeii.h:229
TransitionProxy::setLo
void setLo(int ipLo) const
Definition: transition.h:400
Fe2CPump
static double ** Fe2CPump
Definition: atom_feii.cpp:140
USE_OLD
#define USE_OLD
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
ipIRON
const int ipIRON
Definition: cddefines.h:330
t_FeII::ipRedisFcnResonance
int ipRedisFcnResonance
Definition: atomfeii.h:244
thirdparty.h
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
t_LineSave::ipNormWavL
long int ipNormWavL
Definition: lines.h:81
TransitionList::begin
iterator begin(void)
Definition: transition.h:305
ipoint.h
ProxyIterator
Definition: proxy_iterator.h:58
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
lines_service.h
getrs_wrapper
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
Definition: thirdparty_lapack.cpp:69
EmissionProxy::Pdest
realnum & Pdest() const
Definition: emission.h:543
FeII_OTS
void FeII_OTS(void)
Definition: atom_feii.cpp:2504
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:450
t_FeII::lgFeIILargeOn
bool lgFeIILargeOn
Definition: atomfeii.h:198
t_FeII::feconwlHi
realnum feconwlHi
Definition: atomfeii.h:239
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
FeIIReset
void FeIIReset(void)
Definition: atom_feii.cpp:2112
iso.h
PI4
const UNUSED double PI4
Definition: physconst.h:35
t_FeII::nFeIILevel_local
long int nFeIILevel_local
Definition: atomfeii.h:191
FeII_RT_TauInc
void FeII_RT_TauInc(void)
Definition: atom_feii.cpp:1396
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
TransitionProxy::Coll
CollisionProxy Coll() const
Definition: transition.h:424
t_FeII::for7
double for7
Definition: atomfeii.h:259
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:423
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:16
POW2
#define POW2
Definition: cddefines.h:929
MIN2
#define MIN2
Definition: cddefines.h:761
PhotOccNumLyaCenter
static double PhotOccNumLyaCenter
Definition: atom_feii.cpp:164
TransitionProxy
Definition: transition.h:23
t_FeII::FeIIAul
double FeIIAul[NFE2LEVN][NFE2LEVN]
Definition: atomfeii.h:213
CollisionProxy::heat
double & heat() const
Definition: collision.h:194
lgFeIIMalloc
bool lgFeIIMalloc
Definition: cdinit.cpp:90
cddrive.h
ncs1
int ** ncs1
Definition: atom_feii.cpp:94
t_FeII::nfe2con
long int nfe2con
Definition: atomfeii.h:241
nzone
long int nzone
Definition: cddefines.cpp:14
LineSave
t_LineSave LineSave
Definition: lines.cpp:5
FeIIPunPop
void FeIIPunPop(FILE *ioPUN, bool lgPunchRange, long int ipRangeLo, long int ipRangeHi, bool lgPunchDensity)
Definition: atom_feii.cpp:2189
FeII_InterEnergy
double FeII_InterEnergy(void)
FeII_Cont
realnum ** FeII_Cont
Definition: atom_feii.cpp:113
multi_arr::size
size_type size() const
Definition: container_classes.h:1753
NLEVDEP
#define NLEVDEP
FeIISumBand
double FeIISumBand(realnum wl1, realnum wl2, double *SumBandInward)
Definition: atom_feii.cpp:1355
FeII_LineZero
void FeII_LineZero(void)
Definition: atom_feii.cpp:1825
FeIIPunData
void FeIIPunData(FILE *ioPUN, bool lgDoAll)
Definition: atom_feii.cpp:1954
radius
t_radius radius
Definition: radius.cpp:5
FeIIAddLines
void FeIIAddLines(void)
Definition: atom_feii.cpp:1574
lgFeIIEverCalled
STATIC bool lgFeIIEverCalled
Definition: atom_feii.cpp:717
t_FeII::nFeIILevel_malloc
long int nFeIILevel_malloc
Definition: atomfeii.h:193
TransitionProxy::EnergyRyd
double EnergyRyd() const
Definition: transition.h:83
t_hydro::HLineWidth
realnum HLineWidth
Definition: hydrogenic.h:63
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
FeIIPrint
void FeIIPrint(void)
Definition: atom_feii.cpp:1340
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_FeII::feconwlLo
realnum feconwlLo
Definition: atomfeii.h:239
col_str
static double ** col_str
Definition: species2.cpp:29
t_FeII::lgSlow
bool lgSlow
Definition: atomfeii.h:201
Parser
Definition: parser.h:31
dense.h
FeIIIntenZero
void FeIIIntenZero(void)
Definition: atom_feii.cpp:1851
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
trace
t_trace trace
Definition: trace.cpp:5
EmissionProxy::phots
double & phots() const
Definition: emission.h:503
FeIIPunchLevels
void FeIIPunchLevels(FILE *ioPUN)
Definition: atom_feii.cpp:1625
Fe2LevN
TransitionList Fe2LevN("Fe2LevN", &Fe2LevNStates)
cddefines.h
qList::resize
void resize(size_t i)
Definition: quantumstate.h:83
lgCalled
static bool lgCalled
Definition: cddrive.cpp:425
TransitionList::states
qList *& states()
Definition: transition.h:325
FeIIPunchOpticalDepth
void FeIIPunchOpticalDepth(FILE *ioPUN)
Definition: atom_feii.cpp:1673
t_LineSave::ScaleNormLine
double ScaleNormLine
Definition: lines.h:94
t_LineSave::ipass
long int ipass
Definition: lines.h:75
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_rfield::chLineLabel
char ** chLineLabel
Definition: rfield.h:220
POW3
#define POW3
Definition: cddefines.h:936
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
t_FeII::Fe2_large_heat
double Fe2_large_heat
Definition: atomfeii.h:252
EnerLyaProf4
static double EnerLyaProf4
Definition: atom_feii.cpp:163
multi_arr::alloc
void alloc()
Definition: container_classes.h:1116
t_tag_LineSv::SumLine
double SumLine[4]
Definition: lines.h:125
EmissionProxy::ots
double & ots() const
Definition: emission.h:623
FeIICreate
void FeIICreate(void)
Definition: atom_feii.cpp:219
radius.h
FeIIPunchColden
void FeIIPunchColden(FILE *ioPUN)
Definition: atom_feii.cpp:1647
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
ipTFe56
long ipTFe56
Definition: atmdat_readin.cpp:88
FeII
t_FeII FeII
Definition: atomfeii.cpp:5
ipFe2LevN
multi_arr< int, 2 > ipFe2LevN
Definition: taulines.cpp:34
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:141
t_FeII::FeIINRGs
double FeIINRGs[NFE2LEVN]
Definition: atomfeii.h:207
EmissionProxy::pump
double & pump() const
Definition: emission.h:473
FeIIPoint
void FeIIPoint(void)
Definition: atom_feii.cpp:1456
yVector
static double * yVector
Definition: atom_feii.cpp:167
MAX2
#define MAX2
Definition: cddefines.h:782
t_phycon::te_wn
double te_wn
Definition: phycon.h:20
pressure.h
CollisionProxy::col_str
realnum & col_str() const
Definition: collision.h:167
Fe2Energies
realnum * Fe2Energies
Definition: atom_feii.cpp:143
Save1Line
void Save1Line(const TransitionProxy &t, FILE *io, realnum xLimit, long index, realnum DopplerWidth)
Definition: save_do.cpp:4347
getrf_wrapper
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
Definition: thirdparty_lapack.cpp:47
FeII_RT_tau_reset
void FeII_RT_tau_reset(void)
Definition: atom_feii.cpp:1425
CS2SMALL
realnum CS2SMALL
Definition: atom_feii.cpp:78
sPradDat
static realnum *** sPradDat
Definition: atom_feii.cpp:129
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
t_FeII::lgSimulate
bool lgSimulate
Definition: atomfeii.h:219
FeII_RT_Out
void FeII_RT_Out(void)
Definition: atom_feii.cpp:2542
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1228
save.h
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
TransitionProxy::setHi
void setHi(int ipHi) const
Definition: transition.h:404
Parser::getNumberCheck
double getNumberCheck(const char *chDesc)
Definition: parser.cpp:273
FeIIRadPress
double FeIIRadPress(void)
Definition: atom_feii.cpp:2838
EmissionProxy::ColOvTot
double & ColOvTot() const
Definition: emission.h:573
TauLines
TransitionList TauLines("TauLines", &AnonStates)
abscf
double abscf(double gf, double enercm, double gl)
Definition: lines_service.cpp:122
ipTFe26
long ipTFe26
Definition: atmdat_readin.cpp:87
nnPradDat
static long int * nnPradDat
Definition: atom_feii.cpp:124
hydrogenic.h
nFeIIBands
long int nFeIIBands
Definition: atom_feii.cpp:116
EmissionProxy::Pelec_esc
realnum & Pelec_esc() const
Definition: emission.h:533
doppvel.h
t_FeII::ipRedisFcnSubordinate
int ipRedisFcnSubordinate
Definition: atomfeii.h:245
ipPRD
const int ipPRD
Definition: cddefines.h:290
Fe2Coll
static realnum ** Fe2Coll
Definition: atom_feii.cpp:146
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
ipH2p
const int ipH2p
Definition: iso.h:29
FeIIPunchLineStuff
void FeIIPunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: atom_feii.cpp:2820
ipTFe35
long ipTFe35
Definition: atmdat_readin.cpp:88
NFE2LEVN
#define NFE2LEVN
Definition: atomfeii.h:180
FeIIPunDepart
void FeIIPunDepart(FILE *ioPUN, bool lgDoAll)
Definition: atom_feii.cpp:2029
rt.h
powi
double powi(double, long int)
Definition: service.cpp:604
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:664
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
parser.h
COLL_CONST
const UNUSED double COLL_CONST
Definition: physconst.h:229
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:523
physconst.h
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:499
nLine
static long int nLine
Definition: save_line.cpp:288
FeIILevelPops
void FeIILevelPops(void)
Definition: atom_feii.cpp:718
t_FeII::Fe2_large_cool
double Fe2_large_cool
Definition: atomfeii.h:250
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
RT_line_one
void RT_line_one(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
Definition: rt_line_one.cpp:387
AMAT
#define AMAT(I_, J_)
Fe2DepCoef
static double * Fe2DepCoef
Definition: atom_feii.cpp:153
t_FeII::FeIISTWT
double FeIISTWT[NFE2LEVN]
Definition: atomfeii.h:210
t_radius::drad_x_fillfac
double drad_x_fillfac
Definition: radius.h:71
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
ipT191
long ipT191
Definition: atmdat_readin.cpp:89
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
FeIIBandsCreate
STATIC int FeIIBandsCreate(const char chFile[])
Definition: atom_feii.cpp:2293
TRANS_PROB_CONST
const UNUSED double TRANS_PROB_CONST
Definition: physconst.h:237
taulines.h
FeIIZero
void FeIIZero(void)
Definition: atom_feii.cpp:2127
FeII_Bands
realnum ** FeII_Bands
Definition: atom_feii.cpp:106
phycon.h
EnerLyaProf1
static double EnerLyaProf1
Definition: atom_feii.cpp:162
t_FeII::fe2thresh
realnum fe2thresh
Definition: atomfeii.h:236
NLEVPOP
#define NLEVPOP
EmissionProxy::TauTot
realnum & TauTot() const
Definition: emission.h:433
ipH1s
const int ipH1s
Definition: iso.h:27
RT_OTS_AddLine
void RT_OTS_AddLine(double ots, long int ip)
Definition: rt_ots.cpp:402
t_phycon::sqrte
double sqrte
Definition: phycon.h:48
FeIIContCreate
STATIC void FeIIContCreate(double xLamLow, double xLamHigh, long int ncell)
Definition: atom_feii.cpp:2574
AssertFeIIDep
void AssertFeIIDep(double *pred, double *BigError, double *StdDev)
Definition: atom_feii.cpp:2452
FeIIAccel
void FeIIAccel(double *fe2drive)
Definition: atom_feii.cpp:1509
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
amat
static double * amat
Definition: atom_feii.cpp:173
ParseAtomFeII
void ParseAtomFeII(Parser &p)
Definition: atom_feii.cpp:2619
dsexp
double dsexp(double x)
Definition: service.cpp:953
LineSv
LinSv * LineSv
Definition: cdinit.cpp:70
t_thermal::halfte
double halfte
Definition: thermal.h:123
Fe2ColDen
static double * Fe2ColDen
Definition: atom_feii.cpp:157
EmissionProxy::opacity
realnum & opacity() const
Definition: emission.h:593
nFeIIConBins
long int nFeIIConBins
Definition: atom_feii.cpp:119
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
t_thermal::tsq1
double tsq1
Definition: thermal.h:122
ipCRDW
const int ipCRDW
Definition: cddefines.h:294
ipCRD
const int ipCRD
Definition: cddefines.h:292
t_phycon::te
double te
Definition: phycon.h:11
EmissionProxy::FracInwd
realnum & FracInwd() const
Definition: emission.h:463
FeII_RT_Make
void FeII_RT_Make(void)
Definition: atom_feii.cpp:1538
TransitionProxy::Zero
void Zero(void) const
Definition: transition.cpp:505
FeIIPun1Depart
void FeIIPun1Depart(FILE *ioPUN, long int nPUN)
Definition: atom_feii.cpp:2082
RT_line_one_tau_reset
void RT_line_one_tau_reset(const TransitionProxy &t)
Definition: rt_line_one_tau_reset.cpp:12
t_FeII::ddT_Fe2_large_cool
double ddT_Fe2_large_cool
Definition: atomfeii.h:251
Fe2LevelPop
static double * Fe2LevelPop
Definition: atom_feii.cpp:155
xMatrix
static double ** xMatrix
Definition: atom_feii.cpp:170
atomfeii.h
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
save
t_save save
Definition: save.cpp:5
TransitionProxy::outline_resonance
void outline_resonance() const
Definition: transition.cpp:37
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
Fe2SavN
static double ** Fe2SavN
Definition: atom_feii.cpp:134
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17