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