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

Generated for cloudy by doxygen 1.8.3.1