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