cloudy  trunk
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
21 typedef 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 
31 STATIC void GetStandardHeLines(void);
32 STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements, double Te );
33 STATIC void DoSatelliteLines( long nelem );
34 
35 static bool lgFirstRun = true;
36 static double CaBDensities[NUMDENS];
37 static double CaBTemps[NUMTEMPS];
38 static long NumLines;
39 static double ****CaBIntensity;
40 static stdLines **CaBLines;
41 
42 void 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  {
74  GetStandardHeLines();
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  {
458  t_iso_sp* sp = &iso_sp[ipHE_LIKE][ipHELIUM];
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 */
466  sp->trans(ipHe3s3S,ipHe2p3P2).Emis().xIntensity()/
467  sp->trans(ipHe3s3S,ipHe2p3P2).EnergyErg() +
468  sp->trans(ipHe3d3D,ipHe2p3P2).Emis().xIntensity()/
469  sp->trans(ipHe3d3D,ipHe2p3P2).EnergyErg() +
470  sp->trans(ipHe4s3S,ipHe2p3P2).Emis().xIntensity()/
471  sp->trans(ipHe4s3S,ipHe2p3P2).EnergyErg() +
472  /* to 2p3P1 */
473  sp->trans(ipHe3s3S,ipHe2p3P1).Emis().xIntensity()/
474  sp->trans(ipHe3s3S,ipHe2p3P1).EnergyErg() +
475  sp->trans(ipHe3d3D,ipHe2p3P1).Emis().xIntensity()/
476  sp->trans(ipHe3d3D,ipHe2p3P1).EnergyErg() +
477  sp->trans(ipHe4s3S,ipHe2p3P1).Emis().xIntensity()/
478  sp->trans(ipHe4s3S,ipHe2p3P1).EnergyErg() +
479  /* to 2p3P0 */
480  sp->trans(ipHe3s3S,ipHe2p3P0).Emis().xIntensity()/
481  sp->trans(ipHe3s3S,ipHe2p3P0).EnergyErg() +
482  sp->trans(ipHe3d3D,ipHe2p3P0).Emis().xIntensity()/
483  sp->trans(ipHe3d3D,ipHe2p3P0).EnergyErg() +
484  sp->trans(ipHe4s3S,ipHe2p3P0).Emis().xIntensity()/
485  sp->trans(ipHe4s3S,ipHe2p3P0).EnergyErg() +
486  /* to 2s3S */
487  sp->trans(ipHe3p3P,ipHe2s3S).Emis().xIntensity()/
488  sp->trans(ipHe3p3P,ipHe2s3S).EnergyErg() +
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 
519 STATIC void GetStandardHeLines(void)
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");
539  cdEXIT(EXIT_FAILURE);
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 );
556  cdEXIT(EXIT_FAILURE);
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");
602  cdEXIT(EXIT_FAILURE);
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 );
690  cdEXIT(EXIT_FAILURE);
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 );
706  cdEXIT(EXIT_FAILURE);
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 
726 STATIC 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. */
760 STATIC 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 }
TransitionProxy::EnergyErg
realnum EnergyErg() const
Definition: transition.h:78
helike.h
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
lines.h
CollisionProxy::cool
double & cool() const
Definition: collision.h:190
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
t_dense::eden
double eden
Definition: dense.h:190
dense
t_dense dense
Definition: dense.cpp:24
elementnames.h
stdLines
Definition: prt_lines_helium.cpp:21
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
EmissionProxy::xIntensity
double & xIntensity() const
Definition: emission.h:483
chLine_LENGTH
#define chLine_LENGTH
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
realnum
float realnum
Definition: cddefines.h:103
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
ipHe3s3S
const int ipHe3s3S
Definition: iso.h:52
STATIC
#define STATIC
Definition: cddefines.h:97
NumLines
static long NumLines
Definition: prt_lines_helium.cpp:38
linint
double linint(const double x[], const double y[], long n, double xval)
Definition: thirdparty_interpolate.cpp:456
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
DoSatelliteLines
STATIC void DoSatelliteLines(long nelem)
Definition: prt_lines_helium.cpp:760
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
L_
#define L_(A_)
Definition: iso.h:21
thirdparty.h
t_isoCTRL::lgRandErrGen
bool lgRandErrGen[NISO]
Definition: iso.h:403
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
NUMDENS
#define NUMDENS
Definition: prt_lines_helium.cpp:20
t_LineSave::WavLNorm
realnum WavLNorm
Definition: lines.h:84
t_iso_sp::TwoNu
vector< two_photon > TwoNu
Definition: iso.h:586
lines_service.h
NUMTEMPS
#define NUMTEMPS
Definition: prt_lines_helium.cpp:19
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:450
t_iso_sp
Definition: iso.h:441
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
TransitionProxy::Coll
CollisionProxy Coll() const
Definition: transition.h:424
atmdat.h
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:423
t_iso_sp::ex
multi_arr< extra_tr, 2 > ex
Definition: iso.h:449
t_isoCTRL::lgFSM
bool lgFSM[NISO]
Definition: iso.h:399
TransitionProxy
Definition: transition.h:23
CollisionProxy::heat
double & heat() const
Definition: collision.h:194
CaBLines
static stdLines ** CaBLines
Definition: prt_lines_helium.cpp:40
LineSave
t_LineSave LineSave
Definition: lines.cpp:5
ipHe3d3D
const int ipHe3d3D
Definition: iso.h:55
hunt_bisect
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:270
CaBDensities
static double CaBDensities[NUMDENS]
Definition: prt_lines_helium.cpp:36
stdLines::label
char label[5]
Definition: prt_lines_helium.cpp:27
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
linadd
void linadd(double xInten, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
Definition: lines_service.cpp:316
dense.h
t_isoCTRL::lgDielRecom
bool lgDielRecom[NISO]
Definition: iso.h:365
trace
t_trace trace
Definition: trace.cpp:5
EmissionProxy::phots
double & phots() const
Definition: emission.h:503
prt
t_prt prt
Definition: prt.cpp:10
cddefines.h
chIonLbl
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)
Definition: transition.cpp:195
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
t_LineSave::ipass
long int ipass
Definition: lines.h:75
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
lagrange
double lagrange(const double x[], const double y[], long n, double xval)
Definition: thirdparty_interpolate.cpp:433
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
IPRAD
#define IPRAD
Definition: iso.h:86
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
t_prt::lgPrnIsoCollapsed
bool lgPrnIsoCollapsed
Definition: prt.h:154
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
EmissionProxy::pump
double & pump() const
Definition: emission.h:473
t_atmdat::CaseBWlHeI
vector< realnum > CaseBWlHeI
Definition: atmdat.h:209
MAX2
#define MAX2
Definition: cddefines.h:782
LIMELM
const int LIMELM
Definition: cddefines.h:258
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:487
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
GetStandardHeLines
STATIC void GetStandardHeLines(void)
Definition: prt_lines_helium.cpp:519
prt.h
EmissionProxy::Pelec_esc
realnum & Pelec_esc() const
Definition: emission.h:533
lgFirstRun
static bool lgFirstRun
Definition: prt_lines_helium.cpp:35
stdLines::ipHi
long int ipHi
Definition: prt_lines_helium.cpp:24
PutLine
void PutLine(const TransitionProxy &t, const char *chComment, const char *chLabelTemp)
Definition: transition.cpp:449
CASEBMAGIC
#define CASEBMAGIC
Definition: helike.h:26
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
min
long min(int a, long b)
Definition: cddefines.h:723
ipHe3p3P
const int ipHe3p3P
Definition: iso.h:54
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:523
physconst.h
SatelliteLines
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:38
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
fixit
void fixit(void)
Definition: service.cpp:991
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
TempInterp2
STATIC double TempInterp2(double *TempArray, double *ValueArray, long NumElements, double Te)
Definition: prt_lines_helium.cpp:726
stdLines::ipLo
long int ipLo
Definition: prt_lines_helium.cpp:25
StuffComment
long int StuffComment(const char *chComment)
Definition: prt_final.cpp:1932
ERG1CM
const UNUSED double ERG1CM
Definition: physconst.h:164
phycon.h
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
lines_helium
void lines_helium(void)
Definition: prt_lines_helium.cpp:42
EmissionProxy::TauTot
realnum & TauTot() const
Definition: emission.h:433
CaBIntensity
static double **** CaBIntensity
Definition: prt_lines_helium.cpp:39
CaBTemps
static double CaBTemps[NUMTEMPS]
Definition: prt_lines_helium.cpp:37
ipSatelliteLines
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:37
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
t_phycon::te
double te
Definition: phycon.h:11
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
EN1RYD
const UNUSED double EN1RYD
Definition: physconst.h:179
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
save
t_save save
Definition: save.cpp:5
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
strchr_s
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1439
Generated by   doxygen 1.8.17