cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cont_setintensity.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 /*ContSetIntensity derive intensity of incident continuum */
4 /*extin do extinction of incident continuum as set by extinguish command */
5 /*sumcon sums L and Q for net incident continuum */
6 /*ptrcer show continuum pointers in real time following drive pointers command */
7 /*conorm normalize continuum to proper intensity */
8 /*qintr integrates Q for any continuum between two limits, used for normalization */
9 /*pintr integrates L for any continuum between two limits, used for normalization */
10 #include "cddefines.h"
11 #include "physconst.h"
12 #include "iso.h"
13 #include "noexec.h"
14 #include "ionbal.h"
15 #include "hextra.h"
16 #include "trace.h"
17 #include "dense.h"
18 #include "oxy.h"
19 #include "conv.h"
20 #include "prt.h"
21 #include "heavy.h"
22 #include "rfield.h"
23 #include "phycon.h"
24 #include "called.h"
25 #include "hydrogenic.h"
26 #include "timesc.h"
27 #include "secondaries.h"
28 #include "opacity.h"
29 #include "thermal.h"
30 #include "ipoint.h"
31 #include "atmdat.h"
32 #include "rt.h"
33 #include "radius.h"
34 #include "geometry.h"
35 #include "grainvar.h"
36 #include "continuum.h"
37 #include "taulines.h"
38 
39 /* these are weights used for continuum integration */
40 static const double aweigh[4]={-0.4305682,-0.1699905, 0.1699905, 0.4305682};
41 static const double fweigh[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274};
42 
43 /*conorm normalize continuum to proper intensity */
44 STATIC void conorm();
45 
46 /*pintr integrates L for any continuum between two limits, used for normalization */
47 STATIC double pintr(double penlo,
48  double penhi);
49 
50 /*qintr integrates Q for any continuum between two limits, used for normalization */
51 STATIC double qintr(double *qenlo,
52  double *qenhi);
53 
54 
55 /*sumcon sums L and Q for net incident continuum */
56 STATIC void sumcon(long int il,
57  long int ih,
58  realnum *q,
59  realnum *p,
60  realnum *panu);
61 
62 /*extin do extinction of incident continuum as set by extinguish command */
63 STATIC void extin(realnum *ex1ryd);
64 
65 /*ptrcer show continuum pointers in real time following drive pointers command */
66 STATIC void ptrcer();
67 
68 void IncidentContinuumHere()
69 {
70 
71  DEBUG_ENTRY( "IncidentContinuumHere()" );
72  double frac_beam_time;
73  /* fraction of beamed continuum that is constant */
74  double frac_beam_const;
75  /* fraction of continuum that is isotropic */
76  double frac_isotropic;
77 
78  double BigLog = 0.;
79  for( long i = 0; i<rfield.nflux; ++i )
80  {
81  double flux_org = rfield.flux[0][i];
82  double flux_now = ffun( rfield.anu[i] ,
83  &frac_beam_time ,
84  &frac_beam_const ,
85  &frac_isotropic )*rfield.widflx[i]*rfield.ExtinguishFactor[i];
86 
87  double ratio = 1.;
88  if( flux_org>SMALLFLOAT && flux_now>SMALLFLOAT )
89  {
90  ratio = fabs( log10( flux_now / flux_org ) );
91  BigLog = max( ratio , BigLog );
92  }
93  }
94 
95  if( BigLog > 0.01 )
96  fprintf(ioQQQ , "DEBUG diff continua %.2e\n", BigLog );
97  return;
98 }
99 /* called by Cloudy to set up continuum */
100 void ContSetIntensity()
101 {
102  bool lgCheckOK;
103 
104  long int i,
105  ip,
106  j,
107  n;
108 
109  realnum EdenHeav,
110  ex1ryd,
111  factor,
112  occ1,
113  p,
114  p1,
115  p2,
116  p3,
117  p4,
118  p5,
119  p6,
120  p7,
121  p8,
122  pgn,
123  phe,
124  pheii,
125  qgn;
126 
127  realnum xIoniz;
128 
129  double HCaseBRecCoeff,
130  wanu[4],
131  alf,
132  bet,
133  fntest,
134  fsum,
135  ecrit,
136  tcompr,
137  tcomp,
138  RatioIonizToRecomb,
139  r3ov2;
140 
141  double amean,
142  amean2,
143  amean3,
144  peak,
145  wfun[4];
146 
147  /* fraction of beamed continuum that is varies with time */
148  double frac_beam_time , frac_beam_time1;
149  /* fraction of beamed continuum that is constant */
150  double frac_beam_const , frac_beam_const1;
151  /* fraction of continuum that is isotropic */
152  double frac_isotropic , frac_isotropic1;
153 
154  long int nelem , ion;
155 
156  DEBUG_ENTRY( "ContSetIntensity()" );
157 
158  /* set continuum */
159  if( trace.lgTrace )
160  {
161  fprintf( ioQQQ, " ContSetIntensity called.\n" );
162  }
163 
164  /* find normalization factors for the continua - this decides whether continuum is
165  * per unit area of luminosity, and which is desired final product */
166  conorm();
167 
168  /* define factors to convert rfeld.flux array into photon occupation array OCCNUM
169  * by multiplication */
170  factor = (realnum)(EN1RYD/PI4/FR1RYD/HNU3C2);
171 
172  /*------------------------------------------------------------- */
173  lgCheckOK = true;
174 
175  for( i=0; i < rfield.nupper; i++ )
176  {
177  /* this was original anu array with no continuum averaging */
178  rfield.anu[i] = rfield.AnuOrg[i];
179  rfield.ContBoltz[i] = 0.;
180  fsum = 0.;
181  amean = 0.;
182  amean2 = 0.;
183  amean3 = 0.;
184  frac_beam_time = 0.;
185  frac_beam_const = 0.;
186  frac_isotropic = 0.;
187 
188  for( j=0; j < 4; j++ )
189  {
190  /* aweigh is symmetric about 0.5 widflx */
191  wanu[j] = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
192  /* >>chng 02 jul 16, add test on continuum limits -
193  * this was exceeded when resolution set very large */
194  wanu[j] = MAX2( wanu[j] , rfield.emm );
195  wanu[j] = MIN2( wanu[j] , rfield.egamry );
196  /* >>chng 06 feb 03, the continuum binning can change dramatically
197  * at some energies - make sure that this cell does not overextend the
198  * boundaries of its neighbors */
199  if( i > 0 && i < rfield.nupper-1 )
200  {
201  wanu[j] = MAX2( wanu[j] , rfield.anu[i-1] + 0.5*rfield.widflx[i-1] );
202  wanu[j] = MIN2( wanu[j] , rfield.anu[i+1] - 0.5*rfield.widflx[i+1] );
203  }
204 
205  wfun[j] = fweigh[j]*ffun(wanu[j] ,
206  &frac_beam_time1 ,
207  &frac_beam_const1 ,
208  &frac_isotropic1 );
209  /*if( i==76 )
210  fprintf(ioQQQ,"DEBUG ffun %li %.4e at %.4e R\n", j ,
211  ffun(wanu[j]) , wanu[j] );*/
212  fsum += wfun[j];
213  amean += wanu[j]*wfun[j];
214  amean2 += wanu[j]*wanu[j]*wfun[j];
215  amean3 += wanu[j]*wanu[j]*wanu[j]*wfun[j];
216  frac_beam_time += fweigh[j]*frac_beam_time1;
217  frac_beam_const += fweigh[j]*frac_beam_const1;
218  frac_isotropic += fweigh[j]*frac_isotropic1;
219  }
220 
221  ASSERT( fabs( 1.-frac_beam_time-frac_beam_const-frac_isotropic)<
222  10.*FLT_EPSILON);
223  /* This is a fix for ticket #1 */
224  if( fsum*rfield.widflx[i] > BIGFLOAT )
225  {
226  fprintf( ioQQQ, "\n Cannot continue. The continuum is far too intense.\n" );
227  for( j=0; j < rfield.nShape; j++ )
228  {
229  if( (wfun[j]*rfield.widflx[i] > BIGFLOAT) && ( rfield.nShape > 1 ) )
230  {
231  fprintf( ioQQQ, " Problem is with source number %li\n", j );
232  break;
233  }
234  }
235  cdEXIT(EXIT_FAILURE);
236  }
237 
238  rfield.flux[0][i] = (realnum)(fsum*rfield.widflx[i]);
239 
240  /* save separation into isotropic constant and beamed, and possibly
241  * time-variable beamed continua */
242  rfield.flux_beam_const[i] = rfield.flux[0][i] * (realnum)frac_beam_const;
243  rfield.flux_beam_time[i] = rfield.flux[0][i] * (realnum)frac_beam_time;
244  rfield.flux_isotropic[i] = rfield.flux[0][i] * (realnum)frac_isotropic;
245 
246  if( rfield.flux[0][i] > 0. )
247  {
248  /*fprintf(ioQQQ,"DEBUG %4li %e %e %.3e %.3e\n",
249  i , rfield.anu[i] , (amean2/amean) , amean2 , amean );*/
250  rfield.anu[i] = (realnum)(amean2/amean);
251  rfield.anu2[i] = (realnum)(amean3/amean);
252  /* mesh must be strictly monotonically increasing - make it so */
253  if( i > 0 && rfield.anu[i] <= rfield.anu[i-1] )
254  {
255  /* prevent roundoff from allowing i cell to lie below i-1
256  * cell when continuum mesh is very fine. */
257  /* use 2*epsilon to protect against unusual rounding modes */
258  rfield.anu[i] = rfield.anu[i-1]*(1.f+2.f*FLT_EPSILON);
259  rfield.anu2[i] = pow2(rfield.anu[i]);
260  }
261  ASSERT( i==0 || rfield.anu[i] > rfield.anu[i-1] );
262  /* define array of LOG10( nu(ryd) ) these are not guaranteed to
263  * be monotonically increasing due to loss of precision in log
264  * of float */
265  rfield.anulog[i] = (realnum)log10(rfield.anu[i]);
266  }
267 
268  else if( rfield.flux[0][i] == 0. )
269  {
270  rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
271  rfield.anulog[i] = (realnum)log10(rfield.anu[i]);
272  }
273 
274  else
275  {
276  rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
277  fprintf( ioQQQ, " negative continuum returned at%6ld%10.2e%10.2e\n",
278  i, rfield.anu[i], rfield.flux[0][i] );
279  lgCheckOK = false;
280  }
281  rfield.anu3[i] = rfield.anu2[i]*rfield.anu[i];
282 
283  rfield.ConEmitReflec[0][i] = 0.;
284  rfield.ConEmitOut[0][i] = 0.;
285  rfield.convoc[i] = factor/rfield.widflx[i]/rfield.anu2[i];
286 
287  /* following are Compton exchange factors from Tarter */
288  alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i]));
289  bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/
290  4.;
291  rfield.csigh[i] = (realnum)(alf*rfield.anu2[i]*3.858e-25);
292  rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25);
293  rfield.lgMeshSetUp = true;
294  }
295 
296  /*i = 76;
297  fprintf(ioQQQ,"\nDEBUG %4li %e \n",
298  i , rfield.anu[i] );*/
299 
300  /* >>chng 05 jul 01 add this
301  * finished with stored continua - return these vectors */
302 #if 0
303  /* commented out since we must conserve energy, and continuum was set with old widflx */
304  /* now fix widflx array so that it is correct */
305  for( i=1; i<rfield.nupper-1; ++i )
306  {
307  /*rfield.widflx[i] = rfield.anu[i+1] - rfield.anu[i];*/
308  rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
309  }
310 #endif
311 
312  if( !lgCheckOK )
313  {
314  ShowMe();
315  cdEXIT(EXIT_FAILURE);
316  }
317 
318  if( trace.lgTrace && trace.lgComBug )
319  {
320  fprintf( ioQQQ, "\n\n Compton heating, cooling coefficients \n" );
321  for( i=0; i < rfield.nupper; i += 2 )
322  {
323  fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i],
324  rfield.csigh[i], rfield.csigc[i] );
325  }
326  fprintf( ioQQQ, "\n" );
327  }
328 
329  /* option to check frequencies in real time, drive pointers command,
330  * routine is below, is file static */
331  if( trace.lgPtrace )
332  ptrcer();
333 
334  /* extinguish continuum if set on */
335  extin(&ex1ryd);
336 
337  /* now find peak of hydrogen ionizing continuum - for PDR calculations
338  * this will remain equal to 1 since the loop will not execute */
339  prt.ipeak = 1;
340  peak = 0.;
341 
342  for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1; i < rfield.nupper; i++ )
343  {
344  if( rfield.flux[0][i]*rfield.anu[i]/rfield.widflx[i] > (realnum)peak )
345  {
346  /* prt.ipeak points to largest f_nu at H-ionizing energies
347  * and is passed to other parts of code */
348  /* i+1 to keep ipeak on fortran version of energy array */
349  prt.ipeak = i+1;
350  peak = rfield.flux[0][i]*rfield.anu[i]/rfield.widflx[i];
351  }
352  }
353 
354  /* find highest energy to consider in continuum flux array
355  * peak is the peak product nu*flux */
356  peak = rfield.flux[0][prt.ipeak-1]/rfield.widflx[prt.ipeak-1]*
357  rfield.anu2[prt.ipeak-1];
358 
359  /* say what type of cpu this is, if desired */
360  if( trace.lgTrace )
361  {
362  fprintf( ioQQQ, " ContSetIntensity: The peak of the H-ion continuum is at %.3e Ryd - its value is %.2e\n",
363  rfield.anu[prt.ipeak-1] , peak);
364  }
365 
366  if( peak > 1e38 )
367  {
368  fprintf( ioQQQ, " PROBLEM DISASTER The continuum is too intense to compute. Use a fainter continuum. (This is the nu*f_nu test)\n" );
369  fprintf( ioQQQ, " Sorry.\n" );
370  cdEXIT(EXIT_FAILURE);
371  }
372 
373  /* FluxFaint set in zero.c; normally 1e-10 */
374  /* this will be faintest level of continuum we want to consider.
375  * peak was set above, is peak of hydrogen ionizing radiation field,
376  * and is zero if no H-ionizing radiation */
377  fntest = peak*rfield.FluxFaint;
378  {
379  enum {DEBUG_LOC=false};
380  /* print flux array then quit */
381  if( DEBUG_LOC )
382  {
383  for( i=0; i<rfield.nupper; ++i )
384  {
385  fprintf(ioQQQ," consetintensityBUGGG\t%.2e\t%.2e\n" ,
386  rfield.anu[i] , rfield.flux[0][i]/rfield.widflx[i] );
387  }
388  cdEXIT(EXIT_SUCCESS);
389  }
390  }
391 
392  if( fntest > 0. )
393  {
394  /* this test is not done in pdr conditions where NO H-ionizing radiation,
395  * since fntest is zero*/
396  i = rfield.nupper;
397  /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs
398  * where continuum barely goes into H-ionizing radiation */
399  while( i > prt.ipeak+3 &&
400  rfield.flux[0][i-1]*rfield.anu2[i-1]/rfield.widflx[i-1] < (realnum)fntest )
401  {
402  --i;
403  }
404  }
405  else
406  {
407  /* when no H-ionizing radiation set to Lyman edge */
408  /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs
409  * where continuum barely goes into H-ionizing radiation */
410  i = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+3;
411  }
412 
413  /*
414  * this line of code dates from 1979 and IOA Cambridge. removed july 19 95
415  * I think it was the last line of the original Cambridge Fortran source
416  nflux = MAX( ineon(1)+4 , i )
417  */
418 
419  /* >>chng 99 apr 28, reinstate the rfield.FluxFaint limit with nflux */
420  rfield.nflux = i;
421 
422  /* trim down nflux, was set to rfield.nupper, the dimension of all vectors, in zero.c,
423  * in ContCreatePointers was set to nupper, the number of cells needed to get up to the
424  * high-energy limit of the code */
425  while( rfield.flux[0][rfield.nflux-1] < SMALLFLOAT && rfield.nflux > 1 )
426  {
427  --rfield.nflux;
428  }
429  /* make sure we go high enough, avoid 1-off bugs as in draine field */
430  ++rfield.nflux;
431 
432  if( rfield.nflux == 1 )
433  {
434  fprintf( ioQQQ, " PROBLEM DISASTER This incident continuum appears to have no radiation.\n" );
435  fprintf( ioQQQ, " Sorry.\n" );
436  cdEXIT(EXIT_FAILURE);
437  }
438 
439  /* >>chng 04 oct 10, add this limit - arrays will malloc to nupper, but will add unit
440  * continuum to [nflux] - this must be within array */
441  rfield.nflux = MIN2( rfield.nupper-1 , rfield.nflux );
442 
443  /* >>chng 05 mar 01, nfine should not extend beyond continuum
444  * make sure fine opacity scale does not extend beyond continuum we will use */
445  rfield.nfine = rfield.nfine_malloc;
446  while( rfield.nfine > 0 && rfield.fine_anu[rfield.nfine-1] > rfield.anu[rfield.nflux-1] )
447  {
448  --rfield.nfine;
449  }
450 
451  /* check that continuum defined everywhere - look for zero's and comment if present */
452  continuum.lgCon0 = false;
453  ip = 0;
454  for( i=0; i < rfield.nflux; i++ )
455  {
456  if( rfield.flux[0][i] == 0. )
457  {
458  if( called.lgTalk && !continuum.lgCon0 )
459  {
460  fprintf( ioQQQ, " NOTE Setcon: continuum has zero intensity starting at %11.4e Ryd.\n",
461  rfield.anu[i] );
462  continuum.lgCon0 = true;
463  }
464  ++ip;
465  }
466  }
467 
468  if( continuum.lgCon0 && called.lgTalk )
469  {
470  fprintf( ioQQQ,
471  "%6ld cells in the incident continuum have zero intensity. Problems???\n\n",
472  ip );
473  }
474 
475  /* check for devastating error in the continuum mesh or intensity */
476  lgCheckOK = true;
477  for( i=1; i < rfield.nflux; i++ )
478  {
479  if( rfield.flux[0][i] < 0. )
480  {
481  fprintf( ioQQQ,
482  " PROBLEM DISASTER Continuum has negative intensity at %.4e Ryd=%.2e %4.4s %4.4s\n",
483  rfield.anu[i], rfield.flux[0][i], rfield.chLineLabel[i], rfield.chContLabel[i] );
484  lgCheckOK = false;
485  }
486  else if( rfield.anu[i] <= rfield.anu[i-1] )
487  {
488  fprintf( ioQQQ,
489  " PROBLEM DISASTER cont_setintensity - internal error - continuum energies not in increasing order: energies follow\n" );
490  fprintf( ioQQQ,
491  "%ld %e %ld %e %ld %e\n",
492  i -1 , rfield.anu[i-1], i, rfield.anu[i], i +1, rfield.anu[i+1] );
493  lgCheckOK = false;
494  }
495  }
496 
497  /* either of the ways lgCheckOK would be set true would be a major internal error */
498  if( !lgCheckOK )
499  {
500  ShowMe();
501  cdEXIT(EXIT_FAILURE);
502  }
503 
504  /* turn off recoil ionization if high energy < 190R */
505  if( rfield.anu[rfield.nflux-1] <= 190 )
506  {
507  ionbal.lgCompRecoil = false;
508  }
509 
510  /* sum photons and energy, save mean */
511 
512  /* sum from low energy to Balmer edge */
513  sumcon(1,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon-1,&rfield.qrad,&prt.pradio,&p1);
514 
515  /* sum over Balmer continuum */
516  sumcon(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1,&rfield.qbal,&prt.pbal,&p2);
517 
518  /* sum from Lyman edge to HeI edge */
519  sumcon(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon,
520  iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1,&prt.q,&p,&p3);
521 
522  /* sum from HeI to HeII edges */
523  sumcon(iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon,
524  iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1,&rfield.qhe,&phe,&p4);
525 
526  /* sum from Lyman edge to carbon k-shell */
527  sumcon(iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon,opac.ipCKshell-1,&rfield.qheii,&pheii,&p5);
528 
529  /* sum from c k-shell to gamma ray - where pairs start */
530  sumcon( opac.ipCKshell , rfield.ipEnerGammaRay-1 , &prt.qx,
531  &prt.xpow , &p6);
532 
533  /* complete sum up to high-energy limit */
534  sumcon(rfield.ipEnerGammaRay,rfield.nflux,&prt.qgam,&prt.GammaLumin, &p7);
535 
536  /* find to estimate photoerosion timescale */
537  n = ipoint(7.35e5);
538  sumcon(n,rfield.nflux,&qgn,&pgn,&p8);
539  timesc.TimeErode = qgn;
540 
541  /* find Compton temp */
542  tcompr = (p1 + p2 + p3 + p4 + p5 + p6 + p7)/(prt.pradio + prt.pbal +
543  p + phe + pheii + prt.xpow + prt.GammaLumin);
544 
545  tcomp = tcompr/(4.*6.272e-6);
546 
547  if( trace.lgTrace )
548  {
549  fprintf( ioQQQ,
550  " mean photon energy=%10.3eR =%10.3eK low, high nu=%12.4e%12.4e\n",
551  tcompr, tcomp, rfield.anu[0] - rfield.widflx[0]/2., rfield.anu[rfield.nflux-1] +
552  rfield.widflx[rfield.nflux-1]/2. );
553  }
554 
555  /* this is total power in ionizing radiation, per unit area */
556  prt.powion = p + phe + pheii + prt.xpow + prt.GammaLumin;
557 
558  /* this is the total radiation intensity, erg cm-2 s-1 */
559  continuum.TotalLumin = prt.pradio + prt.powion + prt.pbal;
560 
561  /* this is placed into the line stack on the first zone, then
562  * reset to zero, to end up with luminosity in the emission lines array.
563  * at end of iteration it is reset to TotalLumin */
564  continuum.totlsv = continuum.TotalLumin;
565 
566  /* total H-ionizing photon number, */
567  rfield.qhtot = prt.q + rfield.qhe + rfield.qheii + prt.qx + prt.qgam;
568 
569  /* ftotal photon number, all energies */
570  rfield.qtot = rfield.qhtot + rfield.qbal + rfield.qrad;
571 
572  if( prt.powion <= 0. && called.lgTalk )
573  {
574  rfield.lgHionRad = true;
575  fprintf( ioQQQ, " NOTE There is no hydrogen-ionizing radiation.\n" );
576  fprintf( ioQQQ, " Was this intended?\n\n" );
577  /* check if any Balmer ionizing radiation since even metals will be
578  * totally neutral */
579  if( prt.pbal <=0. && called.lgTalk )
580  {
581  fprintf( ioQQQ, " NOTE There is no Balmer continuum radiation.<<<<\n" );
582  fprintf( ioQQQ, " Was this intended?\n\n" );
583  }
584  }
585 
586  else
587  {
588  rfield.lgHionRad = false;
589  }
590 
591  /* option to add energy deposition due to fast neutrons,
592  * entered as fraction of total photon luminosity */
593  if( hextra.lgNeutrnHeatOn )
594  {
595  /* hextra.totneu is erg cm-2 s-1 in neutrons
596  * hextra.effneu - efficiency default is unity */
597  hextra.totneu = (realnum)(hextra.effneu*continuum.TotalLumin*
598  pow((realnum)10.f,hextra.frcneu));
599  }
600  else
601  {
602  hextra.totneu = (realnum)0.;
603  }
604 
605  /* temp correspond to energy density, printed in STARTR */
606  phycon.TEnerDen = pow(continuum.TotalLumin/SPEEDLIGHT/7.56464e-15,0.25);
607 
608  /* sanity check for single blackbody, that energy density temperature
609  * is smaller than black body temperature */
610  if( rfield.nShape==1 &&
611  strcmp( rfield.chSpType[rfield.nShape-1], "BLACK" )==0 )
612  {
613  /* single black body, now confirm that TEnerDen is less than this temperature,
614  * >>>chng 99 may 02,
615  * in lte these are very very close, factor of 1.00001 allows for numerical
616  * errors, and apparently slightly different atomic coef in different parts
617  * of code. eventaully all mustuse physonst.h and agree exactly */
618  if( phycon.TEnerDen > 1.0001f*rfield.slope[rfield.nShape-1] )
619  {
620  fprintf( ioQQQ,
621  "\n WARNING: The energy density temperature (%g) is greater than the"
622  " black body temperature (%g). This is unphysical.\n\n",
623  phycon.TEnerDen , rfield.slope[rfield.nShape-1]);
624  }
625  }
626 
627  /* incident continuum nu*f_nu at Hbeta and Ly-alpha */
628  continuum.cn4861 = (realnum)(ffun(0.1875)*HPLANCK*FR1RYD*0.1875*0.1875);
629  continuum.cn1216 = (realnum)(ffun(0.75)*HPLANCK*FR1RYD*0.75*0.75);
630  continuum.sv4861 = continuum.cn4861;
631  continuum.sv1216 = continuum.cn1216;
632  /* flux density in V, erg / s / cm2 / hz */
633  continuum.fluxv = (realnum)(ffun(0.1643)*HPLANCK*0.1643);
634  continuum.fbeta = (realnum)(ffun(0.1875)*HPLANCK*0.1875*6.167e14);
635 
636  /* flux density nu*Fnu = erg / s / cm2
637  * EX1RYD is optional extinction factor at 1 ryd */
638  prt.fx1ryd = (realnum)(ffun(1.000)*HPLANCK*ex1ryd*FR1RYD);
639 
640  realnum plsFrqConstant = (realnum)(ELEM_CHARGE_ESU/sqrt(PI*ELECTRON_MASS)/FR1RYD);
641  ASSERT( plsFrqConstant > 2.7e-12 && plsFrqConstant < 2.8e-12 );
642 
643  /* check for plasma frequency - then zero out incident continuum
644  * for energies below this
645  * this is critical electron density, shielding of incident continuum
646  * if electron density is greater than this */
647  ecrit = POW2(rfield.anu[0]/plsFrqConstant);
648 
649  if( dense.gas_phase[ipHYDROGEN]*1.2 > ecrit )
650  {
651  rfield.lgPlasNu = true;
652  // This should be electron density, but we don't know that yet.
653  // We use n_H (with 1.2 factor for He) as an ansatz.
654  rfield.plsfrq = plsFrqConstant*sqrt(dense.gas_phase[ipHYDROGEN]*1.2f);
655  rfield.plsfrqmax = rfield.plsfrq;
656  rfield.ipPlasma = ipoint(rfield.plsfrq);
657 
658  /* save max pointer too */
659  rfield.ipPlasmax = rfield.ipPlasma;
660 
661  /* now loop over all affected energies, setting incident continuum
662  * to zero there, and counting all as reflected */
663  /* >>chng 01 jul 14, from i < ipPlasma to ipPlasma-1 -
664  * when ipPlasma is 1 plasma freq is not on energy scale */
665  for( i=0; i < rfield.ipPlasma-1; i++ )
666  {
667  /* count as reflected incident continuum */
668  rfield.ConRefIncid[0][i] = rfield.flux[0][i];
669  /* set continuum to zero there */
670  rfield.flux_beam_const[i] = 0.;
671  rfield.flux_beam_time[i] = 0.;
672  rfield.flux_isotropic[i] = 0.;
673  rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
674  rfield.flux_isotropic[i];
675  }
676  }
677  else
678  {
679  rfield.lgPlasNu = false;
680  /* >>chng 01 jul 14, from 0 to 1 - 1 is the first array element on the F scale,
681  * ipoint would return this, so rest of code assumes ipPlasma is 1 plus correct index */
682  rfield.ipPlasma = 1;
683  rfield.plsfrqmax = 0.;
684  rfield.plsfrq = 0.;
685  }
686 
687  if( rfield.ipPlasma > 1 && called.lgTalk )
688  {
689  fprintf( ioQQQ,
690  " !The plasma frequency is %.2e Ryd. The incident continuum is set to 0 below this.\n",
691  rfield.plsfrq );
692  }
693 
694  rfield.occmax = 0.;
695  rfield.tbrmax = 0.;
696  for( i=0; i < rfield.nupper; i++ )
697  {
698  /* set up occupation number array */
699  rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
700  if( rfield.OccNumbIncidCont[i] > rfield.occmax )
701  {
702  rfield.occmax = rfield.OccNumbIncidCont[i];
703  rfield.occmnu = rfield.anu[i];
704  }
705  /* following product is continuum brightness temperature */
706  if( rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i] > rfield.tbrmax )
707  {
708  rfield.tbrmax = (realnum)(rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i]);
709  rfield.tbrmnu = rfield.anu[i];
710  }
711  /* save continuum for next iteration */
712  rfield.flux_total_incident[0][i] = rfield.flux[0][i];
713  rfield.flux_beam_const_save[i] = rfield.flux_beam_const[i];
714  rfield.flux_time_beam_save[i] = rfield.flux_beam_time[i];
715  rfield.flux_isotropic_save[i] = rfield.flux_isotropic[i];
716  /*fprintf(ioQQQ,"DEBUG type cont %li\t%.3e\t%.2e\t%.2e\t%.2e\t%.2e\n",
717  i, rfield.anu[i],
718  rfield.flux[0][i],rfield.flux_beam_const[i],rfield.flux_beam_time[i],
719  rfield.flux_isotropic[i]);
720  fflush(ioQQQ);*/
721  }
722 
723  /* if continuum brightness temp is large, where does it fall below
724  * 1e4k??? */
725  if( rfield.tbrmax > 1e4 )
726  {
727  i = ipoint(rfield.tbrmnu)-1;
728  while( i < rfield.nupper-1 && (rfield.OccNumbIncidCont[i]*TE1RYD*
729  rfield.anu[i] > 1e4) )
730  {
731  ++i;
732  }
733  rfield.tbr4nu = rfield.anu[i];
734  }
735  else
736  {
737  rfield.tbr4nu = 0.;
738  }
739 
740  /* if continuum occ num is large, where does it fall below 1? */
741  if( rfield.occmax > 1. )
742  {
743  i = ipoint(rfield.occmnu)-1;
744  while( i < rfield.nupper && (rfield.OccNumbIncidCont[i] > 1.) )
745  {
746  ++i;
747  }
748  rfield.occ1nu = rfield.anu[i];
749  }
750  else
751  {
752  rfield.occ1nu = 0.;
753  }
754 
755  /* remember if incident radiation field is less than 10*Habing ISM */
756  /* >>chng 06 aug 01, change this test from continuum.TotalLumin to
757  * energy in balmer and ionizing continuum, since this is the true habing field
758  * and is the continuum that interacts with gas. When CMB set this
759  * tests on total did not trigger due to cold blackbody, which has little
760  * effect on gas, other than compton */
761  if( (prt.powion + prt.pbal) < 1.8e-2 )
762  {
763  /* thermal.ConstTemp def is zero, set pos when constant temperature is set */
764  rfield.lgHabing = true;
765  /* >>chng 06 aug 01 also print warning if substantially below Habing, this may be a mistake */
766  if( ((prt.powion + prt.pbal) < 1.8e-12) &&
767  /* this is test for not constant temperature */
768  (!thermal.lgTemperatureConstant) )
769  {
770  fprintf( ioQQQ, "\n >>>\n"
771  " >>> NOTE The incident continuum is surprisingly faint.\n" );
772  fprintf( ioQQQ,
773  " >>> The total energy in the Balmer and Lyman continua is %.2e erg cm-2 s-1.\n"
774  ,(prt.powion + prt.pbal));
775  fprintf( ioQQQ, " >>> This is many orders of magnitude fainter than the ISM galactic background.\n" );
776  fprintf( ioQQQ, " >>> This seems unphysical - please check that the continuum intensity has been properly set.\n" );
777  fprintf( ioQQQ, " >>> YOU MAY BE MAKING A BIG MISTAKE!!\n >>>\n\n\n\n" );
778  }
779  }
780 
781  /* fix ionization parameter (per hydrogen) at inner edge */
782  rfield.uh = (realnum)(rfield.qhtot/
783  (dense.gas_phase[ipHYDROGEN]*SPEEDLIGHT));
784  rfield.uheii = (realnum)((rfield.qheii + prt.qx)/
785  (dense.gas_phase[ipHYDROGEN]*SPEEDLIGHT));
786  if( rfield.uh > 1e10 )
787  {
788  fprintf( ioQQQ, "\n\n"
789  " CAUTION The incident radiation field is surprisingly intense.\n" );
790  fprintf( ioQQQ,
791  " The dimensionless hydrogen ionization parameter is %.2e.\n"
792  , rfield.uh );
793  fprintf( ioQQQ, " This is many orders of magnitude brighter than commonly seen.\n" );
794  fprintf( ioQQQ, " This seems unphysical - please check that the radiation field intensity has been properly set.\n" );
795  fprintf( ioQQQ, " YOU MAY BE MAKING A BIG MISTAKE!!\n\n\n\n\n" );
796  }
797 
798  /* guess first temperature and neutral h density */
799  double TeNew;
800  if( thermal.ConstTemp > 0. )
801  {
802  TeNew = thermal.ConstTemp;
803  }
804  else
805  {
806  if( rfield.uh > 0. )
807  {
808  TeNew = (20000.+log10(rfield.uh)*5000.);
809  TeNew = MAX2(8000. , TeNew );
810  }
811  else
812  {
813  TeNew = 1000.;
814  }
815  }
816  TempChange( TeNew );
817 
818  /* this is an option to stop after printing header only */
819  if( noexec.lgNoExec )
820  return;
821 
822  /* estimate secondary ionization rate - probably 0, but possible extra
823  * SetCsupra set with "set csupra" */
824  /* >>>chng 99 apr 29, added cosmic ray ionization since this is used to get
825  * helium ionization fraction, and was zero in pdr, so He turned off at start,
826  * and never turned back on */
827  /* coef on cryden is from highen.c */
828  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
829  {
830  for( ion=0; ion<nelem+1; ++ion )
831  {
832  secondaries.csupra[nelem][ion] =
833  secondaries.SetCsupra + hextra.cryden*2e-9f;
834  }
835  }
836 
837  /*********************************************************************
838  * *
839  * estimate hydrogen's level of ionization *
840  * *
841  *********************************************************************/
842 
843  /* create fake ionization balance, but will conserve number of hydrogens */
844  dense.xIonDense[ipHYDROGEN][0] = 0.;
845  dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN];
846  /* this must be zero since PresTotCurrent will do radiation pressure due to H */
847  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = 0.;
848 
849  /* "extra" electrons from command line, or assumed residual electrons */
850  double EdenExtraLocal = dense.EdenExtra +
851  /* if we are in a molecular cloud the current logic could badly fail
852  * do not let electron density fall below 1e-7 of H density */
853  1e-7*dense.gas_phase[ipHYDROGEN];
854  EdenChange( dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal );
855 
856  /* hydrogen case B recombination coefficient */
857  HCaseBRecCoeff = (-9.9765209 + 0.158607055*phycon.telogn[0] + 0.30112749*
858  phycon.telogn[1] - 0.063969007*phycon.telogn[2] + 0.0012691546*
859  phycon.telogn[3])/(1. + 0.035055422*phycon.telogn[0] -
860  0.037621619*phycon.telogn[1] + 0.0076175048*phycon.telogn[2] -
861  0.00023432613*phycon.telogn[3]);
862  HCaseBRecCoeff = pow(10.,HCaseBRecCoeff)/phycon.te;
863 
864  double CollIoniz = t_ADfA::Inst().coll_ion_wrapper(0,0,phycon.te);
865 
866  // mean H-ionizing photon energy
867  double PhotonEnergy = 1.;
868  if( rfield.qhtot>SMALLFLOAT )
869  PhotonEnergy = prt.powion/rfield.qhtot/EN1RYD;
870 
871  double OtherIonization = rfield.qhtot*6.3e-18/POW3(PhotonEnergy) +
872  secondaries.csupra[ipHYDROGEN][0];
873 
874  double newEden = dense.eden;
875  long loopCount = 0;
876 
877  do
878  {
879  /* update electron density */
880  EdenChange( newEden );
881  double RatioIoniz =
882  (CollIoniz*dense.eden+OtherIonization)/(HCaseBRecCoeff*dense.eden);
883  if( RatioIoniz<1e-3 )
884  {
885  /* very low ionization solution */
886  dense.xIonDense[ipHYDROGEN][1] = (realnum)(
887  dense.gas_phase[ipHYDROGEN]*RatioIoniz);
888  dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
889  dense.xIonDense[ipHYDROGEN][1];
890  ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
891  dense.xIonDense[ipHYDROGEN][0]<=dense.gas_phase[ipHYDROGEN]);
892  ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
893  dense.xIonDense[ipHYDROGEN][1]<dense.gas_phase[ipHYDROGEN]);
894  //fprintf(ioQQQ,"DEBUG br 1 H0 %.2e\n", dense.xIonDense[ipHYDROGEN][0]);
895  }
896  else if( RatioIoniz>1e3 )
897  {
898  /* very high ionization solution */
899  dense.xIonDense[ipHYDROGEN][0] = (realnum)(
900  dense.gas_phase[ipHYDROGEN]/RatioIoniz);
901  dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN] -
902  dense.xIonDense[ipHYDROGEN][0];
903  ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
904  dense.xIonDense[ipHYDROGEN][0]<dense.gas_phase[ipHYDROGEN]);
905  ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
906  dense.xIonDense[ipHYDROGEN][1]<=dense.gas_phase[ipHYDROGEN]);
907  //fprintf(ioQQQ,"DEBUG br 2 H0 %.2e rat %.2e\n", dense.xIonDense[ipHYDROGEN][0],
908  // dense.gas_phase[ipHYDROGEN]/RatioIoniz);
909  }
910  else
911  {
912  /* intermediate ionization - solve quadratic */
913  double alpha = HCaseBRecCoeff + CollIoniz ;
914  double beta = HCaseBRecCoeff*EdenExtraLocal + OtherIonization +
915  (EdenExtraLocal - dense.gas_phase[ipHYDROGEN])*CollIoniz;
916  double gamma = -dense.gas_phase[ipHYDROGEN]*(OtherIonization+EdenExtraLocal*CollIoniz);
917 
918  double discriminant = POW2(beta) - 4.*alpha*gamma;
919  if( discriminant <0 )
920  {
921  /* oops */
922  fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found negative discriminant.\n");
923  TotalInsanity();
924  }
925 
926  dense.xIonDense[ipHYDROGEN][1] = (realnum)((-beta + sqrt(discriminant))/(2.*alpha));
927  if( dense.xIonDense[ipHYDROGEN][1]> dense.gas_phase[ipHYDROGEN] )
928  {
929  /* oops */
930  fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H+)>n(H).\n");
931  TotalInsanity();
932  }
933  dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
934  dense.xIonDense[ipHYDROGEN][1];
935  if( dense.xIonDense[ipHYDROGEN][0]<= 0 )
936  {
937  /* oops */
938  fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H0)<0.\n");
939  TotalInsanity();
940  }
941  //fprintf(ioQQQ,"DEBUG br 3 H0 %.2e\n", dense.xIonDense[ipHYDROGEN][0]);
942  }
943 
944 
945  if( dense.xIonDense[ipHYDROGEN][1] > 1e-30 )
946  {
947  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = dense.xIonDense[ipHYDROGEN][0];
948  }
949  else
950  {
951  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = 0.;
952  }
953 
954  /* now save estimates of whether induced recombination is going
955  * to be important -this is a code pacesetter since GammaBn is slower
956  * than GammaK */
957  hydro.lgHInducImp = false;
958  for( i=ipH1s; i < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max; i++ )
959  {
960  if( rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].ipIsoLevNIonCon-1] > 0.01 )
961  hydro.lgHInducImp = true;
962  }
963 
964  /*******************************************************************
965  * *
966  * estimate helium's level of ionization *
967  * *
968  *******************************************************************/
969 
970  /* only if helium is turned on */
971  if( dense.lgElmtOn[ipHELIUM] )
972  {
973  /* next estimate level of helium singly ionized */
974  xIoniz = (realnum)t_ADfA::Inst().coll_ion_wrapper(1,0,phycon.te);
975  /* >>chng 05 jan 05, add cosmic rays */
976  xIoniz = (realnum)(xIoniz*dense.eden + rfield.qhe*1e-18 + secondaries.csupra[ipHELIUM][1] );
977  double RecTot = HCaseBRecCoeff*dense.eden;
978  RatioIonizToRecomb = xIoniz/RecTot;
979 
980  /* now estimate level of helium doubly ionized */
981  xIoniz = (realnum)t_ADfA::Inst().coll_ion_wrapper(1,1,phycon.te);
982  /* >>chng 05 jan 05, add cosmic rays */
983  xIoniz = (realnum)(xIoniz*dense.eden + rfield.qheii*1e-18 + ionbal.CosRayIonRate );
984 
985  /* rough charge dependence */
986  RecTot *= 4.;
987  r3ov2 = xIoniz/RecTot;
988 
989  /* now set level of helium ionization */
990  if( RatioIonizToRecomb > 0. )
991  {
992  double r1 = dense.gas_phase[ipHELIUM]/(1./RatioIonizToRecomb + 1. + r3ov2);
993  dense.xIonDense[ipHELIUM][1] = (realnum)(r1);
994  dense.xIonDense[ipHELIUM][0] = (realnum)(r1/RatioIonizToRecomb);
995  dense.xIonDense[ipHELIUM][2] = (realnum)(r1*r3ov2);
996  }
997  else
998  {
999  /* no He ionizing radiation */
1000  dense.xIonDense[ipHELIUM][1] = dense.xIonDense[ipHELIUM][2] = 0.;
1001  dense.xIonDense[ipHELIUM][0] = dense.gas_phase[ipHELIUM];
1002  }
1003 
1004  if( dense.xIonDense[ipHELIUM][2] > 1e-30 )
1005  {
1006  iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = dense.xIonDense[ipHELIUM][1];
1007  }
1008  else
1009  {
1010  iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = 0.;
1011  }
1012  }
1013  else
1014  {
1015  /* case where helium is turned off */
1016  dense.xIonDense[ipHELIUM][1] = 0.;
1017  dense.xIonDense[ipHELIUM][0] = 0.;
1018  dense.xIonDense[ipHELIUM][2] = 0.;
1019  iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = 0.;
1020  }
1021 
1022  /* update electron density */
1023  newEden = dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal + dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2];
1024 
1025  loopCount++;
1026  }
1027  /* repeat above until guessed and calculated eden agree to at least 0.1%. */
1028  while( loopCount < 10 && fabs(newEden/dense.eden - 1.) > 0.001 );
1029 
1030  if( dense.xIonDense[ipHYDROGEN][0]<0.)
1031  TotalInsanity();
1032  else if( dense.xIonDense[ipHYDROGEN][0] == 0. )
1033  {
1034  fprintf(ioQQQ,"PROBLEM the derived atomic hydrogen density is zero.\n");
1035  if( dense.gas_phase[ipHYDROGEN]<1e-5 && rfield.uh > 1e10)
1036  {
1037  fprintf(ioQQQ,"This is almost certainly due to floating point "
1038  "limits on this computer.\nThe ionization parameter is very large,"
1039  " the density is very small,\nand the H^0 density cannot be"
1040  " stored in a float.\n");
1041  //cdEXIT( EXIT_FAILURE );
1042  }
1043  }
1044  //ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);
1045 
1046  /* update electron density */
1047  EdenChange( newEden );
1048 
1049  if( dense.eden <= SMALLFLOAT )
1050  {
1051  /* no electrons! */
1052  fprintf(ioQQQ,"\n PROBLEM DISASTER - this simulation has no source"
1053  " of ionization. The electron density is zero. Consider "
1054  "adding a source of ionization such as cosmic rays.\n\n");
1055  cdEXIT(EXIT_FAILURE);
1056  }
1057 
1058  /* fix range of stages of ionization */
1059  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
1060  {
1061  if( dense.lgElmtOn[nelem] )
1062  {
1063  // IonHigh[n] is the highest stage of ionization present
1064  // the IonHigh array index is on the C scale, so [0] is hydrogen
1065  // the value is also on the C scale, so element [nelem] can range
1066  // from 0 to nelem+1
1067  dense.IonHigh[nelem] = nelem + 1;
1068 
1069  dense.IonLow[nelem] = 0;
1070  // for very intense radiation fields very heavy elements, N>Fe, will fail
1071  // in ion_solver due to ill conditioned matrix. all populations are in
1072  // fully stripped state. Start will fully stripped ion distribution in this case.
1073  if( rfield.uh > 1e15 )
1074  {
1075  //trim down highest stage to be within incident radiation field
1076  while ( rfield.anu[Heavy.ipHeavy[nelem][dense.IonHigh[nelem]-1]] >
1077  rfield.anu[rfield.nflux] && dense.IonHigh[nelem]>1 )
1078  --dense.IonHigh[nelem];
1079 
1080  dense.IonLow[nelem] = max( 0 , dense.IonHigh[nelem]-1 );
1081  }
1082 
1083  /* >>chng 04 jan 13, add this test, caught by Orly Gnat */
1084  /* check on actual zero abundances of lower stages - this will only
1085  * happen when ionization is set with element ionization command */
1086  if( dense.lgSetIoniz[nelem] )
1087  {
1088  while( dense.SetIoniz[nelem][dense.IonLow[nelem]] < dense.density_low_limit )
1089  ++dense.IonLow[nelem];
1090  while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] < dense.density_low_limit )
1091  --dense.IonHigh[nelem];
1092  }
1093 
1094  // make low-stage populations zero
1095  for( ion=0; ion<dense.IonLow[nelem]; ++ion )
1096  {
1097  dense.xIonDense[nelem][dense.IonLow[nelem]] += dense.xIonDense[nelem][ion];
1098  dense.xIonDense[nelem][ion] = 0.;
1099  }
1100  for( ion=nelem+1; ion>dense.IonHigh[nelem]; --ion )
1101  {
1102  dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion];
1103  dense.xIonDense[nelem][ion] = 0.;
1104  }
1105  }
1106  else
1107  {
1108  /* this element is turned off, make stages impossible */
1109  dense.IonLow[nelem] = -1;
1110  dense.IonHigh[nelem] = -1;
1111  }
1112  }
1113 
1114  // make first estimate of iso continuum lowering
1115  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1116  {
1117  for( long nelem=ipISO; nelem<LIMELM; ++nelem)
1118  {
1119  if( nelem < 2 || dense.lgElmtOn[nelem] )
1120  {
1121  if( iso_ctrl.lgContinuumLoweringEnabled[ipISO] )
1122  iso_continuum_lower( ipISO, nelem );
1123  }
1124  }
1125  }
1126 
1127  /* estimate electrons from heavies, assuming each at least
1128  * 1 times ionized */
1129  EdenHeav = 0.;
1130  realnum atomFrac = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
1131  realnum firstIonFrac = dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN];
1132  for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
1133  {
1134  if( dense.lgElmtOn[nelem] )
1135  {
1136  if( dense.IonLow[nelem] == 0 && dense.IonHigh[nelem] >=1)
1137  {
1138  realnum low2dens = dense.xIonDense[nelem][0]+dense.xIonDense[nelem][1];
1139  dense.xIonDense[nelem][0] = low2dens * atomFrac;
1140  dense.xIonDense[nelem][1] = low2dens * firstIonFrac;
1141  }
1142  for (long ion=1;ion<dense.IonHigh[nelem];++ion)
1143  {
1144  EdenHeav += ion*dense.xIonDense[nelem][ion];
1145  }
1146  }
1147  }
1148 
1149  /* modify estimate of electron density */
1150  dense.eden += EdenHeav;
1151 
1152  /* >>chng 05 jan 05, insure positive eden */
1153  EdenChange( MAX2( SMALLFLOAT , dense.eden ) );
1154 
1155  if( dense.EdenSet > 0. )
1156  {
1157  EdenChange( dense.EdenSet );
1158  }
1159 
1160  dense.EdenHCorr = dense.eden;
1161  dense.EdenHCorr_f = (realnum)dense.EdenHCorr;
1162 
1163  if( dense.eden < 0. )
1164  {
1165  fprintf( ioQQQ, " PROBLEM DISASTER Negative electron density results in ContSetIntensity.\n" );
1166  fprintf( ioQQQ, "%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n",
1167  dense.eden, dense.xIonDense[ipHYDROGEN][1], dense.xIonDense[ipHELIUM][1],
1168  dense.xIonDense[ipHELIUM][2], dense.gas_phase[ipCARBON], dense.EdenExtra );
1169  TotalInsanity();
1170  }
1171 
1172  if( dense.EdenSet > 0. )
1173  {
1174  dense.EdenTrue = dense.EdenSet;
1175  }
1176  else
1177  dense.EdenTrue = dense.eden;
1178 
1179  if( trace.lgTrace )
1180  {
1181  fprintf( ioQQQ,
1182  " ContSetIntensity sets initial EDEN to %.4e, contributors H+=%.2e He+, ++= %.2e %.2e Heav %.2e extra %.2e\n",
1183  dense.eden ,
1184  dense.xIonDense[ipHYDROGEN][1],
1185  dense.xIonDense[ipHELIUM][1],
1186  2.*dense.xIonDense[ipHELIUM][2],
1187  EdenHeav,
1188  dense.EdenExtra);
1189  }
1190 
1191  // photon occupation number at 1 Ryd - used for printout
1192  occ1 = (realnum)(prt.fx1ryd/HNU3C2/PI4/FR1RYD);
1193  if( occ1 > 1. )
1194  rfield.lgOcc1Hi = true;
1195  else
1196  rfield.lgOcc1Hi = false;
1197 
1198  if( trace.lgTrace && trace.lgConBug )
1199  {
1200  fprintf(ioQQQ,"\ntrace continuum print of %li incident spectral "
1201  "components\n", rfield.nShape);
1202  fprintf(ioQQQ," # type Illum Beam? 1/cos TimeVary?\n");
1203  for( rfield.ipSpec=0; rfield.ipSpec<rfield.nShape; ++rfield.ipSpec )
1204  {
1205  fprintf(ioQQQ,"%3li %6s %5i %c %.3f %c\n",
1206  rfield.ipSpec ,
1207  rfield.chSpType[rfield.ipSpec],
1208  rfield.Illumination[rfield.ipSpec],
1209  TorF(rfield.lgBeamed[rfield.ipSpec]),
1210  rfield.OpticalDepthScaleFactor[rfield.ipSpec],
1211  TorF(rfield.lgTimeVary[rfield.ipSpec]));
1212  }
1213  fprintf(ioQQQ,"\n");
1214 
1215  /* print some useful pointers to ionization edges */
1216  fprintf( ioQQQ, " H2,1=%5ld%5ld NX=%5ld IRC=%5ld\n",
1217  iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon,
1218  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon,
1219  opac.ipCKshell,
1220  ionbal.ipCompRecoil[ipHYDROGEN][0] );
1221  fprintf( ioQQQ, " CARBON" );
1222  for( i=0; i < 6; i++ )
1223  fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipCARBON][i] );
1224  fprintf( ioQQQ, "\n" );
1225 
1226  fprintf( ioQQQ, " OXY" );
1227  for( i=0; i < 8; i++ )
1228  fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipOXYGEN][i] );
1229  fprintf( ioQQQ, "%5ld%5ld%5ld\n", opac.ipo3exc[0],
1230  oxy.i2d, oxy.i2p );
1231 
1232  double sum = 0.;
1233  ASSERT( rfield.ipG0_DB96_lo>0 );
1234  for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; i++ )
1235  sum += rfield.flux[0][i];
1236  fprintf(ioQQQ,"\n Sum DB96 photons %.2e", sum);
1237  sum = 0.;
1238  for( i=(Heavy.ipHeavy[ipHYDROGEN][0] - 1); i<rfield.nflux; i++ )
1239  sum += rfield.flux[0][i];
1240  fprintf(ioQQQ,", sum H-ioniz photons %.2e\n", sum);
1241 
1242 
1243  fprintf( ioQQQ, "\n\n PHOTONS PER CELL (NOT RYD)\n" );
1244  fprintf( ioQQQ, " nu, flux, wid, occ \n" );
1245  for( i=0; i < rfield.nflux; i++ )
1246  {
1247  fprintf( ioQQQ, "%4ld%10.2e%10.2e%10.2e%10.2e\n", i,
1248  rfield.anu[i], rfield.flux[0][i], rfield.widflx[i],
1249  rfield.OccNumbIncidCont[i] );
1250  }
1251  fprintf( ioQQQ, " \n" );
1252  }
1253 
1254  /* zero out some continua related to the ots rates,
1255  * prototype and routine in RT_OTS_Update. This is done here since summed cont will
1256  * be set to rfield */
1257  RT_OTS_Zero();
1258 
1259  if( trace.lgTrace )
1260  {
1261  fprintf( ioQQQ, " ContSetIntensity returns, nflux=%5ld anu(nflux)=%11.4e eden=%10.2e\n",
1262  rfield.nflux, rfield.anu[rfield.nflux-1], dense.eden );
1263  }
1264 
1265  return;
1266 }
1267 
1268 /*sumcon sums L and Q for net incident continuum */
1269 STATIC void sumcon(long int il,
1270  long int ih,
1271  realnum *q,
1272  realnum *p,
1273  realnum *panu)
1274 {
1275  long int i,
1276  iupper; /* used as upper limit to the sum */
1277 
1278  DEBUG_ENTRY( "sumcon()" );
1279 
1280  *q = 0.;
1281  *p = 0.;
1282  *panu = 0.;
1283 
1284  /* soft continua may not go as high as the requested bin */
1285  iupper = MIN2(rfield.nflux,ih);
1286 
1287  /* n.b. - in F77 loop IS NOT executed when IUPPER < IL */
1288  for( i=il-1; i < iupper; i++ )
1289  {
1290  /* sum photon number */
1291  *q += rfield.flux[0][i];
1292  /* en1ryd is needed to stop overflow */
1293  /* sum flux */
1294  *p += (realnum)(rfield.flux[0][i]*(rfield.anu[i]*EN1RYD));
1295  /* this sum needed for means */
1296  *panu += (realnum)(rfield.flux[0][i]*(rfield.anu2[i]*EN1RYD));
1297  }
1298 
1299  return;
1300 }
1301 
1302 /*ptrcer show continuum pointers in real time following drive pointers command */
1303 STATIC void ptrcer()
1304 {
1305  char chCard[INPUT_LINE_LENGTH];
1306  /* in case of checking everything, will write errors to this file */
1307  FILE * ioERRORS=NULL;
1308  bool lgEOL;
1309  char chKey;
1310  long int i,
1311  ipnt,
1312  j;
1313  double pnt,
1314  t1,
1315  t2;
1316 
1317  DEBUG_ENTRY( "ptrcer()" );
1318 
1319  fprintf( ioQQQ, " There are two ways to do this:\n");
1320  fprintf( ioQQQ, " do you want me to test all the pointers (enter y)\n");
1321  fprintf( ioQQQ, " or do you want to enter energies yourself? (enter n)\n" );
1322 
1323  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
1324  {
1325  fprintf( ioQQQ, " error getting input \n" );
1326  cdEXIT(EXIT_FAILURE);
1327  }
1328 
1329  /* this must be either y or n */
1330  chKey = chCard[0];
1331 
1332  if( chKey == 'n' )
1333  {
1334  /* this branch, enter energies by hand, and see what happens */
1335  fprintf( ioQQQ, " Enter energy (Ryd); 0 to stop; negative is log.\n" );
1336  pnt = 1.;
1337  while( pnt!=0. )
1338  {
1339  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
1340  {
1341  fprintf( ioQQQ, " error getting input2 \n" );
1342  cdEXIT(EXIT_FAILURE);
1343  }
1344  /* now get the number off the line */
1345  i = 1;
1346  pnt = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
1347 
1348  /* bail if no number at all, or it is zero*/
1349  if( lgEOL || pnt==0. )
1350  {
1351  break;
1352  }
1353 
1354  /* if number negative then interpret as log */
1355  if( pnt < 0. )
1356  {
1357  pnt = pow(10.,pnt);
1358  }
1359 
1360  /* get pointer to call */
1361  ipnt = ipoint(pnt);
1362  fprintf( ioQQQ, " Cell num%4ld center:%10.2e width:%10.2e low:%10.2e hi:%10.2e convoc:%10.2e\n",
1363  ipnt, rfield.anu[ipnt-1], rfield.widflx[ipnt-1],
1364  rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]/2.,
1365  rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]/2.,
1366  rfield.convoc[ipnt-1] );
1367  }
1368  }
1369 
1370  else if( chKey == 'y' )
1371  {
1372  /* first check that ipoint will not crash due to out of range call*/
1373  if( rfield.anu[0] - rfield.widflx[0]/2.*0.9 < continuum.filbnd[0] )
1374  {
1375  fprintf( ioQQQ," ipoint would crash since lowest desired energy of %e ryd is below limit of %e\n",
1376  rfield.anu[0] - rfield.widflx[0]/2.*0.9 , continuum.filbnd[0] );
1377  fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[0]);
1378  cdEXIT(EXIT_FAILURE);
1379  }
1380 
1381  else if( rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 >
1382  continuum.filbnd[continuum.nrange] )
1383  {
1384  fprintf( ioQQQ," ipoint would crash since highest desired energy of %e ryd is above limit of %e\n",
1385  rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 ,
1386  continuum.filbnd[continuum.nrange-1] );
1387  fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[rfield.nflux]);
1388  fprintf( ioQQQ," this, previous cells are %e %e\n",
1389  rfield.anu[rfield.nflux-1],rfield.anu[rfield.nflux-2]);
1390  cdEXIT(EXIT_FAILURE);
1391  }
1392 
1393  /* this branch check everything, write errors to error file */
1394  fprintf( ioQQQ, " errors output on errors.txt\n");
1395  fprintf( ioQQQ, " IP(cor),IP(fount),nu lower, upper of found, desired cell.\n" );
1396 
1397  /* error file not open, set to null so we can check later */
1398  ioERRORS = NULL;
1399  for( i=0; i < rfield.nflux-1; i++ )
1400  {
1401  t1 = rfield.anu[i] - rfield.widflx[i]/2.*0.9;
1402  t2 = rfield.anu[i] + rfield.widflx[i]/2.*0.9;
1403 
1404  j = ipoint(t1);
1405  if( j != i+1 )
1406  {
1407  /* open file for errors if not already open */
1408  if( ioERRORS == NULL )
1409  {
1410  ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY );
1411  fprintf( ioQQQ," created errors.txt file with error summary\n");
1412  }
1413 
1414  fprintf( ioQQQ, " Pointers do not agree for lower bound of cell%4ld, %e\n",
1415  i, rfield.anu[i]);
1416  fprintf( ioERRORS, " Pointers do not agree for lower bound of cell%4ld, %e\n",
1417  i, rfield.anu[i] );
1418  }
1419 
1420  j = ipoint(t2);
1421  if( j != i+1 )
1422  {
1423  /* open file for errors if not already open */
1424  if( ioERRORS == NULL )
1425  {
1426  ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY );
1427  fprintf( ioQQQ," created errors.txt file with error summary\n");
1428  }
1429  fprintf( ioQQQ, " Pointers do not agree for upper bound of cell%4ld, %e\n",
1430  i , rfield.anu[i]);
1431  fprintf( ioERRORS, " Pointers do not agree for upper bound of cell%4ld, %e\n",
1432  i , rfield.anu[i]);
1433  }
1434 
1435  }
1436  }
1437 
1438  else
1439  {
1440  fprintf( ioQQQ, "I do not understand this key, sorry. %c\n", chKey );
1441  cdEXIT(EXIT_FAILURE);
1442  }
1443 
1444  if( ioERRORS!=NULL )
1445  fclose( ioERRORS );
1446  cdEXIT(EXIT_FAILURE);
1447 }
1448 
1449 /*extin do extinction of incident continuum as set by extinguish command */
1450 STATIC void extin(realnum *ex1ryd)
1451 {
1452 
1453  DEBUG_ENTRY( "extin()" );
1454 
1455  /* modify input continuum by leaky absorber
1456  * power law fit to
1457  * >>refer XUV extinction Cruddace, R., Paresce, F., Bowyer, S., & Lampton, M. 1974, ApJ, 187, 497. */
1458  if( rfield.ExtinguishColumnDensity == 0. )
1459  {
1460  *ex1ryd = 1.;
1461 
1462  for( long i=0; i<rfield.nupper; ++i )
1463  rfield.ExtinguishFactor[i] = 1.;
1464  }
1465  else
1466  {
1467  double absorb = 1. - rfield.ExtinguishLeakage;
1468  double factor = rfield.ExtinguishColumnDensity*
1469  rfield.ExtinguishConvertColDen2OptDepth;
1470  /* extinction at 1 4 Ryd */
1471  *ex1ryd = (realnum)(rfield.ExtinguishLeakage + absorb*sexp(factor));
1472 
1473  // low energy limit of extinction
1474  long low = ipoint(rfield.ExtinguishLowEnergyLimit);
1475 
1476  for( long i=0; i<low-1; ++i )
1477  rfield.ExtinguishFactor[i] = 1.;
1478 
1479  for( long i=low-1; i < rfield.nupper; i++ )
1480  {
1481  realnum extfactor = (realnum)(rfield.ExtinguishLeakage + absorb*
1482  sexp(factor*(pow(rfield.anu[i],(double)rfield.ExtinguishEnergyPowerLow))));
1483 
1484  rfield.ExtinguishFactor[i] = extfactor;
1485  rfield.flux_beam_const[i] *= extfactor;
1486  rfield.flux_beam_time[i] *= extfactor;
1487  rfield.flux_isotropic[i] *= extfactor;
1488 
1489  rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
1490  rfield.flux_isotropic[i];
1491  }
1492  }
1493  return;
1494 }
1495 
1496 /*conorm normalize continuum to proper intensity */
1497 STATIC void conorm()
1498 {
1499  long int i;
1500  double xLog_radius_inner,
1501  diff,
1502  f,
1503  flx1,
1504  flx2,
1505  pentrd,
1506  qentrd;
1507 
1508  DEBUG_ENTRY( "conorm()" );
1509 
1510  xLog_radius_inner = log10(radius.rinner);
1511 
1512  /* this is a sanity check, it can't happen */
1513  for( i=0; i < rfield.nShape; i++ )
1514  {
1515  if( strcmp(rfield.chRSpec[i],"UNKN") == 0 )
1516  {
1517  fprintf( ioQQQ, " UNKN spectral normalization cannot happen.\n" );
1518  fprintf( ioQQQ, " conorm punts.\n" );
1519  cdEXIT(EXIT_FAILURE);
1520  }
1521 
1522  else if( strcmp(rfield.chRSpec[i],"SQCM") != 0 &&
1523  strcmp(rfield.chRSpec[i],"4 PI") != 0 )
1524  {
1525  fprintf( ioQQQ, " chRSpec must be SQCM or 4 PI, and it was %4.4s. This cannot happen.\n",
1526  rfield.chRSpec[i] );
1527  fprintf( ioQQQ, " conorm punts.\n" );
1528  cdEXIT(EXIT_FAILURE);
1529  }
1530 
1531 
1532  /* this sanity check makes sure that atlas.mod or werner.mod grids
1533  * are for the current version of the code */
1534  if( strcmp(rfield.chSpType[i],"VOLK ") == 0 )
1535  {
1536  if( !fp_equal( rfield.RSFCheck[i], continuum.ResolutionScaleFactor ) )
1537  {
1538  fprintf( ioQQQ,"\n\n PROBLEM DISASTER At least one of the compiled stellar atmosphere"
1539  " grids has been compiled with a different energy grid resolution factor.\n" );
1540  fprintf( ioQQQ, " Please recompile this file using the COMPILE STARS command "
1541  "and make sure that you use the correct SET CONTINUUM RESOLUTION factor.\n" );
1542  cdEXIT(EXIT_FAILURE);
1543  }
1544  }
1545  else if( strcmp(rfield.chSpType[i],"READ ") == 0 )
1546  {
1547  if( !fp_equal( rfield.RSFCheck[i], continuum.ResolutionScaleFactor ) )
1548  {
1549  fprintf( ioQQQ,"\n\n PROBLEM DISASTER The file read by the TABLE READ command "
1550  "has been compiled with a different energy grid resolution factor.\n" );
1551  fprintf( ioQQQ, " Please recompile this file using the SAVE TRANSMITTED CONTINUUM "
1552  "command and use the correct SET CONTINUUM RESOLUTION factor.\n" );
1553  cdEXIT(EXIT_FAILURE);
1554  }
1555  }
1556  }
1557 
1558  /* this sanity check is that the grains we have read in from opacity files agree
1559  * with the energy grid in this version of cloudy */
1560  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1561  {
1562  if( !fp_equal( gv.bin[nd]->RSFCheck, continuum.ResolutionScaleFactor ) )
1563  {
1564  fprintf( ioQQQ,"\n\n PROBLEM DISASTER At least one of the grain opacity files "
1565  "has been compiled with a different energy grid resolution factor.\n" );
1566  fprintf( ioQQQ, " Please recompile this file using the COMPILE GRAINS command "
1567  "and make sure that you use the correct SET CONTINUUM RESOLUTION factor.\n" );
1568  cdEXIT(EXIT_FAILURE);
1569  }
1570  }
1571 
1572  /* default is is to predict line intensities,
1573  * but if any continuum specified as luminosity then would override this -
1574  * following two values are correct for intensities */
1575  radius.pirsq = 0.;
1576  radius.lgPredLumin = false;
1577 
1578  /* check whether ANY luminosities are present */
1579  for( i=0; i < rfield.nShape; i++ )
1580  {
1581  if( strcmp(rfield.chRSpec[i],"4 PI") == 0 )
1582  {
1583  radius.pirsq = (realnum)(1.0992099 + 2.*xLog_radius_inner);
1584  radius.lgPredLumin = true;
1585  /* convert down to intensity */
1586  rfield.totpow[i] -= radius.pirsq;
1587 
1588  if( trace.lgTrace )
1589  {
1590  fprintf( ioQQQ,
1591  " conorm converts continuum %ld from luminosity to intensity.\n",
1592  i );
1593  }
1594  }
1595  }
1596 
1597  /* if total luminosities are present, must have specified a starting radius */
1598  if( radius.lgPredLumin && !radius.lgRadiusKnown )
1599  {
1600  fprintf(ioQQQ,"PROBLEM DISASTER conorm: - A continuum source was specified as a luminosity, but the inner radius of the cloud was not set.\n");
1601  fprintf(ioQQQ,"Please set an inner radius.\nSorry.\n");
1602  cdEXIT(EXIT_FAILURE);
1603  }
1604 
1605  /* convert ionization parameters to number of photons, called "q(h)"
1606  * at this stage q(h) and "PHI " are the same */
1607  for( i=0; i < rfield.nShape; i++ )
1608  {
1609  if( strcmp(rfield.chSpNorm[i],"IONI") == 0 )
1610  {
1611  /* the log of the ionization parameter was stored here, this converts
1612  * it to the log of the number of photons per sq cm */
1613  rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) + log10(SPEEDLIGHT);
1614  strcpy( rfield.chSpNorm[i], "Q(H)" );
1615  if( trace.lgTrace )
1616  {
1617  fprintf( ioQQQ,
1618  " conorm converts continuum %ld from ionizat par to q(h).\n",
1619  i );
1620  }
1621  }
1622  }
1623 
1624  /* convert x-ray ionization parameter xi to intensity */
1625  for( i=0; i < rfield.nShape; i++ )
1626  {
1627  if( strcmp(rfield.chSpNorm[i],"IONX") == 0 )
1628  {
1629  /* this converts it to an intensity */
1630  rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) - log10(PI4);
1631  strcpy( rfield.chSpNorm[i], "LUMI" );
1632  if( trace.lgTrace )
1633  {
1634  fprintf( ioQQQ, " conorm converts continuum%3ld from x-ray ionizat par to I.\n",
1635  i );
1636  }
1637  }
1638  }
1639 
1640  /* indicate whether we ended up with luminosity or intensity */
1641  if( trace.lgTrace )
1642  {
1643  if( radius.lgPredLumin )
1644  {
1645  fprintf( ioQQQ, " Cloudy will predict lumin into 4pi\n" );
1646  }
1647  else
1648  {
1649  fprintf( ioQQQ, " Cloudy will do surface flux for lumin\n" );
1650  }
1651  }
1652 
1653  /* if intensity per unit area is predicted then geometric
1654  * covering factor must be unity
1655  * variable can also be set elsewhere */
1656  if( !radius.lgPredLumin )
1657  {
1658  geometry.covgeo = 1.;
1659  }
1660 
1661  /* main loop over all continuum shapes to find continuum normalization
1662  * for each one */
1663  for( i=0; i < rfield.nShape; i++ )
1664  {
1665  rfield.ipSpec = i;
1666 
1667  /* check that, if laser, bounds include laser energy */
1668  if( strcmp(rfield.chSpType[rfield.ipSpec],"LASER") == 0 )
1669  {
1670  if( !( rfield.range[rfield.ipSpec][0] < rfield.slope[rfield.ipSpec] &&
1671  rfield.range[rfield.ipSpec][1] > rfield.slope[rfield.ipSpec]) )
1672  {
1673  fprintf(ioQQQ,"PROBLEM DISASTER, a continuum source is a laser at %f Ryd, but the intensity was specified over a range from %f to %f Ryd.\n",
1674  rfield.slope[rfield.ipSpec],
1675  rfield.range[rfield.ipSpec][0],
1676  rfield.range[rfield.ipSpec][1]);
1677  fprintf(ioQQQ,"Please specify the continuum flux where the laser is active.\n");
1678  cdEXIT(EXIT_FAILURE);
1679  }
1680  }
1681 
1682  if( trace.lgTrace )
1683  {
1684  long int jj;
1685  fprintf( ioQQQ, " conorm continuum number %ld is shape %s range is %.2e %.2e\n",
1686  i,
1687  rfield.chSpType[i],
1688  rfield.range[i][0],
1689  rfield.range[i][1] );
1690  fprintf( ioQQQ, "the continuum points follow\n");
1691  jj = 0;
1692  if( rfield.lgContMalloc[rfield.ipSpec] )
1693  {
1694  while( rfield.tNu[rfield.ipSpec][jj].Ryd() != 0. && jj < 100 )
1695  {
1696  fprintf( ioQQQ, "%li %e %e\n",
1697  jj ,
1698  rfield.tNu[rfield.ipSpec][jj].Ryd(),
1699  rfield.tslop[rfield.ipSpec][jj] );
1700  ++jj;
1701  }
1702  }
1703  }
1704 
1705  if( strcmp(rfield.chSpNorm[i],"RATI") == 0 )
1706  {
1707  /* option to scale relative to previous continua
1708  * this must come first since otherwise may be too late
1709  * BUT ratio cannot be the first continuum source */
1710  if( trace.lgTrace )
1711  {
1712  fprintf( ioQQQ, " conorm this is ratio to 1st con\n" );
1713  }
1714 
1715  /* check that this is not the first continuum source, we must ratio */
1716  if( i == 0 )
1717  {
1718  fprintf( ioQQQ, " I cant form a ratio if continuum is first source.\n" );
1719  cdEXIT(EXIT_FAILURE);
1720  }
1721 
1722  /* first find photon flux and Q of previous continuum */
1723  rfield.ipSpec -= 1;
1724  flx1 = ffun1(rfield.range[i][0])*rfield.spfac[rfield.ipSpec]*
1725  rfield.range[i][0];
1726 
1727  /* check that previous continua were not zero where ratio is formed */
1728  if( flx1 <= 0. )
1729  {
1730  fprintf( ioQQQ, " Previous continua were zero where ratio is desired.\n" );
1731  cdEXIT(EXIT_FAILURE);
1732  }
1733 
1734  /* return pointer to previous (correct) value, find F, Q */
1735  rfield.ipSpec += 1;
1736 
1737  /* we want a continuum totpow as powerful, flx is now desired flx */
1738  flx1 *= rfield.totpow[i];
1739 
1740  /*. find flux of this new continuum at that point */
1741  flx2 = ffun1(rfield.range[i][1])*rfield.range[i][1];
1742 
1743  /* this is ratio of desired to actual */
1744  rfield.spfac[i] = flx1/flx2;
1745  if( trace.lgTrace )
1746  {
1747  fprintf( ioQQQ, " conorm ratio will set scale fac to%10.3e at%10.2e Ryd.\n",
1748  rfield.totpow[i], rfield.range[i][0] );
1749  }
1750  }
1751 
1752  else if( strcmp(rfield.chSpNorm[i],"FLUX") == 0 )
1753  {
1754  /* specify flux density
1755  * option to use arbitrary frequency or range */
1756  f = ffun1(rfield.range[i][0]);
1757 
1758  /* make sure this is positive, could be zero if out of range of table,
1759  * or for various forms of insanity */
1760  if( f<=SMALLDOUBLE )
1761  {
1762  fprintf( ioQQQ, "\n\n PROBLEM DISASTER\n The intensity of continuum source %ld is non-positive at the energy used to normalize it (%.3e Ryd). Something is seriously wrong.\n",
1763  i , rfield.range[i][0]);
1764  /* is this a table? if so, is ffun within its bounds? */
1765  if( strcmp(rfield.chSpType[i],"INTER") == 0 )
1766  fprintf( ioQQQ, " This continuum shape given by a table of points - check that the intensity is specified at an energy within the range of that table.\n");
1767  fprintf( ioQQQ, " Also check that the numbers used to specify the shape and intensity do not under or overflow on this cpu.\n\n");
1768 
1769  cdEXIT(EXIT_FAILURE);
1770  }
1771 
1772  /* now convert to log in continuum units we shall use */
1773  f = log10(f) + log10(rfield.range[i][0]*EN1RYD/FR1RYD);
1774  f = rfield.totpow[i] - f;
1775  rfield.spfac[i] = pow(10.,f);
1776 
1777  if( trace.lgTrace )
1778  {
1779  fprintf( ioQQQ, " conorm will set log fnu to%10.3e at%10.2e Ryd. Factor=%11.4e\n",
1780  rfield.totpow[i], rfield.range[i][0], rfield.spfac[i] );
1781  }
1782  }
1783 
1784  else if( strcmp(rfield.chSpNorm[i],"Q(H)") == 0 ||
1785  strcmp(rfield.chSpNorm[i],"PHI ") == 0 )
1786  {
1787  /* some type of photon density entered */
1788  if( trace.lgTrace )
1789  {
1790  fprintf( ioQQQ, " conorm calling qintr range=%11.3e %11.3e desired val is %11.3e\n",
1791  rfield.range[i][0],
1792  rfield.range[i][1] ,
1793  rfield.totpow[i]);
1794  }
1795 
1796  /* the total number of photons over the specified range in
1797  * the arbitrary system of units that the code save the continuum shape */
1798  ASSERT( rfield.range[i][0] < rfield.range[i][1] );
1799  qentrd = qintr(&rfield.range[i][0],&rfield.range[i][1]);
1800  /* this is the log of the scale factor that must multiply the
1801  * continuum shape to get the final set of numbers */
1802  diff = rfield.totpow[i] - qentrd;
1803 
1804  /* >>chng 03 mar 13, from diff < -25 to <-35,
1805  * tripped for very low U models used for H2 simulations */
1806  /*if( diff < -25. || diff > 35. )*/
1807  if( diff < -35. || diff > 35. )
1808  {
1809  fprintf( ioQQQ, " PROBLEM DISASTER Continuum source specified is too extreme.\n" );
1810  fprintf( ioQQQ,
1811  " The integral over the continuum shape gave (log) %.3e photons, and the command requested (log) %.3e.\n" ,
1812  qentrd , rfield.totpow[i]);
1813  fprintf( ioQQQ,
1814  " The difference in the log is %.3e.\n" ,
1815  diff );
1816  if( diff>0. )
1817  {
1818  fprintf( ioQQQ, " The continuum source is too bright.\n" );
1819  }
1820  else
1821  {
1822  fprintf( ioQQQ, " The continuum source is too faint.\n" );
1823  }
1824  /* explain what happened */
1825  fprintf( ioQQQ, " The usual cause for this problem is an incorrect continuum intensity/luminosity or radius command.\n" );
1826  fprintf( ioQQQ, " There were a total of %li continuum shape commands entered - the problem is with number %li.\n",
1827  rfield.nShape , i+1 );
1828  cdEXIT(EXIT_FAILURE);
1829  }
1830 
1831  else
1832  {
1833  rfield.spfac[i] = pow(10.,diff);
1834  }
1835 
1836  if( trace.lgTrace )
1837  {
1838  fprintf( ioQQQ, " conorm finds Q over range from%11.4e-%11.4e Ryd, integral= %10.4e Factor=%11.4e\n",
1839  rfield.range[i][0],
1840  rfield.range[i][1],
1841  qentrd ,
1842  rfield.spfac[i] );
1843  }
1844  }
1845 
1846  else if( strcmp(rfield.chSpNorm[i],"LUMI") == 0 )
1847  {
1848  /* luminosity entered, special since default is TOTAL lumin */
1849  /*pintr integrates L for any continuum between two limits, used for normalization,
1850  * return units are log of ryd cm-2 s-1, last log conv to ergs */
1851  pentrd = pintr(rfield.range[i][0],rfield.range[i][1]) + log10(EN1RYD);
1852  f = rfield.totpow[i] - pentrd;
1853  rfield.spfac[i] = pow(10.,f);
1854 
1855  if( trace.lgTrace )
1856  {
1857  fprintf( ioQQQ, " conorm finds luminosity range is %10.3e to %9.3e Ryd, factor is %11.4e\n",
1858  rfield.range[i][0], rfield.range[i][1],
1859  rfield.spfac[i] );
1860  }
1861  }
1862 
1863  else
1864  {
1865  fprintf( ioQQQ, "PROBLEM DISASTER What chSpNorm label is this? =%s=\n", rfield.chSpNorm[i]);
1866  TotalInsanity();
1867  }
1868 
1869  /* spfac used to renormalize SED into flux */
1870  if( fabs(rfield.spfac[i]) <=SMALLDOUBLE )
1871  {
1872  fprintf( ioQQQ, "PROBLEM DISASTER conorm finds infinite continuum scale factor.\n" );
1873  fprintf( ioQQQ, "The continuum is too intense to compute with this cpu.\n" );
1874  fprintf( ioQQQ, "Were the intensity and luminosity commands switched?\n" );
1875  fprintf( ioQQQ, "Sorry, but I cannot go on.\n" );
1876  cdEXIT(EXIT_FAILURE);
1877  }
1878  }
1879 
1880  /* this is conversion factor for final units of line intensities or luminosities in printout,
1881  * will be intensities (==0) unless luminosity is to be printed, or flux at Earth
1882  * pirsq is the log of 4 pi r_in^2 */
1883  radius.Conv2PrtInten = radius.pirsq;
1884 
1885  /* >>chng 02 apr 25, add option for slit on aperture command */
1886  if( geometry.iEmissPower == 1 )
1887  {
1888  if( radius.lgPredLumin )
1889  {
1890  /* factor should be divided by 2 r_in (so that 2pi*r_in remains) */
1891  radius.Conv2PrtInten -= (log10(2.) + xLog_radius_inner);
1892  }
1893  else if( !radius.lgPredLumin )
1894  {
1895  /* this is an error - slit requested but radius is not known */
1896  fprintf( ioQQQ, "PROBLEM DISASTER conorm: Aperture slit specified, but not predicting luminosity.\n" );
1897  fprintf( ioQQQ, "conorm: Please specify an inner radius to determine L.\nSorry\n" );
1898  cdEXIT(EXIT_FAILURE);
1899  }
1900  }
1901  if( geometry.iEmissPower == 0 && radius.lgPredLumin )
1902  {
1903  /* leave Conv2PrtInten at zero if not predicting luminosity */
1904  radius.Conv2PrtInten = log10(2.);
1905  }
1906 
1907  /* this is option to give final absolute results as flux observed at Earth */
1908  if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth )
1909  {
1910  /* this implements the conversion from Q_alpha -> F_alpha
1911  * described in the section on the APERTURE command in Hazy 1 */
1912  if( geometry.iEmissPower == 0 )
1913  radius.Conv2PrtInten += log10(double(geometry.size)/SQAS_SKY);
1914  else if( geometry.iEmissPower == 1 )
1915  radius.Conv2PrtInten += log10(double(geometry.size)/(PI4*AS1RAD*radius.distance));
1916  else if( geometry.iEmissPower == 2 )
1917  radius.Conv2PrtInten -= log10( PI4*pow2(radius.distance) );
1918  else
1919  TotalInsanity();
1920  }
1921 
1922  /* normally lines are into 4pi, this is option to do per sr or arcsec^2 */
1923  if( prt.lgSurfaceBrightness )
1924  {
1925  if( radius.pirsq != 0. )
1926  {
1927  /* make sure we are predicting line intensities, not luminosity */
1928  fprintf( ioQQQ, " PROBLEM DISASTER Sorry, but both luminosity and surface brightness have been requested for lines.\n" );
1929  fprintf( ioQQQ, " the PRINT LINE SURFACE BRIGHTNESS command can only be used when lines are predicted per unit cloud area.\n" );
1930  cdEXIT(EXIT_FAILURE);
1931  }
1932  if( prt.lgSurfaceBrightness_SR )
1933  {
1934  /* we want final units to be per sr */
1935  radius.Conv2PrtInten -= log10( PI4 );
1936  }
1937  else
1938  {
1939  /* we want final units to be per square arcsec */
1940  radius.Conv2PrtInten -= log10(SQAS_SKY);
1941  }
1942  }
1943  return;
1944 }
1945 
1946 /*qintr integrates Q for any continuum between two limits, used for normalization */
1947 STATIC double qintr(double *qenlo,
1948  double *qenhi)
1949 {
1950  long int i,
1951  ipHi,
1952  ipLo,
1953  j;
1954  double qintr_v,
1955  sum,
1956  wanu;
1957 
1958  DEBUG_ENTRY( "qintr()" );
1959 
1960  /* returns LOG of number of photons over energy interval */
1961 
1962  /* this is copy of logic that occurs three times across code */
1963  ASSERT(*qenhi > *qenlo);
1964  ipLo = ipoint(*qenlo);
1965  ipHi = ipoint(*qenhi);
1966  /* this is actual sum of photons within band */
1967  sum = 0.;
1968  for( i=ipLo-1; i < (ipHi - 1); i++ )
1969  {
1970  /*sum += ffun1(rfield.anu[i])*rfield.widflx[i];*/
1971  for( j=0; j < 4; j++ )
1972  {
1973  wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
1974  /* >>chng 02 jul 16, add test on continuum limits -
1975  * this was exceeded when resolution set very large */
1976  wanu = MAX2( wanu , rfield.emm );
1977  wanu = MIN2( wanu , rfield.egamry );
1978  sum += fweigh[j]*ffun1(wanu)*rfield.widflx[i];
1979  }
1980  }
1981 
1982  if( sum <= 0. )
1983  {
1984  fprintf( ioQQQ, " PROBLEM DISASTER Photon number sum in QINTR is %.3e\n",
1985  sum );
1986  fprintf( ioQQQ, " This source has no ionizing radiation, and the number of ionizing photons was specified.\n" );
1987  fprintf( ioQQQ, " This was continuum source number%3ld\n",
1988  rfield.ipSpec );
1989  fprintf( ioQQQ, " Sorry, but I cannot go on. ANU and FLUX arrays follow. Enjoy.\n" );
1990  fprintf( ioQQQ, "\n\n This error is also caused by an old table read file whose energy mesh does not agree with the code.\n" );
1991  for( i=0; i < rfield.nupper; i++ )
1992  {
1993  fprintf( ioQQQ, "%.2e\t%.2e\n",
1994  rfield.anu[i],
1995  rfield.flux[0][i] );
1996  }
1997  cdEXIT(EXIT_FAILURE);
1998  }
1999  else
2000  {
2001  qintr_v = log10(sum);
2002  }
2003  return qintr_v;
2004 }
2005 
2006 /*pintr integrates L for any continuum between two limits, used for normalization,
2007  * return units are log of ryd cm-2 s-1 */
2008 STATIC double pintr(double penlo,
2009  double penhi)
2010 {
2011  long int i,
2012  j;
2013  double fsum,
2014  pintr_v,
2015  sum,
2016  wanu,
2017  wfun;
2018  long int ip1 , ip2;
2019 
2020  DEBUG_ENTRY( "pintr()" );
2021 
2022  /* computes log of luminosity in radiation over some integral
2023  * answer is in Ryd per sec */
2024 
2025  sum = 0.;
2026  /* >>chng 02 oct 27, do not call qg32, do same type sum as
2027  * final integration */
2028  /* laser is special since delta function, this is center of laser */
2029  /* >>chng 01 jul 01, was +-21 cells, change to call to ipoint */
2030  ip1 = ipoint( penlo );
2031 
2032  ip2 = ipoint( penhi );
2033 
2034  for( i=ip1-1; i < ip2-1; i++ )
2035  {
2036  fsum = 0.;
2037  for( j=0; j < 4; j++ )
2038  {
2039  wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
2040  /*++iiii;fprintf(ioQQQ,"DEBUG iii %li %e \n",iiii, wanu );*/
2041  wfun = fweigh[j]*ffun1(wanu)*wanu;
2042  fsum += wfun;
2043  }
2044  sum += fsum*rfield.widflx[i];
2045  }
2046 
2047  if( sum > 0. )
2048  {
2049  pintr_v = log10(sum);
2050  }
2051  else
2052  {
2053  pintr_v = -38.;
2054  }
2055 
2056  return pintr_v;
2057 }
MIN2
#define MIN2
Definition: cddefines.h:765
t_phycon::TEnerDen
double TEnerDen
Definition: phycon.h:98
t_rfield::fine_anu
realnum * fine_anu
Definition: rfield.h:412
t_rfield::csigh
realnum * csigh
Definition: rfield.h:288
t_prt::q
realnum q
Definition: prt.h:228
t_rfield::widflx
realnum * widflx
Definition: rfield.h:65
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:616
t_rfield::lgBeamed
bool lgBeamed[LIMSPC]
Definition: rfield.h:310
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_rfield::flux_isotropic
realnum * flux_isotropic
Definition: rfield.h:89
t_isoCTRL::lgContinuumLoweringEnabled
bool lgContinuumLoweringEnabled[NISO]
Definition: iso.h:358
t_timesc::TimeErode
realnum TimeErode
Definition: timesc.h:48
t_geometry::size
realnum size
Definition: geometry.h:67
t_iso_sp::st
qList st
Definition: iso.h:453
ffun
double ffun(double anu)
Definition: cont_ffun.cpp:17
t_continuum::sv4861
realnum sv4861
Definition: continuum.h:101
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
iso_continuum_lower
void iso_continuum_lower(long ipISO, long nelem)
Definition: iso_continuum_lower.cpp:12
t_rfield::qtot
realnum qtot
Definition: rfield.h:356
t_rfield::ipG0_DB96_hi
long int ipG0_DB96_hi
Definition: rfield.h:267
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
t_rfield::qbal
realnum qbal
Definition: rfield.h:356
t_rfield::flux_beam_const_save
realnum * flux_beam_const_save
Definition: rfield.h:210
opac
t_opac opac
Definition: opacity.cpp:5
t_trace::lgPtrace
bool lgPtrace
Definition: trace.h:118
t_rfield::flux
realnum ** flux
Definition: rfield.h:86
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
t_rfield::occ1nu
realnum occ1nu
Definition: rfield.h:472
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:178
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
AS1RAD
const double AS1RAD
Definition: physconst.h:129
NISO
const int NISO
Definition: cddefines.h:265
t_rfield::totpow
double totpow[LIMSPC]
Definition: rfield.h:300
t_rfield::chRSpec
char chRSpec[LIMSPC][5]
Definition: rfield.h:351
t_rfield::flux_total_incident
realnum ** flux_total_incident
Definition: rfield.h:209
t_continuum::sv1216
realnum sv1216
Definition: continuum.h:101
TorF
char TorF(bool l)
Definition: cddefines.h:714
t_rfield::anu3
realnum * anu3
Definition: rfield.h:77
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:316
MAX2
#define MAX2
Definition: cddefines.h:786
t_rfield::lgOcc1Hi
bool lgOcc1Hi
Definition: rfield.h:480
t_geometry::iEmissPower
long int iEmissPower
Definition: geometry.h:61
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:182
t_rfield::FluxFaint
realnum FluxFaint
Definition: rfield.h:55
t_trace::lgComBug
bool lgComBug
Definition: trace.h:37
hextra
t_hextra hextra
Definition: hextra.cpp:5
t_hextra::lgNeutrnHeatOn
bool lgNeutrnHeatOn
Definition: hextra.h:71
t_radius::distance
double distance
Definition: radius.h:65
t_rfield::qhe
realnum qhe
Definition: rfield.h:356
t_rfield::emm
realnum emm
Definition: rfield.h:49
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_oxy::i2d
long int i2d
Definition: oxy.h:30
t_rfield::lgMeshSetUp
bool lgMeshSetUp
Definition: rfield.h:131
t_continuum::fluxv
realnum fluxv
Definition: continuum.h:107
t_hextra::effneu
realnum effneu
Definition: hextra.h:75
dense
t_dense dense
Definition: dense.cpp:24
t_dense::EdenSet
realnum EdenSet
Definition: dense.h:203
t_rfield::RSFCheck
double RSFCheck[LIMSPC]
Definition: rfield.h:339
t_ADfA::coll_ion_wrapper
double coll_ion_wrapper(long int z, long int n, double t)
Definition: atmdat_adfa.cpp:814
t_rfield::tNu
vector< Energy > tNu[LIMSPC]
Definition: rfield.h:330
ffun1
double ffun1(double xnu)
Definition: cont_ffun.cpp:112
t_continuum::cn4861
realnum cn4861
Definition: continuum.h:101
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
t_ionbal::CosRayIonRate
double CosRayIonRate
Definition: ionbal.h:123
t_opac::ipo3exc
long int ipo3exc[3]
Definition: opacity.h:210
noexec
t_noexec noexec
Definition: noexec.cpp:5
t_hydro::lgHInducImp
bool lgHInducImp
Definition: hydrogenic.h:95
t_rfield::ipEnerGammaRay
long int ipEnerGammaRay
Definition: rfield.h:466
t_geometry::covgeo
realnum covgeo
Definition: geometry.h:35
t_prt::pradio
realnum pradio
Definition: prt.h:228
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_rfield::spfac
double spfac[LIMSPC]
Definition: rfield.h:300
t_rfield::tslop
vector< realnum > tslop[LIMSPC]
Definition: rfield.h:331
t_called::lgTalk
bool lgTalk
Definition: called.h:12
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
TempChange
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:51
SPEEDLIGHT
const double SPEEDLIGHT
Definition: physconst.h:100
t_dense::density_low_limit
double density_low_limit
Definition: dense.h:197
t_noexec::lgNoExec
bool lgNoExec
Definition: noexec.h:14
t_rfield::ExtinguishEnergyPowerLow
realnum ExtinguishEnergyPowerLow
Definition: rfield.h:99
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
qintr
STATIC double qintr(double *qenlo, double *qenhi)
Definition: cont_setintensity.cpp:1947
t_hextra::frcneu
realnum frcneu
Definition: hextra.h:73
t_rfield::lgTimeVary
bool lgTimeVary[LIMSPC]
Definition: rfield.h:306
t_prt::fx1ryd
realnum fx1ryd
Definition: prt.h:228
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:179
t_rfield::range
double range[LIMSPC][2]
Definition: rfield.h:347
t_rfield::lgHabing
bool lgHabing
Definition: rfield.h:376
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
PI4
const double PI4
Definition: physconst.h:35
t_rfield::ExtinguishConvertColDen2OptDepth
realnum ExtinguishConvertColDen2OptDepth
Definition: rfield.h:99
t_rfield::egamry
realnum egamry
Definition: rfield.h:52
t_rfield::nupper
long int nupper
Definition: rfield.h:46
SQAS_SKY
const double SQAS_SKY
Definition: physconst.h:135
t_secondaries::SetCsupra
realnum SetCsupra
Definition: secondaries.h:33
t_thermal::lgTemperatureConstant
bool lgTemperatureConstant
Definition: thermal.h:32
trace
t_trace trace
Definition: trace.cpp:5
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:816
RT_OTS_Zero
void RT_OTS_Zero(void)
Definition: rt_ots.cpp:599
t_continuum::nrange
long int nrange
Definition: continuum.h:75
t_dense::SetIoniz
realnum SetIoniz[LIMELM][LIMELM+1]
Definition: dense.h:154
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
t_radius::rinner
double rinner
Definition: radius.h:22
t_rfield::ExtinguishLeakage
realnum ExtinguishLeakage
Definition: rfield.h:99
t_rfield::ipG0_DB96_lo
long int ipG0_DB96_lo
Definition: rfield.h:267
geometry
t_geometry geometry
Definition: geometry.cpp:5
FR1RYD
const double FR1RYD
Definition: physconst.h:195
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:16
t_trace::lgConBug
bool lgConBug
Definition: trace.h:100
t_rfield::slope
double slope[LIMSPC]
Definition: rfield.h:300
aweigh
static const double aweigh[4]
Definition: cont_setintensity.cpp:40
t_prt::pbal
realnum pbal
Definition: prt.h:228
POW2
#define POW2
Definition: cddefines.h:933
t_rfield::tbrmax
realnum tbrmax
Definition: rfield.h:472
t_rfield::chLineLabel
char ** chLineLabel
Definition: rfield.h:220
ipH1s
const int ipH1s
Definition: iso.h:27
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:216
t_rfield::anulog
realnum * anulog
Definition: rfield.h:77
STATIC
#define STATIC
Definition: cddefines.h:101
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
t_rfield::ipPlasma
long int ipPlasma
Definition: rfield.h:453
extin
STATIC void extin(realnum *ex1ryd)
Definition: cont_setintensity.cpp:1450
continuum
t_continuum continuum
Definition: continuum.cpp:5
t_opac::ipCKshell
long int ipCKshell
Definition: opacity.h:291
t_prt::lgSurfaceBrightness
bool lgSurfaceBrightness
Definition: prt.h:138
t_rfield::qhtot
realnum qhtot
Definition: rfield.h:356
t_dense::EdenHCorr_f
realnum EdenHCorr_f
Definition: dense.h:218
t_thermal::ConstTemp
realnum ConstTemp
Definition: thermal.h:44
rfield
t_rfield rfield
Definition: rfield.cpp:8
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
t_rfield::anu2
realnum * anu2
Definition: rfield.h:77
t_rfield::flux_time_beam_save
realnum * flux_time_beam_save
Definition: rfield.h:210
pintr
STATIC double pintr(double penlo, double penhi)
Definition: cont_setintensity.cpp:2008
t_continuum::ResolutionScaleFactor
double ResolutionScaleFactor
Definition: continuum.h:90
t_rfield::convoc
realnum * convoc
Definition: rfield.h:134
realnum
float realnum
Definition: cddefines.h:107
t_rfield::uh
realnum uh
Definition: rfield.h:364
t_continuum::filbnd
realnum * filbnd
Definition: continuum.h:69
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:144
t_prt::GammaLumin
realnum GammaLumin
Definition: prt.h:237
BIGFLOAT
const realnum BIGFLOAT
Definition: cpu.h:176
t_rfield::Illumination
Illuminate::IlluminationType Illumination[LIMSPC]
Definition: rfield.h:316
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:258
t_rfield::uheii
realnum uheii
Definition: rfield.h:367
max
long max(int a, long b)
Definition: cddefines.h:779
t_prt::lgSurfaceBrightness_SR
bool lgSurfaceBrightness_SR
Definition: prt.h:138
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:438
t_radius::lgRadiusKnown
bool lgRadiusKnown
Definition: radius.h:116
t_rfield::ContBoltz
double * ContBoltz
Definition: rfield.h:145
t_rfield::OccNumbIncidCont
realnum * OccNumbIncidCont
Definition: rfield.h:138
EN1RYD
const double EN1RYD
Definition: physconst.h:179
t_rfield::chSpNorm
char chSpNorm[LIMSPC][5]
Definition: rfield.h:351
t_hextra::cryden
realnum cryden
Definition: hextra.h:14
radius
t_radius radius
Definition: radius.cpp:5
timesc
t_timesc timesc
Definition: timesc.cpp:5
prt
t_prt prt
Definition: prt.cpp:10
t_radius::pirsq
realnum pirsq
Definition: radius.h:143
t_rfield::ConEmitOut
realnum ** ConEmitOut
Definition: rfield.h:161
t_hextra::totneu
realnum totneu
Definition: hextra.h:69
t_radius::lgPredLumin
bool lgPredLumin
Definition: radius.h:139
HNU3C2
const double HNU3C2
Definition: physconst.h:198
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
t_radius::Conv2PrtInten
double Conv2PrtInten
Definition: radius.h:147
ipH2p
const int ipH2p
Definition: iso.h:29
t_continuum::TotalLumin
double TotalLumin
Definition: continuum.h:97
t_dense::IonLow
long int IonLow[LIMELM+1]
Definition: dense.h:119
t_prt::xpow
realnum xpow
Definition: prt.h:228
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:582
t_rfield::anu
double * anu
Definition: rfield.h:58
t_continuum::fbeta
realnum fbeta
Definition: continuum.h:107
t_continuum::lgCon0
bool lgCon0
Definition: continuum.h:93
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
LIMELM
const int LIMELM
Definition: cddefines.h:262
pow2
T pow2(T a)
Definition: cddefines.h:935
t_rfield::nfine
long nfine
Definition: rfield.h:402
IncidentContinuumHere
void IncidentContinuumHere()
Definition: cont_setintensity.cpp:68
sumcon
STATIC void sumcon(long int il, long int ih, realnum *q, realnum *p, realnum *panu)
Definition: cont_setintensity.cpp:1269
t_rfield::flux_isotropic_save
realnum * flux_isotropic_save
Definition: rfield.h:210
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:688
ELECTRON_MASS
const double ELECTRON_MASS
Definition: physconst.h:91
ELEM_CHARGE_ESU
const double ELEM_CHARGE_ESU
Definition: physconst.h:147
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:310
t_prt::qgam
realnum qgam
Definition: prt.h:228
t_prt::powion
realnum powion
Definition: prt.h:228
t_rfield::ExtinguishColumnDensity
realnum ExtinguishColumnDensity
Definition: rfield.h:99
t_rfield::qheii
realnum qheii
Definition: rfield.h:356
t_rfield::csigc
realnum * csigc
Definition: rfield.h:288
ptrcer
STATIC void ptrcer()
Definition: cont_setintensity.cpp:1303
t_rfield::lgHionRad
bool lgHionRad
Definition: rfield.h:469
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
fweigh
static const double fweigh[4]
Definition: cont_setintensity.cpp:41
t_dense::eden
double eden
Definition: dense.h:190
t_continuum::totlsv
double totlsv
Definition: continuum.h:97
t_rfield::AnuOrg
double * AnuOrg
Definition: rfield.h:62
t_dense::EdenTrue
double EdenTrue
Definition: dense.h:221
oxy
t_oxy oxy
Definition: oxy.cpp:5
t_prt::lgPrintFluxEarth
bool lgPrintFluxEarth
Definition: prt.h:134
t_ionbal::lgCompRecoil
bool lgCompRecoil
Definition: ionbal.h:149
t_rfield::ExtinguishFactor
realnum * ExtinguishFactor
Definition: rfield.h:98
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
t_rfield::OpticalDepthScaleFactor
realnum OpticalDepthScaleFactor[LIMSPC]
Definition: rfield.h:314
t_rfield::ipPlasmax
long int ipPlasmax
Definition: rfield.h:453
t_rfield::nfine_malloc
long nfine_malloc
Definition: rfield.h:404
conorm
STATIC void conorm()
Definition: cont_setintensity.cpp:1497
t_phycon::telogn
double telogn[7]
Definition: phycon.h:76
ipCARBON
const int ipCARBON
Definition: cddefines.h:314
t_rfield::nShape
long int nShape
Definition: rfield.h:322
t_rfield::ipSpec
long int ipSpec
Definition: rfield.h:322
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
t_rfield::occmnu
realnum occmnu
Definition: rfield.h:472
EdenChange
void EdenChange(double EdenNew)
Definition: eden_change.cpp:12
t_rfield::ConEmitReflec
realnum ** ConEmitReflec
Definition: rfield.h:155
t_rfield::flux_beam_time
realnum * flux_beam_time
Definition: rfield.h:92
gv
GrainVar gv
Definition: grainvar.cpp:5
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
t_dense::EdenExtra
realnum EdenExtra
Definition: dense.h:206
POW3
#define POW3
Definition: cddefines.h:940
ShowMe
void ShowMe(void)
Definition: service.cpp:181
t_rfield::flux_beam_const
realnum * flux_beam_const
Definition: rfield.h:92
t_phycon::te
double te
Definition: phycon.h:11
t_oxy::i2p
long int i2p
Definition: oxy.h:30
t_rfield::plsfrqmax
realnum plsfrqmax
Definition: rfield.h:447
t_rfield::qrad
realnum qrad
Definition: rfield.h:356
t_prt::qx
realnum qx
Definition: prt.h:228
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:309
t_rfield::plsfrq
realnum plsfrq
Definition: rfield.h:447
t_rfield::nflux
long int nflux
Definition: rfield.h:43
ipLITHIUM
const int ipLITHIUM
Definition: cddefines.h:311
t_rfield::lgPlasNu
bool lgPlasNu
Definition: rfield.h:445
t_rfield::ConRefIncid
realnum ** ConRefIncid
Definition: rfield.h:167
t_rfield::tbr4nu
realnum tbr4nu
Definition: rfield.h:472
t_dense::lgSetIoniz
bool lgSetIoniz[LIMELM]
Definition: dense.h:149
t_rfield::chSpType
char chSpType[LIMSPC][6]
Definition: rfield.h:351
called
t_called called
Definition: called.cpp:5
t_rfield::chContLabel
char ** chContLabel
Definition: rfield.h:223
HPLANCK
const double HPLANCK
Definition: physconst.h:103
t_Heavy::ipHeavy
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
t_continuum::cn1216
realnum cn1216
Definition: continuum.h:101
t_rfield::lgContMalloc
bool lgContMalloc[LIMSPC]
Definition: rfield.h:343
PI
const double PI
Definition: physconst.h:29
t_rfield::ExtinguishLowEnergyLimit
realnum ExtinguishLowEnergyLimit
Definition: rfield.h:99
t_ionbal::ipCompRecoil
long int ** ipCompRecoil
Definition: ionbal.h:155
TE1RYD
const double TE1RYD
Definition: physconst.h:183
ioStdin
FILE * ioStdin
Definition: cddefines.cpp:8
t_rfield::occmax
realnum occmax
Definition: rfield.h:472
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:142
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
ContSetIntensity
void ContSetIntensity()
Definition: cont_setintensity.cpp:100
t_prt::ipeak
long int ipeak
Definition: prt.h:236
t_rfield::tbrmnu
realnum tbrmnu
Definition: rfield.h:472
Generated on Sat Oct 19 2013 10:43:09 for cloudy by   doxygen 1.8.5