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