00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "cddefines.h"
00011 #include "physconst.h"
00012 #include "iso.h"
00013 #include "noexec.h"
00014 #include "ionbal.h"
00015 #include "hextra.h"
00016 #include "trace.h"
00017 #include "dense.h"
00018 #include "oxy.h"
00019 #include "prt.h"
00020 #include "heavy.h"
00021 #include "rfield.h"
00022 #include "phycon.h"
00023 #include "called.h"
00024 #include "hydrogenic.h"
00025 #include "timesc.h"
00026 #include "secondaries.h"
00027 #include "opacity.h"
00028 #include "thermal.h"
00029 #include "ipoint.h"
00030 #include "atmdat.h"
00031 #include "rt.h"
00032 #include "radius.h"
00033 #include "geometry.h"
00034 #include "grainvar.h"
00035 #include "continuum.h"
00036 #include "taulines.h"
00037
00038
00039 static const double aweigh[4]={-0.4305682,-0.1699905, 0.1699905, 0.4305682};
00040 static const double fweigh[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274};
00041
00042
00043 STATIC void conorm(void);
00044
00045
00046 STATIC double pintr(double penlo,
00047 double penhi);
00048
00049
00050 STATIC double qintr(double *qenlo,
00051 double *qenhi);
00052
00053
00054
00055 STATIC void sumcon(long int il,
00056 long int ih,
00057 realnum *q,
00058 realnum *p,
00059 realnum *panu);
00060
00061
00062 STATIC void extin(realnum *ex1ryd);
00063
00064
00065 STATIC void ptrcer(void);
00066
00067 void IncidentContinuumHere( void )
00068 {
00069
00070 DEBUG_ENTRY( "IncidentContinuumHere()" );
00071 double frac_beam_time;
00072
00073 double frac_beam_const;
00074
00075 double frac_isotropic;
00076
00077 double BigLog = 0.;
00078 for( long i = 0; i<rfield.nflux; ++i )
00079 {
00080 double flux_org = rfield.flux[0][i];
00081 double flux_now = ffun( rfield.anu[i] ,
00082 &frac_beam_time ,
00083 &frac_beam_const ,
00084 &frac_isotropic )*rfield.widflx[i]*rfield.ExtinguishFactor[i];
00085
00086 double ratio = 1.;
00087 if( flux_org>SMALLFLOAT && flux_now>SMALLFLOAT )
00088 {
00089 ratio = fabs( log10( flux_now / flux_org ) );
00090 BigLog = max( ratio , BigLog );
00091 }
00092 }
00093
00094 if( BigLog > 0.01 )
00095 fprintf(ioQQQ , "DEBUG diff continua %.2e\n", BigLog );
00096 return;
00097 }
00098
00099
00100 int ContSetIntensity(void)
00101 {
00102 bool lgCheckOK;
00103
00104 long int i,
00105 ip,
00106 j,
00107 n;
00108
00109 realnum EdenHeav,
00110 ex1ryd,
00111 factor,
00112 occ1,
00113 p,
00114 p1,
00115 p2,
00116 p3,
00117 p4,
00118 p5,
00119 p6,
00120 p7,
00121 p8,
00122 pgn,
00123 phe,
00124 pheii,
00125 qgn;
00126
00127 realnum xIoniz;
00128
00129 double HCaseBRecCoeff,
00130 wanu[4],
00131 alf,
00132 bet,
00133 fntest,
00134 fsum,
00135 ecrit,
00136 tcompr,
00137 tcomp,
00138 RatioIonizToRecomb,
00139 r3ov2;
00140
00141 double amean,
00142 amean2,
00143 amean3,
00144 peak,
00145 wfun[4];
00146
00147
00148 double frac_beam_time , frac_beam_time1;
00149
00150 double frac_beam_const , frac_beam_const1;
00151
00152 double frac_isotropic , frac_isotropic1;
00153
00154 long int nelem , ion;
00155
00156 DEBUG_ENTRY( "ContSetIntensity()" );
00157
00158
00159 if( trace.lgTrace )
00160 {
00161 fprintf( ioQQQ, " ContSetIntensity called.\n" );
00162 }
00163
00164
00165
00166 conorm();
00167
00168
00169
00170 factor = (realnum)(EN1RYD/PI4/FR1RYD/HNU3C2);
00171
00172
00173 lgCheckOK = true;
00174
00175 for( i=0; i < rfield.nupper; i++ )
00176 {
00177
00178 rfield.anu[i] = rfield.AnuOrg[i];
00179 rfield.ContBoltz[i] = 0.;
00180 fsum = 0.;
00181 amean = 0.;
00182 amean2 = 0.;
00183 amean3 = 0.;
00184 frac_beam_time = 0.;
00185 frac_beam_const = 0.;
00186 frac_isotropic = 0.;
00187
00188 for( j=0; j < 4; j++ )
00189 {
00190
00191 wanu[j] = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
00192
00193
00194 wanu[j] = MAX2( wanu[j] , rfield.emm );
00195 wanu[j] = MIN2( wanu[j] , rfield.egamry );
00196
00197
00198
00199 if( i > 0 && i < rfield.nupper-1 )
00200 {
00201 wanu[j] = MAX2( wanu[j] , rfield.anu[i-1] + 0.5*rfield.widflx[i-1] );
00202 wanu[j] = MIN2( wanu[j] , rfield.anu[i+1] - 0.5*rfield.widflx[i+1] );
00203 }
00204
00205 wfun[j] = fweigh[j]*ffun(wanu[j] ,
00206 &frac_beam_time1 ,
00207 &frac_beam_const1 ,
00208 &frac_isotropic1 );
00209
00210
00211
00212 fsum += wfun[j];
00213 amean += wanu[j]*wfun[j];
00214 amean2 += wanu[j]*wanu[j]*wfun[j];
00215 amean3 += wanu[j]*wanu[j]*wanu[j]*wfun[j];
00216 frac_beam_time += fweigh[j]*frac_beam_time1;
00217 frac_beam_const += fweigh[j]*frac_beam_const1;
00218 frac_isotropic += fweigh[j]*frac_isotropic1;
00219 }
00220
00221 ASSERT( fabs( 1.-frac_beam_time-frac_beam_const-frac_isotropic)<
00222 10.*FLT_EPSILON);
00223
00224 if( fsum*rfield.widflx[i] > BIGFLOAT )
00225 {
00226 fprintf( ioQQQ, "\n Cannot continue. The continuum is far too intense.\n" );
00227 for( j=0; j < rfield.nShape; j++ )
00228 {
00229 if( (wfun[j]*rfield.widflx[i] > BIGFLOAT) && ( rfield.nShape > 1 ) )
00230 {
00231 fprintf( ioQQQ, " Problem is with source number %li\n", j );
00232 break;
00233 }
00234 }
00235 cdEXIT(EXIT_FAILURE);
00236 }
00237
00238 rfield.flux[0][i] = (realnum)(fsum*rfield.widflx[i]);
00239
00240
00241
00242 rfield.flux_beam_const[i] = rfield.flux[0][i] * (realnum)frac_beam_const;
00243 rfield.flux_beam_time[i] = rfield.flux[0][i] * (realnum)frac_beam_time;
00244 rfield.flux_isotropic[i] = rfield.flux[0][i] * (realnum)frac_isotropic;
00245
00246 if( rfield.flux[0][i] > 0. )
00247 {
00248
00249
00250 rfield.anu[i] = (realnum)(amean2/amean);
00251 rfield.anu2[i] = (realnum)(amean3/amean);
00252
00253 if( i > 0 && rfield.anu[i] <= rfield.anu[i-1] )
00254 {
00255
00256
00257
00258 rfield.anu[i] = rfield.anu[i-1]*(1.f+2.f*FLT_EPSILON);
00259 rfield.anu2[i] = pow2(rfield.anu[i]);
00260 }
00261 ASSERT( i==0 || rfield.anu[i] > rfield.anu[i-1] );
00262
00263
00264
00265 rfield.anulog[i] = (realnum)log10(rfield.anu[i]);
00266 }
00267
00268 else if( rfield.flux[0][i] == 0. )
00269 {
00270 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
00271 rfield.anulog[i] = (realnum)log10(rfield.anu[i]);
00272 }
00273
00274 else
00275 {
00276 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
00277 fprintf( ioQQQ, " negative continuum returned at%6ld%10.2e%10.2e\n",
00278 i, rfield.anu[i], rfield.flux[0][i] );
00279 lgCheckOK = false;
00280 }
00281 rfield.anu3[i] = rfield.anu2[i]*rfield.anu[i];
00282
00283 rfield.ConEmitReflec[0][i] = 0.;
00284 rfield.ConEmitOut[0][i] = 0.;
00285 rfield.convoc[i] = factor/rfield.widflx[i]/rfield.anu2[i];
00286
00287
00288 alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i]));
00289 bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/
00290 4.;
00291 rfield.csigh[i] = (realnum)(alf*rfield.anu2[i]*3.858e-25);
00292 rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25);
00293 }
00294
00295
00296
00297
00298
00299
00300
00301 #if 0
00302
00303
00304 for( i=1; i<rfield.nupper-1; ++i )
00305 {
00306
00307 rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
00308 }
00309 #endif
00310
00311 if( !lgCheckOK )
00312 {
00313 ShowMe();
00314 cdEXIT(EXIT_FAILURE);
00315 }
00316
00317 if( trace.lgTrace && trace.lgComBug )
00318 {
00319 fprintf( ioQQQ, "\n\n Compton heating, cooling coefficients \n" );
00320 for( i=0; i < rfield.nupper; i += 2 )
00321 {
00322 fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i],
00323 rfield.csigh[i], rfield.csigc[i] );
00324 }
00325 fprintf( ioQQQ, "\n" );
00326 }
00327
00328
00329
00330 if( trace.lgPtrace )
00331 ptrcer();
00332
00333
00334 extin(&ex1ryd);
00335
00336
00337
00338 prt.ipeak = 1;
00339 peak = 0.;
00340
00341 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
00342 {
00343 if( rfield.flux[0][i]*rfield.anu[i]/rfield.widflx[i] > (realnum)peak )
00344 {
00345
00346
00347
00348 prt.ipeak = i+1;
00349 peak = rfield.flux[0][i]*rfield.anu[i]/rfield.widflx[i];
00350 }
00351 }
00352
00353
00354
00355 peak = rfield.flux[0][prt.ipeak-1]/rfield.widflx[prt.ipeak-1]*
00356 rfield.anu2[prt.ipeak-1];
00357
00358
00359 if( trace.lgTrace )
00360 {
00361 fprintf( ioQQQ, " ContSetIntensity: The peak of the H-ion continuum is at %.3e Ryd - its value is %.2e\n",
00362 rfield.anu[prt.ipeak-1] , peak);
00363 }
00364
00365 if( peak > 1e38 )
00366 {
00367 fprintf( ioQQQ, " PROBLEM DISASTER The continuum is too intense to compute. Use a fainter continuum. (This is the nu*f_nu test)\n" );
00368 fprintf( ioQQQ, " Sorry.\n" );
00369 cdEXIT(EXIT_FAILURE);
00370 }
00371
00372
00373
00374
00375
00376 fntest = peak*rfield.FluxFaint;
00377 {
00378 enum {DEBUG_LOC=false};
00379
00380 if( DEBUG_LOC )
00381 {
00382 for( i=0; i<rfield.nupper; ++i )
00383 {
00384 fprintf(ioQQQ," consetintensityBUGGG\t%.2e\t%.2e\n" ,
00385 rfield.anu[i] , rfield.flux[0][i]/rfield.widflx[i] );
00386 }
00387 cdEXIT(EXIT_SUCCESS);
00388 }
00389 }
00390
00391 if( fntest > 0. )
00392 {
00393
00394
00395 i = rfield.nupper;
00396
00397
00398 while( i > prt.ipeak+3 &&
00399 rfield.flux[0][i-1]*rfield.anu2[i-1]/rfield.widflx[i-1] < (realnum)fntest )
00400 {
00401 --i;
00402 }
00403 }
00404 else
00405 {
00406
00407
00408
00409 i = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]+3;
00410 }
00411
00412
00413
00414
00415
00416
00417
00418
00419 rfield.nflux = i;
00420
00421
00422
00423
00424 while( rfield.flux[0][rfield.nflux-1] < SMALLFLOAT && rfield.nflux > 1 )
00425 {
00426 --rfield.nflux;
00427 }
00428
00429 ++rfield.nflux;
00430
00431 if( rfield.nflux == 1 )
00432 {
00433 fprintf( ioQQQ, " PROBLEM DISASTER This incident continuum appears to have no radiation.\n" );
00434 fprintf( ioQQQ, " Sorry.\n" );
00435 cdEXIT(EXIT_FAILURE);
00436 }
00437
00438
00439
00440 rfield.nflux = MIN2( rfield.nupper-1 , rfield.nflux );
00441
00442
00443
00444 rfield.nfine = rfield.nfine_malloc;
00445 while( rfield.nfine > 0 && rfield.fine_anu[rfield.nfine-1] > rfield.anu[rfield.nflux-1] )
00446 {
00447 --rfield.nfine;
00448 }
00449
00450
00451 continuum.lgCon0 = false;
00452 ip = 0;
00453 for( i=0; i < rfield.nflux; i++ )
00454 {
00455 if( rfield.flux[0][i] == 0. )
00456 {
00457 if( called.lgTalk && !continuum.lgCon0 )
00458 {
00459 fprintf( ioQQQ, " NOTE Setcon: continuum has zero intensity starting at %11.4e Ryd.\n",
00460 rfield.anu[i] );
00461 continuum.lgCon0 = true;
00462 }
00463 ++ip;
00464 }
00465 }
00466
00467 if( continuum.lgCon0 && called.lgTalk )
00468 {
00469 fprintf( ioQQQ,
00470 "%6ld cells in the incident continuum have zero intensity. Problems???\n\n",
00471 ip );
00472 }
00473
00474
00475 lgCheckOK = true;
00476 for( i=1; i < rfield.nflux; i++ )
00477 {
00478 if( rfield.flux[0][i] < 0. )
00479 {
00480 fprintf( ioQQQ,
00481 " PROBLEM DISASTER Continuum has negative intensity at %.4e Ryd=%.2e %4.4s %4.4s\n",
00482 rfield.anu[i], rfield.flux[0][i], rfield.chLineLabel[i], rfield.chContLabel[i] );
00483 lgCheckOK = false;
00484 }
00485 else if( rfield.anu[i] <= rfield.anu[i-1] )
00486 {
00487 fprintf( ioQQQ,
00488 " PROBLEM DISASTER cont_setintensity - internal error - continuum energies not in increasing order: energies follow\n" );
00489 fprintf( ioQQQ,
00490 "%ld %e %ld %e %ld %e\n",
00491 i -1 , rfield.anu[i-1], i, rfield.anu[i], i +1, rfield.anu[i+1] );
00492 lgCheckOK = false;
00493 }
00494 }
00495
00496
00497 if( !lgCheckOK )
00498 {
00499 ShowMe();
00500 cdEXIT(EXIT_FAILURE);
00501 }
00502
00503
00504 if( rfield.anu[rfield.nflux-1] <= 190 )
00505 {
00506 ionbal.lgCompRecoil = false;
00507 }
00508
00509
00510
00511
00512 sumcon(1,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,&rfield.qrad,&prt.pradio,&p1);
00513
00514
00515 sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p],iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,&rfield.qbal,&prt.pbal,&p2);
00516
00517
00518 sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
00519 iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1,&prt.q,&p,&p3);
00520
00521
00522 sumcon(iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0],
00523 iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1,&rfield.qhe,&phe,&p4);
00524
00525
00526 sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s],opac.ipCKshell-1,&rfield.qheii,&pheii,&p5);
00527
00528
00529 sumcon( opac.ipCKshell , rfield.ipEnerGammaRay-1 , &prt.qx,
00530 &prt.xpow , &p6);
00531
00532
00533 sumcon(rfield.ipEnerGammaRay,rfield.nflux,&prt.qgam,&prt.GammaLumin, &p7);
00534
00535
00536 n = ipoint(7.35e5);
00537 sumcon(n,rfield.nflux,&qgn,&pgn,&p8);
00538 timesc.TimeErode = qgn;
00539
00540
00541 tcompr = (p1 + p2 + p3 + p4 + p5 + p6 + p7)/(prt.pradio + prt.pbal +
00542 p + phe + pheii + prt.xpow + prt.GammaLumin);
00543
00544 tcomp = tcompr/(4.*6.272e-6);
00545
00546 if( trace.lgTrace )
00547 {
00548 fprintf( ioQQQ,
00549 " mean photon energy=%10.3eR =%10.3eK low, high nu=%12.4e%12.4e\n",
00550 tcompr, tcomp, rfield.anu[0] - rfield.widflx[0]/2., rfield.anu[rfield.nflux-1] +
00551 rfield.widflx[rfield.nflux-1]/2. );
00552 }
00553
00554
00555 prt.powion = p + phe + pheii + prt.xpow + prt.GammaLumin;
00556
00557
00558 continuum.TotalLumin = prt.pradio + prt.powion + prt.pbal;
00559
00560
00561
00562
00563 continuum.totlsv = continuum.TotalLumin;
00564
00565
00566 rfield.qhtot = prt.q + rfield.qhe + rfield.qheii + prt.qx + prt.qgam;
00567
00568
00569 rfield.qtot = rfield.qhtot + rfield.qbal + rfield.qrad;
00570
00571 if( prt.powion <= 0. && called.lgTalk )
00572 {
00573 rfield.lgHionRad = true;
00574 fprintf( ioQQQ, " NOTE There is no hydrogen-ionizing radiation.\n" );
00575 fprintf( ioQQQ, " Was this intended?\n\n" );
00576
00577
00578 if( prt.pbal <=0. && called.lgTalk )
00579 {
00580 fprintf( ioQQQ, " NOTE There is no Balmer continuum radiation.<<<<\n" );
00581 fprintf( ioQQQ, " Was this intended?\n\n" );
00582 }
00583 }
00584
00585 else
00586 {
00587 rfield.lgHionRad = false;
00588 }
00589
00590
00591
00592 if( hextra.lgNeutrnHeatOn )
00593 {
00594
00595
00596 hextra.totneu = (realnum)(hextra.effneu*continuum.TotalLumin*
00597 pow((realnum)10.f,hextra.frcneu));
00598 }
00599 else
00600 {
00601 hextra.totneu = (realnum)0.;
00602 }
00603
00604
00605 phycon.TEnerDen = pow(continuum.TotalLumin/SPEEDLIGHT/7.56464e-15,0.25);
00606
00607
00608
00609 if( rfield.nShape==1 &&
00610 strcmp( rfield.chSpType[rfield.nShape-1], "BLACK" )==0 )
00611 {
00612
00613
00614
00615
00616
00617 if( phycon.TEnerDen > 1.0001f*rfield.slope[rfield.nShape-1] )
00618 {
00619 fprintf( ioQQQ,
00620 "\n WARNING: The energy density temperature (%g) is greater than the"
00621 " black body temperature (%g). This is unphysical.\n\n",
00622 phycon.TEnerDen , rfield.slope[rfield.nShape-1]);
00623 }
00624 }
00625
00626
00627 continuum.cn4861 = (realnum)(ffun(0.1875)*HPLANCK*FR1RYD*0.1875*0.1875);
00628 continuum.cn1216 = (realnum)(ffun(0.75)*HPLANCK*FR1RYD*0.75*0.75);
00629 continuum.sv4861 = continuum.cn4861;
00630 continuum.sv1216 = continuum.cn1216;
00631
00632 continuum.fluxv = (realnum)(ffun(0.1643)*HPLANCK*0.1643);
00633 continuum.fbeta = (realnum)(ffun(0.1875)*HPLANCK*0.1875*6.167e14);
00634
00635
00636
00637 prt.fx1ryd = (realnum)(ffun(1.000)*HPLANCK*ex1ryd*FR1RYD);
00638
00639 realnum plsFrqConstant = (realnum)(ELEM_CHARGE_ESU/sqrt(PI*ELECTRON_MASS)/FR1RYD);
00640 ASSERT( plsFrqConstant > 2.7e-12 && plsFrqConstant < 2.8e-12 );
00641
00642
00643
00644
00645
00646 ecrit = POW2(rfield.anu[0]/plsFrqConstant);
00647
00648 if( dense.gas_phase[ipHYDROGEN]*1.2 > ecrit )
00649 {
00650 rfield.lgPlasNu = true;
00651
00652
00653 rfield.plsfrq = plsFrqConstant*sqrt(dense.gas_phase[ipHYDROGEN]*1.2f);
00654 rfield.plsfrqmax = rfield.plsfrq;
00655 rfield.ipPlasma = ipoint(rfield.plsfrq);
00656
00657
00658 rfield.ipPlasmax = rfield.ipPlasma;
00659
00660
00661
00662
00663
00664 for( i=0; i < rfield.ipPlasma-1; i++ )
00665 {
00666
00667 rfield.ConRefIncid[0][i] = rfield.flux[0][i];
00668
00669 rfield.flux_beam_const[i] = 0.;
00670 rfield.flux_beam_time[i] = 0.;
00671 rfield.flux_isotropic[i] = 0.;
00672 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
00673 rfield.flux_isotropic[i];
00674 }
00675 }
00676 else
00677 {
00678 rfield.lgPlasNu = false;
00679
00680
00681 rfield.ipPlasma = 1;
00682 rfield.plsfrqmax = 0.;
00683 rfield.plsfrq = 0.;
00684 }
00685
00686 if( rfield.ipPlasma > 1 && called.lgTalk )
00687 {
00688 fprintf( ioQQQ,
00689 " !The plasma frequency is %.2e Ryd. The incident continuum is set to 0 below this.\n",
00690 rfield.plsfrq );
00691 }
00692
00693 rfield.occmax = 0.;
00694 rfield.tbrmax = 0.;
00695 for( i=0; i < rfield.nupper; i++ )
00696 {
00697
00698 rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
00699 if( rfield.OccNumbIncidCont[i] > rfield.occmax )
00700 {
00701 rfield.occmax = rfield.OccNumbIncidCont[i];
00702 rfield.occmnu = rfield.anu[i];
00703 }
00704
00705 if( rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i] > rfield.tbrmax )
00706 {
00707 rfield.tbrmax = (realnum)(rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i]);
00708 rfield.tbrmnu = rfield.anu[i];
00709 }
00710
00711 rfield.flux_total_incident[0][i] = rfield.flux[0][i];
00712 rfield.flux_beam_const_save[i] = rfield.flux_beam_const[i];
00713 rfield.flux_time_beam_save[i] = rfield.flux_beam_time[i];
00714 rfield.flux_isotropic_save[i] = rfield.flux_isotropic[i];
00715
00716
00717
00718
00719
00720 }
00721
00722
00723
00724 if( rfield.tbrmax > 1e4 )
00725 {
00726 i = ipoint(rfield.tbrmnu)-1;
00727 while( i < rfield.nupper-1 && (rfield.OccNumbIncidCont[i]*TE1RYD*
00728 rfield.anu[i] > 1e4) )
00729 {
00730 ++i;
00731 }
00732 rfield.tbr4nu = rfield.anu[i];
00733 }
00734 else
00735 {
00736 rfield.tbr4nu = 0.;
00737 }
00738
00739
00740 if( rfield.occmax > 1. )
00741 {
00742 i = ipoint(rfield.occmnu)-1;
00743 while( i < rfield.nupper && (rfield.OccNumbIncidCont[i] > 1.) )
00744 {
00745 ++i;
00746 }
00747 rfield.occ1nu = rfield.anu[i];
00748 }
00749 else
00750 {
00751 rfield.occ1nu = 0.;
00752 }
00753
00754
00755
00756
00757
00758
00759
00760 if( (prt.powion + prt.pbal) < 1.8e-2 )
00761 {
00762
00763 rfield.lgHabing = true;
00764
00765 if( ((prt.powion + prt.pbal) < 1.8e-12) &&
00766
00767 (!thermal.lgTemperatureConstant) )
00768 {
00769 fprintf( ioQQQ, "\n >>>\n"
00770 " >>> NOTE The incident continuum is surprisingly faint.\n" );
00771 fprintf( ioQQQ,
00772 " >>> The total energy in the Balmer and Lyman continua is %.2e erg cm-2 s-1.\n"
00773 ,(prt.powion + prt.pbal));
00774 fprintf( ioQQQ, " >>> This is many orders of magnitude fainter than the ISM galactic background.\n" );
00775 fprintf( ioQQQ, " >>> This seems unphysical - please check that the continuum intensity has been properly set.\n" );
00776 fprintf( ioQQQ, " >>> YOU MAY BE MAKING A BIG MISTAKE!!\n >>>\n\n\n\n" );
00777 }
00778 }
00779
00780
00781 rfield.uh = (realnum)(rfield.qhtot/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT);
00782 rfield.uheii = (realnum)((rfield.qheii + prt.qx)/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT);
00783 if( rfield.uh > 1e10 )
00784 {
00785 fprintf( ioQQQ, "\n\n"
00786 " CAUTION The incident radiation field is surprisingly intense.\n" );
00787 fprintf( ioQQQ,
00788 " The dimensionless hydrogen ionization parameter is %.2e.\n"
00789 , rfield.uh );
00790 fprintf( ioQQQ, " This is many orders of magnitude brighter than commonly seen.\n" );
00791 fprintf( ioQQQ, " This seems unphysical - please check that the radiation field intensity has been properly set.\n" );
00792 fprintf( ioQQQ, " YOU MAY BE MAKING A BIG MISTAKE!!\n\n\n\n\n" );
00793 }
00794
00795
00796 if( thermal.ConstTemp > 0. )
00797 {
00798 phycon.te = thermal.ConstTemp;
00799 }
00800 else
00801 {
00802 if( rfield.uh > 0. )
00803 {
00804 phycon.te = (20000.+log10(rfield.uh)*5000.);
00805 phycon.te = MAX2(8000. , phycon.te );
00806 }
00807 else
00808 {
00809 phycon.te = 1000.;
00810 }
00811 }
00812
00813
00814 if( noexec.lgNoExec )
00815 return(0);
00816
00817
00818
00819
00820
00821
00822
00823 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00824 {
00825 for( ion=0; ion<nelem+1; ++ion )
00826 {
00827 secondaries.csupra[nelem][ion] =
00828 secondaries.SetCsupra + hextra.cryden*2e-9f;
00829 }
00830 }
00831
00832
00833
00834
00835
00836
00837
00838
00839 dense.xIonDense[ipHYDROGEN][0] = 0.;
00840 dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN];
00841
00842 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop = 0.;
00843
00844
00845 double EdenExtraLocal = dense.EdenExtra +
00846
00847
00848 1e-7*dense.gas_phase[ipHYDROGEN];
00849 dense.eden = dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal;
00850
00851
00852 HCaseBRecCoeff = (-9.9765209 + 0.158607055*phycon.telogn[0] + 0.30112749*
00853 phycon.telogn[1] - 0.063969007*phycon.telogn[2] + 0.0012691546*
00854 phycon.telogn[3])/(1. + 0.035055422*phycon.telogn[0] -
00855 0.037621619*phycon.telogn[1] + 0.0076175048*phycon.telogn[2] -
00856 0.00023432613*phycon.telogn[3]);
00857 HCaseBRecCoeff = pow(10.,HCaseBRecCoeff)/phycon.te;
00858
00859 double CollIoniz = t_ADfA::Inst().coll_ion(1,1,phycon.te);
00860 double OtherIonization = rfield.qhtot*2e-18 + secondaries.csupra[ipHYDROGEN][0];
00861 double newEden = dense.eden;
00862 long loopCount = 0;
00863
00864 do
00865 {
00866
00867 dense.eden = newEden;
00868 double RatioIoniz =
00869 (CollIoniz*dense.eden+OtherIonization)/(HCaseBRecCoeff*dense.eden);
00870 if( RatioIoniz<1e-3 )
00871 {
00872
00873 dense.xIonDense[ipHYDROGEN][1] = (realnum)(
00874 dense.gas_phase[ipHYDROGEN]*RatioIoniz);
00875 dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
00876 dense.xIonDense[ipHYDROGEN][1];
00877 ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
00878 dense.xIonDense[ipHYDROGEN][0]<=dense.gas_phase[ipHYDROGEN]);
00879 ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
00880 dense.xIonDense[ipHYDROGEN][1]<dense.gas_phase[ipHYDROGEN]);
00881 }
00882 else if( RatioIoniz>1e3 )
00883 {
00884
00885 dense.xIonDense[ipHYDROGEN][0] = (realnum)(
00886 dense.gas_phase[ipHYDROGEN]/RatioIoniz);
00887 dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN] -
00888 dense.xIonDense[ipHYDROGEN][0];
00889 ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. &&
00890 dense.xIonDense[ipHYDROGEN][0]<dense.gas_phase[ipHYDROGEN]);
00891 ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. &&
00892 dense.xIonDense[ipHYDROGEN][1]<=dense.gas_phase[ipHYDROGEN]);
00893 }
00894 else
00895 {
00896
00897 double alpha = HCaseBRecCoeff + CollIoniz ;
00898 double beta = HCaseBRecCoeff*EdenExtraLocal + OtherIonization +
00899 (EdenExtraLocal - dense.gas_phase[ipHYDROGEN])*CollIoniz;
00900 double gamma = -dense.gas_phase[ipHYDROGEN]*(OtherIonization+EdenExtraLocal*CollIoniz);
00901
00902 double discriminant = POW2(beta) - 4.*alpha*gamma;
00903 if( discriminant <0 )
00904 {
00905
00906 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found negative discriminant.\n");
00907 TotalInsanity();
00908 }
00909
00910 dense.xIonDense[ipHYDROGEN][1] = (realnum)((-beta + sqrt(discriminant))/(2.*alpha));
00911 if( dense.xIonDense[ipHYDROGEN][1]> dense.gas_phase[ipHYDROGEN] )
00912 {
00913
00914 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H+)>n(H).\n");
00915 TotalInsanity();
00916 }
00917 dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] -
00918 dense.xIonDense[ipHYDROGEN][1];
00919 if( dense.xIonDense[ipHYDROGEN][0]<= 0 )
00920 {
00921
00922 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H0)<0.\n");
00923 TotalInsanity();
00924 }
00925 }
00926
00927
00928 if( dense.xIonDense[ipHYDROGEN][1] > 1e-30 )
00929 {
00930 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop = dense.xIonDense[ipHYDROGEN][0];
00931 }
00932 else
00933 {
00934 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop = 0.;
00935 }
00936
00937
00938
00939
00940 hydro.lgHInducImp = false;
00941 for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; i++ )
00942 {
00943 if( rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i]-1] > 0.01 )
00944 hydro.lgHInducImp = true;
00945 }
00946
00947
00948
00949
00950
00951
00952
00953
00954 if( dense.lgElmtOn[ipHELIUM] )
00955 {
00956
00957 xIoniz = (realnum)t_ADfA::Inst().coll_ion(2,2,phycon.te);
00958
00959 xIoniz = (realnum)(xIoniz*dense.eden + rfield.qhe*1e-18 + secondaries.csupra[ipHELIUM][1] );
00960 double RecTot = HCaseBRecCoeff*dense.eden;
00961 RatioIonizToRecomb = xIoniz/RecTot;
00962
00963
00964 xIoniz = (realnum)t_ADfA::Inst().coll_ion(2,1,phycon.te);
00965
00966 xIoniz = (realnum)(xIoniz*dense.eden + rfield.qheii*1e-18 + ionbal.CosRayIonRate );
00967
00968
00969 RecTot *= 4.;
00970 r3ov2 = xIoniz/RecTot;
00971
00972
00973 if( RatioIonizToRecomb > 0. )
00974 {
00975 dense.xIonDense[ipHELIUM][1] = (realnum)(dense.gas_phase[ipHELIUM]/(1./RatioIonizToRecomb + 1. + r3ov2));
00976 dense.xIonDense[ipHELIUM][0] = (realnum)(dense.xIonDense[ipHELIUM][1]/RatioIonizToRecomb);
00977 }
00978 else
00979 {
00980
00981 dense.xIonDense[ipHELIUM][1] = 0.;
00982 dense.xIonDense[ipHELIUM][0] = dense.gas_phase[ipHELIUM];
00983 }
00984
00985 dense.xIonDense[ipHELIUM][2] = (realnum)(dense.xIonDense[ipHELIUM][1]*r3ov2);
00986
00987 if( dense.xIonDense[ipHELIUM][2] > 1e-30 )
00988 {
00989 StatesElemNEW[ipHELIUM][ipHELIUM-ipH_LIKE][ipH1s].Pop = dense.xIonDense[ipHELIUM][1];
00990 }
00991 else
00992 {
00993 StatesElemNEW[ipHELIUM][ipHELIUM-ipH_LIKE][ipH1s].Pop = 0.;
00994 }
00995 }
00996 else
00997 {
00998
00999 dense.xIonDense[ipHELIUM][1] = 0.;
01000 dense.xIonDense[ipHELIUM][0] = 0.;
01001 dense.xIonDense[ipHELIUM][2] = 0.;
01002 StatesElemNEW[ipHELIUM][ipHELIUM-ipH_LIKE][ipH1s].Pop = 0.;
01003 }
01004
01005
01006 newEden = dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal + dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2];
01007
01008 loopCount++;
01009 }
01010
01011 while( loopCount < 10 && fabs(newEden/dense.eden - 1.) > 0.001 );
01012
01013 if( dense.xIonDense[ipHYDROGEN][0]<0.)
01014 TotalInsanity();
01015 else if( dense.xIonDense[ipHYDROGEN][0] == 0. )
01016 {
01017 fprintf(ioQQQ,"PROBLEM the derived atomic hydrogen density is zero.\n");
01018 if( dense.gas_phase[ipHYDROGEN]<1e-5 && rfield.uh > 1e10)
01019 {
01020 fprintf(ioQQQ,"This is almost certainly due to floating point "
01021 "limits on this computer.\nThe ionization parameter is very large,"
01022 " the density is very small,\nand the H^0 density cannot be"
01023 " stored in a float.\n");
01024 cdEXIT( EXIT_FAILURE );
01025 }
01026 }
01027 ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);
01028
01029
01030 dense.eden = newEden;
01031
01032 if( dense.eden <= SMALLFLOAT )
01033 {
01034
01035 fprintf(ioQQQ,"\n PROBLEM DISASTER - this simulation has no source"
01036 " of ionization. The electron density is zero. Consider "
01037 "adding a source of ionization such as cosmic rays.\n\n");
01038 cdEXIT(EXIT_FAILURE);
01039 }
01040
01041
01042 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
01043 {
01044 if( dense.lgElmtOn[nelem] )
01045 {
01046
01047
01048
01049
01050 dense.IonHigh[nelem] = nelem + 1;
01051
01052 dense.IonLow[nelem] = 0;
01053
01054
01055
01056 if( rfield.uh > 1e15 )
01057 {
01058
01059 while ( rfield.anu[Heavy.ipHeavy[nelem][dense.IonHigh[nelem]-1]] >
01060 rfield.anu[rfield.nflux] && dense.IonHigh[nelem]>1 )
01061 --dense.IonHigh[nelem];
01062 dense.IonLow[nelem] = max( 0 , dense.IonHigh[nelem]-1 );
01063 }
01064
01065
01066 for( ion=0; ion<dense.IonLow[nelem]; ++ion )
01067 dense.xIonDense[nelem][ion] = 0.;
01068
01069
01070
01071
01072 if( dense.lgSetIoniz[nelem] )
01073 {
01074 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] < dense.density_low_limit )
01075 ++dense.IonLow[nelem];
01076 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] < dense.density_low_limit )
01077 --dense.IonHigh[nelem];
01078 }
01079 }
01080 else
01081 {
01082
01083 dense.IonLow[nelem] = -1;
01084 dense.IonHigh[nelem] = -1;
01085 }
01086 }
01087
01088
01089
01090 EdenHeav = 0.;
01091 realnum atomFrac = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
01092 realnum firstIonFrac = dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN];
01093 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
01094 {
01095
01096
01097 if( dense.lgElmtOn[nelem] )
01098 {
01099 if( dense.IonLow[nelem] == 0)
01100 dense.xIonDense[nelem][0] = dense.gas_phase[nelem] * atomFrac;
01101 if( dense.IonLow[nelem] <= 1)
01102 {
01103 dense.xIonDense[nelem][1] = dense.gas_phase[nelem] * firstIonFrac;
01104 EdenHeav += dense.xIonDense[nelem][1];
01105 }
01106 }
01107 }
01108
01109
01110 dense.eden += EdenHeav;
01111
01112
01113 dense.eden = MAX2( SMALLFLOAT , dense.eden );
01114
01115 if( dense.EdenSet > 0. )
01116 {
01117 dense.eden = dense.EdenSet;
01118 }
01119
01120 dense.EdenHCorr = dense.eden;
01121
01122 if( dense.eden < 0. )
01123 {
01124 fprintf( ioQQQ, " PROBLEM DISASTER Negative electron density results in ContSetIntensity.\n" );
01125 fprintf( ioQQQ, "%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n",
01126 dense.eden, dense.xIonDense[ipHYDROGEN][1], dense.xIonDense[ipHELIUM][1],
01127 dense.xIonDense[ipHELIUM][2], dense.gas_phase[ipCARBON], dense.EdenExtra );
01128 TotalInsanity();
01129 }
01130
01131 if( dense.EdenSet > 0. )
01132 {
01133 dense.EdenTrue = dense.EdenSet;
01134 }
01135 else
01136 dense.EdenTrue = dense.eden;
01137
01138 if( trace.lgTrace )
01139 {
01140 fprintf( ioQQQ,
01141 " ContSetIntensity sets initial EDEN to %.4e, contributors H+=%.2e He+, ++= %.2e %.2e Heav %.2e extra %.2e\n",
01142 dense.eden ,
01143 dense.xIonDense[ipHYDROGEN][1],
01144 dense.xIonDense[ipHELIUM][1],
01145 2.*dense.xIonDense[ipHELIUM][2],
01146 EdenHeav,
01147 dense.EdenExtra);
01148 }
01149
01150
01151 occ1 = (realnum)(prt.fx1ryd/HNU3C2/PI4/FR1RYD);
01152 if( occ1 > 1. )
01153 rfield.lgOcc1Hi = true;
01154 else
01155 rfield.lgOcc1Hi = false;
01156
01157 if( trace.lgTrace && trace.lgConBug )
01158 {
01159 fprintf(ioQQQ,"\ntrace continuum print of %li incident spectral "
01160 "components\n", rfield.nShape);
01161 fprintf(ioQQQ," # type Illum Beam? 1/cos TimeVary?\n");
01162 for( rfield.ipSpec=0; rfield.ipSpec<rfield.nShape; ++rfield.ipSpec )
01163 {
01164 fprintf(ioQQQ,"%3li %6s %5i %c %.3f %c\n",
01165 rfield.ipSpec ,
01166 rfield.chSpType[rfield.ipSpec],
01167 rfield.Illumination[rfield.ipSpec],
01168 TorF(rfield.lgBeamed[rfield.ipSpec]),
01169 rfield.OpticalDepthScaleFactor[rfield.ipSpec],
01170 TorF(rfield.lgTimeVary[rfield.ipSpec]));
01171 }
01172 fprintf(ioQQQ,"\n");
01173
01174
01175 fprintf( ioQQQ, " H2,1=%5ld%5ld NX=%5ld IRC=%5ld\n",
01176 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2],
01177 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
01178 opac.ipCKshell,
01179 ionbal.ipCompRecoil[ipHYDROGEN][0] );
01180
01181 fprintf( ioQQQ, " CARBON" );
01182 for( i=0; i < 6; i++ )
01183 fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipCARBON][i] );
01184 fprintf( ioQQQ, "\n" );
01185
01186 fprintf( ioQQQ, " OXY" );
01187 for( i=0; i < 8; i++ )
01188 fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipOXYGEN][i] );
01189 fprintf( ioQQQ, "%5ld%5ld%5ld\n", opac.ipo3exc[0],
01190 oxy.i2d, oxy.i2p );
01191
01192 fprintf( ioQQQ,
01193 "\n\n PHOTONS PER CELL (NOT RYD)\n" );
01194 fprintf( ioQQQ,
01195 "\n\n nu, flux, wid, occ \n" );
01196 fprintf( ioQQQ,
01197 " " );
01198
01199 for( i=0; i < rfield.nflux; i++ )
01200 {
01201 fprintf( ioQQQ, "%4ld%10.2e%10.2e%10.2e%10.2e", i,
01202 rfield.anu[i], rfield.flux[0][i], rfield.widflx[i],
01203 rfield.OccNumbIncidCont[i] );
01204 }
01205 fprintf( ioQQQ, " \n" );
01206 }
01207
01208
01209
01210
01211 RT_OTS_Zero();
01212
01213 if( trace.lgTrace )
01214 {
01215 fprintf( ioQQQ, " ContSetIntensity returns, nflux=%5ld anu(nflux)=%11.4e eden=%10.2e\n",
01216 rfield.nflux, rfield.anu[rfield.nflux-1], dense.eden );
01217 }
01218
01219 return(0);
01220 }
01221
01222
01223 STATIC void sumcon(long int il,
01224 long int ih,
01225 realnum *q,
01226 realnum *p,
01227 realnum *panu)
01228 {
01229 long int i,
01230 iupper;
01231
01232 DEBUG_ENTRY( "sumcon()" );
01233
01234 *q = 0.;
01235 *p = 0.;
01236 *panu = 0.;
01237
01238
01239 iupper = MIN2(rfield.nflux,ih);
01240
01241
01242 for( i=il-1; i < iupper; i++ )
01243 {
01244
01245 *q += rfield.flux[0][i];
01246
01247
01248 *p += (realnum)(rfield.flux[0][i]*(rfield.anu[i]*EN1RYD));
01249
01250 *panu += (realnum)(rfield.flux[0][i]*(rfield.anu2[i]*EN1RYD));
01251 }
01252
01253 return;
01254 }
01255
01256
01257 STATIC void ptrcer(void)
01258 {
01259 char chCard[INPUT_LINE_LENGTH];
01260
01261 FILE * ioERRORS=NULL;
01262 bool lgEOL;
01263 char chKey;
01264 long int i,
01265 ipnt,
01266 j;
01267 double pnt,
01268 t1,
01269 t2;
01270
01271 DEBUG_ENTRY( "ptrcer()" );
01272
01273 fprintf( ioQQQ, " There are two ways to do this:\n");
01274 fprintf( ioQQQ, " do you want me to test all the pointers (enter y)\n");
01275 fprintf( ioQQQ, " or do you want to enter energies yourself? (enter n)\n" );
01276
01277 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
01278 {
01279 fprintf( ioQQQ, " error getting input \n" );
01280 cdEXIT(EXIT_FAILURE);
01281 }
01282
01283
01284 chKey = chCard[0];
01285
01286 if( chKey == 'n' )
01287 {
01288
01289 fprintf( ioQQQ, " Enter energy (Ryd); 0 to stop; negative is log.\n" );
01290 pnt = 1.;
01291 while( pnt!=0. )
01292 {
01293 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
01294 {
01295 fprintf( ioQQQ, " error getting input2 \n" );
01296 cdEXIT(EXIT_FAILURE);
01297 }
01298
01299 i = 1;
01300 pnt = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
01301
01302
01303 if( lgEOL || pnt==0. )
01304 {
01305 break;
01306 }
01307
01308
01309 if( pnt < 0. )
01310 {
01311 pnt = pow(10.,pnt);
01312 }
01313
01314
01315 ipnt = ipoint(pnt);
01316 fprintf( ioQQQ, " Cell num%4ld center:%10.2e width:%10.2e low:%10.2e hi:%10.2e convoc:%10.2e\n",
01317 ipnt, rfield.anu[ipnt-1], rfield.widflx[ipnt-1],
01318 rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]/2.,
01319 rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]/2.,
01320 rfield.convoc[ipnt-1] );
01321 }
01322 }
01323
01324 else if( chKey == 'y' )
01325 {
01326
01327 if( rfield.anu[0] - rfield.widflx[0]/2.*0.9 < continuum.filbnd[0] )
01328 {
01329 fprintf( ioQQQ," ipoint would crash since lowest desired energy of %e ryd is below limit of %e\n",
01330 rfield.anu[0] - rfield.widflx[0]/2.*0.9 , continuum.filbnd[0] );
01331 fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[0]);
01332 cdEXIT(EXIT_FAILURE);
01333 }
01334
01335 else if( rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 >
01336 continuum.filbnd[continuum.nrange] )
01337 {
01338 fprintf( ioQQQ," ipoint would crash since highest desired energy of %e ryd is above limit of %e\n",
01339 rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 ,
01340 continuum.filbnd[continuum.nrange-1] );
01341 fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[rfield.nflux]);
01342 fprintf( ioQQQ," this, previous cells are %e %e\n",
01343 rfield.anu[rfield.nflux-1],rfield.anu[rfield.nflux-2]);
01344 cdEXIT(EXIT_FAILURE);
01345 }
01346
01347
01348 fprintf( ioQQQ, " errors output on errors.txt\n");
01349 fprintf( ioQQQ, " IP(cor),IP(fount),nu lower, upper of found, desired cell.\n" );
01350
01351
01352 ioERRORS = NULL;
01353 for( i=0; i < rfield.nflux-1; i++ )
01354 {
01355 t1 = rfield.anu[i] - rfield.widflx[i]/2.*0.9;
01356 t2 = rfield.anu[i] + rfield.widflx[i]/2.*0.9;
01357
01358 j = ipoint(t1);
01359 if( j != i+1 )
01360 {
01361
01362 if( ioERRORS == NULL )
01363 {
01364 ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY );
01365 fprintf( ioQQQ," created errors.txt file with error summary\n");
01366 }
01367
01368 fprintf( ioQQQ, " Pointers do not agree for lower bound of cell%4ld, %e\n",
01369 i, rfield.anu[i]);
01370 fprintf( ioERRORS, " Pointers do not agree for lower bound of cell%4ld, %e\n",
01371 i, rfield.anu[i] );
01372 }
01373
01374 j = ipoint(t2);
01375 if( j != i+1 )
01376 {
01377
01378 if( ioERRORS == NULL )
01379 {
01380 ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY );
01381 fprintf( ioQQQ," created errors.txt file with error summary\n");
01382 }
01383 fprintf( ioQQQ, " Pointers do not agree for upper bound of cell%4ld, %e\n",
01384 i , rfield.anu[i]);
01385 fprintf( ioERRORS, " Pointers do not agree for upper bound of cell%4ld, %e\n",
01386 i , rfield.anu[i]);
01387 }
01388
01389 }
01390 }
01391
01392 else
01393 {
01394 fprintf( ioQQQ, "I do not understand this key, sorry. %c\n", chKey );
01395 cdEXIT(EXIT_FAILURE);
01396 }
01397
01398 if( ioERRORS!=NULL )
01399 fclose( ioERRORS );
01400 cdEXIT(EXIT_FAILURE);
01401 }
01402
01403
01404 STATIC void extin(realnum *ex1ryd)
01405 {
01406
01407 DEBUG_ENTRY( "extin()" );
01408
01409
01410
01411
01412 if( rfield.ExtinguishColumnDensity == 0. )
01413 {
01414 *ex1ryd = 1.;
01415
01416 for( long i=0; i<rfield.nupper; ++i )
01417 rfield.ExtinguishFactor[i] = 1.;
01418 }
01419 else
01420 {
01421 double absorb = 1. - rfield.ExtinguishLeakage;
01422 double factor = rfield.ExtinguishColumnDensity*
01423 rfield.ExtinguishConvertColDen2OptDepth;
01424
01425 *ex1ryd = (realnum)(rfield.ExtinguishLeakage + absorb*sexp(factor));
01426
01427
01428 long low = ipoint(rfield.ExtinguishLowEnergyLimit);
01429
01430 for( long i=0; i<low-1; ++i )
01431 rfield.ExtinguishFactor[i] = 1.;
01432
01433 for( long i=low-1; i < rfield.nupper; i++ )
01434 {
01435 realnum extfactor = (realnum)(rfield.ExtinguishLeakage + absorb*
01436 sexp(factor*(pow(rfield.anu[i],rfield.ExtinguishEnergyPowerLow))));
01437
01438 rfield.ExtinguishFactor[i] = extfactor;
01439 rfield.flux_beam_const[i] *= extfactor;
01440 rfield.flux_beam_time[i] *= extfactor;
01441 rfield.flux_isotropic[i] *= extfactor;
01442
01443 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
01444 rfield.flux_isotropic[i];
01445 }
01446 }
01447 return;
01448 }
01449
01450
01451 STATIC void conorm(void)
01452 {
01453 long int i;
01454 double xLog_radius_inner,
01455 diff,
01456 f,
01457 flx1,
01458 flx2,
01459 pentrd,
01460 qentrd;
01461
01462 DEBUG_ENTRY( "conorm()" );
01463
01464 xLog_radius_inner = log10(radius.Radius);
01465
01466
01467 for( i=0; i < rfield.nShape; i++ )
01468 {
01469 if( strcmp(rfield.chRSpec[i],"UNKN") == 0 )
01470 {
01471 fprintf( ioQQQ, " UNKN spectral normalization cannot happen.\n" );
01472 fprintf( ioQQQ, " conorm punts.\n" );
01473 cdEXIT(EXIT_FAILURE);
01474 }
01475
01476 else if( strcmp(rfield.chRSpec[i],"SQCM") != 0 &&
01477 strcmp(rfield.chRSpec[i],"4 PI") != 0 )
01478 {
01479 fprintf( ioQQQ, " chRSpec must be SQCM or 4 PI, and it was %4.4s. This cannot happen.\n",
01480 rfield.chRSpec[i] );
01481 fprintf( ioQQQ, " conorm punts.\n" );
01482 cdEXIT(EXIT_FAILURE);
01483 }
01484
01485
01486
01487
01488 if( strcmp(rfield.chSpType[i],"VOLK ") == 0 )
01489 {
01490 if( !fp_equal( rfield.RSFCheck[i], continuum.ResolutionScaleFactor ) )
01491 {
01492 fprintf( ioQQQ,"\n\n PROBLEM DISASTER At least one of the compiled stellar atmosphere"
01493 " grids has been compiled with a different energy grid resolution factor.\n" );
01494 fprintf( ioQQQ, " Please recompile this file using the COMPILE STARS command "
01495 "and make sure that you use the correct SET CONTINUUM RESOLUTION factor.\n" );
01496 cdEXIT(EXIT_FAILURE);
01497 }
01498 }
01499 else if( strcmp(rfield.chSpType[i],"READ ") == 0 )
01500 {
01501 if( !fp_equal( rfield.RSFCheck[i], continuum.ResolutionScaleFactor ) )
01502 {
01503 fprintf( ioQQQ,"\n\n PROBLEM DISASTER The file read by the TABLE READ command "
01504 "has been compiled with a different energy grid resolution factor.\n" );
01505 fprintf( ioQQQ, " Please recompile this file using the SAVE TRANSMITTED CONTINUUM "
01506 "command and use the correct SET CONTINUUM RESOLUTION factor.\n" );
01507 cdEXIT(EXIT_FAILURE);
01508 }
01509 }
01510 }
01511
01512
01513
01514 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01515 {
01516 if( !fp_equal( gv.bin[nd]->RSFCheck, continuum.ResolutionScaleFactor ) )
01517 {
01518 fprintf( ioQQQ,"\n\n PROBLEM DISASTER At least one of the grain opacity files "
01519 "has been compiled with a different energy grid resolution factor.\n" );
01520 fprintf( ioQQQ, " Please recompile this file using the COMPILE GRAINS command "
01521 "and make sure that you use the correct SET CONTINUUM RESOLUTION factor.\n" );
01522 cdEXIT(EXIT_FAILURE);
01523 }
01524 }
01525
01526
01527
01528
01529 radius.pirsq = 0.;
01530 radius.lgPredLumin = false;
01531
01532
01533 for( i=0; i < rfield.nShape; i++ )
01534 {
01535 if( strcmp(rfield.chRSpec[i],"4 PI") == 0 )
01536 {
01537 radius.pirsq = (realnum)(1.0992099 + 2.*xLog_radius_inner);
01538 radius.lgPredLumin = true;
01539
01540 rfield.totpow[i] -= radius.pirsq;
01541
01542 if( trace.lgTrace )
01543 {
01544 fprintf( ioQQQ,
01545 " conorm converts continuum %ld from luminosity to intensity.\n",
01546 i );
01547 }
01548 }
01549 }
01550
01551
01552 if( radius.lgPredLumin && !radius.lgRadiusKnown )
01553 {
01554 fprintf(ioQQQ,"PROBLEM DISASTER conorm: - A continuum source was specified as a luminosity, but the inner radius of the cloud was not set.\n");
01555 fprintf(ioQQQ,"Please set an inner radius.\nSorry.\n");
01556 cdEXIT(EXIT_FAILURE);
01557 }
01558
01559
01560
01561 for( i=0; i < rfield.nShape; i++ )
01562 {
01563 if( strcmp(rfield.chSpNorm[i],"IONI") == 0 )
01564 {
01565
01566
01567 rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) + log10(SPEEDLIGHT);
01568 strcpy( rfield.chSpNorm[i], "Q(H)" );
01569 if( trace.lgTrace )
01570 {
01571 fprintf( ioQQQ,
01572 " conorm converts continuum %ld from ionizat par to q(h).\n",
01573 i );
01574 }
01575 }
01576 }
01577
01578
01579 for( i=0; i < rfield.nShape; i++ )
01580 {
01581 if( strcmp(rfield.chSpNorm[i],"IONX") == 0 )
01582 {
01583
01584 rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) - log10(PI4);
01585 strcpy( rfield.chSpNorm[i], "LUMI" );
01586 if( trace.lgTrace )
01587 {
01588 fprintf( ioQQQ, " conorm converts continuum%3ld from x-ray ionizat par to I.\n",
01589 i );
01590 }
01591 }
01592 }
01593
01594
01595 if( trace.lgTrace )
01596 {
01597 if( radius.lgPredLumin )
01598 {
01599 fprintf( ioQQQ, " Cloudy will predict lumin into 4pi\n" );
01600 }
01601 else
01602 {
01603 fprintf( ioQQQ, " Cloudy will do surface flux for lumin\n" );
01604 }
01605 }
01606
01607
01608
01609
01610 if( !radius.lgPredLumin )
01611 {
01612 geometry.covgeo = 1.;
01613 }
01614
01615
01616
01617 for( i=0; i < rfield.nShape; i++ )
01618 {
01619 rfield.ipSpec = i;
01620
01621
01622 if( strcmp(rfield.chSpType[rfield.ipSpec],"LASER") == 0 )
01623 {
01624 if( !( rfield.range[rfield.ipSpec][0] < rfield.slope[rfield.ipSpec] &&
01625 rfield.range[rfield.ipSpec][1] > rfield.slope[rfield.ipSpec]) )
01626 {
01627 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",
01628 rfield.slope[rfield.ipSpec],
01629 rfield.range[rfield.ipSpec][0],
01630 rfield.range[rfield.ipSpec][1]);
01631 fprintf(ioQQQ,"Please specify the continuum flux where the laser is active.\n");
01632 cdEXIT(EXIT_FAILURE);
01633 }
01634 }
01635
01636 if( trace.lgTrace )
01637 {
01638 long int jj;
01639 fprintf( ioQQQ, " conorm continuum number %ld is shape %s range is %.2e %.2e\n",
01640 i,
01641 rfield.chSpType[i],
01642 rfield.range[i][0],
01643 rfield.range[i][1] );
01644 fprintf( ioQQQ, "the continuum points follow\n");
01645 jj = 0;
01646 if( rfield.lgContMalloc[rfield.ipSpec] )
01647 {
01648 while( rfield.tNu[rfield.ipSpec][jj].Ryd() != 0. && jj < 100 )
01649 {
01650 fprintf( ioQQQ, "%li %e %e\n",
01651 jj ,
01652 rfield.tNu[rfield.ipSpec][jj].Ryd(),
01653 rfield.tslop[rfield.ipSpec][jj] );
01654 ++jj;
01655 }
01656 }
01657 }
01658
01659 if( strcmp(rfield.chSpNorm[i],"RATI") == 0 )
01660 {
01661
01662
01663
01664 if( trace.lgTrace )
01665 {
01666 fprintf( ioQQQ, " conorm this is ratio to 1st con\n" );
01667 }
01668
01669
01670 if( i == 0 )
01671 {
01672 fprintf( ioQQQ, " I cant form a ratio if continuum is first source.\n" );
01673 cdEXIT(EXIT_FAILURE);
01674 }
01675
01676
01677 rfield.ipSpec -= 1;
01678 flx1 = ffun1(rfield.range[i][0])*rfield.spfac[rfield.ipSpec]*
01679 rfield.range[i][0];
01680
01681
01682 if( flx1 <= 0. )
01683 {
01684 fprintf( ioQQQ, " Previous continua were zero where ratio is desired.\n" );
01685 cdEXIT(EXIT_FAILURE);
01686 }
01687
01688
01689 rfield.ipSpec += 1;
01690
01691
01692 flx1 *= rfield.totpow[i];
01693
01694
01695 flx2 = ffun1(rfield.range[i][1])*rfield.range[i][1];
01696
01697
01698 rfield.spfac[i] = flx1/flx2;
01699 if( trace.lgTrace )
01700 {
01701 fprintf( ioQQQ, " conorm ratio will set scale fac to%10.3e at%10.2e Ryd.\n",
01702 rfield.totpow[i], rfield.range[i][0] );
01703 }
01704 }
01705
01706 else if( strcmp(rfield.chSpNorm[i],"FLUX") == 0 )
01707 {
01708
01709
01710 f = ffun1(rfield.range[i][0]);
01711
01712
01713
01714 if( f<=0. )
01715 {
01716 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",
01717 i , rfield.range[i][0]);
01718
01719 if( strcmp(rfield.chSpType[i],"INTER") == 0 )
01720 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");
01721 fprintf( ioQQQ, " Also check that the numbers used to specify the shape and intensity do not under or overflow on this cpu.\n\n");
01722
01723 cdEXIT(EXIT_FAILURE);
01724 }
01725
01726
01727 f = MAX2(1e-37,f);
01728 f = log10(f) + log10(rfield.range[i][0]*EN1RYD/FR1RYD);
01729
01730 f = rfield.totpow[i] - f;
01731
01732 if( f > 35. )
01733 {
01734 fprintf( ioQQQ, " PROBLEM DISASTER Continuum source %ld is too intense for this cpu - is it normalized correctly?\n",
01735 i );
01736 cdEXIT(EXIT_FAILURE);
01737 }
01738
01739 rfield.spfac[i] = pow(10.,f);
01740 if( trace.lgTrace )
01741 {
01742 fprintf( ioQQQ, " conorm will set log fnu to%10.3e at%10.2e Ryd. Factor=%11.4e\n",
01743 rfield.totpow[i], rfield.range[i][0], rfield.spfac[i] );
01744 }
01745 }
01746
01747 else if( strcmp(rfield.chSpNorm[i],"Q(H)") == 0 ||
01748 strcmp(rfield.chSpNorm[i],"PHI ") == 0 )
01749 {
01750
01751 if( trace.lgTrace )
01752 {
01753 fprintf( ioQQQ, " conorm calling qintr range=%11.3e %11.3e desired val is %11.3e\n",
01754 rfield.range[i][0],
01755 rfield.range[i][1] ,
01756 rfield.totpow[i]);
01757 }
01758
01759
01760
01761 ASSERT( rfield.range[i][0] < rfield.range[i][1] );
01762 qentrd = qintr(&rfield.range[i][0],&rfield.range[i][1]);
01763
01764
01765 diff = rfield.totpow[i] - qentrd;
01766
01767
01768
01769
01770 if( diff < -35. || diff > 35. )
01771 {
01772 fprintf( ioQQQ, " PROBLEM DISASTER Continuum source specified is too extreme.\n" );
01773 fprintf( ioQQQ,
01774 " The integral over the continuum shape gave (log) %.3e photons, and the command requested (log) %.3e.\n" ,
01775 qentrd , rfield.totpow[i]);
01776 fprintf( ioQQQ,
01777 " The difference in the log is %.3e.\n" ,
01778 diff );
01779 if( diff>0. )
01780 {
01781 fprintf( ioQQQ, " The continuum source is too bright.\n" );
01782 }
01783 else
01784 {
01785 fprintf( ioQQQ, " The continuum source is too faint.\n" );
01786 }
01787
01788 fprintf( ioQQQ, " The usual cause for this problem is an incorrect continuum intensity/luminosity or radius command.\n" );
01789 fprintf( ioQQQ, " There were a total of %li continuum shape commands entered - the problem is with number %li.\n",
01790 rfield.nShape , i+1 );
01791 cdEXIT(EXIT_FAILURE);
01792 }
01793
01794 else
01795 {
01796 rfield.spfac[i] = pow(10.,diff);
01797 }
01798
01799 if( trace.lgTrace )
01800 {
01801 fprintf( ioQQQ, " conorm finds Q over range from%11.4e-%11.4e Ryd, integral= %10.4e Factor=%11.4e\n",
01802 rfield.range[i][0],
01803 rfield.range[i][1],
01804 qentrd ,
01805 rfield.spfac[i] );
01806 }
01807 }
01808
01809 else if( strcmp(rfield.chSpNorm[i],"LUMI") == 0 )
01810 {
01811
01812
01813
01814 pentrd = pintr(rfield.range[i][0],rfield.range[i][1]) +
01815 log10(EN1RYD);
01816 f = rfield.totpow[i] - pentrd;
01817
01818
01819 if( f > 35. )
01820 {
01821 fprintf( ioQQQ, " PROBLEM DISASTER Continuum source %ld is too intense for this cpu - is it normalized correctly?\n",
01822 i );
01823 cdEXIT(EXIT_FAILURE);
01824 }
01825
01826 rfield.spfac[i] = pow(10.,f);
01827 if( trace.lgTrace )
01828 {
01829 fprintf( ioQQQ, " conorm finds luminosity range is %10.3e to %9.3e Ryd, factor is %11.4e\n",
01830 rfield.range[i][0], rfield.range[i][1],
01831 rfield.spfac[i] );
01832 }
01833 }
01834
01835 else
01836 {
01837 fprintf( ioQQQ, "PROBLEM DISASTER What chSpNorm label is this? =%s=\n", rfield.chSpNorm[i]);
01838 TotalInsanity();
01839 }
01840
01841
01842
01843 if( 1./rfield.spfac[i] == 0. || rfield.spfac[i] == 0. )
01844 {
01845 fprintf( ioQQQ, "PROBLEM DISASTER conorm finds infinite continuum scale factor.\n" );
01846 fprintf( ioQQQ, "The continuum is too intense to compute with this cpu.\n" );
01847 fprintf( ioQQQ, "Were the intensity and luminosity commands switched?\n" );
01848 fprintf( ioQQQ, "Sorry, but I cannot go on.\n" );
01849 cdEXIT(EXIT_FAILURE);
01850 }
01851 }
01852
01853
01854
01855
01856 radius.Conv2PrtInten = radius.pirsq;
01857
01858
01859 if( geometry.iEmissPower == 1 )
01860 {
01861 if( radius.lgPredLumin )
01862 {
01863
01864 radius.Conv2PrtInten -= (log10(2.) + xLog_radius_inner);
01865 }
01866 else if( !radius.lgPredLumin )
01867 {
01868
01869 fprintf( ioQQQ, "PROBLEM DISASTER conorm: Aperture slit specified, but not predicting luminosity.\n" );
01870 fprintf( ioQQQ, "conorm: Please specify an inner radius to determine L.\nSorry\n" );
01871 cdEXIT(EXIT_FAILURE);
01872 }
01873 }
01874 if( geometry.iEmissPower == 0 && radius.lgPredLumin )
01875 {
01876
01877 radius.Conv2PrtInten = log10(2.);
01878 }
01879
01880
01881 if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth )
01882 {
01883
01884
01885 if( geometry.iEmissPower == 0 )
01886 radius.Conv2PrtInten += log10(double(geometry.size)/SQAS_SKY);
01887 else if( geometry.iEmissPower == 1 )
01888 radius.Conv2PrtInten += log10(double(geometry.size)/(PI4*AS1RAD*radius.distance));
01889 else if( geometry.iEmissPower == 2 )
01890 radius.Conv2PrtInten -= log10( PI4*pow2(radius.distance) );
01891 else
01892 TotalInsanity();
01893 }
01894
01895
01896 if( prt.lgSurfaceBrightness )
01897 {
01898 if( radius.pirsq != 0. )
01899 {
01900
01901 fprintf( ioQQQ, " PROBLEM DISASTER Sorry, but both luminosity and surface brightness have been requested for lines.\n" );
01902 fprintf( ioQQQ, " the PRINT LINE SURFACE BRIGHTNESS command can only be used when lines are predicted per unit cloud area.\n" );
01903 cdEXIT(EXIT_FAILURE);
01904 }
01905 if( prt.lgSurfaceBrightness_SR )
01906 {
01907
01908 radius.Conv2PrtInten -= log10( PI4 );
01909 }
01910 else
01911 {
01912
01913 radius.Conv2PrtInten -= log10(SQAS_SKY);
01914 }
01915 }
01916 return;
01917 }
01918
01919
01920 STATIC double qintr(double *qenlo,
01921 double *qenhi)
01922 {
01923 long int i,
01924 ipHi,
01925 ipLo,
01926 j;
01927 double qintr_v,
01928 sum,
01929 wanu;
01930
01931 DEBUG_ENTRY( "qintr()" );
01932
01933
01934
01935
01936 ASSERT(*qenhi > *qenlo);
01937 ipLo = ipoint(*qenlo);
01938 ipHi = ipoint(*qenhi);
01939
01940 sum = 0.;
01941 for( i=ipLo-1; i < (ipHi - 1); i++ )
01942 {
01943
01944 for( j=0; j < 4; j++ )
01945 {
01946 wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
01947
01948
01949 wanu = MAX2( wanu , rfield.emm );
01950 wanu = MIN2( wanu , rfield.egamry );
01951 sum += fweigh[j]*ffun1(wanu)*rfield.widflx[i];
01952 }
01953 }
01954
01955 if( sum <= 0. )
01956 {
01957 fprintf( ioQQQ, " PROBLEM DISASTER Photon number sum in QINTR is %.3e\n",
01958 sum );
01959 fprintf( ioQQQ, " This source has no ionizing radiation, and the number of ionizing photons was specified.\n" );
01960 fprintf( ioQQQ, " This was continuum source number%3ld\n",
01961 rfield.ipSpec );
01962 fprintf( ioQQQ, " Sorry, but I cannot go on. ANU and FLUX arrays follow. Enjoy.\n" );
01963 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" );
01964 for( i=0; i < rfield.nupper; i++ )
01965 {
01966 fprintf( ioQQQ, "%.2e\t%.2e\n",
01967 rfield.anu[i],
01968 rfield.flux[0][i] );
01969 }
01970 cdEXIT(EXIT_FAILURE);
01971 }
01972 else
01973 {
01974 qintr_v = log10(sum);
01975 }
01976 return( qintr_v );
01977 }
01978
01979
01980
01981 STATIC double pintr(double penlo,
01982 double penhi)
01983 {
01984 long int i,
01985 j;
01986 double fsum,
01987 pintr_v,
01988 sum,
01989 wanu,
01990 wfun;
01991 long int ip1 , ip2;
01992
01993 DEBUG_ENTRY( "pintr()" );
01994
01995
01996
01997
01998 sum = 0.;
01999
02000
02001
02002
02003 ip1 = ipoint( penlo );
02004
02005 ip2 = ipoint( penhi );
02006
02007 for( i=ip1-1; i < ip2-1; i++ )
02008 {
02009 fsum = 0.;
02010 for( j=0; j < 4; j++ )
02011 {
02012 wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
02013
02014 wfun = fweigh[j]*ffun1(wanu)*wanu;
02015 fsum += wfun;
02016 }
02017 sum += fsum*rfield.widflx[i];
02018 }
02019
02020 if( sum > 0. )
02021 {
02022 pintr_v = log10(sum);
02023 }
02024 else
02025 {
02026 pintr_v = -38.;
02027 }
02028
02029 return( pintr_v );
02030 }