cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
prt_lines_helium.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*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 "helike.h"
8 #include "iso.h"
9 #include "atmdat.h"
10 #include "lines.h"
11 #include "lines_service.h"
12 #include "phycon.h"
13 #include "physconst.h"
14 #include "taulines.h"
15 #include "thirdparty.h"
16 #include "trace.h"
17 
18 #define NUMTEMPS 22
19 
20 typedef struct
21 {
22  /* index for upper and lower levels of line */
23  long int ipHi;
24  long int ipLo;
25 
26  char label[5];
27 
28 } stdLines;
29 
30 STATIC void GetStandardHeLines(void);
31 STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements );
32 STATIC void DoSatelliteLines( long nelem );
33 
34 static bool lgFirstRun = true;
35 static double CaABTemps[NUMTEMPS];
36 static long NumLines;
37 static double ***CaABIntensity;
39 
40 void lines_helium(void)
41 {
42  long ipISO = ipHE_LIKE;
43  long int i, nelem, ipHi, ipLo;
44  char chLabel[5]=" ";
45 
46  long int j;
47 
48  double
49  sum,
50  Pop2_3S,
51  photons_3889_plus_7065 = 0.;
52 
53  DEBUG_ENTRY( "lines_helium()" );
54 
55  if( trace.lgTrace )
56  fprintf( ioQQQ, " prt_lines_helium called\n" );
57 
58  i = StuffComment( "He-like iso-sequence" );
59  linadd( 0., (realnum)i , "####", 'i',
60  " start He-like iso sequence");
61 
62  /* read in Case A and B lines from data file */
63  if( lgFirstRun )
64  {
66  lgFirstRun = false;
67  }
68 
69  /* store labels for all case b HeI lines in case we assert case b
70  * ipass == -1 only counting number of lines, =0, malloc then set wl */
71  static bool lgMustMalloc=true;
72  if( LineSave.ipass == 0 && atmdat.nCaseBHeI>0 && lgMustMalloc )
73  {
74  /* second time through - on ipass=-1 we counted number of lines
75  * atmdat.nCaseBHeI, now create space but only if there are He I lines
76  * this is not done if He is turned off */
78  lgMustMalloc=false;
79  }
80  atmdat.nCaseBHeI = 0;
81 
82  /* this is the main printout, where line intensities are entered into the stack */
83  for( nelem=ipISO; nelem < LIMELM; nelem++ )
84  {
85  if( dense.lgElmtOn[nelem] )
86  {
87  if( nelem == ipHELIUM )
88  {
89  double *qTotEff;
90 
91  /* >>chng 06 aug 17, all of these from _max to _local */
92  /* >>chng 06 dec 21, mistake - change back to _max */
93  qTotEff = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) );
94 
95  qTotEff[0] = 0.;
96  qTotEff[1] = 0.;
97 
98  for( i=ipHe2s3S+1; i<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; i++ )
99  {
100  qTotEff[i] = 0.;
101  for( j = i; j<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; j++ )
102  {
103  /*if( StatesElem[ipHE_LIKE][nelem][i].S == 3 )
104  {*/
105  qTotEff[i] +=
106  Transitions[ipHE_LIKE][nelem][j][ipHe2s3S].Coll.ColUL*dense.EdenHCorr*
107  iso.Boltzmann[ipHE_LIKE][nelem][j][ipHe2s3S] *
108  (double)Transitions[ipHE_LIKE][nelem][j][ipHe2s3S].Hi->g /
109  (double)Transitions[ipHE_LIKE][nelem][j][ipHe2s3S].Lo->g*
110  iso.CascadeProb[ipISO][nelem][j][i];
111  /*}*/
112  }
113  }
114 
115  /* get simple 2^3S pop, assume recombinations in are just 0.75 * case B */
116  Pop2_3S = dense.eden*(0.75*iso.RadRec_caseB[ipHE_LIKE][nelem])/
117  ( Transitions[ipHE_LIKE][nelem][ipHe2s3S][ipHe1s1S].Emis->Aul+ dense.eden*iso.qTot2S[ipISO][nelem]);
118 
119  for( i=0; i< NumLines; i++ )
120  {
121  ipHi = CaABLines[nelem][i].ipHi;
122  ipLo = CaABLines[nelem][i].ipLo;
123 
124  /* >>chng 06 aug 17, from _max to _local */
125  /* >>chng 06 dec 21, mistake - change back to _max */
126  if( ipHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem]*(iso.n_HighestResolved_max[ipHE_LIKE][nelem]+1))
127  {
128  double intens = TempInterp2( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
129  dense.xIonDense[nelem][nelem+1-ipISO]*dense.eden;
130 
131  linadd( intens,
132  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
133  CaABLines[nelem][i].label,'i',
134  "Case B intensity ");
135 
136  if( nMatch("Ca B",CaABLines[nelem][i].label) )
137  {
138  /* all lines to/from 2^3Pj are stored as lines to/from 2^3P1, so make sure this loop never tries to
139  * explicitly consider 2^3P0 or 2^3P2 */
140  ASSERT( ipLo!=ipHe2p3P0 && ipLo!=ipHe2p3P2 );
141  ASSERT( ipHi!=ipHe2p3P0 && ipHi!=ipHe2p3P2 );
142 
143  double totBranch = iso.BranchRatio[ipISO][nelem][ipHi][ipLo];
144  if( ipLo==4 )
145  totBranch += iso.BranchRatio[ipISO][nelem][ipHi][3] + iso.BranchRatio[ipISO][nelem][ipHi][5];
146 
147  if( LineSave.ipass < 0 )
148  ++atmdat.nCaseBHeI;
149  else if( LineSave.ipass == 0 )
150  {
151  /* save wavelengths */
153  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng;
154  ++atmdat.nCaseBHeI;
155  }
156 
157  if( ipHi==4 )
158  {
159  linadd( intens +
160  Pop2_3S*dense.xIonDense[nelem][nelem+1-ipISO]*
161  (
162  qTotEff[ipHe2p3P0]*iso.BranchRatio[ipISO][nelem][ipHe2p3P0][ipLo]+
163  qTotEff[ipHe2p3P1]*iso.BranchRatio[ipISO][nelem][ipHe2p3P1][ipLo]+
164  qTotEff[ipHe2p3P2]*iso.BranchRatio[ipISO][nelem][ipHe2p3P2][ipLo]
165  )*
166  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg,
167  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
168  "+Col",'i',
169  "Case B intensity with collisions included");
170 
171  }
172  else
173  {
174  /* chng 05 dec 14, branching ratio was missing here!
175  * not a big effect because lines with biggest collision
176  * enhancements tend to be dominant decay route from upper level. */
177  linadd( intens +
178  Pop2_3S*qTotEff[ipHi]*dense.xIonDense[nelem][nelem+1-ipISO]*totBranch*
179  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg,
180  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
181  "+Col",'i',
182  "Case B intensity with collisions included");
183  }
184  }
185  }
186  /* Make sure to at least do 4471 */
187  else if( ipLo==ipHe2p3P1 && ipHi==ipHe4d3D && nMatch("Ca B",CaABLines[nelem][i].label) )
188  {
189  double intens = TempInterp2( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
190  dense.xIonDense[nelem][nelem+1-ipISO]*dense.eden;
191 
192  linadd( intens, 4471, CaABLines[nelem][i].label, 'i',
193  "Case B intensity ");
194  }
195 
196  }
197  free( qTotEff );
198  }
199 
200  /* NB NB - low and high must be in this order so that all balmer, paschen,
201  * etc series line up correctly in final printout */
202  /* >>chng 01 jun 13, bring 23P lines back together */
203  /* two photon is special, not a line and no ipCont array index, add here */
204  Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->phots =
205  Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->Aul*
206  StatesElem[ipHE_LIKE][nelem][ipHe2s1S].Pop*
207  Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->Pesc*
208  dense.xIonDense[nelem][nelem+1-ipISO];
209 
210  Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->xIntensity =
211  Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->phots*
212  Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg;
213 
214  if( LineSave.ipass == 0 )
215  {
216  /* chIonLbl is function that generates a null terminated 4 char string, of form "C 2"
217  * the result, chLable, is only used when ipass == 0, can be undefined otherwise */
218  /* total two photon emission */
219  chIonLbl(chLabel, &Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S]);
220  }
221 
222  linadd( Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->xIntensity , 0,chLabel,'r',
223  " two photon continuum ");
224 
225  linadd(
226  StatesElem[ipHE_LIKE][nelem][ipHe2s1S].Pop*
227  dense.xIonDense[nelem][nelem+1-ipISO]*
228  iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*
229  Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg,
230  22, chLabel ,'i',
231  " induced two photon emission ");
232 
233  /* here we will create an entry for the three lines
234  * coming from 2 3P to 1 1S - the entry called TOTL will
235  * appear before the lines of the multiplet */
236  sum = 0.;
237  for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
238  {
239  if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 )
240  continue;
241 
242  sum +=
243  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul*
244  StatesElem[ipHE_LIKE][nelem][i].Pop*
245  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Pesc*
246  dense.xIonDense[nelem][nelem+1-ipISO]*
247  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].EnergyErg;
248  }
249 
250  linadd(sum,Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe1s1S].WLAng,"TOTL",'i' ,
251  " total emission in He-like intercombination lines from 2p3P to ground ");
252 
253  /* now do real permitted lines */
254  for( ipLo=0; ipLo < ipHe2p3P0; ipLo++ )
255  {
256  for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
257  {
258  /* >>chng 01 may 30, do not add fake he-like lines (majority) to line stack */
259  /* >>chng 01 dec 11, use variable for smallest A */
260  if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 )
261  continue;
262 
263  /* recombine fine-structure lines since the energies are
264  * not resolved anyway. */
265  if( iso.lgFSM[ipISO] && ( abs(StatesElem[ipHE_LIKE][nelem][ipHi].l -
266  StatesElem[ipHE_LIKE][nelem][ipLo].l)==1 )
267  && (StatesElem[ipHE_LIKE][nelem][ipLo].l>1)
268  && (StatesElem[ipHE_LIKE][nelem][ipHi].l>1)
269  && ( StatesElem[ipHE_LIKE][nelem][ipHi].n ==
270  StatesElem[ipHE_LIKE][nelem][ipLo].n ) )
271  {
272  /* skip if both singlets. */
273  if( (StatesElem[ipHE_LIKE][nelem][ipHi].S==1)
274  && (StatesElem[ipHE_LIKE][nelem][ipLo].S==1) )
275  {
276  continue;
277  }
278  else if( (StatesElem[ipHE_LIKE][nelem][ipHi].S==3)
279  && (StatesElem[ipHE_LIKE][nelem][ipLo].S==3) )
280  {
281 
282  /* singlet to singlet*/
283  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->phots =
284  ( Transitions[ipHE_LIKE][nelem][ipHi][ipLo+1].Emis->Aul*
285  StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
286  Transitions[ipHE_LIKE][nelem][ipHi][ipLo+1].Emis->Pesc +
287  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->Aul*
288  StatesElem[ipHE_LIKE][nelem][ipHi+1].Pop*
289  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->Pesc )*
290  dense.xIonDense[nelem][nelem+1-ipISO];
291 
292  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->xIntensity =
293  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->phots *
294  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].EnergyErg;
295 
296  PutLine(&Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1],
297  " ");
298 
299  /* triplet to triplet */
300  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots =
301  ( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
302  StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
303  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc +
304  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo].Emis->Aul*
305  StatesElem[ipHE_LIKE][nelem][ipHi+1].Pop*
306  Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo].Emis->Pesc )*
307  dense.xIonDense[nelem][nelem+1-ipISO];
308 
309  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity =
310  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *
311  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
312 
313  PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
314  " ");
315  }
316  }
317 
318  else if( ipLo==ipHe2s3S && ipHi == ipHe2p3P0 )
319  {
320  /* here we will create an entry for the three lines
321  * coming from 2 3P to 2 3S - the entry called TOTL will
322  * appear before the lines of the multiplet
323  * for He I this is 10830 */
324 
325  realnum av_wl = 0.;
326  sum = 0.;
327  for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
328  {
329  sum +=
330  Transitions[ipHE_LIKE][nelem][i][ipLo].Emis->Aul*
331  StatesElem[ipHE_LIKE][nelem][i].Pop*
332  Transitions[ipHE_LIKE][nelem][i][ipLo].Emis->Pesc*
333  dense.xIonDense[nelem][nelem+1-ipISO]*
334  Transitions[ipHE_LIKE][nelem][i][ipLo].EnergyErg;
335  av_wl += Transitions[ipHE_LIKE][nelem][i][ipLo].WLAng;
336  }
337  av_wl /= 3.;
338 # if 0
339  {
340 # include "elementnames.h"
341 # include "prt.h"
342  fprintf(ioQQQ,"DEBUG 2P - 2S multiplet wl %s ",
343  elementnames.chElementSym[nelem] );
344  prt_wl( ioQQQ , av_wl );
345  fprintf(ioQQQ,"\n" );
346  }
347 # endif
348 
349  linadd(sum,av_wl,"TOTL",'i',
350  "total emission in He-like lines, use average of three line wavelengths " );
351 
352  /* also add this with the regular label, so it is correctly picked up by assert case-b command */
353  linadd(sum,av_wl,chLabel,'i',
354  "total emission in He-like lines, use average of three line wavelengths " );
355 
356  /*>>chng 05 sep 8, added the following so that the component
357  * from ipHe2p3P0 is printed, in addition to the total. */
358 
359  /* find number of photons escaping cloud */
360  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots =
361  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
362  StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
363  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
364  dense.xIonDense[nelem][nelem+1-ipISO];
365 
366  /* now find line intensity */
367  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity =
368  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
369  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
370 
371  if( iso.lgRandErrGen[ipISO] )
372  {
373  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
374  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
375  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *=
376  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
377  }
378  PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
379  " ");
380  }
381 
382  else
383  {
384 
385  /* find number of photons escaping cloud */
386  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots =
387  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
388  StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
389  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
390  dense.xIonDense[nelem][nelem+1-ipISO];
391 
392  /* now find line intensity */
393  /* >>chng 01 jan 15, put cast double to force double evaluation */
394  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity =
395  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
396  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
397 
398  if( iso.lgRandErrGen[ipISO] )
399  {
400  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
401  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
402  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *=
403  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
404  }
405 
406  /*
407  fprintf(ioQQQ,"1 loop %li %li %.1f\n", ipLo,ipHi,
408  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng ); */
409  PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
410  "total intensity of He-like line");
411  {
412  /* option to print particulars of some line when called
413  * a prettier print statement is near where chSpin is defined below*/
414  enum {DEBUG_LOC=false};
415  if( DEBUG_LOC )
416  {
417  if( nelem==1 && ipLo==0 && ipHi==1 )
418  {
419  fprintf(ioQQQ,"he1 626 %.2e %.2e \n",
420  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauIn,
421  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauTot
422  );
423  }
424  }
425  }
426  }
427  }
428  }
429 
430  /* this sum will bring together the three lines going to J levels within 23P */
431  for( ipHi=ipHe2p3P2+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
432  {
433  double sumcool , sumheat ,
434  save , savecool , saveheat;
435 
436  sum = 0;
437  sumcool = 0.;
438  sumheat = 0.;
439  for( ipLo=ipHe2p3P0; ipLo <= ipHe2p3P2; ++ipLo )
440  {
441  if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
442  continue;
443 
444  /* find number of photons escaping cloud */
445  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots =
446  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
447  StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
448  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
449  dense.xIonDense[nelem][nelem+1-ipISO];
450 
451  /* now find line intensity */
452  /* >>chng 01 jan 15, put cast double to force double evaluation */
453  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity =
454  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
455  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
456 
457  if( iso.lgRandErrGen[ipISO] )
458  {
459  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
460  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
461  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *=
462  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
463  }
464 
465  sumcool += Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.cool;
466  sumheat += Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.heat;
467  sum += Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity;
468  }
469 
470  /* skip non-radiative lines */
471  if( Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].ipCont < 1 )
472  continue;
473 
474  /* this will enter .xIntensity into the line stack */
475  save = Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Emis->xIntensity;
476  savecool = Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.cool;
477  saveheat = Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.heat;
478 
479  Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Emis->xIntensity = sum;
480  Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.cool = sumcool;
481  Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.heat = sumheat;
482 
483  /*fprintf(ioQQQ,"2 loop %li %li %.1f\n", ipHe2p3P2,ipHi,
484  Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].WLAng );*/
485  PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2],
486  "predicted line, all processes included");
487 
488  Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Emis->xIntensity = save;
489  Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.cool = savecool;
490  Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.heat = saveheat;
491  }
492  for( ipLo=ipHe2p3P2+1; ipLo < iso.numPrintLevels[ipHE_LIKE][nelem]-1; ipLo++ )
493  {
494  for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
495  {
496  /* skip non-radiative lines */
497  if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 )
498  continue;
499 
500  /* find number of photons escaping cloud */
501  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots =
502  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
503  StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
504  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
505  dense.xIonDense[nelem][nelem+1-ipISO];
506 
507  /* now find line intensity */
508  /* >>chng 01 jan 15, put cast double to force double evaluation */
509  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity =
510  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
511  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
512 
513  if( iso.lgRandErrGen[ipISO] )
514  {
515  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
516  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
517  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *=
518  iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
519  }
520 
521  /* this will enter .xIntensity into the line stack */
522  PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
523  "predicted line, all processes included");
524  }
525  }
526 
527  /* Now put the satellite lines in */
528  if( iso.lgDielRecom[ipISO] )
529  DoSatelliteLines(nelem);
530  }
531  }
532 
535  {
536  /* For info only, add the total photon flux in 3889 and 7065,
537  * and 3188, 4713, and 5876. */
538  photons_3889_plus_7065 =
539  /* to 2p3P2 */
540  Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P2].Emis->xIntensity/
542  Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P2].Emis->xIntensity/
544  Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P2].Emis->xIntensity/
546  /* to 2p3P1 */
547  Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P1].Emis->xIntensity/
549  Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P1].Emis->xIntensity/
551  Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P1].Emis->xIntensity/
553  /* to 2p3P0 */
554  Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P0].Emis->xIntensity/
556  Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P0].Emis->xIntensity/
558  Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P0].Emis->xIntensity/
560  /* to 2s3S */
561  Transitions[ipHE_LIKE][ipHELIUM][ipHe3p3P][ipHe2s3S].Emis->xIntensity/
563  Transitions[ipHE_LIKE][ipHELIUM][ipHe4p3P][ipHe2s3S].Emis->xIntensity/
565 
566  long upperIndexofH8 = iso.QuantumNumbers2Index[ipH_LIKE][ipHYDROGEN][8][1][2];
567 
568  /* Add in photon flux of H8 3889 */
569  photons_3889_plus_7065 +=
570  Transitions[ipH_LIKE][ipHYDROGEN][upperIndexofH8][1].Emis->xIntensity/
571  Transitions[ipH_LIKE][ipHYDROGEN][upperIndexofH8][1].EnergyErg;
572 
573  /* now multiply by ergs of normalization line, so that relative flux of
574  * this line will be ratio of photon fluxes. */
575  photons_3889_plus_7065 *= (ERG1CM*1.e8)/LineSave.WavLNorm;
576  linadd( photons_3889_plus_7065, 3889., "Pho+", 'i',
577  "photon sum given in Porter et al. 2007 (astro-ph/0611579)");
578  }
579 
580  /* ====================================================
581  * end helium
582  * ====================================================*/
583 
584  if( trace.lgTrace )
585  {
586  fprintf( ioQQQ, " lines_helium returns\n" );
587  }
588  return;
589 }
590 
591 
593 {
594  FILE *ioDATA;
595  bool lgEOL, lgHIT;
596  long i, i1, i2, j, nelem;
597 
598 # define chLine_LENGTH 1000
599  char chLine[chLine_LENGTH];
600 
601  DEBUG_ENTRY( "GetStandardHeLines()" );
602 
603  if( trace.lgTrace )
604  fprintf( ioQQQ," lines_helium opening he1_case_ab.dat\n");
605 
606  ioDATA = open_data( "he1_case_ab.dat", "r" );
607 
608  /* check that magic number is ok */
609  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
610  {
611  fprintf( ioQQQ, " lines_helium could not read first line of he1_case_ab.dat.\n");
612  cdEXIT(EXIT_FAILURE);
613  }
614  i = 1;
615  /* first number is magic number, second is number of lines in file */
616  i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
617  i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
618  NumLines = i2;
619 
620  /* the following is to check the numbers that appear at the start of he1_case_ab.dat */
621  if( i1 !=CASEABMAGIC )
622  {
623  fprintf( ioQQQ,
624  " lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
625  fprintf( ioQQQ,
626  " lines_helium: I expected to find the number %i and got %li instead.\n" ,
627  CASEABMAGIC, i1 );
628  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
629  cdEXIT(EXIT_FAILURE);
630  }
631 
632  /* get the array of temperatures */
633  lgHIT = false;
634  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
635  {
636  /* only look at lines without '#' in first col */
637  if( chLine[0] != '#')
638  {
639  lgHIT = true;
640  j = 0;
641  lgEOL = false;
642  i = 1;
643  while( !lgEOL && j < NUMTEMPS)
644  {
645  CaABTemps[j] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
646  ++j;
647  }
648  break;
649  }
650  }
651 
652  if( !lgHIT )
653  {
654  fprintf( ioQQQ, " lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
655  cdEXIT(EXIT_FAILURE);
656  }
657 
658  /* create space for array of values, if not already done */
659  CaABIntensity = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
660  /* create space for array of values, if not already done */
661  CaABLines = (stdLines **)MALLOC(sizeof(stdLines *)*(unsigned)LIMELM );
662 
663  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
664  {
668  if( nelem != ipHELIUM )
669  continue;
670 
671  /* only grab core for elements that are turned on */
672  if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
673  {
674  /* create space for array of values, if not already done */
675  CaABIntensity[nelem] = (double **)MALLOC(sizeof(double *)*(unsigned)(i2) );
676  CaABLines[nelem] = (stdLines *)MALLOC(sizeof(stdLines )*(unsigned)(i2) );
677 
678  /* avoid allocating 0 bytes, some OS return NULL pointer, PvH
679  CaABIntensity[nelem][0] = NULL;*/
680  for( j = 0; j < i2; ++j )
681  {
682  CaABIntensity[nelem][j] = (double *)MALLOC(sizeof(double)*(unsigned)NUMTEMPS );
683 
684  CaABLines[nelem][j].ipHi = -1;
685  CaABLines[nelem][j].ipLo = -1;
686  strcpy( CaABLines[nelem][j].label , " " );
687 
688  for( i=0; i<NUMTEMPS; ++i )
689  {
690  CaABIntensity[nelem][j][i] = 0.;
691  }
692  }
693  }
694  }
695 
696  /* now read in the case A and B line data */
697  lgHIT = false;
698  nelem = ipHELIUM;
699  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
700  {
701  static long line = 0;
702  char *chTemp;
703 
704  /* only look at lines without '#' in first col */
705  if( (chLine[0] == ' ') || (chLine[0]=='\n') )
706  break;
707  if( chLine[0] != '#')
708  {
709  lgHIT = true;
710 
711  /* get lower and upper level index */
712  j = 1;
713  /* the first number is the wavelength, which is not used
714  * for anything, but must skip over it. */
715  i1 = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
716  CaABLines[nelem][line].ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
717  CaABLines[nelem][line].ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
718 
719  ASSERT( CaABLines[nelem][line].ipLo < CaABLines[nelem][line].ipHi );
720  /*if( CaABLines[nelem][line].ipHi >= iso.numLevels_max[ipHE_LIKE][nelem] )
721  continue;*/
722 
723  chTemp = chLine;
724  /* skip over 4 tabs to start of data */
725  for( i=0; i<3; ++i )
726  {
727  if( (chTemp = strchr( chTemp, '\t' )) == NULL )
728  {
729  fprintf( ioQQQ, " lines_helium no init case A and B\n" );
730  cdEXIT(EXIT_FAILURE);
731  }
732  ++chTemp;
733  }
734 
735  strncpy( CaABLines[nelem][line].label, chTemp , 4 );
736  CaABLines[nelem][line].label[4] = 0;
737 
738  for( i=0; i<NUMTEMPS; ++i )
739  {
740  double b;
741  if( (chTemp = strchr( chTemp, '\t' )) == NULL )
742  {
743  fprintf( ioQQQ, " lines_helium could not scan case A and B lines, current indices: %li %li\n",
744  CaABLines[nelem][line].ipHi,
745  CaABLines[nelem][line].ipLo );
746  cdEXIT(EXIT_FAILURE);
747  }
748  ++chTemp;
749  sscanf( chTemp , "%le" , &b );
750  CaABIntensity[nelem][line][i] = pow(10.,b);
751  }
752  line++;
753  }
754  }
755 
756  /* close the data file */
757  fclose( ioDATA );
758  return;
759 }
760 
762 STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements )
763 {
764  long int ipTe=-1;
765  double rate = 0.;
766  long i0;
767 
768  DEBUG_ENTRY( "TempInterp2()" );
769 
770  ASSERT( fabs( 1. - (double)phycon.alogte/log10(phycon.te) ) < 0.0001 );
771 
772  /* te totally unknown */
773  if( phycon.alogte <= TempArray[0] )
774  {
775  return ValueArray[0];
776  }
777  else if( phycon.alogte >= TempArray[NumElements-1] )
778  {
779  return ValueArray[NumElements-1];
780  }
781 
782  /* now search for temperature */
783  ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte );
784 
785  ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
786 
787  /* Do a four-point interpolation */
788  const int ORDER = 3; /* order of the fitting polynomial */
789  i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
790  rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte );
791 
792  return rate;
793 }
794 
796 /* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002 */
797 /* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001. */
798 STATIC void DoSatelliteLines( long nelem )
799 {
800  long ipISO = ipHE_LIKE;
801 
802  DEBUG_ENTRY( "DoSatelliteLines()" );
803 
804  ASSERT( iso.lgDielRecom[ipISO] && dense.lgElmtOn[nelem] );
805 
806  for( long i=0; i<iso.numPrintLevels[ipISO][nelem]; i++ )
807  {
808  double dr_rate = iso.DielecRecomb[ipISO][nelem][i];
809 
810  SatelliteLines[ipISO][nelem][i].Emis->phots =
811  dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
812 
813  SatelliteLines[ipISO][nelem][i].Emis->xIntensity =
814  SatelliteLines[ipISO][nelem][i].Emis->phots * ERG1CM * SatelliteLines[ipISO][nelem][i].EnergyWN;
815 
816  PutLine( &SatelliteLines[ipISO][nelem][i], "iso satellite line", "Sate" );
817 
818  //linadd( SatelliteLines[ipISO][nelem][i].Emis->xIntensity ,
819  // SatelliteLines[ipISO][nelem][i].WLAng, "Sate",'i', "iso satellite line" );
820  }
821 
822  return;
823 }

Generated for cloudy by doxygen 1.8.1.2