cloudy trunk
Loading...
Searching...
No Matches
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 */
40static const double aweigh[4]={-0.4305682,-0.1699905, 0.1699905, 0.4305682};
41static const double fweigh[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274};
42
43/*conorm normalize continuum to proper intensity */
44STATIC void conorm();
45
46/*pintr integrates L for any continuum between two limits, used for normalization */
47STATIC double pintr(double penlo,
48 double penhi);
49
50/*qintr integrates Q for any continuum between two limits, used for normalization */
51STATIC double qintr(double *qenlo,
52 double *qenhi);
53
54
55/*sumcon sums L and Q for net incident continuum */
56STATIC 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 */
63STATIC void extin(realnum *ex1ryd);
64
65/*ptrcer show continuum pointers in real time following drive pointers command */
66STATIC void ptrcer();
67
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 */
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 }
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();
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" );
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 }
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" );
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();
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] = dense.gas_phase[ipHYDROGEN]*RatioIoniz;
887 dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
888 dense.xIonDense[ipHYDROGEN][1];
889 ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
890 dense.xIonDense[ipHYDROGEN][0]<=dense.gas_phase[ipHYDROGEN]);
891 ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
892 dense.xIonDense[ipHYDROGEN][1]<dense.gas_phase[ipHYDROGEN]);
893 //fprintf(ioQQQ,"DEBUG br 1 H0 %.2e\n", dense.xIonDense[ipHYDROGEN][0]);
894 }
895 else if( RatioIoniz>1e3 )
896 {
897 /* very high ionization solution */
898 dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN]/RatioIoniz;
899 dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN] -
900 dense.xIonDense[ipHYDROGEN][0];
901 ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
902 dense.xIonDense[ipHYDROGEN][0]<dense.gas_phase[ipHYDROGEN]);
903 ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
904 dense.xIonDense[ipHYDROGEN][1]<=dense.gas_phase[ipHYDROGEN]);
905 //fprintf(ioQQQ,"DEBUG br 2 H0 %.2e rat %.2e\n", dense.xIonDense[ipHYDROGEN][0],
906 // dense.gas_phase[ipHYDROGEN]/RatioIoniz);
907 }
908 else
909 {
910 /* intermediate ionization - solve quadratic */
911 double alpha = HCaseBRecCoeff + CollIoniz ;
912 double beta = HCaseBRecCoeff*EdenExtraLocal + OtherIonization +
913 (EdenExtraLocal - dense.gas_phase[ipHYDROGEN])*CollIoniz;
914 double gamma = -dense.gas_phase[ipHYDROGEN]*(OtherIonization+EdenExtraLocal*CollIoniz);
915
916 double discriminant = POW2(beta) - 4.*alpha*gamma;
917 if( discriminant <0 )
918 {
919 /* oops */
920 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found negative discriminant.\n");
922 }
923
924 dense.xIonDense[ipHYDROGEN][1] = (realnum)((-beta + sqrt(discriminant))/(2.*alpha));
925 if( dense.xIonDense[ipHYDROGEN][1]> dense.gas_phase[ipHYDROGEN] )
926 {
927 /* oops */
928 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H+)>n(H).\n");
930 }
931 dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
932 dense.xIonDense[ipHYDROGEN][1];
933 if( dense.xIonDense[ipHYDROGEN][0]<= 0 )
934 {
935 /* oops */
936 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H0)<0.\n");
938 }
939 //fprintf(ioQQQ,"DEBUG br 3 H0 %.2e\n", dense.xIonDense[ipHYDROGEN][0]);
940 }
941
942
943 if( dense.xIonDense[ipHYDROGEN][1] > 1e-30 )
944 {
945 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = dense.xIonDense[ipHYDROGEN][0];
946 }
947 else
948 {
949 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = 0.;
950 }
951
952 /* now save estimates of whether induced recombination is going
953 * to be important -this is a code pacesetter since GammaBn is slower
954 * than GammaK */
955 hydro.lgHInducImp = false;
956 for( i=ipH1s; i < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max; i++ )
957 {
958 if( rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].ipIsoLevNIonCon-1] > 0.01 )
959 hydro.lgHInducImp = true;
960 }
961
962 /*******************************************************************
963 * *
964 * estimate helium's level of ionization *
965 * *
966 *******************************************************************/
967
968 /* only if helium is turned on */
969 if( dense.lgElmtOn[ipHELIUM] )
970 {
971 /* next estimate level of helium singly ionized */
972 xIoniz = (realnum)t_ADfA::Inst().coll_ion_wrapper(1,0,phycon.te);
973 /* >>chng 05 jan 05, add cosmic rays */
974 xIoniz = (realnum)(xIoniz*dense.eden + rfield.qhe*1e-18 + secondaries.csupra[ipHELIUM][1] );
975 double RecTot = HCaseBRecCoeff*dense.eden;
976 RatioIonizToRecomb = xIoniz/RecTot;
977
978 /* now estimate level of helium doubly ionized */
979 xIoniz = (realnum)t_ADfA::Inst().coll_ion_wrapper(1,1,phycon.te);
980 /* >>chng 05 jan 05, add cosmic rays */
981 xIoniz = (realnum)(xIoniz*dense.eden + rfield.qheii*1e-18 + ionbal.CosRayIonRate );
982
983 /* rough charge dependence */
984 RecTot *= 4.;
985 r3ov2 = xIoniz/RecTot;
986
987 /* now set level of helium ionization */
988 if( RatioIonizToRecomb > 0. )
989 {
990 double r1 = dense.gas_phase[ipHELIUM]/(1./RatioIonizToRecomb + 1. + r3ov2);
991 dense.xIonDense[ipHELIUM][1] = (realnum)(r1);
992 dense.xIonDense[ipHELIUM][0] = (realnum)(r1/RatioIonizToRecomb);
993 dense.xIonDense[ipHELIUM][2] = (realnum)(r1*r3ov2);
994 }
995 else
996 {
997 /* no He ionizing radiation */
998 dense.xIonDense[ipHELIUM][1] = dense.xIonDense[ipHELIUM][2] = 0.;
999 dense.xIonDense[ipHELIUM][0] = dense.gas_phase[ipHELIUM];
1000 }
1001
1002 if( dense.xIonDense[ipHELIUM][2] > 1e-30 )
1003 {
1004 iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = dense.xIonDense[ipHELIUM][1];
1005 }
1006 else
1007 {
1008 iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = 0.;
1009 }
1010 }
1011 else
1012 {
1013 /* case where helium is turned off */
1014 dense.xIonDense[ipHELIUM][1] = 0.;
1015 dense.xIonDense[ipHELIUM][0] = 0.;
1016 dense.xIonDense[ipHELIUM][2] = 0.;
1017 iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = 0.;
1018 }
1019
1020 /* update electron density */
1021 newEden = dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal + dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2];
1022
1023 loopCount++;
1024 }
1025 /* repeat above until guessed and calculated eden agree to at least 0.1%. */
1026 while( loopCount < 10 && fabs(newEden/dense.eden - 1.) > 0.001 );
1027
1028 if( dense.xIonDense[ipHYDROGEN][0]<0.)
1029 TotalInsanity();
1030 else if( dense.xIonDense[ipHYDROGEN][0] == 0. )
1031 {
1032 fprintf(ioQQQ,"PROBLEM the derived atomic hydrogen density is zero.\n");
1033 if( dense.gas_phase[ipHYDROGEN]<1e-5 && rfield.uh > 1e10)
1034 {
1035 fprintf(ioQQQ,"This is almost certainly due to floating point "
1036 "limits on this computer.\nThe ionization parameter is very large,"
1037 " the density is very small,\nand the H^0 density cannot be"
1038 " stored in a double.\n");
1039 //cdEXIT( EXIT_FAILURE );
1040 }
1041 }
1042 //ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);
1043
1044 /* update electron density */
1045 EdenChange( newEden );
1046
1047 if( dense.eden <= SMALLFLOAT )
1048 {
1049 /* no electrons! */
1050 fprintf(ioQQQ,"\n PROBLEM DISASTER - this simulation has no source"
1051 " of ionization. The electron density is zero. Consider "
1052 "adding a source of ionization such as cosmic rays.\n\n");
1054 }
1055
1056 /* fix range of stages of ionization */
1057 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
1058 {
1059 if( dense.lgElmtOn[nelem] )
1060 {
1061 // IonHigh[n] is the highest stage of ionization present
1062 // the IonHigh array index is on the C scale, so [0] is hydrogen
1063 // the value is also on the C scale, so element [nelem] can range
1064 // from 0 to nelem+1
1065 dense.IonHigh[nelem] = nelem + 1;
1066
1067 dense.IonLow[nelem] = 0;
1068 // for very intense radiation fields very heavy elements, N>Fe, will fail
1069 // in ion_solver due to ill conditioned matrix. all populations are in
1070 // fully stripped state. Start will fully stripped ion distribution in this case.
1071 if( rfield.uh > 1e15 )
1072 {
1073 //trim down highest stage to be within incident radiation field
1074 while ( rfield.anu[Heavy.ipHeavy[nelem][dense.IonHigh[nelem]-1]] >
1075 rfield.anu[rfield.nflux] && dense.IonHigh[nelem]>1 )
1076 --dense.IonHigh[nelem];
1077
1078 dense.IonLow[nelem] = max( 0 , dense.IonHigh[nelem]-1 );
1079 }
1080
1081 /* >>chng 04 jan 13, add this test, caught by Orly Gnat */
1082 /* check on actual zero abundances of lower stages - this will only
1083 * happen when ionization is set with element ionization command */
1084 if( dense.lgSetIoniz[nelem] )
1085 {
1086 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] < dense.density_low_limit )
1087 ++dense.IonLow[nelem];
1088 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] < dense.density_low_limit )
1089 --dense.IonHigh[nelem];
1090 }
1091
1092 // make low-stage populations zero
1093 for( ion=0; ion<dense.IonLow[nelem]; ++ion )
1094 {
1095 dense.xIonDense[nelem][dense.IonLow[nelem]] += dense.xIonDense[nelem][ion];
1096 dense.xIonDense[nelem][ion] = 0.;
1097 }
1098 for( ion=nelem+1; ion>dense.IonHigh[nelem]; --ion )
1099 {
1100 dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion];
1101 dense.xIonDense[nelem][ion] = 0.;
1102 }
1103 }
1104 else
1105 {
1106 /* this element is turned off, make stages impossible */
1107 dense.IonLow[nelem] = -1;
1108 dense.IonHigh[nelem] = -1;
1109 }
1110 }
1111
1112 // make first estimate of iso continuum lowering
1113 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1114 {
1115 for( long nelem=ipISO; nelem<LIMELM; ++nelem)
1116 {
1117 if( nelem < 2 || dense.lgElmtOn[nelem] )
1118 {
1119 if( iso_ctrl.lgContinuumLoweringEnabled[ipISO] )
1120 iso_continuum_lower( ipISO, nelem );
1121 }
1122 }
1123 }
1124
1125 /* estimate electrons from heavies, assuming each at least
1126 * 1 times ionized */
1127 EdenHeav = 0.;
1128 realnum atomFrac = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
1129 realnum firstIonFrac = dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN];
1130 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
1131 {
1132 if( dense.lgElmtOn[nelem] )
1133 {
1134 if( dense.IonLow[nelem] == 0 && dense.IonHigh[nelem] >=1)
1135 {
1136 realnum low2dens = dense.xIonDense[nelem][0]+dense.xIonDense[nelem][1];
1137 dense.xIonDense[nelem][0] = low2dens * atomFrac;
1138 dense.xIonDense[nelem][1] = low2dens * firstIonFrac;
1139 }
1140 for (long ion=1;ion<dense.IonHigh[nelem];++ion)
1141 {
1142 EdenHeav += ion*dense.xIonDense[nelem][ion];
1143 }
1144 }
1145 }
1146
1147 /* modify estimate of electron density */
1148 dense.eden += EdenHeav;
1149
1150 /* >>chng 05 jan 05, insure positive eden */
1151 EdenChange( MAX2( SMALLFLOAT , dense.eden ) );
1152
1153 if( dense.EdenSet > 0. )
1154 {
1155 EdenChange( dense.EdenSet );
1156 }
1157
1158 dense.EdenHCorr = dense.eden;
1159 dense.EdenHCorr_f = (realnum)dense.EdenHCorr;
1160
1161 if( dense.eden < 0. )
1162 {
1163 fprintf( ioQQQ, " PROBLEM DISASTER Negative electron density results in ContSetIntensity.\n" );
1164 fprintf( ioQQQ, "%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n",
1165 dense.eden, dense.xIonDense[ipHYDROGEN][1], dense.xIonDense[ipHELIUM][1],
1166 dense.xIonDense[ipHELIUM][2], dense.gas_phase[ipCARBON], dense.EdenExtra );
1167 TotalInsanity();
1168 }
1169
1170 if( dense.EdenSet > 0. )
1171 {
1172 dense.EdenTrue = dense.EdenSet;
1173 }
1174 else
1175 dense.EdenTrue = dense.eden;
1176
1177 if( trace.lgTrace )
1178 {
1179 fprintf( ioQQQ,
1180 " ContSetIntensity sets initial EDEN to %.4e, contributors H+=%.2e He+, ++= %.2e %.2e Heav %.2e extra %.2e\n",
1181 dense.eden ,
1182 dense.xIonDense[ipHYDROGEN][1],
1183 dense.xIonDense[ipHELIUM][1],
1184 2.*dense.xIonDense[ipHELIUM][2],
1185 EdenHeav,
1186 dense.EdenExtra);
1187 }
1188
1189 // photon occupation number at 1 Ryd - used for printout
1190 occ1 = (realnum)(prt.fx1ryd/HNU3C2/PI4/FR1RYD);
1191 if( occ1 > 1. )
1192 rfield.lgOcc1Hi = true;
1193 else
1194 rfield.lgOcc1Hi = false;
1195
1196 if( trace.lgTrace && trace.lgConBug )
1197 {
1198 fprintf(ioQQQ,"\ntrace continuum print of %li incident spectral "
1199 "components\n", rfield.nShape);
1200 fprintf(ioQQQ," # type Illum Beam? 1/cos TimeVary?\n");
1201 for( rfield.ipSpec=0; rfield.ipSpec<rfield.nShape; ++rfield.ipSpec )
1202 {
1203 fprintf(ioQQQ,"%3li %6s %5i %c %.3f %c\n",
1204 rfield.ipSpec ,
1205 rfield.chSpType[rfield.ipSpec],
1206 rfield.Illumination[rfield.ipSpec],
1207 TorF(rfield.lgBeamed[rfield.ipSpec]),
1208 rfield.OpticalDepthScaleFactor[rfield.ipSpec],
1209 TorF(rfield.lgTimeVary[rfield.ipSpec]));
1210 }
1211 fprintf(ioQQQ,"\n");
1212
1213 /* print some useful pointers to ionization edges */
1214 fprintf( ioQQQ, " H2,1=%5ld%5ld NX=%5ld IRC=%5ld\n",
1215 iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon,
1216 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon,
1217 opac.ipCKshell,
1218 ionbal.ipCompRecoil[ipHYDROGEN][0] );
1219 fprintf( ioQQQ, " CARBON" );
1220 for( i=0; i < 6; i++ )
1221 fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipCARBON][i] );
1222 fprintf( ioQQQ, "\n" );
1223
1224 fprintf( ioQQQ, " OXY" );
1225 for( i=0; i < 8; i++ )
1226 fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipOXYGEN][i] );
1227 fprintf( ioQQQ, "%5ld%5ld%5ld\n", opac.ipo3exc[0],
1228 oxy.i2d, oxy.i2p );
1229
1230 double sum = 0.;
1231 ASSERT( rfield.ipG0_DB96_lo>0 );
1232 for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; i++ )
1233 sum += rfield.flux[0][i];
1234 fprintf(ioQQQ,"\n Sum DB96 photons %.2e", sum);
1235 sum = 0.;
1236 for( i=(Heavy.ipHeavy[ipHYDROGEN][0] - 1); i<rfield.nflux; i++ )
1237 sum += rfield.flux[0][i];
1238 fprintf(ioQQQ,", sum H-ioniz photons %.2e\n", sum);
1239
1240
1241 fprintf( ioQQQ, "\n\n PHOTONS PER CELL (NOT RYD)\n" );
1242 fprintf( ioQQQ, " nu, flux, wid, occ \n" );
1243 for( i=0; i < rfield.nflux; i++ )
1244 {
1245 fprintf( ioQQQ, "%4ld%10.2e%10.2e%10.2e%10.2e\n", i,
1246 rfield.anu[i], rfield.flux[0][i], rfield.widflx[i],
1247 rfield.OccNumbIncidCont[i] );
1248 }
1249 fprintf( ioQQQ, " \n" );
1250 }
1251
1252 /* zero out some continua related to the ots rates,
1253 * prototype and routine in RT_OTS_Update. This is done here since summed cont will
1254 * be set to rfield */
1255 RT_OTS_Zero();
1256
1257 if( trace.lgTrace )
1258 {
1259 fprintf( ioQQQ, " ContSetIntensity returns, nflux=%5ld anu(nflux)=%11.4e eden=%10.2e\n",
1260 rfield.nflux, rfield.anu[rfield.nflux-1], dense.eden );
1261 }
1262
1263 return;
1264}
1265
1266/*sumcon sums L and Q for net incident continuum */
1267STATIC void sumcon(long int il,
1268 long int ih,
1269 realnum *q,
1270 realnum *p,
1271 realnum *panu)
1272{
1273 long int i,
1274 iupper; /* used as upper limit to the sum */
1275
1276 DEBUG_ENTRY( "sumcon()" );
1277
1278 *q = 0.;
1279 *p = 0.;
1280 *panu = 0.;
1281
1282 /* soft continua may not go as high as the requested bin */
1283 iupper = MIN2(rfield.nflux,ih);
1284
1285 /* n.b. - in F77 loop IS NOT executed when IUPPER < IL */
1286 for( i=il-1; i < iupper; i++ )
1287 {
1288 /* sum photon number */
1289 *q += rfield.flux[0][i];
1290 /* en1ryd is needed to stop overflow */
1291 /* sum flux */
1292 *p += (realnum)(rfield.flux[0][i]*(rfield.anu[i]*EN1RYD));
1293 /* this sum needed for means */
1294 *panu += (realnum)(rfield.flux[0][i]*(rfield.anu2[i]*EN1RYD));
1295 }
1296
1297 return;
1298}
1299
1300/*ptrcer show continuum pointers in real time following drive pointers command */
1302{
1303 char chCard[INPUT_LINE_LENGTH];
1304 /* in case of checking everything, will write errors to this file */
1305 FILE * ioERRORS=NULL;
1306 bool lgEOL;
1307 char chKey;
1308 long int i,
1309 ipnt,
1310 j;
1311 double pnt,
1312 t1,
1313 t2;
1314
1315 DEBUG_ENTRY( "ptrcer()" );
1316
1317 fprintf( ioQQQ, " There are two ways to do this:\n");
1318 fprintf( ioQQQ, " do you want me to test all the pointers (enter y)\n");
1319 fprintf( ioQQQ, " or do you want to enter energies yourself? (enter n)\n" );
1320
1321 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
1322 {
1323 fprintf( ioQQQ, " error getting input \n" );
1325 }
1326
1327 /* this must be either y or n */
1328 chKey = chCard[0];
1329
1330 if( chKey == 'n' )
1331 {
1332 /* this branch, enter energies by hand, and see what happens */
1333 fprintf( ioQQQ, " Enter energy (Ryd); 0 to stop; negative is log.\n" );
1334 pnt = 1.;
1335 while( pnt!=0. )
1336 {
1337 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
1338 {
1339 fprintf( ioQQQ, " error getting input2 \n" );
1341 }
1342 /* now get the number off the line */
1343 i = 1;
1344 pnt = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
1345
1346 /* bail if no number at all, or it is zero*/
1347 if( lgEOL || pnt==0. )
1348 {
1349 break;
1350 }
1351
1352 /* if number negative then interpret as log */
1353 if( pnt < 0. )
1354 {
1355 pnt = pow(10.,pnt);
1356 }
1357
1358 /* get pointer to call */
1359 ipnt = ipoint(pnt);
1360 fprintf( ioQQQ, " Cell num%4ld center:%10.2e width:%10.2e low:%10.2e hi:%10.2e convoc:%10.2e\n",
1361 ipnt, rfield.anu[ipnt-1], rfield.widflx[ipnt-1],
1362 rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]/2.,
1363 rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]/2.,
1364 rfield.convoc[ipnt-1] );
1365 }
1366 }
1367
1368 else if( chKey == 'y' )
1369 {
1370 /* first check that ipoint will not crash due to out of range call*/
1371 if( rfield.anu[0] - rfield.widflx[0]/2.*0.9 < continuum.filbnd[0] )
1372 {
1373 fprintf( ioQQQ," ipoint would crash since lowest desired energy of %e ryd is below limit of %e\n",
1374 rfield.anu[0] - rfield.widflx[0]/2.*0.9 , continuum.filbnd[0] );
1375 fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[0]);
1377 }
1378
1379 else if( rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 >
1380 continuum.filbnd[continuum.nrange] )
1381 {
1382 fprintf( ioQQQ," ipoint would crash since highest desired energy of %e ryd is above limit of %e\n",
1383 rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 ,
1384 continuum.filbnd[continuum.nrange-1] );
1385 fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[rfield.nflux]);
1386 fprintf( ioQQQ," this, previous cells are %e %e\n",
1387 rfield.anu[rfield.nflux-1],rfield.anu[rfield.nflux-2]);
1389 }
1390
1391 /* this branch check everything, write errors to error file */
1392 fprintf( ioQQQ, " errors output on errors.txt\n");
1393 fprintf( ioQQQ, " IP(cor),IP(fount),nu lower, upper of found, desired cell.\n" );
1394
1395 /* error file not open, set to null so we can check later */
1396 ioERRORS = NULL;
1397 for( i=0; i < rfield.nflux-1; i++ )
1398 {
1399 t1 = rfield.anu[i] - rfield.widflx[i]/2.*0.9;
1400 t2 = rfield.anu[i] + rfield.widflx[i]/2.*0.9;
1401
1402 j = ipoint(t1);
1403 if( j != i+1 )
1404 {
1405 /* open file for errors if not already open */
1406 if( ioERRORS == NULL )
1407 {
1408 ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY );
1409 fprintf( ioQQQ," created errors.txt file with error summary\n");
1410 }
1411
1412 fprintf( ioQQQ, " Pointers do not agree for lower bound of cell%4ld, %e\n",
1413 i, rfield.anu[i]);
1414 fprintf( ioERRORS, " Pointers do not agree for lower bound of cell%4ld, %e\n",
1415 i, rfield.anu[i] );
1416 }
1417
1418 j = ipoint(t2);
1419 if( j != i+1 )
1420 {
1421 /* open file for errors if not already open */
1422 if( ioERRORS == NULL )
1423 {
1424 ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY );
1425 fprintf( ioQQQ," created errors.txt file with error summary\n");
1426 }
1427 fprintf( ioQQQ, " Pointers do not agree for upper bound of cell%4ld, %e\n",
1428 i , rfield.anu[i]);
1429 fprintf( ioERRORS, " Pointers do not agree for upper bound of cell%4ld, %e\n",
1430 i , rfield.anu[i]);
1431 }
1432
1433 }
1434 }
1435
1436 else
1437 {
1438 fprintf( ioQQQ, "I do not understand this key, sorry. %c\n", chKey );
1440 }
1441
1442 if( ioERRORS!=NULL )
1443 fclose( ioERRORS );
1445}
1446
1447/*extin do extinction of incident continuum as set by extinguish command */
1448STATIC void extin(realnum *ex1ryd)
1449{
1450
1451 DEBUG_ENTRY( "extin()" );
1452
1453 /* modify input continuum by leaky absorber
1454 * power law fit to
1455 * >>refer XUV extinction Cruddace, R., Paresce, F., Bowyer, S., & Lampton, M. 1974, ApJ, 187, 497. */
1456 if( rfield.ExtinguishColumnDensity == 0. )
1457 {
1458 *ex1ryd = 1.;
1459
1460 for( long i=0; i<rfield.nupper; ++i )
1461 rfield.ExtinguishFactor[i] = 1.;
1462 }
1463 else
1464 {
1465 double absorb = 1. - rfield.ExtinguishLeakage;
1466 double factor = rfield.ExtinguishColumnDensity*
1467 rfield.ExtinguishConvertColDen2OptDepth;
1468 /* extinction at 1 4 Ryd */
1469 *ex1ryd = (realnum)(rfield.ExtinguishLeakage + absorb*sexp(factor));
1470
1471 // low energy limit of extinction
1472 long low = ipoint(rfield.ExtinguishLowEnergyLimit);
1473
1474 for( long i=0; i<low-1; ++i )
1475 rfield.ExtinguishFactor[i] = 1.;
1476
1477 for( long i=low-1; i < rfield.nupper; i++ )
1478 {
1479 realnum extfactor = (realnum)(rfield.ExtinguishLeakage + absorb*
1480 sexp(factor*(pow(rfield.anu[i],(double)rfield.ExtinguishEnergyPowerLow))));
1481
1482 rfield.ExtinguishFactor[i] = extfactor;
1483 rfield.flux_beam_const[i] *= extfactor;
1484 rfield.flux_beam_time[i] *= extfactor;
1485 rfield.flux_isotropic[i] *= extfactor;
1486
1487 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
1488 rfield.flux_isotropic[i];
1489 }
1490 }
1491 return;
1492}
1493
1494/*conorm normalize continuum to proper intensity */
1496{
1497 long int i;
1498 double xLog_radius_inner,
1499 diff,
1500 f,
1501 flx1,
1502 flx2,
1503 pentrd,
1504 qentrd;
1505
1506 DEBUG_ENTRY( "conorm()" );
1507
1508 xLog_radius_inner = log10(radius.rinner);
1509
1510 /* this is a sanity check, it can't happen */
1511 for( i=0; i < rfield.nShape; i++ )
1512 {
1513 if( strcmp(rfield.chRSpec[i],"UNKN") == 0 )
1514 {
1515 fprintf( ioQQQ, " UNKN spectral normalization cannot happen.\n" );
1516 fprintf( ioQQQ, " conorm punts.\n" );
1518 }
1519
1520 else if( strcmp(rfield.chRSpec[i],"SQCM") != 0 &&
1521 strcmp(rfield.chRSpec[i],"4 PI") != 0 )
1522 {
1523 fprintf( ioQQQ, " chRSpec must be SQCM or 4 PI, and it was %4.4s. This cannot happen.\n",
1524 rfield.chRSpec[i] );
1525 fprintf( ioQQQ, " conorm punts.\n" );
1527 }
1528
1529
1530 /* this sanity check makes sure that atlas.mod or werner.mod grids
1531 * are for the current version of the code */
1532 if( strcmp(rfield.chSpType[i],"VOLK ") == 0 )
1533 {
1534 if( !fp_equal( rfield.RSFCheck[i], continuum.ResolutionScaleFactor ) )
1535 {
1536 fprintf( ioQQQ,"\n\n PROBLEM DISASTER At least one of the compiled stellar atmosphere"
1537 " grids has been compiled with a different energy grid resolution factor.\n" );
1538 fprintf( ioQQQ, " Please recompile this file using the COMPILE STARS command "
1539 "and make sure that you use the correct SET CONTINUUM RESOLUTION factor.\n" );
1541 }
1542 }
1543 else if( strcmp(rfield.chSpType[i],"READ ") == 0 )
1544 {
1545 if( !fp_equal( rfield.RSFCheck[i], continuum.ResolutionScaleFactor ) )
1546 {
1547 fprintf( ioQQQ,"\n\n PROBLEM DISASTER The file read by the TABLE READ command "
1548 "has been compiled with a different energy grid resolution factor.\n" );
1549 fprintf( ioQQQ, " Please recompile this file using the SAVE TRANSMITTED CONTINUUM "
1550 "command and use the correct SET CONTINUUM RESOLUTION factor.\n" );
1552 }
1553 }
1554 }
1555
1556 /* this sanity check is that the grains we have read in from opacity files agree
1557 * with the energy grid in this version of cloudy */
1558 for( size_t nd=0; nd < gv.bin.size(); nd++ )
1559 {
1560 if( !fp_equal( gv.bin[nd]->RSFCheck, continuum.ResolutionScaleFactor ) )
1561 {
1562 fprintf( ioQQQ,"\n\n PROBLEM DISASTER At least one of the grain opacity files "
1563 "has been compiled with a different energy grid resolution factor.\n" );
1564 fprintf( ioQQQ, " Please recompile this file using the COMPILE GRAINS command "
1565 "and make sure that you use the correct SET CONTINUUM RESOLUTION factor.\n" );
1567 }
1568 }
1569
1570 /* default is is to predict line intensities,
1571 * but if any continuum specified as luminosity then would override this -
1572 * following two values are correct for intensities */
1573 radius.pirsq = 0.;
1574 radius.lgPredLumin = false;
1575
1576 /* check whether ANY luminosities are present */
1577 for( i=0; i < rfield.nShape; i++ )
1578 {
1579 if( strcmp(rfield.chRSpec[i],"4 PI") == 0 )
1580 {
1581 radius.pirsq = (realnum)(1.0992099 + 2.*xLog_radius_inner);
1582 radius.lgPredLumin = true;
1583 /* convert down to intensity */
1584 rfield.totpow[i] -= radius.pirsq;
1585
1586 if( trace.lgTrace )
1587 {
1588 fprintf( ioQQQ,
1589 " conorm converts continuum %ld from luminosity to intensity.\n",
1590 i );
1591 }
1592 }
1593 }
1594
1595 /* if total luminosities are present, must have specified a starting radius */
1596 if( radius.lgPredLumin && !radius.lgRadiusKnown )
1597 {
1598 fprintf(ioQQQ,"PROBLEM DISASTER conorm: - A continuum source was specified as a luminosity, but the inner radius of the cloud was not set.\n");
1599 fprintf(ioQQQ,"Please set an inner radius.\nSorry.\n");
1601 }
1602
1603 /* convert ionization parameters to number of photons, called "q(h)"
1604 * at this stage q(h) and "PHI " are the same */
1605 for( i=0; i < rfield.nShape; i++ )
1606 {
1607 if( strcmp(rfield.chSpNorm[i],"IONI") == 0 )
1608 {
1609 /* the log of the ionization parameter was stored here, this converts
1610 * it to the log of the number of photons per sq cm */
1611 rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) + log10(SPEEDLIGHT);
1612 strcpy( rfield.chSpNorm[i], "Q(H)" );
1613 if( trace.lgTrace )
1614 {
1615 fprintf( ioQQQ,
1616 " conorm converts continuum %ld from ionizat par to q(h).\n",
1617 i );
1618 }
1619 }
1620 }
1621
1622 /* convert x-ray ionization parameter xi to intensity */
1623 for( i=0; i < rfield.nShape; i++ )
1624 {
1625 if( strcmp(rfield.chSpNorm[i],"IONX") == 0 )
1626 {
1627 /* this converts it to an intensity */
1628 rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) - log10(PI4);
1629 strcpy( rfield.chSpNorm[i], "LUMI" );
1630 if( trace.lgTrace )
1631 {
1632 fprintf( ioQQQ, " conorm converts continuum%3ld from x-ray ionizat par to I.\n",
1633 i );
1634 }
1635 }
1636 }
1637
1638 /* indicate whether we ended up with luminosity or intensity */
1639 if( trace.lgTrace )
1640 {
1641 if( radius.lgPredLumin )
1642 {
1643 fprintf( ioQQQ, " Cloudy will predict lumin into 4pi\n" );
1644 }
1645 else
1646 {
1647 fprintf( ioQQQ, " Cloudy will do surface flux for lumin\n" );
1648 }
1649 }
1650
1651 /* if intensity per unit area is predicted then geometric
1652 * covering factor must be unity
1653 * variable can also be set elsewhere */
1654 if( !radius.lgPredLumin )
1655 {
1656 geometry.covgeo = 1.;
1657 }
1658
1659 /* main loop over all continuum shapes to find continuum normalization
1660 * for each one */
1661 for( i=0; i < rfield.nShape; i++ )
1662 {
1663 rfield.ipSpec = i;
1664
1665 /* check that, if laser, bounds include laser energy */
1666 if( strcmp(rfield.chSpType[rfield.ipSpec],"LASER") == 0 )
1667 {
1668 if( !( rfield.range[rfield.ipSpec][0] < rfield.slope[rfield.ipSpec] &&
1669 rfield.range[rfield.ipSpec][1] > rfield.slope[rfield.ipSpec]) )
1670 {
1671 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",
1672 rfield.slope[rfield.ipSpec],
1673 rfield.range[rfield.ipSpec][0],
1674 rfield.range[rfield.ipSpec][1]);
1675 fprintf(ioQQQ,"Please specify the continuum flux where the laser is active.\n");
1677 }
1678 }
1679
1680 if( trace.lgTrace )
1681 {
1682 long int jj;
1683 fprintf( ioQQQ, " conorm continuum number %ld is shape %s range is %.2e %.2e\n",
1684 i,
1685 rfield.chSpType[i],
1686 rfield.range[i][0],
1687 rfield.range[i][1] );
1688 fprintf( ioQQQ, "the continuum points follow\n");
1689 jj = 0;
1690 if( rfield.lgContMalloc[rfield.ipSpec] )
1691 {
1692 while( rfield.tNu[rfield.ipSpec][jj].Ryd() != 0. && jj < 100 )
1693 {
1694 fprintf( ioQQQ, "%li %e %e\n",
1695 jj ,
1696 rfield.tNu[rfield.ipSpec][jj].Ryd(),
1697 rfield.tslop[rfield.ipSpec][jj] );
1698 ++jj;
1699 }
1700 }
1701 }
1702
1703 if( strcmp(rfield.chSpNorm[i],"RATI") == 0 )
1704 {
1705 /* option to scale relative to previous continua
1706 * this must come first since otherwise may be too late
1707 * BUT ratio cannot be the first continuum source */
1708 if( trace.lgTrace )
1709 {
1710 fprintf( ioQQQ, " conorm this is ratio to 1st con\n" );
1711 }
1712
1713 /* check that this is not the first continuum source, we must ratio */
1714 if( i == 0 )
1715 {
1716 fprintf( ioQQQ, " I cant form a ratio if continuum is first source.\n" );
1718 }
1719
1720 /* first find photon flux and Q of previous continuum */
1721 rfield.ipSpec -= 1;
1722 flx1 = ffun1(rfield.range[i][0])*rfield.spfac[rfield.ipSpec]*
1723 rfield.range[i][0];
1724
1725 /* check that previous continua were not zero where ratio is formed */
1726 if( flx1 <= 0. )
1727 {
1728 fprintf( ioQQQ, " Previous continua were zero where ratio is desired.\n" );
1730 }
1731
1732 /* return pointer to previous (correct) value, find F, Q */
1733 rfield.ipSpec += 1;
1734
1735 /* we want a continuum totpow as powerful, flx is now desired flx */
1736 flx1 *= rfield.totpow[i];
1737
1738 /*. find flux of this new continuum at that point */
1739 flx2 = ffun1(rfield.range[i][1])*rfield.range[i][1];
1740
1741 /* this is ratio of desired to actual */
1742 rfield.spfac[i] = flx1/flx2;
1743 if( trace.lgTrace )
1744 {
1745 fprintf( ioQQQ, " conorm ratio will set scale fac to%10.3e at%10.2e Ryd.\n",
1746 rfield.totpow[i], rfield.range[i][0] );
1747 }
1748 }
1749
1750 else if( strcmp(rfield.chSpNorm[i],"FLUX") == 0 )
1751 {
1752 /* specify flux density
1753 * option to use arbitrary frequency or range */
1754 f = ffun1(rfield.range[i][0]);
1755
1756 /* make sure this is positive, could be zero if out of range of table,
1757 * or for various forms of insanity */
1758 if( f<=SMALLDOUBLE )
1759 {
1760 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",
1761 i , rfield.range[i][0]);
1762 /* is this a table? if so, is ffun within its bounds? */
1763 if( strcmp(rfield.chSpType[i],"INTER") == 0 )
1764 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");
1765 fprintf( ioQQQ, " Also check that the numbers used to specify the shape and intensity do not under or overflow on this cpu.\n\n");
1766
1768 }
1769
1770 /* now convert to log in continuum units we shall use */
1771 f = log10(f) + log10(rfield.range[i][0]*EN1RYD/FR1RYD);
1772 f = rfield.totpow[i] - f;
1773 rfield.spfac[i] = pow(10.,f);
1774
1775 if( trace.lgTrace )
1776 {
1777 fprintf( ioQQQ, " conorm will set log fnu to%10.3e at%10.2e Ryd. Factor=%11.4e\n",
1778 rfield.totpow[i], rfield.range[i][0], rfield.spfac[i] );
1779 }
1780 }
1781
1782 else if( strcmp(rfield.chSpNorm[i],"Q(H)") == 0 ||
1783 strcmp(rfield.chSpNorm[i],"PHI ") == 0 )
1784 {
1785 /* some type of photon density entered */
1786 if( trace.lgTrace )
1787 {
1788 fprintf( ioQQQ, " conorm calling qintr range=%11.3e %11.3e desired val is %11.3e\n",
1789 rfield.range[i][0],
1790 rfield.range[i][1] ,
1791 rfield.totpow[i]);
1792 }
1793
1794 /* the total number of photons over the specified range in
1795 * the arbitrary system of units that the code save the continuum shape */
1796 ASSERT( rfield.range[i][0] < rfield.range[i][1] );
1797 qentrd = qintr(&rfield.range[i][0],&rfield.range[i][1]);
1798 /* this is the log of the scale factor that must multiply the
1799 * continuum shape to get the final set of numbers */
1800 diff = rfield.totpow[i] - qentrd;
1801
1802 /* >>chng 03 mar 13, from diff < -25 to <-35,
1803 * tripped for very low U models used for H2 simulations */
1804 /*if( diff < -25. || diff > 35. )*/
1805 if( diff < -35. || diff > 35. )
1806 {
1807 fprintf( ioQQQ, " PROBLEM DISASTER Continuum source specified is too extreme.\n" );
1808 fprintf( ioQQQ,
1809 " The integral over the continuum shape gave (log) %.3e photons, and the command requested (log) %.3e.\n" ,
1810 qentrd , rfield.totpow[i]);
1811 fprintf( ioQQQ,
1812 " The difference in the log is %.3e.\n" ,
1813 diff );
1814 if( diff>0. )
1815 {
1816 fprintf( ioQQQ, " The continuum source is too bright.\n" );
1817 }
1818 else
1819 {
1820 fprintf( ioQQQ, " The continuum source is too faint.\n" );
1821 }
1822 /* explain what happened */
1823 fprintf( ioQQQ, " The usual cause for this problem is an incorrect continuum intensity/luminosity or radius command.\n" );
1824 fprintf( ioQQQ, " There were a total of %li continuum shape commands entered - the problem is with number %li.\n",
1825 rfield.nShape , i+1 );
1827 }
1828
1829 else
1830 {
1831 rfield.spfac[i] = pow(10.,diff);
1832 }
1833
1834 if( trace.lgTrace )
1835 {
1836 fprintf( ioQQQ, " conorm finds Q over range from%11.4e-%11.4e Ryd, integral= %10.4e Factor=%11.4e\n",
1837 rfield.range[i][0],
1838 rfield.range[i][1],
1839 qentrd ,
1840 rfield.spfac[i] );
1841 }
1842 }
1843
1844 else if( strcmp(rfield.chSpNorm[i],"LUMI") == 0 )
1845 {
1846 /* luminosity entered, special since default is TOTAL lumin */
1847 /*pintr integrates L for any continuum between two limits, used for normalization,
1848 * return units are log of ryd cm-2 s-1, last log conv to ergs */
1849 pentrd = pintr(rfield.range[i][0],rfield.range[i][1]) + log10(EN1RYD);
1850 f = rfield.totpow[i] - pentrd;
1851 rfield.spfac[i] = pow(10.,f);
1852
1853 if( trace.lgTrace )
1854 {
1855 fprintf( ioQQQ, " conorm finds luminosity range is %10.3e to %9.3e Ryd, factor is %11.4e\n",
1856 rfield.range[i][0], rfield.range[i][1],
1857 rfield.spfac[i] );
1858 }
1859 }
1860
1861 else
1862 {
1863 fprintf( ioQQQ, "PROBLEM DISASTER What chSpNorm label is this? =%s=\n", rfield.chSpNorm[i]);
1864 TotalInsanity();
1865 }
1866
1867 /* spfac used to renormalize SED into flux */
1868 if( fabs(rfield.spfac[i]) <=SMALLDOUBLE )
1869 {
1870 fprintf( ioQQQ, "PROBLEM DISASTER conorm finds infinite continuum scale factor.\n" );
1871 fprintf( ioQQQ, "The continuum is too intense to compute with this cpu.\n" );
1872 fprintf( ioQQQ, "Were the intensity and luminosity commands switched?\n" );
1873 fprintf( ioQQQ, "Sorry, but I cannot go on.\n" );
1875 }
1876 }
1877
1878 /* this is conversion factor for final units of line intensities or luminosities in printout,
1879 * will be intensities (==0) unless luminosity is to be printed, or flux at Earth
1880 * pirsq is the log of 4 pi r_in^2 */
1881 radius.Conv2PrtInten = radius.pirsq;
1882
1883 /* >>chng 02 apr 25, add option for slit on aperture command */
1884 if( geometry.iEmissPower == 1 )
1885 {
1886 if( radius.lgPredLumin )
1887 {
1888 /* factor should be divided by 2 r_in (so that 2pi*r_in remains) */
1889 radius.Conv2PrtInten -= (log10(2.) + xLog_radius_inner);
1890 }
1891 else if( !radius.lgPredLumin )
1892 {
1893 /* this is an error - slit requested but radius is not known */
1894 fprintf( ioQQQ, "PROBLEM DISASTER conorm: Aperture slit specified, but not predicting luminosity.\n" );
1895 fprintf( ioQQQ, "conorm: Please specify an inner radius to determine L.\nSorry\n" );
1897 }
1898 }
1899 if( geometry.iEmissPower == 0 && radius.lgPredLumin )
1900 {
1901 /* leave Conv2PrtInten at zero if not predicting luminosity */
1902 radius.Conv2PrtInten = log10(2.);
1903 }
1904
1905 /* this is option to give final absolute results as flux observed at Earth */
1906 if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth )
1907 {
1908 /* this implements the conversion from Q_alpha -> F_alpha
1909 * described in the section on the APERTURE command in Hazy 1 */
1910 if( geometry.iEmissPower == 0 )
1911 radius.Conv2PrtInten += log10(double(geometry.size)/SQAS_SKY);
1912 else if( geometry.iEmissPower == 1 )
1913 radius.Conv2PrtInten += log10(double(geometry.size)/(PI4*AS1RAD*radius.distance));
1914 else if( geometry.iEmissPower == 2 )
1915 radius.Conv2PrtInten -= log10( PI4*pow2(radius.distance) );
1916 else
1917 TotalInsanity();
1918 }
1919
1920 /* normally lines are into 4pi, this is option to do per sr or arcsec^2 */
1921 if( prt.lgSurfaceBrightness )
1922 {
1923 if( radius.pirsq != 0. )
1924 {
1925 /* make sure we are predicting line intensities, not luminosity */
1926 fprintf( ioQQQ, " PROBLEM DISASTER Sorry, but both luminosity and surface brightness have been requested for lines.\n" );
1927 fprintf( ioQQQ, " the PRINT LINE SURFACE BRIGHTNESS command can only be used when lines are predicted per unit cloud area.\n" );
1929 }
1930 if( prt.lgSurfaceBrightness_SR )
1931 {
1932 /* we want final units to be per sr */
1933 radius.Conv2PrtInten -= log10( PI4 );
1934 }
1935 else
1936 {
1937 /* we want final units to be per square arcsec */
1938 radius.Conv2PrtInten -= log10(SQAS_SKY);
1939 }
1940 }
1941 return;
1942}
1943
1944/*qintr integrates Q for any continuum between two limits, used for normalization */
1945STATIC double qintr(double *qenlo,
1946 double *qenhi)
1947{
1948 long int i,
1949 ipHi,
1950 ipLo,
1951 j;
1952 double qintr_v,
1953 sum,
1954 wanu;
1955
1956 DEBUG_ENTRY( "qintr()" );
1957
1958 /* returns LOG of number of photons over energy interval */
1959
1960 /* this is copy of logic that occurs three times across code */
1961 ASSERT(*qenhi > *qenlo);
1962 ipLo = ipoint(*qenlo);
1963 ipHi = ipoint(*qenhi);
1964 /* this is actual sum of photons within band */
1965 sum = 0.;
1966 for( i=ipLo-1; i < (ipHi - 1); i++ )
1967 {
1968 /*sum += ffun1(rfield.anu[i])*rfield.widflx[i];*/
1969 for( j=0; j < 4; j++ )
1970 {
1971 wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
1972 /* >>chng 02 jul 16, add test on continuum limits -
1973 * this was exceeded when resolution set very large */
1974 wanu = MAX2( wanu , rfield.emm );
1975 wanu = MIN2( wanu , rfield.egamry );
1976 sum += fweigh[j]*ffun1(wanu)*rfield.widflx[i];
1977 }
1978 }
1979
1980 if( sum <= 0. )
1981 {
1982 fprintf( ioQQQ, " PROBLEM DISASTER Photon number sum in QINTR is %.3e\n",
1983 sum );
1984 fprintf( ioQQQ, " This source has no ionizing radiation, and the number of ionizing photons was specified.\n" );
1985 fprintf( ioQQQ, " This was continuum source number%3ld\n",
1986 rfield.ipSpec );
1987 fprintf( ioQQQ, " Sorry, but I cannot go on. ANU and FLUX arrays follow. Enjoy.\n" );
1988 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" );
1989 for( i=0; i < rfield.nupper; i++ )
1990 {
1991 fprintf( ioQQQ, "%.2e\t%.2e\n",
1992 rfield.anu[i],
1993 rfield.flux[0][i] );
1994 }
1996 }
1997 else
1998 {
1999 qintr_v = log10(sum);
2000 }
2001 return qintr_v;
2002}
2003
2004/*pintr integrates L for any continuum between two limits, used for normalization,
2005 * return units are log of ryd cm-2 s-1 */
2006STATIC double pintr(double penlo,
2007 double penhi)
2008{
2009 long int i,
2010 j;
2011 double fsum,
2012 pintr_v,
2013 sum,
2014 wanu,
2015 wfun;
2016 long int ip1 , ip2;
2017
2018 DEBUG_ENTRY( "pintr()" );
2019
2020 /* computes log of luminosity in radiation over some integral
2021 * answer is in Ryd per sec */
2022
2023 sum = 0.;
2024 /* >>chng 02 oct 27, do not call qg32, do same type sum as
2025 * final integration */
2026 /* laser is special since delta function, this is center of laser */
2027 /* >>chng 01 jul 01, was +-21 cells, change to call to ipoint */
2028 ip1 = ipoint( penlo );
2029
2030 ip2 = ipoint( penhi );
2031
2032 for( i=ip1-1; i < ip2-1; i++ )
2033 {
2034 fsum = 0.;
2035 for( j=0; j < 4; j++ )
2036 {
2037 wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
2038 /*++iiii;fprintf(ioQQQ,"DEBUG iii %li %e \n",iiii, wanu );*/
2039 wfun = fweigh[j]*ffun1(wanu)*wanu;
2040 fsum += wfun;
2041 }
2042 sum += fsum*rfield.widflx[i];
2043 }
2044
2045 if( sum > 0. )
2046 {
2047 pintr_v = log10(sum);
2048 }
2049 else
2050 {
2051 pintr_v = -38.;
2052 }
2053
2054 return pintr_v;
2055}
t_called called
Definition called.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
FILE * ioStdin
Definition cddefines.cpp:8
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
#define POW3
Definition cddefines.h:936
const int ipCARBON
Definition cddefines.h:310
T pow2(T a)
Definition cddefines.h:931
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
const int ipLITHIUM
Definition cddefines.h:307
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
char TorF(bool l)
Definition cddefines.h:710
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long max(int a, long b)
Definition cddefines.h:775
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
static t_ADfA & Inst()
Definition cddefines.h:175
double coll_ion_wrapper(long int z, long int n, double t)
double ffun1(double xnu)
double ffun(double anu)
Definition cont_ffun.cpp:17
long ipoint(double energy_ryd)
STATIC void ptrcer()
static const double aweigh[4]
STATIC void sumcon(long int il, long int ih, realnum *q, realnum *p, realnum *panu)
static const double fweigh[4]
STATIC double qintr(double *qenlo, double *qenhi)
STATIC void conorm()
STATIC double pintr(double penlo, double penhi)
void ContSetIntensity()
STATIC void extin(realnum *ex1ryd)
void IncidentContinuumHere()
t_continuum continuum
Definition continuum.cpp:5
void EdenChange(double EdenNew)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
@ AS_LOCAL_ONLY
Definition cpu.h:208
const double SMALLDOUBLE
Definition cpu.h:195
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_geometry geometry
Definition geometry.cpp:5
GrainVar gv
Definition grainvar.cpp:5
t_Heavy Heavy
Definition heavy.cpp:5
t_hextra hextra
Definition hextra.cpp:5
t_hydro hydro
Definition hydrogenic.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH1s
Definition iso.h:27
void iso_continuum_lower(long ipISO, long nelem)
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
t_noexec noexec
Definition noexec.cpp:5
t_opac opac
Definition opacity.cpp:5
t_oxy oxy
Definition oxy.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double AS1RAD
Definition physconst.h:129
UNUSED const double FR1RYD
Definition physconst.h:195
UNUSED const double PI
Definition physconst.h:29
UNUSED const double HPLANCK
Definition physconst.h:103
UNUSED const double PI4
Definition physconst.h:35
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double SQAS_SKY
Definition physconst.h:135
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double HNU3C2
Definition physconst.h:198
UNUSED const double ELEM_CHARGE_ESU
Definition physconst.h:147
UNUSED const double TE1RYD
Definition physconst.h:183
t_prt prt
Definition prt.cpp:10
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
void RT_OTS_Zero(void)
Definition rt_ots.cpp:599
t_secondaries secondaries
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
t_timesc timesc
Definition timesc.cpp:5
t_trace trace
Definition trace.cpp:5