root/branches/newmole/source/cont_setintensity.cpp

Revision 2601, 61.9 kB (checked in by rporter, 8 days ago)

branches/newmole - bring in r2595:2600

  • Property svn:eol-style set to native
Line 
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 */
39static double aweigh[4]={-0.4305682,-0.1699905, 0.1699905, 0.4305682};
40static double fweigh[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274};
41
42/*conorm normalize continuum to proper intensity */
43STATIC void conorm(void);
44
45/*pintr integrates L for any continuum between two limits, used for normalization */
46STATIC double pintr(double penlo,
47  double penhi);
48
49/*qintr integrates Q for any continuum between two limits, used for normalization */
50STATIC double qintr(double *qenlo,
51  double *qenhi);
52
53
54/*sumcon sums L and Q for net incident continuum */
55STATIC 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 */
62STATIC void extin(realnum *ex1ryd);
63
64/*ptrcer show continuum pointers in real time following drive pointers command */
65STATIC 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 */
69int ContSetIntensity(void)
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[0][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[0][i] * (realnum)frac_beam_const;
211                rfield.flux_beam_time[i] = rfield.flux[0][i] * (realnum)frac_beam_time;
212                rfield.flux_isotropic[i] = rfield.flux[0][i] * (realnum)frac_isotropic;
213
214                if( rfield.flux[0][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) ) these are not guaranteed to
231                         * be monotonically increasing due to loss of precision in log
232                         * of float */
233                        rfield.anulog[i] = (realnum)log10(rfield.anu[i]);
234                }
235
236                else if( rfield.flux[0][i] == 0. )
237                {
238                        rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
239                        rfield.anulog[i] = (realnum)log10(rfield.anu[i]);
240                }
241
242                else
243                {
244                        rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
245                        fprintf( ioQQQ, " negative continuum returned at%6ld%10.2e%10.2e\n",
246                          i, rfield.anu[i], rfield.flux[0][i] );
247                        lgCheckOK = false;
248                }
249                rfield.anu3[i] = rfield.anu2[i]*rfield.anu[i];
250
251                rfield.ConEmitReflec[0][i] = 0.;
252                rfield.ConEmitOut[0][i] = 0.;
253                rfield.convoc[i] = factor/rfield.widflx[i]/rfield.anu2[i];
254
255                /* following are Compton exchange factors from Tarter */
256                alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i]));
257                bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/
258                  4.;
259                rfield.csigh[i] = (realnum)(alf*rfield.anu2[i]*3.858e-25);
260                rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25);
261        }
262
263        /*i = 76;
264        fprintf(ioQQQ,"\nDEBUG %4li %e \n",
265                i , rfield.anu[i] );*/
266
267        /* >>chng 05 jul 01 add this
268         * finished with stored continua - return these vectors */
269#if 0
270        /* commented out since we must conserve energy, and continuum was set with old widflx */
271        /* now fix widflx array so that it is correct */
272        for( i=1; i<rfield.nupper-1; ++i )
273        {
274                /*rfield.widflx[i] = rfield.anu[i+1] - rfield.anu[i];*/
275                rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
276        }
277#endif
278
279        if( !lgCheckOK )
280        {
281                ShowMe();
282                cdEXIT(EXIT_FAILURE);
283        }
284
285        if( trace.lgTrace && trace.lgComBug )
286        {
287                fprintf( ioQQQ, "\n\n Compton heating, cooling coefficients \n" );
288                for( i=0; i < rfield.nupper; i += 2 )
289                {
290                        fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i],
291                          rfield.csigh[i], rfield.csigc[i] );
292                }
293                fprintf( ioQQQ, "\n" );
294        }
295
296        /* option to check frequencies in real time, drive pointers command,
297         * routine is below, is file static */
298        if( trace.lgPtrace )
299                ptrcer();
300
301        /* extinguish continuum if set on */
302        extin(&ex1ryd);
303
304        /* now find peak of hydrogen ionizing continuum - for PDR calculations
305         * this will remain equal to 1 since the loop will not execute */
306        prt.ipeak = 1;
307        peak = 0.;
308
309        for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
310        {
311                if( rfield.flux[0][i]*rfield.anu[i]/rfield.widflx[i] > (realnum)peak )
312                {
313                        /* prt.ipeak points to largest f_nu at H-ionizing energies
314                         * and is passed to other parts of code */
315                        /* i+1 to keep ipeak on fortran version of energy array */
316                        prt.ipeak = i+1;
317                        peak = rfield.flux[0][i]*rfield.anu[i]/rfield.widflx[i];
318                }
319        }
320
321        /* find highest energy to consider in continuum flux array
322         * peak is the peak product nu*flux */
323        peak = rfield.flux[0][prt.ipeak-1]/rfield.widflx[prt.ipeak-1]*
324          rfield.anu2[prt.ipeak-1];
325
326        /* say what type of cpu this is, if desired */
327        if( trace.lgTrace )
328        {
329                fprintf( ioQQQ, " ContSetIntensity: The peak of the H-ion continuum is at %.3e Ryd - its value is %.2e\n",
330                  rfield.anu[prt.ipeak-1] , peak);
331        }
332
333        if( peak > 1e38 )
334        {
335                fprintf( ioQQQ, " PROBLEM DISASTER The continuum is too intense to compute. Use a fainter continuum. (This is the nu*f_nu test)\n" );
336                fprintf( ioQQQ, " Sorry.\n" );
337                cdEXIT(EXIT_FAILURE);
338        }
339
340        /* FluxFaint set in zero.c; normally 1e-10 */
341        /* this will be faintest level of continuum we want to consider.
342         * peak was set above, is peak of hydrogen ionizing radiation field,
343         * and is zero if no H-ionizing radiation */
344        fntest = peak*rfield.FluxFaint;
345        {
346                enum {DEBUG_LOC=false};
347                /* print flux array then quit */
348                if( DEBUG_LOC )
349                {
350                        for( i=0; i<rfield.nupper; ++i )
351                        {
352                                fprintf(ioQQQ," consetintensityBUGGG\t%.2e\t%.2e\n" ,
353                                        rfield.anu[i] , rfield.flux[0][i]/rfield.widflx[i] );
354                        }
355                        cdEXIT(EXIT_SUCCESS);
356                }
357        }
358
359        if( fntest > 0. )
360        {
361                /* this test is not done in pdr conditions where NO H-ionizing radiation,
362                 * since fntest is zero*/
363                i = rfield.nupper;
364                /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs
365                 * where continuum barely goes into H-ionizing radiation */
366                while( i > prt.ipeak+3 &&
367                        rfield.flux[0][i-1]*rfield.anu2[i-1]/rfield.widflx[i-1] < (realnum)fntest )
368                {
369                        --i;
370                }
371        }
372        else
373        {
374                /* when no H-ionizing radiation set to Lyman edge */
375                /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs
376                 * where continuum barely goes into H-ionizing radiation */
377                i = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]+3;
378        }
379
380        /*
381         * this line of code dates from 1979 and IOA Cambridge.  removed july 19 95
382         * I think it was the last line of the original Cambridge Fortran source
383           nflux = MAX( ineon(1)+4 , i )
384         */
385
386        /* >>chng 99 apr 28, reinstate the rfield.FluxFaint limit with nflux */
387        rfield.nflux = i;
388
389        /* trim down nflux, was set to rfield.nupper, the dimension of all vectors, in zero.c,
390         * in ContCreatePointers was set to nupper, the number of cells needed to get up to the
391         * high-energy limit of the code */
392        while( rfield.flux[0][rfield.nflux-1] < SMALLFLOAT && rfield.nflux > 1 )
393        {
394                --rfield.nflux;
395        }
396        /* make sure we go high enough, avoid 1-off bugs as in draine field */
397        ++rfield.nflux;
398
399        if( rfield.nflux == 1 )
400        {
401                fprintf( ioQQQ, " PROBLEM DISASTER This incident continuum appears to have no radiation.\n" );
402                fprintf( ioQQQ, " Sorry.\n" );
403                cdEXIT(EXIT_FAILURE);
404        }
405
406        /* >>chng 04 oct 10, add this limit - arrays will malloc to nupper, but will add unit
407         * continuum to [nflux] - this must be within array */
408        rfield.nflux = MIN2( rfield.nupper-1 , rfield.nflux );
409
410        /* >>chng 05 mar 01, nfine should not extend beyond continuum
411         * make sure fine opacity scale does not extend beyond continuum we will use */
412        rfield.nfine = rfield.nfine_malloc;
413        while( rfield.nfine > 0 && rfield.fine_anu[rfield.nfine-1] > rfield.anu[rfield.nflux-1] )
414        {
415                --rfield.nfine;
416        }
417
418        /* check that continuum defined everywhere - look for zero's and comment if present */
419        continuum.lgCon0 = false;
420        ip = 0;
421        for( i=0; i < rfield.nflux; i++ )
422        {
423                if( rfield.flux[0][i] == 0. )
424                {
425                        if( called.lgTalk && !continuum.lgCon0 )
426                        {
427                                fprintf( ioQQQ, " NOTE Setcon: continuum has zero intensity starting at %11.4e Ryd.\n",
428                                  rfield.anu[i] );
429                                continuum.lgCon0 = true;
430                        }
431                        ++ip;
432                }
433        }
434
435        if( continuum.lgCon0 && called.lgTalk )
436        {
437                fprintf( ioQQQ,
438                        "%6ld cells in the incident continuum have zero intensity.  Problems???\n\n",
439                  ip );
440        }
441
442        /* check for devastating error in the continuum mesh or intensity */
443        lgCheckOK = true;
444        for( i=1; i < rfield.nflux; i++ )
445        {
446                if( rfield.flux[0][i] < 0. )
447                {
448                        fprintf( ioQQQ,
449                                " PROBLEM DISASTER Continuum has negative intensity at %.4e Ryd=%.2e %4.4s %4.4s\n",
450                          rfield.anu[i], rfield.flux[0][i], rfield.chLineLabel[i], rfield.chContLabel[i] );
451                        lgCheckOK = false;
452                }
453                else if( rfield.anu[i] <= rfield.anu[i-1] )
454                {
455                        fprintf( ioQQQ,
456                                " PROBLEM DISASTER cont_setintensity - internal error - continuum energies not in increasing order: energies follow\n" );
457                        fprintf( ioQQQ,
458                                "%ld %e %ld %e %ld %e\n",
459                          i -1 , rfield.anu[i-1], i, rfield.anu[i], i +1, rfield.anu[i+1] );
460                        lgCheckOK = false;
461                }
462        }
463
464        /* either of the ways lgCheckOK would be set true would be a major internal error */
465        if( !lgCheckOK )
466        {
467                ShowMe();
468                cdEXIT(EXIT_FAILURE);
469        }
470
471        /* turn off recoil ionization if high energy < 190R */
472        if( rfield.anu[rfield.nflux-1] <= 190 )
473        {
474                ionbal.lgCompRecoil = false;
475        }
476
477        /* sum photons and energy, save mean */
478
479        /* sum from low energy to Balmer edge */
480        sumcon(1,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,&rfield.qrad,&prt.pradio,&p1);
481
482        /* sum over Balmer continuum */
483        sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2],iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,&rfield.qbal,&prt.pbal,&p1);
484
485        /* sum from Lyman edge to HeI edge */
486        sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
487                iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1,&prt.q,&p,&p2);
488
489        /* sum from HeI to HeII edges */
490        sumcon(iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0],
491                iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s]-1,&rfield.qhe,&phe,&p3);
492
493        /* sum from Lyman edge to carbon k-shell */
494        sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s],opac.ipCKshell-1,&rfield.qheii,&pheii,&p4);
495
496        /* sum from c k-shell to gamma ray - where pairs start */
497        sumcon( opac.ipCKshell , rfield.ipEnerGammaRay-1 , &prt.qx,
498                &prt.xpow ,  &p5);
499
500        /* complete sum up to high-energy limit */
501        sumcon(rfield.ipEnerGammaRay,rfield.nflux,&prt.qgam,&prt.GammaLumin, &p6);
502
503        /* find to estimate photoerosion timescale */
504        n = ipoint(7.35e5);
505        sumcon(n,rfield.nflux,&qgn,&pgn,&p7);
506        timesc.TimeErode = qgn;
507
508        /* find Compton temp */
509        tcompr = (p1 + p2 + p3 + p4 + p5 + p6)/(prt.pradio + prt.pbal +
510          p + phe + pheii + prt.xpow + prt.GammaLumin);
511
512        tcomp = tcompr/(4.*6.272e-6);
513
514        if( trace.lgTrace )
515        {
516                fprintf( ioQQQ,
517                        " mean photon energy=%10.3eR =%10.3eK low, high nu=%12.4e%12.4e\n",
518                  tcompr, tcomp, rfield.anu[0] - rfield.widflx[0]/2., rfield.anu[rfield.nflux-1] +
519                  rfield.widflx[rfield.nflux-1]/2. );
520        }
521
522        /* this is total power in ionizing radiation, per unit area */
523        prt.powion = p + phe + pheii + prt.xpow + prt.GammaLumin;
524
525        /* this is the total radiation intensity, erg cm-2 s-1 */
526        continuum.TotalLumin = prt.pradio + prt.powion + prt.pbal;
527
528        /* this is placed into the line stack on the first zone, then
529         * reset to zero, to end up with luminosity in the emission lines array.
530         * at end of iteration it is reset to TotalLumin */
531        continuum.totlsv = continuum.TotalLumin;
532
533        /* total H-ionizing photon number, */
534        rfield.qhtot = prt.q + rfield.qhe + rfield.qheii + prt.qx + prt.qgam;
535
536        /* ftotal photon number, all energies  */
537        rfield.qtot = rfield.qhtot + rfield.qbal + rfield.qrad;
538
539        if( prt.powion <= 0. && called.lgTalk )
540        {
541                rfield.lgHionRad = true;
542                fprintf( ioQQQ, " NOTE There is no hydrogen-ionizing radiation.\n" );
543                fprintf( ioQQQ, " Was this intended?\n\n" );
544                /* check if any Balmer ionizing radiation since even metals will be
545                 * totally neutral */
546                if( prt.pbal <=0. && called.lgTalk  )
547                {
548                        fprintf( ioQQQ, " NOTE There is no Balmer continuum radiation.<<<<\n" );
549                        fprintf( ioQQQ, " Was this intended?\n\n" );
550                }
551        }
552
553        else
554        {
555                rfield.lgHionRad = false;
556        }
557
558        /* option to add energy deposition due to fast neutrons,
559         * entered as fraction of total photon luminosity */
560        if( hextra.lgNeutrnHeatOn )
561        {
562                /* hextra.totneu is erg cm-2 s-1 in neutrons
563                 * hextra.effneu - efficiency default is unity */
564                hextra.totneu = (realnum)(hextra.effneu*continuum.TotalLumin*
565                        pow((realnum)10.f,hextra.frcneu));
566        }
567        else
568        {
569                hextra.totneu = (realnum)0.;
570        }
571
572        /* temp correspond to energy density, printed in STARTR */
573        phycon.TEnerDen = pow(continuum.TotalLumin/SPEEDLIGHT/7.56464e-15,0.25);
574
575        /* sanity check for single blackbody, that energy density temperature
576         * is smaller than black body temperature */
577        if( rfield.nspec==1 &&
578                strcmp( rfield.chSpType[rfield.nspec-1], "BLACK" )==0 )
579        {
580                /* single black body, now confirm that TEnerDen is less than this temperature,
581                 * >>>chng 99 may 02,
582                 * in lte these are very very close, factor of 1.00001 allows for numerical
583                 * errors, and apparently slightly different atomic coef in different parts
584                 * of code.  eventaully all mustuse physonst.h and agree exactly */
585                if( phycon.TEnerDen > 1.0001f*rfield.slope[rfield.nspec-1] )
586                {
587                        fprintf( ioQQQ,
588                                "\n WARNING:  The energy density temperature (%g) is greater than the"
589                                " black body temperature (%g).  This is unphysical.\n\n",
590                                phycon.TEnerDen , rfield.slope[rfield.nspec-1]);
591                }
592        }
593
594        /* incident continuum nu*f_nu at Hbeta and Ly-alpha */
595        continuum.cn4861 = (realnum)(ffun(0.1875)*HPLANCK*FR1RYD*0.1875*0.1875);
596        continuum.cn1216 = (realnum)(ffun(0.75)*HPLANCK*FR1RYD*0.75*0.75);
597        continuum.sv4861 = continuum.cn4861;
598        continuum.sv1216 = continuum.cn1216;
599        /* flux density in V, erg / s / cm2 / hz */
600        continuum.fluxv = (realnum)(ffun(0.1643)*HPLANCK*0.1643);
601        continuum.fbeta = (realnum)(ffun(0.1875)*HPLANCK*0.1875*6.167e14);
602
603        /* flux density nu*Fnu = erg / s / cm2
604         * EX1RYD is optional extinction factor at 1 ryd */
605        prt.fx1ryd = (realnum)(ffun(1.000)*HPLANCK*ex1ryd*FR1RYD);
606
607        /* check for plasma frequency - then zero out incident continuum
608         * for energies below this
609         * this is critical electron density, shielding of incident continuum
610         * if electron density is greater than this */
611        ecrit = POW2(rfield.anu[0]/2.729e-12);
612
613        if( dense.gas_phase[ipHYDROGEN]*1.2 > ecrit )
614        {
615                rfield.lgPlasNu = true;
616                rfield.plsfrq = (realnum)(2.729e-12*sqrt(dense.gas_phase[ipHYDROGEN]*1.2));
617                rfield.plsfrqmax = rfield.plsfrq;
618                rfield.ipPlasma = ipoint(rfield.plsfrq);
619
620                /* save max pointer too */
621                rfield.ipPlasmax = rfield.ipPlasma;
622
623                /* now loop over all affected energies, setting incident continuum
624                 * to zero there, and counting all as reflected */
625                /* >>chng 01 jul 14, from i < ipPlasma to ipPlasma-1 -
626                 * when ipPlasma is 1 plasma freq is not on energy scale */
627                for( i=0; i < rfield.ipPlasma-1; i++ )
628                {
629                        /* count as reflected incident continuum */
630                        rfield.ConRefIncid[0][i] = rfield.flux[0][i];
631                        /* set continuum to zero there */
632                        rfield.flux_beam_const[i] = 0.;
633                        rfield.flux_beam_time[i] = 0.;
634                        rfield.flux_isotropic[i] = 0.;
635                        rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
636                                rfield.flux_isotropic[i];
637                }
638        }
639        else
640        {
641                rfield.lgPlasNu = false;
642                /* >>chng 01 jul 14, from 0 to 1 - 1 is the first array element on the F scale,
643                 * ipoint would return this, so rest of code assumes ipPlasma is 1 plus correct index */
644                rfield.ipPlasma = 1;
645                rfield.plsfrqmax = 0.;
646                rfield.plsfrq = 0.;
647        }
648
649        if( rfield.ipPlasma > 1 && called.lgTalk )
650        {
651                fprintf( ioQQQ,
652                        "           !The plasma frequency is %.2e Ryd.  The incident continuum is set to 0 below this.\n",
653                  rfield.plsfrq );
654        }
655
656        rfield.occmax = 0.;
657        rfield.tbrmax = 0.;
658        for( i=0; i < rfield.nupper; i++ )
659        {
660                /* set up occupation number array */
661                rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
662                if( rfield.OccNumbIncidCont[i] > rfield.occmax )
663                {
664                        rfield.occmax = rfield.OccNumbIncidCont[i];
665                        rfield.occmnu = rfield.anu[i];
666                }
667                /* following product is continuum brightness temperature */
668                if( rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i] > rfield.tbrmax )
669                {
670                        rfield.tbrmax = (realnum)(rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i]);
671                        rfield.tbrmnu = rfield.anu[i];
672                }
673                /* save continuum for next iteration */
674                rfield.flux_total_incident[0][i] = rfield.flux[0][i];
675                rfield.flux_beam_const_save[i] = rfield.flux_beam_const[i];
676                rfield.flux_time_beam_save[i] = rfield.flux_beam_time[i];
677                rfield.flux_isotropic_save[i] = rfield.flux_isotropic[i];
678                /*fprintf(ioQQQ,"DEBUG type cont %li\t%.3e\t%.2e\t%.2e\t%.2e\t%.2e\n",
679                        i, rfield.anu[i],
680                        rfield.flux[0][i],rfield.flux_beam_const[i],rfield.flux_beam_time[i],
681                        rfield.flux_isotropic[i]);
682                fflush(ioQQQ);*/
683        }
684
685        /* if continuum brightness temp is large, where does it fall below
686         * 1e4k??? */
687        if( rfield.tbrmax > 1e4 )
688        {
689                i = ipoint(rfield.tbrmnu)-1;
690                while( i < rfield.nupper-1 && (rfield.OccNumbIncidCont[i]*TE1RYD*
691                  rfield.anu[i] > 1e4) )
692                {
693                        ++i;
694                }
695                rfield.tbr4nu = rfield.anu[i];
696        }
697        else
698        {
699                rfield.tbr4nu = 0.;
700        }
701
702        /* if continuum occ num is large, where does it fall below 1? */
703        if( rfield.occmax > 1. )
704        {
705                i = ipoint(rfield.occmnu)-1;
706                while( i < rfield.nupper && (rfield.OccNumbIncidCont[i] > 1.) )
707                {
708                        ++i;
709                }
710                rfield.occ1nu = rfield.anu[i];
711        }
712        else
713        {
714                rfield.occ1nu = 0.;
715        }
716
717        /* remember if incident radiation field is less than 10*Habing ISM */
718        /* >>chng 06 aug 01, change this test from continuum.TotalLumin to
719         * energy in balmer and ionizing continuum, since this is the true habing field
720         * and is the continuum that interacts with gas.  When CMB set this
721         * tests on total did not trigger due to cold blackbody, which has little
722         * effect on gas, other than compton */
723        if( (prt.powion + prt.pbal) < 1.8e-2 )
724        {
725                /* thermal.ConstTemp def is zero, set pos when constant temperature is set */
726                rfield.lgHabing = true;
727                /* >>chng 06 aug 01 also print warning if substantially below Habing, this may be a mistake */
728                if( ((prt.powion + prt.pbal) < 1.8e-12) &&
729                        /* this is test for not constant temperature */
730                        (!thermal.lgTemperatureConstant) )
731                {
732                        fprintf( ioQQQ, "\n >>>\n"
733                                                        " >>> NOTE The incident continuum is surprisingly faint.\n" );
734                        fprintf( ioQQQ,
735                                " >>> The total energy in the Balmer and Lyman continua is %.2e erg cm-2 s-1.\n"
736                                ,(prt.powion + prt.pbal));
737                        fprintf( ioQQQ, " >>> This is many orders of magnitude fainter than the ISM galactic background.\n" );
738                        fprintf( ioQQQ, " >>> This seems unphysical - please check that the continuum intensity has been properly set.\n" );
739                        fprintf( ioQQQ, " >>> YOU MAY BE MAKING A BIG MISTAKE!!\n >>>\n\n\n\n" );
740                }
741        }
742
743        /* fix ionization parameter (per hydrogen) at inner edge */
744        rfield.uh = (realnum)(rfield.qhtot/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT);
745        rfield.uheii = (realnum)((rfield.qheii + prt.qx)/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT);
746        if( rfield.uh > 1.e10 )
747        {
748                fprintf( ioQQQ, "\n >>>\n"
749                                                " NOTE The incident continuum is surprisingly intense.\n" );
750                fprintf( ioQQQ,
751                        " >>> The hydrogen ionization parameter is %.2e [dimensionless].\n"
752                        , rfield.uh );
753                fprintf( ioQQQ, " This is many orders of magnitude brighter than commonly seen.\n" );
754                fprintf( ioQQQ, " This seems unphysical - please check that the continuum intensity has been properly set.\n" );
755                fprintf( ioQQQ, " YOU MAY BE MAKING A BIG MISTAKE!!\n >>>\n\n\n\n" );
756        }
757
758        /* guess first temperature and neutral h density */
759        if( thermal.ConstTemp > 0. )
760        {
761                phycon.te = thermal.ConstTemp;
762        }
763        else
764        {
765                if( rfield.uh > 0. )
766                {
767                        phycon.te = (20000.+log10(rfield.uh)*5000.);
768                        phycon.te = MAX2(8000. , phycon.te );
769                }
770                else
771                {
772                        phycon.te = 1000.;
773                }
774        }
775
776        /* this is an option to stop after printing header only */
777        if( noexec.lgNoExec )
778                return(0);
779
780        /* estimate secondary ionization rate - probably 0, but possible extra
781         * SetCsupra set with "set csupra" */
782        /* >>>chng 99 apr 29, added cosmic ray ionization since this is used to get
783         * helium ionization fraction, and was zero in pdr, so He turned off at start,
784         * and never turned back on */
785        /* coef on cryden is from highen.c */
786        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
787        {
788                for( ion=0; ion<nelem+1; ++ion )
789                {
790                        secondaries.csupra[nelem][ion] =
791                                secondaries.SetCsupra + hextra.cryden*2e-9f;
792                }
793        }
794
795        /*********************************************************************
796         *                                                                   *
797         * estimate hydrogen's level of ionization                           *
798         *                                                                   *
799         *********************************************************************/
800
801        /* create fake ionization balance, but will conserve number of hydrogens */
802        dense.xIonDense[ipHYDROGEN][0] = 0.;
803        dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN];
804        /* this must be zero since PresTotCurrent will do radiation pressure due to H */
805        StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop = 0.;
806
807        /* "extra" electrons from command line, or assumed residual electrons */
808        double EdenExtraLocal = dense.EdenExtra +
809                /* if we are in a molecular cloud the current logic could badly fail
810                * do not let electron density fall below 1e-7 of H density */
811                1e-7*dense.gas_phase[ipHYDROGEN];
812        dense.eden = dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal;
813
814        /* hydrogen case B recombination coefficient */
815        HCaseBRecCoeff = (-9.9765209 + 0.158607055*phycon.telogn[0] + 0.30112749*
816          phycon.telogn[1] - 0.063969007*phycon.telogn[2] + 0.0012691546*
817          phycon.telogn[3])/(1. + 0.035055422*phycon.telogn[0] -
818          0.037621619*phycon.telogn[1] + 0.0076175048*phycon.telogn[2] -
819          0.00023432613*phycon.telogn[3]);
820        HCaseBRecCoeff = pow(10.,HCaseBRecCoeff)/phycon.te;
821
822        double CollIoniz = t_ADfA::Inst().coll_ion(1,1,phycon.te);
823        double OtherIonization = rfield.qhtot*2e-18 + secondaries.csupra[ipHYDROGEN][0];
824        double newEden = dense.eden;
825        long loopCount = 0;
826
827        do
828        {
829                /* update electron density */
830                dense.eden = newEden;
831                double RatioIoniz =
832                        (CollIoniz*dense.eden+OtherIonization)/(HCaseBRecCoeff*dense.eden);
833                if( RatioIoniz<1e-3 )
834                {
835                        /* very low ionization solution */
836                        dense.xIonDense[ipHYDROGEN][1] = (realnum)(
837                                dense.gas_phase[ipHYDROGEN]*RatioIoniz);
838                        dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
839                                dense.xIonDense[ipHYDROGEN][1];
840                        ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
841                                dense.xIonDense[ipHYDROGEN][0]<=dense.gas_phase[ipHYDROGEN]);
842                        ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
843                                dense.xIonDense[ipHYDROGEN][1]<dense.gas_phase[ipHYDROGEN]);
844                }
845                else if( RatioIoniz>1e3 )
846                {
847                        /* very high ionization solution */
848                        dense.xIonDense[ipHYDROGEN][0] = (realnum)(
849                                dense.gas_phase[ipHYDROGEN]/RatioIoniz);
850                        dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN] -
851                                dense.xIonDense[ipHYDROGEN][0];
852                        ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
853                                dense.xIonDense[ipHYDROGEN][0]<dense.gas_phase[ipHYDROGEN]);
854                        ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
855                                dense.xIonDense[ipHYDROGEN][1]<=dense.gas_phase[ipHYDROGEN]);
856                }
857                else
858                {
859                        /* intermediate ionization - solve quadratic */
860                        double alpha = HCaseBRecCoeff + CollIoniz ;
861                        double beta = HCaseBRecCoeff*EdenExtraLocal + OtherIonization +
862                                (EdenExtraLocal - dense.gas_phase[ipHYDROGEN])*CollIoniz;
863                        double gamma = -dense.gas_phase[ipHYDROGEN]*(OtherIonization+EdenExtraLocal*CollIoniz);
864
865                        double discriminant = POW2(beta) - 4.*alpha*gamma;
866                        if( discriminant <0 )
867                        {
868                                /* oops */
869                                fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found negative discriminant.\n");
870                                TotalInsanity();
871                        }
872
873                        dense.xIonDense[ipHYDROGEN][1] = (realnum)((-beta + sqrt(discriminant))/(2.*alpha));
874                        if( dense.xIonDense[ipHYDROGEN][1]> dense.gas_phase[ipHYDROGEN] )
875                        {
876                                /* oops */
877                                fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H+)>n(H).\n");
878                                TotalInsanity();
879                        }
880                        dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
881                                dense.xIonDense[ipHYDROGEN][1];
882                        if( dense.xIonDense[ipHYDROGEN][0]<= 0 )
883                        {
884                                /* oops */
885                                fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H0)<0.\n");
886                                TotalInsanity();
887                        }
888                }
889
890
891                if( dense.xIonDense[ipHYDROGEN][1] > 1e-30 )
892                {
893                        StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop = dense.xIonDense[ipHYDROGEN][0]/dense.xIonDense[ipHYDROGEN][1];
894                }
895                else
896                {
897                        StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop = 0.;
898                }
899
900                /* now save estimates of whether induced recombination is going
901                 * to be important -this is a code pacesetter since GammaBn is slower
902                 * than GammaK */
903                hydro.lgHInducImp = false;
904                for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; i++ )
905                {
906                        if( rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i]-1] > 0.01 )
907                                hydro.lgHInducImp = true;
908                }
909
910                /*******************************************************************
911                 *                                                                 *
912                 * estimate helium's level of ionization                           *
913                 *                                                                 *
914                 *******************************************************************/
915
916                /* only if helium is turned on */
917                if( dense.lgElmtOn[ipHELIUM] )
918                {
919                        /* next estimate level of helium singly ionized */
920                        xIoniz = (realnum)t_ADfA::Inst().coll_ion(2,2,phycon.te);
921                        /* >>chng 05 jan 05, add cosmic rays */
922                        xIoniz = (realnum)(xIoniz*dense.eden + rfield.qhe*1e-18 +  secondaries.csupra[ipHELIUM][1] );
923                        double RecTot = HCaseBRecCoeff*dense.eden;
924                        RatioIonizToRecomb = xIoniz/RecTot;
925
926                        /* now estimate level of helium doubly ionized */
927                        xIoniz = (realnum)t_ADfA::Inst().coll_ion(2,1,phycon.te);
928                        /* >>chng 05 jan 05, add cosmic rays */
929                        xIoniz = (realnum)(xIoniz*dense.eden + rfield.qheii*1e-18 +  ionbal.CosRayIonRate );
930
931                        /* rough charge dependence */
932                        RecTot *= 4.;
933                        r3ov2 = xIoniz/RecTot;
934
935                        /* now set level of helium ionization */
936                        if( RatioIonizToRecomb > 0. )
937                        {
938                                dense.xIonDense[ipHELIUM][1] = (realnum)(dense.gas_phase[ipHELIUM]/(1./RatioIonizToRecomb + 1. + r3ov2));
939                                dense.xIonDense[ipHELIUM][0] = (realnum)(dense.xIonDense[ipHELIUM][1]/RatioIonizToRecomb);
940                        }
941                        else
942                        {
943                                /* no He ionizing radiation */
944                                dense.xIonDense[ipHELIUM][1] = 0.;
945                                dense.xIonDense[ipHELIUM][0] = dense.gas_phase[ipHELIUM];
946                        }
947
948                        dense.xIonDense[ipHELIUM][2] = (realnum)(dense.xIonDense[ipHELIUM][1]*r3ov2);
949
950                        if( dense.xIonDense[ipHELIUM][2] > 1e-30 )
951                        {
952                                StatesElem[ipH_LIKE][1][ipH1s].Pop = dense.xIonDense[ipHELIUM][1]/dense.xIonDense[ipHELIUM][2];
953                        }
954