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

Generated for cloudy by doxygen 1.8.4