cloudy trunk
Loading...
Searching...
No Matches
prt_lines_helium.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/*lines_helium put He-like iso sequence into line intensity stack */
4/*TempInterp interpolates on a grid of values to produce predicted value at current Te.*/
5#include "cddefines.h"
6#include "dense.h"
7#include "prt.h"
8#include "helike.h"
9#include "iso.h"
10#include "atmdat.h"
11#include "lines.h"
12#include "lines_service.h"
13#include "phycon.h"
14#include "physconst.h"
15#include "taulines.h"
16#include "thirdparty.h"
17#include "trace.h"
18
19#define NUMTEMPS 21
20#define NUMDENS 14
21typedef struct
22{
23 /* index for upper and lower levels of line */
24 long int ipHi;
25 long int ipLo;
26
27 char label[5];
28
29} stdLines;
30
31STATIC void GetStandardHeLines(void);
32STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements, double Te );
33STATIC void DoSatelliteLines( long nelem );
34
35static bool lgFirstRun = true;
36static double CaBDensities[NUMDENS];
37static double CaBTemps[NUMTEMPS];
38static long NumLines;
39static double ****CaBIntensity;
41
42void lines_helium(void)
43{
44 long ipISO = ipHE_LIKE;
45 long int i, nelem, ipHi, ipLo;
46 char chLabel[5]=" ";
47
48 double
49 sum,
50 photons_3889_plus_7065 = 0.;
51
52 DEBUG_ENTRY( "lines_helium()" );
53
54 if( trace.lgTrace )
55 fprintf( ioQQQ, " prt_lines_helium called\n" );
56
57 // this can be changed with the atom levels command but must be at
58 // least 3.
59 ASSERT( iso_sp[ipHE_LIKE][ipHELIUM].n_HighestResolved_max >= 3 );
60
61 i = StuffComment( "He-like iso-sequence" );
62 linadd( 0., (realnum)i , "####", 'i',
63 " start He-like iso sequence");
64
65 linadd(MAX2(0.,iso_sp[ipHE_LIKE][ipHELIUM].xLineTotCool),0,"He1c",'c',
66 " total collisional cooling due to all HeI lines ");
67
68 linadd(MAX2(0.,-iso_sp[ipHE_LIKE][ipHELIUM].xLineTotCool),0,"He1h",'h' ,
69 " total collisional heating due to all HeI lines ");
70
71 /* read in Case A and B lines from data file */
72 if( lgFirstRun )
73 {
75 lgFirstRun = false;
76 }
77
78 /* this is the main printout, where line intensities are entered into the stack */
79 for( nelem=ipISO; nelem < LIMELM; nelem++ )
80 {
81 if( dense.lgElmtOn[nelem] )
82 {
83 t_iso_sp* sp = &iso_sp[ipHE_LIKE][nelem];
84
85 ASSERT( sp->n_HighestResolved_max >= 3 );
86
87 if( nelem == ipHELIUM )
88 {
89 for( i=0; i< NumLines; i++ )
90 {
91 ipHi = CaBLines[nelem][i].ipHi;
92 ipLo = CaBLines[nelem][i].ipLo;
93 double intens_at_Te[NUMDENS];
94 for( long ipDens = 0; ipDens < NUMDENS; ++ipDens )
95 intens_at_Te[ipDens] = TempInterp2( CaBTemps, CaBIntensity[nelem][i][ipDens], NUMTEMPS, phycon.te );
96 double intens = linint( CaBDensities, intens_at_Te, NUMDENS, log10(dense.eden) );
97 intens = pow( 10., intens ) * dense.xIonDense[nelem][nelem+1-ipISO]*dense.eden;
98 ASSERT( intens >= 0. );
99 linadd( intens, atmdat.CaseBWlHeI[i], CaBLines[nelem][i].label,'i', "Case B intensity ");
100 }
101 }
102
103 // add two-photon details here
104 if( LineSave.ipass == 0 )
105 {
106 /* chIonLbl is function that generates a null terminated 4 char string, of form "C 2"
107 * the result, chLable, is only used when ipass == 0, can be undefined otherwise */
108 chIonLbl(chLabel, nelem+1, nelem+1-ipISO);
109 }
110 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
111 {
112 fixit(); // This was multiplied by Pesc when treated as a line, now what? Only used for printout?
113 fixit(); // below should be 'i' instead of 'r' ?
114 linadd( tnu->AulTotal * tnu->E2nu * EN1RYD * (*tnu->Pop),
115 0, chLabel, 'r',
116 " two photon continuum ");
117
118 linadd( tnu->induc_dn * tnu->E2nu * EN1RYD * (*tnu->Pop),
119 22, chLabel ,'i',
120 " induced two photon emission ");
121 }
122
123 /* here we will create an entry for the three lines
124 * coming from 2 3P to 1 1S - the entry called TOTL will
125 * appear before the lines of the multiplet */
126 sum = 0.;
127 for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
128 {
129 if( sp->trans(i,ipHe1s1S).ipCont() <= 0 )
130 continue;
131
132 sum +=
133 sp->trans(i,ipHe1s1S).Emis().Aul()*
134 sp->st[i].Pop()*
135 (sp->trans(i,ipHe1s1S).Emis().Pesc() +
136 sp->trans(i,ipHe1s1S).Emis().Pelec_esc() ) *
137 sp->trans(i,ipHe1s1S).EnergyErg();
138 }
139
140 linadd(sum,sp->trans(ipHe2p3P1,ipHe1s1S).WLAng(),"TOTL",'i' ,
141 " total emission in He-like intercombination lines from 2p3P to ground ");
142
143 /* set number of levels we want to print, first is default,
144 * only print real levels, second is set with "print line
145 * iso collapsed" command */
146 long int nLoop = sp->numLevels_max - sp->nCollapsed_max;
147 if( prt.lgPrnIsoCollapsed )
148 nLoop = sp->numLevels_max;
149
150 /* now do real permitted lines */
151 /* NB NB - low and high must be in this order so that all balmer, paschen,
152 * etc series line up correctly in final printout */
153 /* >>chng 01 jun 13, bring 23P lines back together */
154 for( ipLo=0; ipLo < ipHe2p3P0; ipLo++ )
155 {
156 vector<long> EnterTheseLast;
157 for( ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
158 {
159 /* >>chng 01 may 30, do not add fake he-like lines (majority) to line stack */
160 /* >>chng 01 dec 11, use variable for smallest A */
161 if( sp->trans(ipHi,ipLo).ipCont() < 1 )
162 continue;
163
164 /* recombine fine-structure lines since the energies are
165 * not resolved anyway. */
166 if( iso_ctrl.lgFSM[ipISO] && ( abs(sp->st[ipHi].l() -
167 sp->st[ipLo].l())==1 )
168 && (sp->st[ipLo].l()>1)
169 && (sp->st[ipHi].l()>1)
170 && ( sp->st[ipHi].n() ==
171 sp->st[ipLo].n() ) )
172 {
173 /* skip if both singlets. */
174 if( (sp->st[ipHi].S()==1)
175 && (sp->st[ipLo].S()==1) )
176 {
177 continue;
178 }
179 else if( (sp->st[ipHi].S()==3)
180 && (sp->st[ipLo].S()==3) )
181 {
182
183 /* singlet to singlet*/
184 sp->trans(ipHi+1,ipLo+1).Emis().phots() =
185 sp->trans(ipHi,ipLo+1).Emis().Aul()*
186 sp->st[ipHi].Pop()*
187 (sp->trans(ipHi,ipLo+1).Emis().Pesc() +
188 sp->trans(ipHi,ipLo+1).Emis().Pelec_esc() ) +
189 sp->trans(ipHi+1,ipLo+1).Emis().Aul()*
190 sp->st[ipHi+1].Pop()*
191 (sp->trans(ipHi+1,ipLo+1).Emis().Pesc() +
192 sp->trans(ipHi+1,ipLo+1).Emis().Pelec_esc() );
193
194 sp->trans(ipHi+1,ipLo+1).Emis().xIntensity() =
195 sp->trans(ipHi+1,ipLo+1).Emis().phots() *
196 sp->trans(ipHi+1,ipLo+1).EnergyErg();
197
198 PutLine(sp->trans(ipHi+1,ipLo+1), " ");
199
200 /* triplet to triplet */
201 sp->trans(ipHi,ipLo).Emis().phots() =
202 sp->trans(ipHi,ipLo).Emis().Aul()*
203 sp->st[ipHi].Pop()*
204 (sp->trans(ipHi,ipLo).Emis().Pesc() +
205 sp->trans(ipHi,ipLo).Emis().Pelec_esc() ) +
206 sp->trans(ipHi+1,ipLo).Emis().Aul()*
207 sp->st[ipHi+1].Pop()*
208 (sp->trans(ipHi+1,ipLo).Emis().Pesc() +
209 sp->trans(ipHi+1,ipLo).Emis().Pelec_esc() );
210
211 sp->trans(ipHi,ipLo).Emis().xIntensity() =
212 sp->trans(ipHi,ipLo).Emis().phots() *
213 sp->trans(ipHi,ipLo).EnergyErg();
214
215 PutLine(sp->trans(ipHi,ipLo), " ");
216 }
217 }
218
219 else if( ipLo==ipHe2s3S && ipHi == ipHe2p3P0 )
220 {
221 /* here we will create an entry for the three lines
222 * coming from 2 3P to 2 3S - the entry called TOTL will
223 * appear before the lines of the multiplet
224 * for He I this is 10830 */
225
226 realnum av_wl = 0.;
227 sum = 0.;
228 for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
229 {
230 sum +=
231 sp->trans(i,ipLo).Emis().Aul()*
232 sp->st[i].Pop()*
233 (sp->trans(i,ipLo).Emis().Pesc() +
234 sp->trans(i,ipLo).Emis().Pelec_esc() )*
235 sp->trans(i,ipLo).EnergyErg();
236 av_wl += sp->trans(i,ipLo).WLAng();
237 }
238 av_wl /= 3.;
239# if 0
240 {
241# include "elementnames.h"
242# include "prt.h"
243 fprintf(ioQQQ,"DEBUG 2P - 2S multiplet wl %s ",
244 elementnames.chElementSym[nelem] );
245 prt_wl( ioQQQ , av_wl );
246 fprintf(ioQQQ,"\n" );
247 }
248# endif
249
250 linadd(sum,av_wl,"TOTL",'i',
251 "total emission in He-like lines, use average of three line wavelengths " );
252
253 /* also add this with the regular label, so it is correctly picked up by assert case-b command */
254 linadd(sum,av_wl,chLabel,'i',
255 "total emission in He-like lines, use average of three line wavelengths " );
256
257 /*>>chng 05 sep 8, added the following so that the component
258 * from ipHe2p3P0 is printed, in addition to the total. */
259
260 /* find number of photons escaping cloud */
261 sp->trans(ipHi,ipLo).Emis().phots() =
262 sp->trans(ipHi,ipLo).Emis().Aul()*
263 sp->st[ipHi].Pop()*
264 (sp->trans(ipHi,ipLo).Emis().Pesc() +
265 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
266
267 /* now find line intensity */
268 sp->trans(ipHi,ipLo).Emis().xIntensity() =
269 sp->trans(ipHi,ipLo).Emis().phots()*
270 sp->trans(ipHi,ipLo).EnergyErg();
271
272 if( iso_ctrl.lgRandErrGen[ipISO] )
273 {
274 sp->trans(ipHi,ipLo).Emis().phots() *=
275 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
276 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
277 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
278 }
279 PutLine(sp->trans(ipHi,ipLo), " ");
280 }
281
282 else
283 {
284 /* find number of photons escaping cloud */
285 sp->trans(ipHi,ipLo).Emis().phots() =
286 sp->trans(ipHi,ipLo).Emis().Aul()*
287 sp->st[ipHi].Pop()*
288 (sp->trans(ipHi,ipLo).Emis().Pesc() +
289 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
290
291 /* now find line intensity */
292 /* >>chng 01 jan 15, put cast double to force double evaluation */
293 sp->trans(ipHi,ipLo).Emis().xIntensity() =
294 sp->trans(ipHi,ipLo).Emis().phots()*
295 sp->trans(ipHi,ipLo).EnergyErg();
296
297 if( iso_ctrl.lgRandErrGen[ipISO] )
298 {
299 sp->trans(ipHi,ipLo).Emis().phots() *=
300 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
301 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
302 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
303 }
304
305 if( abs( L_(ipHi) - L_(ipLo) ) != 1 )
306 {
307 EnterTheseLast.push_back( ipHi );
308 continue;
309 }
310
311 /*
312 fprintf(ioQQQ,"1 loop %li %li %.1f\n", ipLo,ipHi,
313 sp->trans(ipHi,ipLo).WLAng() ); */
314 PutLine(sp->trans(ipHi,ipLo),
315 "total intensity of He-like line");
316 {
317 /* option to print particulars of some line when called
318 * a prettier print statement is near where chSpin is defined below*/
319 enum {DEBUG_LOC=false};
320 if( DEBUG_LOC )
321 {
322 if( nelem==1 && ipLo==0 && ipHi==1 )
323 {
324 fprintf(ioQQQ,"he1 626 %.2e %.2e \n",
325 sp->trans(ipHi,ipLo).Emis().TauIn(),
326 sp->trans(ipHi,ipLo).Emis().TauTot()
327 );
328 }
329 }
330 }
331 }
332 }
333
334 // enter these lines last because they are generally weaker quadrupole transitions.
335 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
336 PutLine( sp->trans(*it,ipLo),
337 "predicted line, all processes included");
338 }
339
340 /* this sum will bring together the three lines going to J levels within 23P */
341 for( ipHi=ipHe2p3P2+1; ipHi < nLoop; ipHi++ )
342 {
343 double sumcool , sumheat ,
344 save , savecool , saveheat;
345
346 sum = 0;
347 sumcool = 0.;
348 sumheat = 0.;
349 for( ipLo=ipHe2p3P0; ipLo <= ipHe2p3P2; ++ipLo )
350 {
351 if( sp->trans(ipHi,ipLo).ipCont() <= 0 )
352 continue;
353
354 /* find number of photons escaping cloud */
355 sp->trans(ipHi,ipLo).Emis().phots() =
356 sp->trans(ipHi,ipLo).Emis().Aul()*
357 sp->st[ipHi].Pop()*
358 (sp->trans(ipHi,ipLo).Emis().Pesc() +
359 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
360
361 /* now find line intensity */
362 /* >>chng 01 jan 15, put cast double to force double evaluation */
363 sp->trans(ipHi,ipLo).Emis().xIntensity() =
364 sp->trans(ipHi,ipLo).Emis().phots()*
365 sp->trans(ipHi,ipLo).EnergyErg();
366
367 if( iso_ctrl.lgRandErrGen[ipISO] )
368 {
369 sp->trans(ipHi,ipLo).Emis().phots() *=
370 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
371 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
372 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
373 }
374
375 sumcool += sp->trans(ipHi,ipLo).Coll().cool();
376 sumheat += sp->trans(ipHi,ipLo).Coll().heat();
377 sum += sp->trans(ipHi,ipLo).Emis().xIntensity();
378 }
379
380 /* skip non-radiative lines */
381 if( sp->trans(ipHi,ipHe2p3P2).ipCont() < 1 )
382 continue;
383
384 /* this will enter .xIntensity() into the line stack */
385 save = sp->trans(ipHi,ipHe2p3P2).Emis().xIntensity();
386 savecool = sp->trans(ipHi,ipHe2p3P2).Coll().cool();
387 saveheat = sp->trans(ipHi,ipHe2p3P2).Coll().heat();
388
389 sp->trans(ipHi,ipHe2p3P2).Emis().xIntensity() = sum;
390 sp->trans(ipHi,ipHe2p3P2).Coll().cool() = sumcool;
391 sp->trans(ipHi,ipHe2p3P2).Coll().heat() = sumheat;
392
393 /*fprintf(ioQQQ,"2 loop %li %li %.1f\n", ipHe2p3P2,ipHi,
394 sp->trans(ipHi,ipHe2p3P2).WLAng() );*/
395 PutLine(sp->trans(ipHi,ipHe2p3P2),
396 "predicted line, all processes included");
397
398 sp->trans(ipHi,ipHe2p3P2).Emis().xIntensity() = save;
399 sp->trans(ipHi,ipHe2p3P2).Coll().cool() = savecool;
400 sp->trans(ipHi,ipHe2p3P2).Coll().heat() = saveheat;
401 }
402 for( ipLo=ipHe2p3P2+1; ipLo < nLoop-1; ipLo++ )
403 {
404 vector<long> EnterTheseLast;
405 for( ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
406 {
407 /* skip non-radiative lines */
408 if( sp->trans(ipHi,ipLo).ipCont() < 1 )
409 continue;
410
411 /* find number of photons escaping cloud */
412 sp->trans(ipHi,ipLo).Emis().phots() =
413 sp->trans(ipHi,ipLo).Emis().Aul()*
414 sp->st[ipHi].Pop()*
415 (sp->trans(ipHi,ipLo).Emis().Pesc() +
416 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
417
418 /* now find line intensity */
419 /* >>chng 01 jan 15, put cast double to force double evaluation */
420 sp->trans(ipHi,ipLo).Emis().xIntensity() =
421 sp->trans(ipHi,ipLo).Emis().phots()*
422 sp->trans(ipHi,ipLo).EnergyErg();
423
424 if( iso_ctrl.lgRandErrGen[ipISO] )
425 {
426 sp->trans(ipHi,ipLo).Emis().phots() *=
427 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
428 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
429 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
430 }
431
432 if( abs( L_(ipHi) - L_(ipLo) ) != 1 )
433 {
434 EnterTheseLast.push_back( ipHi );
435 continue;
436 }
437
438 /* this will enter .xIntensity() into the line stack */
439 PutLine(sp->trans(ipHi,ipLo),
440 "predicted line, all processes included");
441 }
442
443 // enter these lines last because they are generally weaker quadrupole transitions.
444 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
445 PutLine(sp->trans(*it,ipLo),
446 "predicted line, all processes included");
447 }
448
449 /* Now put the satellite lines in */
450 if( iso_ctrl.lgDielRecom[ipISO] )
451 DoSatelliteLines(nelem);
452 }
453 }
454
455 if( iso_sp[ipHE_LIKE][ipHELIUM].n_HighestResolved_max >= 4 &&
456 ( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_max + iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_max ) >=8 )
457 {
459 const long ipHe4s3S = iso_sp[ipHE_LIKE][ipHELIUM].QuantumNumbers2Index[4][0][3];
460 const long ipHe4p3P = iso_sp[ipHE_LIKE][ipHELIUM].QuantumNumbers2Index[4][1][3];
461
462 /* For info only, add the total photon flux in 3889 and 7065,
463 * and 3188, 4713, and 5876. */
464 photons_3889_plus_7065 =
465 /* to 2p3P2 */
470 sp->trans(ipHe4s3S,ipHe2p3P2).Emis().xIntensity()/
471 sp->trans(ipHe4s3S,ipHe2p3P2).EnergyErg() +
472 /* to 2p3P1 */
477 sp->trans(ipHe4s3S,ipHe2p3P1).Emis().xIntensity()/
478 sp->trans(ipHe4s3S,ipHe2p3P1).EnergyErg() +
479 /* to 2p3P0 */
484 sp->trans(ipHe4s3S,ipHe2p3P0).Emis().xIntensity()/
485 sp->trans(ipHe4s3S,ipHe2p3P0).EnergyErg() +
486 /* to 2s3S */
489 sp->trans(ipHe4p3P,ipHe2s3S).Emis().xIntensity()/
490 sp->trans(ipHe4p3P,ipHe2s3S).EnergyErg();
491
492 long upperIndexofH8 = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[8][1][2];
493
494 /* Add in photon flux of H8 3889 */
495 photons_3889_plus_7065 +=
496 iso_sp[ipH_LIKE][ipHYDROGEN].trans(upperIndexofH8,1).Emis().xIntensity()/
497 iso_sp[ipH_LIKE][ipHYDROGEN].trans(upperIndexofH8,1).EnergyErg();
498
499 /* now multiply by ergs of normalization line, so that relative flux of
500 * this line will be ratio of photon fluxes. */
501 if( LineSave.WavLNorm > 0 )
502 photons_3889_plus_7065 *= (ERG1CM*1.e8)/LineSave.WavLNorm;
503 linadd( photons_3889_plus_7065, 3889., "Pho+", 'i',
504 "photon sum given in Porter et al. 2007 (astro-ph/0611579)");
505 }
506
507 /* ====================================================
508 * end helium
509 * ====================================================*/
510
511 if( trace.lgTrace )
512 {
513 fprintf( ioQQQ, " lines_helium returns\n" );
514 }
515 return;
516}
517
518
520{
521 FILE *ioDATA;
522 bool lgEOL, lgHIT;
523 long i, i1, i2;
524
525# define chLine_LENGTH 1000
526 char chLine[chLine_LENGTH];
527
528 DEBUG_ENTRY( "GetStandardHeLines()" );
529
530 if( trace.lgTrace )
531 fprintf( ioQQQ," lines_helium opening he1_case_b.dat\n");
532
533 ioDATA = open_data( "he1_case_b.dat", "r" );
534
535 /* check that magic number is ok */
536 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
537 {
538 fprintf( ioQQQ, " lines_helium could not read first line of he1_case_b.dat.\n");
540 }
541 i = 1;
542 /* first number is magic number, second is number of lines in file */
543 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
544 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
545 NumLines = i2;
546
547 /* the following is to check the numbers that appear at the start of he1_case_b.dat */
548 if( i1 !=CASEBMAGIC )
549 {
550 fprintf( ioQQQ,
551 " lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
552 fprintf( ioQQQ,
553 " lines_helium: I expected to find the number %i and got %li instead.\n" ,
554 CASEBMAGIC, i1 );
555 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
557 }
558
559 /* get the array of densities */
560 lgHIT = false;
561 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
562 {
563 /* only look at lines without '#' in first col */
564 if( chLine[0] != '#')
565 {
566 lgHIT = true;
567 long j = 0;
568 lgEOL = false;
569 i = 1;
570 while( !lgEOL && j < NUMDENS)
571 {
572 CaBDensities[j] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
573 ++j;
574 }
575 break;
576 }
577 }
578
579 /* get the array of temperatures */
580 lgHIT = false;
581 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
582 {
583 /* only look at lines without '#' in first col */
584 if( chLine[0] != '#')
585 {
586 lgHIT = true;
587 long j = 0;
588 lgEOL = false;
589 i = 1;
590 while( !lgEOL && j < NUMTEMPS)
591 {
592 CaBTemps[j] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
593 ++j;
594 }
595 break;
596 }
597 }
598
599 if( !lgHIT )
600 {
601 fprintf( ioQQQ, " lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
603 }
604
605 /* create space for array of values, if not already done */
606 CaBIntensity = (double ****)MALLOC(sizeof(double ***)*(unsigned)LIMELM );
607 /* create space for array of values, if not already done */
608 CaBLines = (stdLines **)MALLOC(sizeof(stdLines *)*(unsigned)LIMELM );
609
610 for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem )
611 {
615 if( nelem != ipHELIUM )
616 continue;
617
618 /* only grab core for elements that are turned on */
619 if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
620 {
621 /* create space for array of values, if not already done */
622 CaBIntensity[nelem] = (double ***)MALLOC(sizeof(double **)*(unsigned)(i2) );
623 CaBLines[nelem] = (stdLines *)MALLOC(sizeof(stdLines )*(unsigned)(i2) );
624
625 /* avoid allocating 0 bytes, some OS return NULL pointer, PvH
626 CaBIntensity[nelem][0] = NULL;*/
627 for( long j = 0; j < i2; ++j )
628 {
629 CaBIntensity[nelem][j] = (double **)MALLOC(sizeof(double*)*(unsigned)NUMDENS );
630
631 CaBLines[nelem][j].ipHi = -1;
632 CaBLines[nelem][j].ipLo = -1;
633 strcpy( CaBLines[nelem][j].label , " " );
634
635 for( long k = 0; k < NUMDENS; ++k )
636 {
637 CaBIntensity[nelem][j][k] = (double *)MALLOC(sizeof(double)*(unsigned)NUMTEMPS );
638 for( long l = 0; l < NUMTEMPS; ++l )
639 {
640 CaBIntensity[nelem][j][k][l] = 0.;
641 }
642 }
643 }
644 }
645 }
646
647 /* now read in the case B line data */
648 lgHIT = false;
649 long nelem = ipHELIUM;
650 int line = 0;
651 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
652 {
653 char *chTemp = NULL;
654
655 /* only look at lines without '#' in first col */
656 if( (chLine[0] == ' ') || (chLine[0]=='\n') )
657 break;
658 if( chLine[0] != '#')
659 {
660 lgHIT = true;
661
662 /* get lower and upper level index */
663 long j = 1;
664 // the first number is the wavelength, which is only used if the
665 // model atom is too small to include this transition
666 realnum wl = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
667 long ipLo = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
668 long ipHi = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
669 CaBLines[nelem][line].ipLo = ipLo;
670 CaBLines[nelem][line].ipHi = ipHi;
671
672 ASSERT( CaBLines[nelem][line].ipLo < CaBLines[nelem][line].ipHi );
673
674 strncpy( CaBLines[nelem][line].label, "Ca B" , 4 );
675 CaBLines[nelem][line].label[4] = 0;
676
677 t_iso_sp* sp = &iso_sp[ipHE_LIKE][nelem];
678 if( ipHi < sp->numLevels_max - sp->nCollapsed_max )
679 atmdat.CaseBWlHeI.push_back( sp->trans(ipHi,ipLo).WLAng() );
680 else
681 atmdat.CaseBWlHeI.push_back( wl );
682
683 for( long ipDens = 0; ipDens < NUMDENS; ++ipDens )
684 {
685 if( read_whole_line( chLine, (int)sizeof(chLine) , ioDATA ) == NULL )
686 {
687 fprintf( ioQQQ, " lines_helium could not scan case B lines, current indices: %li %li\n",
688 CaBLines[nelem][line].ipHi,
689 CaBLines[nelem][line].ipLo );
691 }
692
693 chTemp = chLine;
694 j = 1;
695 long den = (long)FFmtRead(chTemp,&j,sizeof(chTemp),&lgEOL);
696 ASSERT( den == ipDens + 1 );
697
698 for( long ipTe = 0; ipTe < NUMTEMPS; ++ipTe )
699 {
700 double b;
701 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
702 {
703 fprintf( ioQQQ, " lines_helium could not scan case B lines, current indices: %li %li\n",
704 CaBLines[nelem][line].ipHi,
705 CaBLines[nelem][line].ipLo );
707 }
708 ++chTemp;
709 sscanf( chTemp, "%le" , &b );
710 CaBIntensity[nelem][line][ipDens][ipTe] = b;
711 }
712 }
713 line++;
714 }
715 }
716
717 ASSERT( line == NumLines );
718 ASSERT( atmdat.CaseBWlHeI.size() == (unsigned)line );
719
720 /* close the data file */
721 fclose( ioDATA );
722 return;
723}
724
726STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements, double Te )
727{
728 long int ipTe=-1;
729 double rate = 0.;
730 long i0;
731
732 DEBUG_ENTRY( "TempInterp2()" );
733
734 /* te totally unknown */
735 if( Te <= TempArray[0] )
736 {
737 return ValueArray[0] + log10( TempArray[0] ) - log10( Te );
738 }
739 else if( Te >= TempArray[NumElements-1] )
740 {
741 return ValueArray[NumElements-1];
742 }
743
744 /* now search for temperature */
745 ipTe = hunt_bisect( TempArray, NumElements, Te );
746
747 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
748
749 /* Do a four-point interpolation */
750 const int ORDER = 3; /* order of the fitting polynomial */
751 i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
752 rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, Te );
753
754 return rate;
755}
756
758/* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002 */
759/* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001. */
760STATIC void DoSatelliteLines( long nelem )
761{
762 long ipISO = ipHE_LIKE;
763
764 DEBUG_ENTRY( "DoSatelliteLines()" );
765
766 ASSERT( iso_ctrl.lgDielRecom[ipISO] && dense.lgElmtOn[nelem] );
767
768 for( long i=0; i < iso_sp[ipISO][nelem].numLevels_max; i++ )
769 {
770 double dr_rate = iso_sp[ipISO][nelem].fb[i].DielecRecomb;
771 TransitionProxy tr = SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][i]];
772
773 tr.Emis().phots() =
774 dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
775
776 tr.Emis().xIntensity() =
777 tr.Emis().phots() * ERG1CM * tr.EnergyWN();
778 tr.Emis().pump() = 0.;
779
780 PutLine( tr, "iso satellite line" );
781 }
782
783 return;
784}
t_atmdat atmdat
Definition atmdat.cpp:6
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MALLOC(exp)
Definition cddefines.h:501
const int LIMELM
Definition cddefines.h:258
#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
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const char * strchr_s(const char *s, int c)
Definition cddefines.h:1439
long max(int a, long b)
Definition cddefines.h:775
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
void fixit(void)
Definition service.cpp:991
double & cool() const
Definition collision.h:190
double & heat() const
Definition collision.h:194
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
double & pump() const
Definition emission.h:473
double & phots() const
Definition emission.h:503
realnum & TauTot() const
Definition emission.h:433
CollisionProxy Coll() const
Definition transition.h:424
realnum & WLAng() const
Definition transition.h:429
realnum EnergyErg() const
Definition transition.h:78
realnum & EnergyWN() const
Definition transition.h:438
long & ipCont() const
Definition transition.h:450
EmissionList::reference Emis() const
Definition transition.h:408
long int numLevels_max
Definition iso.h:493
multi_arr< extra_tr, 2 > ex
Definition iso.h:449
TransitionProxy trans(const long ipHi, const long ipLo)
Definition iso.h:444
long int n_HighestResolved_max
Definition iso.h:505
long int nCollapsed_max
Definition iso.h:487
vector< two_photon > TwoNu
Definition iso.h:586
qList st
Definition iso.h:453
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
#define CASEBMAGIC
Definition helike.h:26
#define chLine_LENGTH
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
#define IPRAD
Definition iso.h:86
const int ipHe2s3S
Definition iso.h:44
const int ipHe3p3P
Definition iso.h:54
const int ipHE_LIKE
Definition iso.h:63
const int ipHe2p3P1
Definition iso.h:47
const int ipHe3d3D
Definition iso.h:55
const int ipHe2p3P0
Definition iso.h:46
const int ipH_LIKE
Definition iso.h:62
const int ipHe3s3S
Definition iso.h:52
const int ipHe2p3P2
Definition iso.h:48
#define L_(A_)
Definition iso.h:21
t_LineSave LineSave
Definition lines.cpp:5
long int StuffComment(const char *chComment)
void linadd(double xInten, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double ERG1CM
Definition physconst.h:164
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
t_prt prt
Definition prt.cpp:10
static double **** CaBIntensity
#define NUMDENS
void lines_helium(void)
static bool lgFirstRun
static double CaBTemps[NUMTEMPS]
static stdLines ** CaBLines
STATIC void DoSatelliteLines(long nelem)
static double CaBDensities[NUMDENS]
STATIC double TempInterp2(double *TempArray, double *ValueArray, long NumElements, double Te)
STATIC void GetStandardHeLines(void)
#define NUMTEMPS
static long NumLines
t_save save
Definition save.cpp:5
vector< vector< TransitionList > > SatelliteLines
Definition taulines.cpp:38
multi_arr< int, 3 > ipSatelliteLines
Definition taulines.cpp:37
double linint(const double x[], const double y[], long n, double xval)
double lagrange(const double x[], const double y[], long n, double xval)
long hunt_bisect(const T x[], long n, T xval)
Definition thirdparty.h:270
t_trace trace
Definition trace.cpp:5
void PutLine(const TransitionProxy &t, const char *chComment, const char *chLabelTemp)
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)