cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cont_createpointers.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 /*ContCreatePointers set up pointers for lines and continua called by cloudy after input read in
4  * and continuum mesh has been set */
5 /*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
6 /*ipShells assign continuum energy pointers to shells for all atoms,
7  * called by ContCreatePointers */
8 /*LimitSh sets upper energy limit to subshell integrations */
9 /*ContBandsCreate - read set of continuum bands to enter total emission into line stack*/
10 #include "cddefines.h"
11 #include "physconst.h"
12 #include "lines_service.h"
13 #include "iso.h"
14 #include "secondaries.h"
15 #include "taulines.h"
16 #include "elementnames.h"
17 #include "ionbal.h"
18 #include "rt.h"
19 #include "opacity.h"
20 #include "yield.h"
21 #include "dense.h"
22 #include "he.h"
23 #include "fe.h"
24 #include "rfield.h"
25 #include "oxy.h"
26 #include "atomfeii.h"
27 #include "atoms.h"
28 #include "trace.h"
29 #include "hmi.h"
30 #include "heavy.h"
31 #include "predcont.h"
32 #include "atmdat.h"
33 #include "ipoint.h"
34 #include "h2.h"
35 #include "continuum.h"
36 
37 /*LimitSh sets upper energy limit to subshell integrations */
38 STATIC long LimitSh(long int ion,
39  long int nshell,
40  long int nelem);
41 
42 STATIC void ipShells(
43  /* nelem is the atomic number on the C scale, Li is 2 */
44  long int nelem);
45 
46 /*ContBandsCreate - read set of continuum bands to enter total emission into line*/
48  /* chFile is optional filename, if void then use default bands,
49  * if not void then use file specified,
50  * return value is 0 for success, 1 for failure */
51  const char chFile[] );
52 
53 /* upper level for two-photon emission in H and He iso sequences */
54 #define TwoS (1+ipISO)
55 
56 /*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
57 STATIC void fiddle(long int ipnt,
58  double exact);
59 
61 {
62  char chLab[5];
63  long int
64  i,
65  ion,
66  ipHi,
67  ipLo,
68  ipISO,
69  iWL_Ang,
70  j,
71  nelem,
72  nshells;
73  double energy,
74  xnew;
75  /* counter to say whether pointers have ever been evaluated */
76  static int nCalled = 0;
77 
78  DEBUG_ENTRY( "ContCreatePointers()" );
79 
80  /* create the hydrogen atom for this core load, routine creates space then zeros it out
81  * on first call, on second and later calls it only zeros things out */
82  iso_create();
83 
84  /* create internal static variables needed to do the H2 molecule */
85  H2_Create();
86 
87  /* nCalled is local static variable defined 0 when defined.
88  * it is incremented below, so that space only allocated one time per coreload. */
89  if( nCalled > 0 )
90  {
91  if( trace.lgTrace )
92  fprintf( ioQQQ, " ContCreatePointers called, not evaluating.\n" );
93  return;
94  }
95  else
96  {
97  if( trace.lgTrace )
98  fprintf( ioQQQ, " ContCreatePointers called first time.\n" );
99  ++nCalled;
100  }
101 
102  for( i=0; i < rfield.nupper; i++ )
103  {
104  /* this is array of labels for lines and continua, set to blanks at first */
105  strcpy( rfield.chContLabel[i], " ");
106  strcpy( rfield.chLineLabel[i], " ");
107  }
108 
109  /* we will generate a set of array indices to ionization edges for
110  * the first thirty elements. First set all array indices to
111  * totally bogus values so we will crash if misused */
112  for( nelem=0; nelem<LIMELM; ++nelem )
113  {
114  if( dense.lgElmtOn[nelem] )
115  {
116  for( ion=0; ion<LIMELM; ++ion )
117  {
118  for( nshells=0; nshells<7; ++nshells )
119  {
120  for( j=0; j<3; ++j )
121  {
122  opac.ipElement[nelem][ion][nshells][j] = INT_MIN;
123  }
124  }
125  }
126  }
127  }
128 
129  /* pointer to excited state of O+2 */
130  opac.ipo3exc[0] = ipContEnergy(3.855,"O3ex");
131 
132  /* main hydrogenic arrays - THIS OCCURS TWICE!! H and He here, then the
133  * remaining hydrogenic species near the bottom. This is so that H and HeII get
134  * their labels stuffed into the arrays, and the rest of the hydrogenic series
135  * get whatever is left over after the level 1 lines.
136  * to find second block, search on "ipZ=2" */
137  /* NB note that no test for H or He being on exists here - we will always
138  * define the H and He arrays even when He is off, since we need to
139  * know where the he edges are for the bookkeeping that occurs in continuum
140  * binning routines */
141  /* this loop is over H, He-like only - DO NOT CHANGE TO NISO */
142  for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
143  {
144  /* this will be over HI, HeII, then HeI only */
145  for( nelem=ipISO; nelem < 2; nelem++ )
146  {
147  /* generate label for this ion */
148  sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
149 
150  /* array index for continuum edges for ground */
151  iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
152  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
153  {
154  /* array index for continuum edges for excited levels */
155  iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
156 
157  /* define all line array indices */
158  for( ipLo=0; ipLo < ipHi; ipLo++ )
159  {
160  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
161  continue;
162 
163  /* some lines have negative or zero energy */
164  /* >>chng 03 apr 22, this check was if less than or equal to zero,
165  * changed to lowest energy point so that ultra low energy transitions are
166  * not considered. */
167  if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
168  continue;
169 
170  /* some energies are negative for inverted levels */
171  Transitions[ipISO][nelem][ipHi][ipLo].ipCont =
172  ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
173  iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
174  Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine =
175  ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
176  /* check that energy scales are the same, to within energy resolution of arrays */
177 # ifndef NDEBUG
178  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine > 0 )
179  {
180  realnum anuCoarse = rfield.anu[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
181  realnum anuFine = rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine];
182  realnum widCoarse = rfield.widflx[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
183  realnum widFine = anuFine - rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine-1];
184  realnum width = MAX2( widFine , widCoarse );
185  /* NB - can only assert more than width of coarse cell
186  * since close to ionization limit, coarse index is
187  * kept below next ionization edge
188  * >>chng 05 mar 02, pre factor below had been 1.5, chng to 2
189  * tripped near H grnd photo threshold */
190  ASSERT( fabs(anuCoarse - anuFine) / anuCoarse <
191  2.*width/anuCoarse);
192  }
193 # endif
194  }
195  }/* ipHi loop */
196  }/* nelem loop */
197  }/* ipISO */
198 
199  /* need to increase the cell for the HeII balmer continuum by one, so that
200  * hydrogen Lyman continuum will pick it up. */
201  nelem = ipHELIUM;
202  /* copy label over to next higher cell */
203  if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "He 2" ) == 0)
204  {
205  strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]],
206  rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] );
207  /* set previous spot to blank */
208  strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , " ");
209  }
210  else if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "H 1" ) == 0)
211  {
212  /* set previous spot to blank */
213  strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]] , "He 2");
214  }
215  else
216  {
217  fprintf(ioQQQ," insanity heii pointer fix contcreatepointers\n");
218  }
219  /* finally increment the two HeII pointers so that they are above the Lyman continuum */
220  ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
221  ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2p];
222 
223  /* we will define an array of either 1 or 0 to show whether photooelectron
224  * energy is large enough to produce secondary ionizations
225  * below 100eV no secondary ionization - this is the threshold*/
226  secondaries.ipSecIon = ipoint(7.353);
227 
228  /* this is highest energy where k-shell opacities are counted
229  * can be adjusted with "set kshell" command */
231 
232  /* pointers for molecules
233  * H2+ dissociation energy 2.647 eV but cs small below 0.638 Ryd */
234  opac.ih2pnt[0] = ipContEnergy(0.502,"H2+ ");
235  opac.ih2pnt[1] = ipoint(1.03);
236 
237  /* pointers for most prominent PAH features
238  * energies given to ipContEnergy are only to put lave in the right place
239  * wavelengths are rough observed values of blends
240  * 7.6 microns */
241  i = ipContEnergy(0.0117, "PAH " );
242 
243  /* feature near 6.2 microns */
244  i = ipContEnergy(0.0147, "PAH " );
245 
246  /* 3.3 microns */
247  i = ipContEnergy(0.028, "PAH " );
248 
249  /* 11.2 microns */
250  i = ipContEnergy(0.0080, "PAH " );
251 
252  /* 12.3 microns */
253  i = ipContEnergy(0.0077, "PAH " );
254 
255  /* 13.5 microns */
256  i = ipContEnergy(0.0069, "PAH " );
257 
258 
259  /* fix pointers for hydrogen and helium */
260 
261  /* pointer to Bowen 374A resonance line */
262  he.ip374 = ipLineEnergy(1.92,"He 2",0);
263  he.ip660 = ipLineEnergy(1.38,"He 2",0);
264 
265  /* set up series of continuum pointers to be used in routine lines to
266  * enter continuum points into line array */
267  for( i=0; i < NPREDCONT; i++ )
268  {
269  /* EnrPredCont contains array of energy points, as set in zerologic.c */
270  ipPredCont[i] = ipoint(EnrPredCont[i]) - 1;
271  }
272 
273  /* pointer to energy defining effect x-ray column */
274  rt.ipxry = ipoint(73.5);
275 
276  /* pointer to Hminus edge at 0.754eV */
277  hmi.iphmin = ipContEnergy(0.05544,"H- ");
278 
279  /* pointer to threshold for H2 photoionization at 15.4 eV */
280  opac.ipH2_photo_thresh = ipContEnergy( 15.4/EVRYD , "H2 ");
281 
282  hmi.iheh1 = ipoint(1.6);
283  hmi.iheh2 = ipoint(2.3);
284 
285  /* pointer to carbon k-shell ionization */
286  opac.ipCKshell = ipoint(20.6);
287 
288  /* pointer to threshold for pair production */
289  opac.ippr = ipoint(7.51155e4) + 1;
290 
291  /* pointer to x-ray - gamma ray bound; 100 kev */
293 
295 
296  /* make low energy edges of some cells exact,
297  * index is on c scale
298  * 0.99946 is correct energy of hydrogen in inf mass ryd */
301  /* confirm that labels are in correct location */
302  ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1], "H 1" ) ==0 );
303  ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1], "H 1" ) ==0 );
304 
306  ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1], "He 2" ) ==0 );
307 
308  /* pointer to excited state of O+2 */
309  fiddle(opac.ipo3exc[0]-1,3.855);
310 
311  /* now fix widflx array so that it is correct */
312  for( i=1; i<rfield.nupper-1; ++i )
313  {
314  rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
315  }
316 
317  /* these are indices for centers of B and V filters,
318  * taken from table on page 202 of Allen, AQ, 3rd ed */
319  /* the B filter array offset */
321  /* the V filter array offset */
323 
324  /* these are the lower and upper bounds for the G0 radiation field
325  * used by Tielens & Hollenbach in their PDR work */
326  rfield.ipG0_TH85_lo = ipoint( 6.0 / EVRYD );
327  rfield.ipG0_TH85_hi = ipoint( 13.6 / EVRYD );
328 
329  /* this is the limits for Draine & Bertoldi Habing field */
330  rfield.ipG0_DB96_lo = ipoint( RYDLAM / 1110. );
331  rfield.ipG0_DB96_hi = ipoint( RYDLAM / 912. );
332 
333  /* this is special form of G0 that could be used in some future version, for now,
334  * use default, TH85 */
335  rfield.ipG0_spec_lo = ipoint( 6.0 / EVRYD );
336  rfield.ipG0_spec_hi = ipoint( RYDLAM / 912. );
337 
338  /* this index is to 1000A to obtain the extinction at 1000A */
339  rfield.ip1000A = ipoint( RYDLAM / 1000. );
340 
341  /* now save current form of array */
342  for( i=0; i < rfield.nupper; i++ )
343  {
344  rfield.AnuOrg[i] = rfield.anu[i];
345  rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
346  }
347 
348  /* following order of elements is in roughly decreasing abundance
349  * when ipShells gets a cell for the valence IP it does so through
350  * ipContEnergy, which makes sure that no two ions get the same cell
351  * earliest elements have most precise ip mapping */
352 
353  /* set up shell pointers for hydrogen */
354  nelem = ipHYDROGEN;
355  ion = 0;
356 
357  /* the number of shells */
358  Heavy.nsShells[nelem][0] = 1;
359 
360  /*pointer to ionization threshold in energy array*/
361  Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
362  opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
363 
364  /* upper limit to energy integration */
365  opac.ipElement[nelem][ion][0][1] = rfield.nupper;
366 
367  if( dense.lgElmtOn[ipHELIUM] )
368  {
369  /* set up shell pointers for helium */
370  nelem = ipHELIUM;
371  ion = 0;
372 
373  /* the number of shells */
374  Heavy.nsShells[nelem][0] = 1;
375 
376  /*pointer to ionization threshold in energy array*/
377  Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
378  opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
379 
380  /* upper limit to energy integration */
381  opac.ipElement[nelem][ion][0][1] = rfield.nupper;
382 
383  /* (hydrogenic) helium ion */
384  ion = 1;
385  /* the number of shells */
386  Heavy.nsShells[nelem][1] = 1;
387 
388  /*pointer to ionization threshold in energy array*/
389  Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
390  opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
391 
392  /* upper limit to energy integration */
393  opac.ipElement[nelem][ion][0][1] = rfield.nupper;
394  }
395 
396  /* check that ionization potential of neutral carbon valence shell is
397  * positive */
398  ASSERT( t_ADfA::Inst().ph1(2,5,5,0) > 0. );
399 
400  /* now fill in all sub-shell ionization array indices for elements heavier than He,
401  * this must be done after previous loop on iso.ipIsoLevNIonCon[ipH_LIKE] since hydrogenic species use
402  * iso.ipIsoLevNIonCon[ipH_LIKE] rather than ipoint in getting array index within continuum array */
403  for( i=NISO; i<LIMELM; ++i )
404  {
405  if( dense.lgElmtOn[i])
406  {
407  /* i is the atomic number on the c scale, 2 for Li */
408  ipShells(i);
409  }
410  }
411 
412  /* most of these are set in ipShells, but not h-like or he-like, so do these here */
413  Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = t_ADfA::Inst().ph1(0,0,ipHYDROGEN,0)/EVRYD* 0.9998787;
414  Heavy.Valence_IP_Ryd[ipHELIUM][0] = t_ADfA::Inst().ph1(0,1,ipHELIUM,0)/EVRYD* 0.9998787;
415  Heavy.Valence_IP_Ryd[ipHELIUM][1] = t_ADfA::Inst().ph1(0,0,ipHELIUM,0)/EVRYD* 0.9998787;
416  for( nelem=2; nelem<LIMELM; ++nelem )
417  {
418  Heavy.Valence_IP_Ryd[nelem][nelem-1] = t_ADfA::Inst().ph1(0,1,nelem,0)/EVRYD* 0.9998787;
419  Heavy.Valence_IP_Ryd[nelem][nelem] = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
420  if( dense.lgElmtOn[nelem])
421  {
422  /* now confirm that all are properly set */
423  for( j=0; j<=nelem; ++j )
424  {
425  ASSERT( Heavy.Valence_IP_Ryd[nelem][j] > 0.05 );
426  }
427  for( j=0; j<nelem; ++j )
428  {
429  ASSERT( Heavy.Valence_IP_Ryd[nelem][j] < Heavy.Valence_IP_Ryd[nelem][j+1]);
430  }
431  }
432  }
433 
434  /* array indices to bound Compton electron recoil ionization of all elements */
435  for( nelem=0; nelem<LIMELM; ++nelem )
436  {
437  if( dense.lgElmtOn[nelem])
438  {
439  for( ion=0; ion<nelem+1; ++ion )
440  {
441  /* this is the threshold energy to Compton ionize valence shell electrons */
442  energy = sqrt( Heavy.Valence_IP_Ryd[nelem][ion] * EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT ) / EN1RYD;
443  /* the array index for this energy */
444  ionbal.ipCompRecoil[nelem][ion] = ipoint( energy );
445  }
446  }
447  }
448 
449  /* oxygen pointers for excited states
450  * IO3 is pointer to O++ exc state, is set above */
451  oxy.i2d = ipoint(1.242);
452  oxy.i2p = ipoint(1.367);
453  opac.ipo1exc[0] = ipContEnergy(0.856,"O1ex");
454  opac.ipo1exc[1] = ipoint(2.0);
455 
456  /* upper limit for excited state photoionization
457  * do not use ipContEnergy since only upper limit */
458  opac.ipo3exc[1] = ipoint(5.0);
459 
460  /* upper level of 4363 */
461  opac.ipo3exc3[0] = ipContEnergy(3.646,"O3ex");
462  opac.ipo3exc3[1] = ipoint(5.0);
463 
464  /* following are various pointers for OI - Ly-beta pump problem
465  * these are delta energies for Boltzmann factors */
470  atoms.ipoiex[0] = ipoint(0.7005);
471  atoms.ipoiex[1] = ipoint(0.10791);
472  atoms.ipoiex[2] = ipoint(0.06925);
473  atoms.ipoiex[3] = ipoint(0.01151);
474  atoms.ipoiex[4] = ipoint(0.01999);
475 
476  /* >>chng 97 jan 27, move nitrogen after oxygen so that O gets the
477  * most accurate pointers
478  * Nitrogen
479  * in1(1) is thresh for photo from excited state */
480  opac.in1[0] = ipContEnergy(0.893,"N1ex");
481 
482  /* upper limit */
483  opac.in1[1] = ipoint(2.);
484 
486  {
487  fprintf( ioQQQ, " ContCreatePointers:%ld energy cells used. N(1R):%4ld N(1.8):%4ld N(4Ryd):%4ld N(O3)%4ld N(x-ray):%5ld N(rcoil)%5ld\n",
488  rfield.nupper,
492  opac.ipo3exc[0],
493  opac.ipCKshell,
495 
496 
497  fprintf( ioQQQ, " ContCreatePointers: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n",
499 
500  fprintf( ioQQQ, " ContCreatePointers: H pointers;" );
501  for( i=0; i <= 6; i++ )
502  {
503  fprintf( ioQQQ, "%4ld%4ld", i, iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i] );
504  }
505  fprintf( ioQQQ, "\n" );
506 
507  fprintf( ioQQQ, " ContCreatePointers: Oxy pnters;" );
508 
509  for( i=1; i <= 8; i++ )
510  {
511  fprintf( ioQQQ, "%4ld%4ld", i, Heavy.ipHeavy[ipOXYGEN][i-1] );
512  }
513  fprintf( ioQQQ, "\n" );
514  }
515 
516  /* Magnesium
517  * following is energy for phot of MG+ from exc state producing 2798 */
518  opac.ipmgex = ipoint(0.779);
519 
520  /* lower, upper edges of Ca+ excited term photoionizaiton */
521  opac.ica2ex[0] = ipContEnergy(0.72,"Ca2x");
522  opac.ica2ex[1] = ipoint(1.);
523 
524  /* set up factors and pointers for Fe continuum fluorescence */
525  fe.ipfe10 = ipoint(2.605);
526 
527  /* following is WL(CM)**2/(8PI) * conv fac for RYD to NU *A21 */
528  fe.pfe10 = (realnum)(2.00e-18/rfield.widflx[fe.ipfe10-1]);
529 
530  /* this is 353 pump, f=0.032 */
531  fe.pfe11a = (realnum)(4.54e-19/rfield.widflx[fe.ipfe10-1]);
532 
533  /* this is 341.1 f=0.012 */
534  fe.pfe11b = (realnum)(2.75e-19/rfield.widflx[fe.ipfe10-1]);
535  fe.pfe14 = (realnum)(1.15e-18/rfield.widflx[fe.ipfe10-1]);
536 
537  /* set up energy pointers for line optical depth arrays
538  * this also increments flux, sets other parameters for lines */
539 
540  /* level 1 heavy elements line array */
541  for( i=1; i <= nLevel1; i++ )
542  {
543  /* put null terminated line label into chLab using line array*/
544  chIonLbl(chLab, &TauLines[i]);
545 
546  TauLines[i].ipCont =
547  ipLineEnergy(TauLines[i].EnergyWN * WAVNRYD, chLab ,0);
548  TauLines[i].Emis->ipFine =
549  ipFineCont(TauLines[i].EnergyWN * WAVNRYD );
550  /* for debugging pointer - index on f not c scale,
551  * this will find all lines that entered a given cell
552  if( TauLines[i].ipCont==603 )
553  fprintf(ioQQQ,"( level1 %s\n", chLab);*/
554 
555  if( TauLines[i].Emis->gf > 0. )
556  {
557  /* the gf value was entered
558  * derive the A, call to function is gf, wl (A), g_up */
559  TauLines[i].Emis->Aul =
560  (realnum)(eina(TauLines[i].Emis->gf,
561  TauLines[i].EnergyWN,TauLines[i].Hi->g));
562  }
563  else if( TauLines[i].Emis->Aul > 0. )
564  {
565  /* the Einstein A value was entered
566  * derive the gf, call to function is A, wl (A), g_up */
567  TauLines[i].Emis->gf =
568  (realnum)(GetGF(TauLines[i].Emis->Aul,
569  TauLines[i].EnergyWN,TauLines[i].Hi->g));
570  }
571  else
572  {
573  fprintf( ioQQQ, " level 1 line does not have valid gf or A\n" );
574  fprintf( ioQQQ, " This is ContCreatePointers\n" );
575  cdEXIT(EXIT_FAILURE);
576  }
577 
578  /* used to get damping constant */
579  TauLines[i].Emis->dampXvel =
580  (realnum)(TauLines[i].Emis->Aul / TauLines[i].EnergyWN/PI4);
581 
582  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
583  TauLines[i].Emis->opacity =
584  (realnum)(abscf(TauLines[i].Emis->gf,
585  TauLines[i].EnergyWN,TauLines[i].Lo->g));
586  /*fprintf(ioQQQ,"TauLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
587  i,TauLines[i].Emis->opacity , TauLines[i].Emis->gf , TauLines[i].Emis->Aul ,TauLines[i].EnergyWN,TauLines[i].Lo->g);*/
588 
589  /* excitation energy of transition in degrees kelvin */
590  TauLines[i].EnergyK =
592 
593  /* energy of photon in ergs */
594  TauLines[i].EnergyErg =
596 
597  /*line wavelength must be gt 0 */
598  ASSERT( TauLines[i].WLAng > 0 );
599  }
600 
601  /*Beginning of the atmolEmis*/
602  for( i=0; i <linesAdded2; i++)
603  {
604  atmolEmis[i].gf = (realnum)GetGF(atmolEmis[i].Aul,
605  atmolEmis[i].tran->EnergyWN,atmolEmis[i].tran->Hi->g);
606  atmolEmis[i].dampXvel = (realnum)(atmolEmis[i].Aul/
607  atmolEmis[i].tran->EnergyWN/PI4);
608  atmolEmis[i].damp = -1000.0;
609  /*Put null terminated line label into chLab*/
610  strncpy(chLab,atmolEmis[i].tran->Hi->chLabel,4);
611  chLab[4]='\0';
612 
613  atmolEmis[i].tran->ipCont =
614  ipLineEnergy(atmolEmis[i].tran->EnergyWN * WAVNRYD, chLab ,0);
615  atmolEmis[i].ipFine = ipFineCont(atmolEmis[i].tran->EnergyWN * WAVNRYD );
616  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
617  atmolEmis[i].opacity =
618  (realnum)(abscf(atmolEmis[i].gf,atmolEmis[i].tran->EnergyWN,
619  atmolEmis[i].tran->Lo->g));
620  /* excitation energy of in degrees kelvin */
621  atmolEmis[i].tran->EnergyK =
623  /* energy of photon in ergs */
624  atmolEmis[i].tran->EnergyErg =
626  }
627  /*end of the atmolEmis*/
628 
629  /* set the ipCont struc element for the H2 molecule */
630  H2_ContPoint();
631 
632  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
633  {
634  /* do remaining part of the he iso sequence */
635  for( nelem=2; nelem < LIMELM; nelem++ )
636  {
637  if( dense.lgElmtOn[nelem])
638  {
639  /* generate label for this ion */
640  sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
641  /* Lya itself is the only transition below n=3 - we explicitly do not
642  * want to generate pointers for 2s-1s or 2p-2s */
643  /* array index for continuum edges for levels in He-like ions */
644  iso.ipIsoLevNIonCon[ipISO][nelem][0] =
645  ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
646 
647  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
648  {
649  /* array index for continuum edges for levels in He-like ions */
650  iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
651 
652  /* define all he-like line pointers */
653  for( ipLo=0; ipLo < ipHi; ipLo++ )
654  {
655 
656  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
657  continue;
658 
659  /* some lines have negative or zero energy */
660  /* >>chng 03 apr 22, this check was if less than or equal to zero,
661  * changed to lowest energy point so that ultra low energy transitions are
662  * not considered. */
663  if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
664  continue;
665 
666  Transitions[ipISO][nelem][ipHi][ipLo].ipCont =
667  ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
668  iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
669  Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine =
670  ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
671  }
672  }
673  iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
674  }
675  }
676  }
677  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
678  {
679  /* this will be over HI, HeII, then HeI only */
680  for( nelem=ipISO; nelem < LIMELM; nelem++ )
681  {
682  if( dense.lgElmtOn[nelem])
683  {
684  ipLo = 0;
685  /* these are the extra Lyman lines */
686  for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )
687  {
688  /* some energies are negative for inverted levels */
689  /*lint -e772 chLab not initialized */
690  ExtraLymanLines[ipISO][nelem][ipHi].ipCont =
691  ipLineEnergy(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
692  iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
693 
694  ExtraLymanLines[ipISO][nelem][ipHi].Emis->ipFine =
695  ipFineCont(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
696  /*lint +e772 chLab not initialized */
697  }
698 
699  if( iso.lgDielRecom[ipISO] )
700  {
701  for( ipHi=0; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
702  {
703 
704  SatelliteLines[ipISO][nelem][ipHi].ipCont = ipLineEnergy(
705  SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
706  iso.ipIsoLevNIonCon[ipISO][nelem][0]);
707 
708  SatelliteLines[ipISO][nelem][ipHi].Emis->ipFine =
709  ipFineCont(SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
710  }
711  }
712  }
713  }
714  }
715 
716  /* some lines do not exist, the flag indicating this is ipCont == -1 */
717  /* for H-like sequence, these are 2p-2s (energies degenerate) and 2s-1s, two photon */
718  ipISO = ipHYDROGEN;
719  /* do remaining part of the he iso sequence */
720  for( nelem=ipISO; nelem < LIMELM; nelem++ )
721  {
722  if( dense.lgElmtOn[nelem])
723  {
724  /* for H-like sequence want 2p-2s (energies degenerate) and 2s-1s, two photon */
725  Transitions[ipISO][nelem][ipH2p][ipH2s].ipCont = -1;
726  //Transitions[ipISO][nelem][ipH2p][ipH2s].Emis->ipFine = -1;
727  Transitions[ipISO][nelem][ipH2s][ipH1s].ipCont = -1;
728  //Transitions[ipISO][nelem][ipH2s][ipH1s].Emis->ipFine = -1;
729  }
730  }
731 
732  fixit(); /* is this redundant? */
733  /* for He-like sequence the majority of the transitions are bogus - A set to special value in this case */
734  ipISO = ipHELIUM;
735  /* do remaining part of the he iso sequence */
736  for( nelem=ipISO; nelem < LIMELM; nelem++ )
737  {
738  if( dense.lgElmtOn[nelem])
739  {
740  /* this is the two photon transition in the singlets */
741  Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].ipCont = -1;
742  Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].Emis->ipFine = -1;
743 
744  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
745  {
746  for( ipLo=0; ipLo < ipHi; ipLo++ )
747  {
748  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
749  continue;
750 
751  if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul - iso.SmallA) < SMALLFLOAT )
752  {
753  /* iso.SmallA is value assigned to bogus transitions */
754  Transitions[ipISO][nelem][ipHi][ipLo].ipCont = -1;
755  Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine = -1;
756  }
757  }
758  }
759  }
760  }
761 
762  /* co carbon monoxide line array */
763  for( i=0; i < nCORotate; i++ )
764  {
765  /* excitation energy of transition in degrees kelvin */
766  C12O16Rotate[i].EnergyK =
768  C13O16Rotate[i].EnergyK =
770 
771  /* energy of photon in ergs */
772  C12O16Rotate[i].EnergyErg =
774  C13O16Rotate[i].EnergyErg =
776 
777  /* put null terminated line label into chLab using line array*/
778  chIonLbl(chLab, &C12O16Rotate[i]);
779  chIonLbl(chLab, &C13O16Rotate[i]);
780 
781  C12O16Rotate[i].ipCont =
782  ipLineEnergy(C12O16Rotate[i].EnergyWN * WAVNRYD, "12CO" ,0);
783  C12O16Rotate[i].Emis->ipFine =
784  ipFineCont(C12O16Rotate[i].EnergyWN * WAVNRYD );
785  C13O16Rotate[i].ipCont =
786  ipLineEnergy(C13O16Rotate[i].EnergyWN * WAVNRYD, "13CO" ,0);
787  C13O16Rotate[i].Emis->ipFine =
788  ipFineCont(C13O16Rotate[i].EnergyWN * WAVNRYD );
789 
790  if( C12O16Rotate[i].Emis->gf > 0. )
791  {
792  /* the gf value was entered
793  * derive the A, call to function is gf, wl (A), g_up */
794  C12O16Rotate[i].Emis->Aul =
795  (realnum)(eina(C12O16Rotate[i].Emis->gf,
797  }
798  else if( C12O16Rotate[i].Emis->Aul > 0. )
799  {
800  /* the Einstein A value was entered
801  * derive the gf, call to function is A, wl (A), g_up */
802  C12O16Rotate[i].Emis->gf =
803  (realnum)(GetGF(C12O16Rotate[i].Emis->Aul,
805  }
806  else
807  {
808  fprintf( ioQQQ, " 12CO line does not have valid gf or A\n" );
809  fprintf( ioQQQ, " This is ContCreatePointers\n" );
810  cdEXIT(EXIT_FAILURE);
811  }
812  if( C13O16Rotate[i].Emis->gf > 0. )
813  {
814  /* the gf value was entered
815  * derive the A, call to function is gf, wl (A), g_up */
816  C13O16Rotate[i].Emis->Aul =
817  (realnum)(eina(C13O16Rotate[i].Emis->gf,
819  }
820  else if( C13O16Rotate[i].Emis->Aul > 0. )
821  {
822  /* the Einstein A value was entered
823  * derive the gf, call to function is A, wl (A), g_up */
824  C13O16Rotate[i].Emis->gf =
825  (realnum)(GetGF(C13O16Rotate[i].Emis->Aul,
827  }
828  else
829  {
830  fprintf( ioQQQ, " 13CO line does not have valid gf or A\n" );
831  fprintf( ioQQQ, " This is ContCreatePointers\n" );
832  cdEXIT(EXIT_FAILURE);
833  }
834 
835  /*line wavelength must be gt 0 */
836  ASSERT( C12O16Rotate[i].WLAng > 0 );
837  ASSERT( C13O16Rotate[i].WLAng > 0 );
838 
839  C12O16Rotate[i].Emis->dampXvel =
840  (realnum)(C12O16Rotate[i].Emis->Aul/
842 
843  C13O16Rotate[i].Emis->dampXvel =
844  (realnum)(C13O16Rotate[i].Emis->Aul/
846 
847  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
848  C12O16Rotate[i].Emis->opacity =
849  (realnum)(abscf(C12O16Rotate[i].Emis->gf,
851  C13O16Rotate[i].Emis->opacity =
852  (realnum)(abscf(C13O16Rotate[i].Emis->gf,
854  }
855 
856  /* inner shell transitions */
857  for( i=0; i<nUTA; ++i )
858  {
859  if( UTALines[i].Emis->Aul > 0. )
860  {
861  UTALines[i].Emis->dampXvel =
862  (realnum)(UTALines[i].Emis->Aul/ UTALines[i].EnergyWN/PI4);
863 
864  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
865  UTALines[i].Emis->opacity =
866  (realnum)(abscf( UTALines[i].Emis->gf, UTALines[i].EnergyWN, UTALines[i].Lo->g));
867 
868  /* excitation energy of transition in degrees kelvin */
869  UTALines[i].EnergyK =
870  (realnum)(T1CM*UTALines[i].EnergyWN);
871 
872  /* energy of photon in ergs */
873  UTALines[i].EnergyErg =
874  (realnum)(ERG1CM*UTALines[i].EnergyWN);
875 
876  /* put null terminated line label into chLab using line array*/
877  chIonLbl(chLab, &UTALines[i]);
878 
879  /* get pointer to energy in continuum mesh */
880  UTALines[i].ipCont = ipLineEnergy(UTALines[i].EnergyWN * WAVNRYD , chLab,0 );
881  UTALines[i].Emis->ipFine = ipFineCont(UTALines[i].EnergyWN * WAVNRYD );
882 
883  /* find heating per absorption,
884  * first find threshold for this shell in ergs */
885  /* ionization threshold in erg */
886  double thresh = Heavy.Valence_IP_Ryd[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] *EN1RYD;
887  UTALines[i].Coll.heat = (UTALines[i].EnergyErg-thresh);
888  ASSERT( UTALines[i].Coll.heat> 0. );
889  }
890  }
891 
892  /* level 2 heavy element lines */
893  for( i=0; i < nWindLine; i++ )
894  {
895 
896  /* derive the A, call to function is gf, wl (A), g_up */
897  TauLine2[i].Emis->Aul =
898  (realnum)(eina(TauLine2[i].Emis->gf,
899  TauLine2[i].EnergyWN,TauLine2[i].Hi->g));
900 
901  /* coefficient needed for damping constant - units cm s-1 */
902  TauLine2[i].Emis->dampXvel =
903  (realnum)(TauLine2[i].Emis->Aul/
904  TauLine2[i].EnergyWN/PI4);
905 
906  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
907  TauLine2[i].Emis->opacity =
908  (realnum)(abscf(TauLine2[i].Emis->gf,
909  TauLine2[i].EnergyWN,TauLine2[i].Lo->g));
910 
911  /* excitation energy of transition in degrees kelvin */
912  TauLine2[i].EnergyK =
913  (realnum)(T1CM*TauLine2[i].EnergyWN);
914 
915  /* energy of photon in ergs */
916  TauLine2[i].EnergyErg =
917  (realnum)(ERG1CM*TauLine2[i].EnergyWN);
918 
919  /* put null terminated line label into chLab using line array*/
920  chIonLbl(chLab, &TauLine2[i]);
921 
922  /* get pointer to energy in continuum mesh */
923  TauLine2[i].ipCont = ipLineEnergy(TauLine2[i].EnergyWN * WAVNRYD , chLab,0 );
924  TauLine2[i].Emis->ipFine = ipFineCont(TauLine2[i].EnergyWN * WAVNRYD );
925  /*if( TauLine2[i].ipCont==751 )
926  fprintf(ioQQQ,"( atom_level2 %s\n", chLab);*/
927  }
928 
929  /* hyperfine structure lines */
930  for( i=0; i < nHFLines; i++ )
931  {
932  /* derive the A, call to function is gf, wl (A), g_up
933  HFLines[i].Emis->Aul =
934  (realnum)(eina(HFLines[i].Emis->gf,
935  HFLines[i].EnergyWN,HFLines[i].Hi->g));*/
936 
937  /* coefficient needed for damping constant */
938  HFLines[i].Emis->dampXvel =
939  (realnum)(HFLines[i].Emis->Aul/
940  HFLines[i].EnergyWN/PI4);
941  HFLines[i].Emis->damp = 1e-20f;
942 
943  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
944  HFLines[i].Emis->opacity =
945  (realnum)(abscf(HFLines[i].Emis->gf,
946  HFLines[i].EnergyWN,
947  HFLines[i].Lo->g));
948  /* gf from this and 21 cm do not agree, A for HFS is 10x larger than level1 dat */
949  /*fprintf(ioQQQ,"HFLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
950  i,HFLines[i].Emis->opacity , HFLines[i].Emis->gf , HFLines[i].Emis->Aul , HFLines[i].EnergyWN,HFLines[i].Lo->g);*/
951 
952  /* excitation energy of transition in degrees kelvin */
953  HFLines[i].EnergyK =
954  (realnum)(T1CM*HFLines[i].EnergyWN);
955 
956  /* energy of photon in ergs */
957  HFLines[i].EnergyErg =
958  (realnum)(ERG1CM*HFLines[i].EnergyWN);
959 
960  /* put null terminated line label into chLab using line array*/
961  chIonLbl(chLab, &HFLines[i]);
962 
963  /* get pointer to energy in continuum mesh */
964  HFLines[i].ipCont = ipLineEnergy(HFLines[i].EnergyWN * WAVNRYD , chLab,0 );
965  HFLines[i].Emis->ipFine = ipFineCont(HFLines[i].EnergyWN * WAVNRYD );
966  }
967 
968  /* Verner's FeII lines - do first so following labels will over write this
969  * only call if large FeII atom is turned on */
970  FeIIPoint();
971 
972  /* the group of inner shell fluorescent lines */
973  for( i=0; i < t_yield::Inst().nlines(); ++i )
974  {
975  strcpy( chLab , elementnames.chElementSym[t_yield::Inst().nelem(i)] );
976  strcat( chLab , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] );
977 
978  t_yield::Inst().set_ipoint( i, ipLineEnergy( t_yield::Inst().energy(i) , chLab , 0 ) );
979  }
980 
981  /* ================================================================================== */
982  /* two photon two-photon 2-nu 2 nu 2 photon 2-photon */
983 
984  /* for induced two photon emission we need mirror image set
985  * of continuum indices for continuum between Lya and half Lya */
986  /* first MALLOC LIMELM dimension of space */
987  /* >>chng 02 jun 28, malloc this NISO instead of 2. */
988  iso.ipSym2nu.reserve( NISO );
989 
990  /* now loop over the two iso-sequences with two photon continua */
991  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
992  {
993  iso.ipSym2nu.reserve( ipISO, LIMELM );
994 
995  /* set up two photon emission */
996  for( nelem=ipISO; nelem<LIMELM; ++nelem )
997  {
998  if( dense.lgElmtOn[nelem] )
999  {
1000  double E2nu = Transitions[ipISO][nelem][TwoS][0].EnergyWN * WAVNRYD;
1001 
1002  /* pointer to the Lya energy */
1003  iso.ipTwoPhoE[ipISO][nelem] = ipoint(E2nu);
1004  /* this half-energy is only used to get induced rates in two photon */
1005  iso.ipHalfTwoPhoE[ipISO][nelem] = ipoint(E2nu / 2.);
1006  while( rfield.anu[iso.ipTwoPhoE[ipISO][nelem]] > E2nu )
1007  {
1008  --iso.ipTwoPhoE[ipISO][nelem];
1009  }
1010  while( rfield.anu[iso.ipHalfTwoPhoE[ipISO][nelem]] > E2nu / 2. )
1011  {
1012  --iso.ipHalfTwoPhoE[ipISO][nelem];
1013  }
1014 
1015  iso.ipSym2nu.reserve( ipISO, nelem, iso.ipTwoPhoE[ipISO][nelem] );
1016  }
1017  }
1018  }
1019 
1020  iso.ipSym2nu.alloc();
1022 
1023  /* now loop over the two iso-sequences with two photon continua */
1024  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1025  {
1026  /* set up two photon emission */
1027  for( nelem=ipISO; nelem<LIMELM; ++nelem )
1028  {
1029  if( dense.lgElmtOn[nelem] )
1030  {
1031  double E2nu = Transitions[ipISO][nelem][TwoS][0].EnergyWN * WAVNRYD;
1032  double Aul = Transitions[ipISO][nelem][TwoS][0].Emis->Aul;
1033  double SumShapeFunction = 0., Renorm= 0.;
1034 
1035  /* >>chng 02 aug 14, change upper limit to full Lya energy */
1036  for( i=0; i < iso.ipTwoPhoE[ipISO][nelem]; i++ )
1037  {
1038  /* energy is symmetric energy, the other side of half E2nu */
1039  energy = E2nu - rfield.anu[i];
1040  /* this is needed since mirror image of cell next to two-nu energy
1041  * may be slightly negative */
1042  energy = MAX2( energy, rfield.anu[0] + rfield.widflx[0]/2. );
1043  /* find index for this symmetric energy */
1044  iso.ipSym2nu[ipISO][nelem][i] = ipoint(energy);
1045  while( rfield.anu[iso.ipSym2nu[ipISO][nelem][i]] > E2nu ||
1046  iso.ipSym2nu[ipISO][nelem][i] >= iso.ipTwoPhoE[ipISO][nelem])
1047  {
1048  --iso.ipSym2nu[ipISO][nelem][i];
1049  }
1050  ASSERT( iso.ipSym2nu[ipISO][nelem][i] >= 0 );
1051  }
1052 
1053  /* ipTwoPhoE is the cell holding the 2nu energy itself, and we do not want
1054  * to include that in the following sum */
1055  for( i=0; i<iso.ipTwoPhoE[ipISO][nelem]; i++ )
1056  {
1057  double ShapeFunction;
1058 
1059  ASSERT( rfield.anu[i]<=E2nu );
1060 
1061  ShapeFunction = atmdat_2phot_shapefunction( rfield.AnuOrg[i]/E2nu, ipISO, nelem )*rfield.widflx[i]/E2nu;
1062  SumShapeFunction += ShapeFunction;
1063 
1064  /* >>refer HI 2nu Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
1065  /* As2nu will add up to the A, so its units are s-1 */
1066  iso.As2nu[ipISO][nelem][i] = (realnum)( Aul * ShapeFunction );
1067  }
1068 
1069  /* The spline function in twophoton.c causes a bit of an error in the integral of the
1070  * shape function. So we renormalize the integral to 1. */
1071  Renorm = 1./SumShapeFunction;
1072 
1073  for( i=0; i<iso.ipTwoPhoE[ipISO][nelem]; i++ )
1074  {
1075  iso.As2nu[ipISO][nelem][i] *= (realnum)Renorm;
1076  }
1077 
1078  /* The result should be VERY close to 1. */
1079  ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );
1080  }
1081  }
1082  }
1083 
1084  {
1085  /* this is an option to print out one of the two photon continua */
1086  enum {DEBUG_LOC=false};
1087  if( DEBUG_LOC )
1088  {
1089 # define NCRS 21
1090  double limit,ener[NCRS]={
1091  0., 0.03738, 0.07506, 0.1124, 0.1498, 0.1875,
1092  0.225, 0.263, 0.300, 0.3373, 0.375, 0.4127,
1093  0.4500, 0.487, 0.525, 0.5625, 0.6002, 0.6376,
1094  0.6749, 0.7126, 0.75};
1095 
1096  nelem = ipHYDROGEN;
1097  ipISO = ipHYDROGEN;
1098 
1099  limit = iso.ipTwoPhoE[ipISO][nelem];
1100 
1101  for( i=0; i < NCRS; i++ )
1102  {
1103  fprintf(ioQQQ,"%.3e\t%.3e\n", ener[i] ,
1104  atmdat_2phot_shapefunction( ener[i]/0.75, ipISO, nelem ) );
1105  }
1106 
1107  xnew = 0.;
1109  for( i=0; i < limit; i++ )
1110  {
1111  double fach = iso.As2nu[ipISO][nelem][i]/2.*rfield.anu2[i]/rfield.widflx[i]*EN1RYD;
1112  fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\n",
1113  rfield.anu[i] ,
1114  iso.As2nu[ipISO][nelem][i] / rfield.widflx[i] ,
1115  fach );
1116  xnew += iso.As2nu[ipISO][nelem][i];
1117  }
1118  fprintf(ioQQQ," sum is %.3e\n", xnew );
1119  cdEXIT(EXIT_FAILURE);
1120  }
1121  }
1122 
1123  {
1124  /* this is an option to print out one of the two photon continua */
1125  enum {DEBUG_LOC=false};
1126  if( DEBUG_LOC )
1127  {
1128  for( i=0; i<11; ++i )
1129  {
1130  char chLsav[10];
1131  TauDummy.WLAng = (realnum)(PI * pow(10.,(double)i));
1132  strcpy( chLsav, chLineLbl(&TauDummy) );
1133  fprintf(ioQQQ,"%.2f\t%s\n", TauDummy.WLAng , chLsav );
1134  }
1135  cdEXIT(EXIT_FAILURE);
1136  }
1137  }
1138 
1139  /* option to print out whole thing with "trace lines" command */
1140  if( trace.lgTrLine )
1141  {
1142  fprintf( ioQQQ, " WL(Ang) E(RYD) IP gl gu gf A damp abs K\n" );
1143  for( i=1; i <= nLevel1; i++ )
1144  {
1145  strcpy( chLab, chLineLbl(&TauLines[i]) );
1146  iWL_Ang = (long)TauLines[i].WLAng;
1147  if( iWL_Ang > 1000000 )
1148  {
1149  iWL_Ang /= 10000;
1150  }
1151  else if( iWL_Ang > 10000 )
1152  {
1153  iWL_Ang /= 1000;
1154  }
1155 
1156  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1157  chLab, iWL_Ang, RYDLAM/TauLines[i].WLAng,
1158  TauLines[i].ipCont, (long)(TauLines[i].Lo->g),
1159  (long)(TauLines[i].Hi->g), TauLines[i].Emis->gf,
1160  TauLines[i].Emis->Aul, TauLines[i].Emis->dampXvel,
1161  TauLines[i].Emis->opacity );
1162  }
1163 
1164  /*Atomic Or Molecular lines*/
1165  for( i=0; i <linesAdded2; i++)
1166  {
1167  strcpy( chLab, chLineLbl(atmolEmis[i].tran) );
1168 
1169  iWL_Ang = (long)atmolEmis[i].tran->WLAng;
1170 
1171  if( iWL_Ang > 1000000 )
1172  {
1173  iWL_Ang /= 10000;
1174  }
1175  else if( iWL_Ang > 10000 )
1176  {
1177  iWL_Ang /= 1000;
1178  }
1179  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1180  chLab, iWL_Ang, RYDLAM/atmolEmis[i].tran->WLAng,
1181  atmolEmis[i].tran->ipCont, (long)(atmolEmis[i].tran->Lo->g),
1182  (long)(atmolEmis[i].tran->Hi->g),atmolEmis[i].gf,
1184  atmolEmis[i].opacity);
1185  }
1186 
1187  /* C12O16 lines */
1188  for( i=0; i < nCORotate; i++ )
1189  {
1190  strcpy( chLab, chLineLbl(&C12O16Rotate[i]) );
1191 
1192  iWL_Ang = (long)C12O16Rotate[i].WLAng;
1193 
1194  if( iWL_Ang > 1000000 )
1195  {
1196  iWL_Ang /= 10000;
1197  }
1198  else if( iWL_Ang > 10000 )
1199  {
1200  iWL_Ang /= 1000;
1201  }
1202  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1203  chLab, iWL_Ang, RYDLAM/C12O16Rotate[i].WLAng,
1204  C12O16Rotate[i].ipCont, (long)(C12O16Rotate[i].Lo->g),
1205  (long)(C12O16Rotate[i].Hi->g), C12O16Rotate[i].Emis->gf,
1207  C12O16Rotate[i].Emis->opacity );
1208  }
1209 
1210  /* C13O16 lines */
1211  for( i=0; i < nCORotate; i++ )
1212  {
1213  strcpy( chLab, chLineLbl(&C13O16Rotate[i]) );
1214 
1215  iWL_Ang = (long)C13O16Rotate[i].WLAng;
1216 
1217  if( iWL_Ang > 1000000 )
1218  {
1219  iWL_Ang /= 10000;
1220  }
1221  else if( iWL_Ang > 10000 )
1222  {
1223  iWL_Ang /= 1000;
1224  }
1225  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1226  chLab, iWL_Ang, RYDLAM/C13O16Rotate[i].WLAng,
1227  C13O16Rotate[i].ipCont, (long)(C13O16Rotate[i].Lo->g),
1228  (long)(C13O16Rotate[i].Hi->g), C13O16Rotate[i].Emis->gf,
1230  C13O16Rotate[i].Emis->opacity );
1231  }
1232 
1233  for( i=0; i < nWindLine; i++ )
1234  {
1235  strcpy( chLab, chLineLbl(&TauLine2[i]) );
1236 
1237  iWL_Ang = (long)TauLine2[i].WLAng;
1238 
1239  if( iWL_Ang > 1000000 )
1240  {
1241  iWL_Ang /= 10000;
1242  }
1243  else if( iWL_Ang > 10000 )
1244  {
1245  iWL_Ang /= 1000;
1246  }
1247  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1248  chLab, iWL_Ang, RYDLAM/TauLine2[i].WLAng,
1249  TauLine2[i].ipCont, (long)(TauLine2[i].Lo->g),
1250  (long)(TauLine2[i].Hi->g), TauLine2[i].Emis->gf,
1251  TauLine2[i].Emis->Aul, TauLine2[i].Emis->dampXvel,
1252  TauLine2[i].Emis->opacity );
1253  }
1254  for( i=0; i < nHFLines; i++ )
1255  {
1256  strcpy( chLab, chLineLbl(&HFLines[i]) );
1257 
1258  iWL_Ang = (long)HFLines[i].WLAng;
1259 
1260  if( iWL_Ang > 1000000 )
1261  {
1262  iWL_Ang /= 10000;
1263  }
1264  else if( iWL_Ang > 10000 )
1265  {
1266  iWL_Ang /= 1000;
1267  }
1268  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
1269  chLab, iWL_Ang, RYDLAM/HFLines[i].WLAng,
1270  HFLines[i].ipCont, (long)(HFLines[i].Lo->g),
1271  (long)(HFLines[i].Hi->g), HFLines[i].Emis->gf,
1272  HFLines[i].Emis->Aul, HFLines[i].Emis->dampXvel,
1273  HFLines[i].Emis->opacity );
1274  }
1275  }
1276 
1277  /* this is an option to kill fine structure line optical depths */
1278  if( !rt.lgFstOn )
1279  {
1280  for( i=1; i <= nLevel1; i++ )
1281  {
1282  if( TauLines[i].EnergyWN < 10000. )
1283  {
1284  TauLines[i].Emis->opacity = 0.;
1285  }
1286  }
1287 
1288  /*Atomic or Molecular Lines-Humeshkar Nemala*/
1289  for( i=0; i <linesAdded2; i++)
1290  {
1291  if(atmolEmis[i].tran->EnergyWN < 10000. )
1292  {
1293  atmolEmis[i].opacity = 0.;
1294  }
1295  }
1296  }
1297 
1298  /* read in continuum bands data set */
1299  ContBandsCreate( "" );
1300 
1301  /* we're done adding lines and states to the stacks.
1302  * This flag is used to make sure we never add them again in this coreload. */
1303  lgLinesAdded = true;
1304  lgStatesAdded = true;
1305 
1306  return;
1307 }
1308 
1309 /*fiddle adjust energy bounds of cell with index ipnt so that lower energy
1310  * matches ionization edges exactly, called by createpoint */
1311 /* ipnt is index on c scale */
1312 STATIC void fiddle(long int ipnt,
1313  double exact)
1314 {
1315  realnum Ehi,
1316  Elo,
1317  OldEner;
1318 
1319  DEBUG_ENTRY( "fiddle()" );
1320 
1321  /* make low edge of cell exact for photo integrals */
1322  ASSERT( ipnt >= 0 );
1323  ASSERT( ipnt < rfield.nupper-1 );
1324 
1325  /* upper edge of higher cell*/
1326  Ehi = rfield.anu[ipnt] + rfield.widflx[ipnt]*0.5f;
1327  /* lower edge of lower cell */
1328  Elo = rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5f;
1329 
1330  /* >>chng 02 nov 11, do nothing if already very close,
1331  * because VERY high res models had negative widths for some very close edges */
1332  if( fabs( Elo/exact - 1. ) < 0.001 )
1333  return;
1334 
1335  ASSERT( Elo <= exact );
1336 
1337  OldEner = rfield.anu[ipnt];
1338 
1339  /* centroid of desired lower energy and upper edge */
1340  rfield.anu[ipnt] = (realnum)((Ehi + exact)/2.);
1341  /* same for lower */
1342  rfield.anu[ipnt-1] = (realnum)((exact + Elo)/2.);
1343 
1344  rfield.widflx[ipnt] = (realnum)(Ehi - exact);
1345  rfield.widflx[ipnt-1] = (realnum)(exact - Elo);
1346 
1347  /* bring upper cell down a bit too, since we dont want large change in widflx */
1348  rfield.anu[ipnt+1] -= (OldEner - rfield.anu[ipnt])/2.f;
1349 
1350  /* sanity check */
1351  ASSERT( rfield.widflx[ipnt-1] > 0. );
1352  ASSERT( rfield.widflx[ipnt] > 0. );
1353  return;
1354 }
1355 
1356 /*ipShells assign continuum energy pointers to shells for all atoms,
1357  * called by ContCreatePointers */
1359  /* nelem is the atomic number on the C scale, Li is 2 */
1360  long int nelem)
1361 {
1362  char chLab[5];
1363  long int
1364  imax,
1365  ion,
1366  nelec,
1367  ns,
1368  nshell;
1369  /* following value cannot be used - will be set to proper threshold */
1370  double thresh=-DBL_MAX;
1371 
1372  DEBUG_ENTRY( "ipShells()" );
1373 
1374  ASSERT( nelem >= NISO);
1375  ASSERT( nelem < LIMELM );
1376 
1377  /* fills in pointers to valence shell ionization threshold
1378  * PH1(a,b,c,d)
1379  * a=1 => thresh, others fitting parameters
1380  * b atomic number
1381  * c number of electrons
1382  * d shell number 7-1 */
1383 
1384  /* threshold in Ryd
1385  * ion=0 for atom, up to nelem-1 for helium like, hydrogenic is elsewhere */
1386  for( ion=0; ion < nelem; ion++ )
1387  {
1388  /* generate label for ionization stage */
1389  /* this is short form of element name */
1390  strcpy( chLab, elementnames.chElementSym[nelem] );
1391 
1392  /* this is a number between 1 and 31 */
1393  strcat( chLab, elementnames.chIonStage[ion] );
1394 
1395  /* this is the iso sequence - must not redo sequence if done as iso */
1396  long int ipISO = nelem-ion;
1397 
1398  /* number of bound electrons */
1399  nelec = ipISO+1;
1400 
1401  /* nsShells(nelem,ion) is the number of shells for ion with nelec electrons,
1402  * physical not c scale */
1403  imax = Heavy.nsShells[nelem][ion];
1404 
1405  /* loop on all inner shells, valence shell */
1406  for( nshell=0; nshell < imax; nshell++ )
1407  {
1408  /* ionization potential of this shell in rydbergs */
1409  thresh = (double)(t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
1410  if( thresh <= 0.1 )
1411  {
1412  /* negative ip shell does not exist, set upper limit
1413  * to less than lower limit so this never looped upon
1414  * these are used as flags by LimitSh to check whether
1415  * this is a real shell - if 1 or 2 is changed - change LimitSh!! */
1416  opac.ipElement[nelem][ion][nshell][0] = 2;
1417  opac.ipElement[nelem][ion][nshell][1] = 1;
1418  }
1419  else
1420  {
1421  /* this is lower bound to energy range for this shell */
1422  /* >>chng 02 may 27, change to version of ip with label, so that
1423  * inner shell edges will appear */
1424  /*opac.ipElement[nelem][ion][nshell][0] = ipoint(thresh);*/
1425  opac.ipElement[nelem][ion][nshell][0] =
1426  ipContEnergy( thresh , chLab );
1427 
1428  /* this is upper bound to energy range for this shell
1429  * LimitSh is an integer function, returns pointer
1430  * to threshold of next major shell. For k-shell it
1431  * returns the values KshellLimit, default=7.35e4
1432  * >>chng 96 sep 26, had been below, result zero cross sec at
1433  * many energies where opacity project did not produce state specific
1434  * cross section */
1435  opac.ipElement[nelem][ion][nshell][1] =
1436  LimitSh(ion+1, nshell+1,nelem+1);
1437  ASSERT( opac.ipElement[nelem][ion][nshell][1] > 0);
1438  }
1439  }
1440 
1441  ASSERT( imax > 0 && imax <= 7 );
1442 
1443  /* this will be index pointing to valence edge */
1444  /* [0] is pointer to threshold in energy array */
1445  opac.ipElement[nelem][ion][imax-1][0] =
1446  ipContEnergy(thresh, chLab);
1447 
1448  /* pointer to valence electron ionization potential */
1449  Heavy.ipHeavy[nelem][ion] = opac.ipElement[nelem][ion][imax-1][0];
1450 
1451  /* ionization potential of valence shell in Ryd
1452  * thresh was evaluated above, now has last value, the valence shell */
1453  Heavy.Valence_IP_Ryd[nelem][ion] = thresh;
1454 
1455  Heavy.xLyaHeavy[nelem][ion] = 0.;
1456  if( ipISO >= NISO )
1457  {
1458  /* this is set of 3/4 of valence shell IP, this is important
1459  * source of ots deep in cloud */
1460  Heavy.ipLyHeavy[nelem][ion] =
1461  ipLineEnergy(thresh*0.75,chLab , 0);
1462 
1463 
1464  Heavy.ipBalHeavy[nelem][ion] =
1465  ipLineEnergy(thresh*0.25,chLab , 0);
1466  }
1467  else
1468  {
1469  /* do not treat this simple way since done exactly with iso
1470  * sequences */
1471  Heavy.ipLyHeavy[nelem][ion] = -1;
1472  Heavy.ipBalHeavy[nelem][ion] = -1;
1473  }
1474  }
1475 
1476  /* above loop did up to hydrogenic, now do hydrogenic -
1477  * hydrogenic is special since arrays already set up */
1478  Heavy.nsShells[nelem][nelem] = 1;
1479 
1480  /* this is lower limit to range */
1481  /* hydrogenic photoionization set to special hydro array
1482  * this is pointer to threshold energy */
1483  /* this statement is in ContCreatePointers but has not been done when this routine called */
1484  /*iso.ipIsoLevNIonCon[ipH_LIKE][ipZ][ipLo] = ipContEnergy(iso.xIsoLevNIonRyd[ipH_LIKE][ipZ][ipLo],chLab);*/
1485  /*opac.ipElement[nelem][nelem][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];*/
1486  opac.ipElement[nelem][nelem][0][0] = ipoint( iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH1s] );
1487  ASSERT( opac.ipElement[nelem][nelem][0][0] > 0 );
1488 
1489  /* this is the high-energy limit */
1490  opac.ipElement[nelem][nelem][0][1] = continuum.KshellLimit;
1491 
1492  Heavy.ipHeavy[nelem][nelem] = opac.ipElement[nelem][nelem][0][0];
1493 
1494  /* this is for backwards computability with Cambridge code */
1495  if( trace.lgTrace && trace.lgPointBug )
1496  {
1497  for( ion=0; ion < (nelem+1); ion++ )
1498  {
1499  fprintf( ioQQQ, "Ion:%3ld%3ld %2.2s%2.2s total shells:%3ld\n",
1500  nelem, ion+1, elementnames.chElementSym[nelem], elementnames.chIonStage[ion]
1501  , Heavy.nsShells[nelem][ion] );
1502  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1503  {
1504  fprintf( ioQQQ, " shell%3ld %2.2s range eV%10.2e-%8.2e\n",
1505  ns+1, Heavy.chShell[ns], rfield.anu[opac.ipElement[nelem][ion][ns][0]-1]*
1506  EVRYD, rfield.anu[opac.ipElement[nelem][ion][ns][1]-1]*EVRYD );
1507  }
1508  }
1509  }
1510  return;
1511 }
1512 
1513 /*LimitSh sets upper energy limit to subshell integrations */
1514 STATIC long LimitSh(long int ion,
1515  long int nshell,
1516  long int nelem)
1517 {
1518  long int LimitSh_v;
1519 
1520  DEBUG_ENTRY( "LimitSh()" );
1521 
1522  /* this routine returns the high-energy limit to the energy range
1523  * for photoionization of a given shell
1524  * */
1525  if( nshell == 1 )
1526  {
1527  /* this limit is high-energy limit to code unless changed with set kshell */
1528  LimitSh_v = continuum.KshellLimit;
1529 
1530  }
1531  else if( nshell == 2 )
1532  {
1533  /* this is 2s shell, upper limit is 1s
1534  * >>chng 96 oct 08, up to high-energy limit
1535  * LimitSh = ipElement(nelem,ion , 1,1)-1 */
1536  LimitSh_v = continuum.KshellLimit;
1537 
1538  }
1539  else if( nshell == 3 )
1540  {
1541  /* this is 2p shell, upper limit is 1s
1542  * >>chng 96 oct 08, up to high-energy limit
1543  * LimitSh = ipElement(nelem,ion , 1,1)-1 */
1544  LimitSh_v = continuum.KshellLimit;
1545 
1546  }
1547  else if( nshell == 4 )
1548  {
1549  /* this is 3s shell, upper limit is 2p
1550  * >>chng 96 oct 08, up to K-shell edge
1551  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1552  LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1553 
1554  }
1555  else if( nshell == 5 )
1556  {
1557  /* this is 3p shell, upper limit is 2p
1558  * >>chng 96 oct 08, up to K-shell edge
1559  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1560  LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1561 
1562  }
1563  else if( nshell == 6 )
1564  {
1565  /* this is 3d shell, upper limit is 2p
1566  * >>chng 96 oct 08, up to K-shell edge
1567  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1568  LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1569 
1570  }
1571  else if( nshell == 7 )
1572  {
1573  /* this is 4s shell, upper limit is 3d */
1574  if( opac.ipElement[nelem-1][ion-1][5][0] < 3 )
1575  {
1576  /* this is check for empty shell 6, 3d
1577  * if so then set to 3p instead */
1578  LimitSh_v = opac.ipElement[nelem-1][ion-1][4][0] -
1579  1;
1580  }
1581  else
1582  {
1583  LimitSh_v = opac.ipElement[nelem-1][ion-1][5][0] -
1584  1;
1585  }
1586  /* >>chng 96 sep 26, set upper limit down to 2s */
1587  LimitSh_v = opac.ipElement[nelem-1][ion-1][2][0] - 1;
1588 
1589  }
1590  else
1591  {
1592  fprintf( ioQQQ, " LimitSh cannot handle nshell as large as%4ld\n",
1593  nshell );
1594  cdEXIT(EXIT_FAILURE);
1595  }
1596  return( LimitSh_v );
1597 }
1598 
1599 /*ContBandsCreate - read set of continuum bands to enter total emission into line*/
1601  /* chFile is optional filename, if void then use default bands,
1602  * if not void then use file specified,
1603  * return value is 0 for success, 1 for failure */
1604  const char chFile[] )
1605 {
1606  char chLine[FILENAME_PATH_LENGTH_2];
1607  const char* chFilename;
1608  FILE *ioDATA;
1609 
1610  bool lgEOL;
1611  long int i,k;
1612 
1613  /* keep track of whether we have been called - want to be
1614  * called a total of one time */
1615  static bool lgCalled=false;
1616 
1617  DEBUG_ENTRY( "ContBandsCreate()" );
1618 
1619  /* do nothing if second or later call*/
1620  if( lgCalled )
1621  {
1622  /* success */
1623  return;
1624  }
1625  lgCalled = true;
1626 
1627  /* use default filename if void string, else use file specified */
1628  chFilename = ( strlen(chFile) == 0 ) ? "bands_continuum.dat" : chFile;
1629 
1630  /* get continuum band data */
1631  if( trace.lgTrace )
1632  {
1633  fprintf( ioQQQ, " ContBandsCreate opening %s:", chFilename );
1634  }
1635 
1636  ioDATA = open_data( chFilename, "r" );
1637 
1638  /* now count how many bands are in the file */
1639  continuum.nContBand = 0;
1640 
1641  /* first line is a magic number and does not count as a band*/
1642  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1643  {
1644  fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
1645  cdEXIT(EXIT_FAILURE);
1646  }
1647  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1648  {
1649  /* we want to count the lines that do not start with #
1650  * since these contain data */
1651  if( chLine[0] != '#')
1652  ++continuum.nContBand;
1653  }
1654 
1655  /* now rewind the file so we can read it a second time*/
1656  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
1657  {
1658  fprintf( ioQQQ, " ContBandsCreate could not rewind %s.\n", chFilename );
1659  cdEXIT(EXIT_FAILURE);
1660  }
1661 
1663  continuum.chContBandLabels = (char **)MALLOC(sizeof(char *)*(unsigned)(continuum.nContBand) );
1664  continuum.ipContBandLow = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
1665  continuum.ipContBandHi = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
1666 
1667  /* now make second dim, id wavelength, and lower and upper bounds */
1668  for( i=0; i<continuum.nContBand; ++i )
1669  {
1670  /* array of labels, each 4 long plus 0 at [4] */
1671  continuum.chContBandLabels[i] = (char *)MALLOC(sizeof(char)*(unsigned)(5) );
1672  }
1673 
1674  /* first line is a versioning magic number - now confirm that it is valid */
1675  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1676  {
1677  fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
1678  cdEXIT(EXIT_FAILURE);
1679  }
1680  /* bands_continuum magic number here <- this string is in band_continuum.dat
1681  * with comment to search for this to find magic number
1682  * >>chng 06 aug 07, update bands and magic number */
1683  {
1684  long int m1 , m2 , m3,
1685  myr = 6,
1686  mmo = 8,
1687  mdy = 7;
1688  i = 1;
1689  m1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1690  m2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1691  m3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1692  if( ( m1 != myr ) ||
1693  ( m2 != mmo ) ||
1694  ( m3 != mdy ) )
1695  {
1696  fprintf( ioQQQ,
1697  " ContBandsCreate: the version of the data file %s I found (%li %li %li)is not the current version (%li %li %li).\n",
1698  chFilename ,
1699  m1 , m2 , m3 ,
1700  myr , mmo , mdy );
1701  fprintf( ioQQQ,
1702  " ContBandsCreate: you need to update this file.\n");
1703  cdEXIT(EXIT_FAILURE);
1704  }
1705  }
1706 
1707  /* now read in data again, but save it this time */
1708  k = 0;
1709  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1710  {
1711  /* we want to count the lines that do not start with #
1712  * since these contain data */
1713  if( chLine[0] != '#')
1714  {
1715  double xLow , xHi;
1716  /* copy 4 char label plus termination */
1717  strncpy( continuum.chContBandLabels[k] , chLine , 4 );
1718  continuum.chContBandLabels[k][4] = 0;
1719 
1720  /* now get central band wavelength
1721  * >>chng 06 aug 11 from 4 to 6, the first 4 char are labels and
1722  * these can contain numbers, next comes a space, then the number */
1723  i = 6;
1725  /* >>chng 06 feb 21, multiply by 1e4 to convert micron wavelength into Angstroms,
1726  * which is assumed by the code. before this correction the band centroid
1727  * wavelength was given in the output incorrectly listed as Angstroms.
1728  * results were correct just label was wrong */
1729  continuum.ContBandWavelength[k] *= 1e4f;
1730 
1731  /* these are short and long wave limits, which are high and
1732  * low energy limits - these are now wl in microns */
1733  xHi = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1734  xLow = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1735  if( lgEOL )
1736  {
1737  fprintf( ioQQQ, " There should have been 3 numbers on this band line. Sorry.\n" );
1738  fprintf( ioQQQ, " string==%s==\n" ,chLine );
1739  cdEXIT(EXIT_FAILURE);
1740  }
1741 
1742  /* make sure bands bounds are in correct order, shorter - longer wavelength*/
1743  if( xHi >= xLow )
1744  {
1745  fprintf( ioQQQ, " ContBandWavelength band %li has improper bounds.\n" ,i);
1746  cdEXIT(EXIT_FAILURE);
1747  }
1748 
1749  /* get continuum index - RYDLAM is 911.6A = 1 Ryd so 1e4 converts
1750  * micron to Angstrom */
1751  continuum.ipContBandHi[k] = ipoint( RYDLAM / (xHi*1e4) );
1752  continuum.ipContBandLow[k] = ipoint( RYDLAM / (xLow*1e4) );
1753  /*fprintf(ioQQQ,"DEBUG bands_continuum %s %.3e %li %li \n",
1754  continuum.chContBandLabels[k],
1755  continuum.ContBandWavelength[k],
1756  continuum.ipContBandHi[k] ,
1757  continuum.ipContBandLow[k]);*/
1758 
1759  if( trace.lgTrace && trace.lgConBug )
1760  {
1761  if( k==0 )
1762  fprintf( ioQQQ, " ContCreatePointer trace bands\n");
1763  fprintf( ioQQQ,
1764  " band %ld label %s low wl= %.3e low ipnt= %li "
1765  " hi wl= %.3e hi ipnt= %li \n",
1766  k,
1768  xLow,
1770  xHi,
1771  continuum.ipContBandHi[k] );
1772  }
1773  /*fprintf(ioQQQ,
1774  "DEBUG band data %s %f %f %f %f %f \n",
1775  continuum.chContBandLabels[k],
1776  continuum.ContBandWavelength[k],
1777  xHi,
1778  rfield.anu[continuum.ipContBandHi[k]-1],
1779  xLow,
1780  rfield.anu[continuum.ipContBandLow[k]-1]);*/
1781  ++k;
1782  }
1783  }
1784  /* now validate this incoming data */
1785  for( i=0; i<continuum.nContBand; ++i )
1786  {
1787  /* make sure all are positive */
1788  if( continuum.ContBandWavelength[i] <=0. )
1789  {
1790  fprintf( ioQQQ, " ContBandWavelength band %li has non-positive entry.\n",i );
1791  cdEXIT(EXIT_FAILURE);
1792  }
1793  }
1794 
1795  fclose(ioDATA);
1796  return;
1797 }

Generated for cloudy by doxygen 1.8.1.2