00001
00002
00003
00004 #include "cddefines.h"
00005 #include "taulines.h"
00006 #include "physconst.h"
00007 #include "iso.h"
00008 #include "geometry.h"
00009 #include "heavy.h"
00010 #include "dense.h"
00011 #include "prt.h"
00012 #include "opacity.h"
00013 #include "coolheavy.h"
00014 #include "optimize.h"
00015 #include "phycon.h"
00016 #include "rfield.h"
00017 #include "predcont.h"
00018 #include "lines_service.h"
00019 #include "radius.h"
00020 #include "continuum.h"
00021 #include "lines.h"
00022 #include "elementnames.h"
00023
00024 void lines_continuum(void)
00025 {
00026
00027 double f1,
00028 f2 ,
00029 bac ,
00030 flow;
00031 long i,nBand;
00032
00033 DEBUG_ENTRY( "lines_continuum()" );
00034
00035
00036
00037 const bool KILL_CONT = false;
00038
00039 i = StuffComment( "continua" );
00040 linadd( 0., (realnum)i , "####", 'i',
00041 " start continua");
00042
00043
00044 if( geometry.iEmissPower == 2 )
00045 {
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00057
00058
00059
00060
00061
00062
00063
00064
00065 f1 = (rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] +
00066 rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/radius.r1r0sq )/
00067 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1];
00068
00069
00070
00071 f2 = (rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]+
00072 rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/radius.r1r0sq )/
00073 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2];
00074
00075
00076 f1 = f1*0.250*0.250*EN1RYD*radius.r1r0sq;
00077 f2 = f2*0.250*0.250*EN1RYD*radius.r1r0sq;
00078 bac = (f1 - f2);
00079
00080
00081
00082
00083
00084
00085 if( LineSave.ipass > 0 )
00086 {
00087 LineSv[LineSave.nsum].SumLine[0] = 0.;
00088 LineSv[LineSave.nsum].SumLine[1] = 0.;
00089 }
00090
00091 linadd(MAX2(0.,bac)/radius.dVeffAper,3646,"Bac ",'i',
00092 "residual flux at head of Balmer continuum, nuFnu ");
00093
00094
00095
00096
00097 if( KILL_CONT && LineSave.ipass > 0 )
00098 {
00099 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00100 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00101 }
00102
00103
00104 if( LineSave.ipass > 0 )
00105 {
00106 LineSv[LineSave.nsum].SumLine[0] = 0.;
00107 LineSv[LineSave.nsum].SumLine[1] = 0.;
00108 }
00109
00110 linadd(f1/radius.dVeffAper,3645,"nFnu",'i',
00111 "total flux above head of Balmer continuum, nuFnu ");
00112
00113
00114
00115
00116 if( KILL_CONT && LineSave.ipass > 0 )
00117 {
00118 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00119 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00120 }
00121
00122
00123 if( LineSave.ipass > 0 )
00124 {
00125 LineSv[LineSave.nsum].SumLine[0] = 0.;
00126 LineSv[LineSave.nsum].SumLine[1] = 0.;
00127 }
00128
00129 linadd(f2/radius.dVeffAper,3647,"nFnu",'i',
00130 "total flux above head of Balmer continuum, nuFnu ");
00131
00132
00133
00134
00135 if( KILL_CONT && LineSave.ipass > 0 )
00136 {
00137 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00138 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00139 }
00140
00141
00142
00143
00144
00145
00146
00147 f1 = rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/
00148 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1];
00149
00150
00151 f2 = rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/
00152 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2];
00153
00154
00155 bac = (f1 - f2)*0.250*0.250*EN1RYD*radius.r1r0sq;
00156
00157
00158 if( LineSave.ipass > 0 )
00159 {
00160 LineSv[LineSave.nsum].SumLine[0] = 0.;
00161 LineSv[LineSave.nsum].SumLine[1] = 0.;
00162 }
00163
00164 linadd(MAX2(0.,bac)/radius.dVeffAper,3646,"cout",'i',
00165 "residual flux in Balmer continuum, nuFnu ");
00166
00167
00168
00169
00170 if( KILL_CONT && LineSave.ipass > 0 )
00171 {
00172 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00173 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00174 }
00175
00176
00177
00178
00179
00180
00181
00182 f1 = rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/
00183 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1];
00184
00185 f2 = rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/
00186 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2];
00187
00188
00189 bac = (f1 - f2)*0.250*0.250*EN1RYD;
00190
00191
00192 if( LineSave.ipass > 0 )
00193 {
00194 LineSv[LineSave.nsum].SumLine[0] = 0.;
00195 LineSv[LineSave.nsum].SumLine[1] = 0.;
00196 }
00197
00198 linadd(MAX2(0.,bac)/radius.dVeffAper,3646,"cref",'i',
00199 "residual flux in Balmer continuum, nuFnu ");
00200
00201
00202
00203
00204 if( KILL_CONT && LineSave.ipass > 0 )
00205 {
00206 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00207 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00208 }
00209
00210
00211
00212 if( nzone > 0 )
00213 {
00214
00215 f1 = rfield.ConEmitLocal[nzone][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/
00216 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1];
00217 f2 = rfield.ConEmitLocal[nzone][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/
00218 rfield.widflx[iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2];
00219 bac = (f1 - f2)*0.250*0.250*EN1RYD;
00220 }
00221 else
00222 {
00223 f1 = 0.;
00224 f2 = 0.;
00225 bac = 0.;
00226 }
00227
00228 linadd(MAX2(0.,bac),3646,"thin",'i',
00229 "residual flux in Balmer continuum, nuFnu ");
00230
00231 linadd(continuum.cn4861/radius.dVeffAper,4860,"Inci",'i',
00232 "incident continuum nu*f_nu at H-beta, at illuminated face of cloud ");
00233
00234
00235
00236
00237
00238
00239
00240 if( LineSave.ipass > 0 )
00241 {
00242 LineSv[LineSave.nsum-1].SumLine[1] = LineSv[LineSave.nsum-1].SumLine[0];
00243 }
00244
00245 linadd(continuum.cn1216/radius.dVeffAper,1215,"Inci",'i',
00246 "incident continuum nu*f_nu near Ly-alpha, at illuminated face of cloud");
00247
00248 if( LineSave.ipass > 0 )
00249 {
00250 LineSv[LineSave.nsum-1].SumLine[1] = LineSv[LineSave.nsum-1].SumLine[0];
00251 }
00252
00253 if( LineSave.ipass > 0 )
00254 {
00255 continuum.cn4861 = 0.;
00256 continuum.cn1216 = 0.;
00257 }
00258 }
00259
00260 flow = (iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].RadRecomb[ipRecRad] +
00261 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].RadRecomb[ipRecRad])*
00262 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].RadRecomb[ipRecEsc]*
00263 dense.eden*dense.xIonDense[ipHYDROGEN][1]* 5.45e-12;
00264 linadd(flow,0,"Ba C",'i',
00265 "integrated Balmer continuum emission");
00266
00267 if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_max >= 3 )
00268 {
00269 flow = ( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3s].RadRecomb[ipRecRad]*
00270 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3s].RadRecomb[ipRecEsc] +
00271 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3p].RadRecomb[ipRecRad]*
00272 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3p].RadRecomb[ipRecEsc] +
00273 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3d].RadRecomb[ipRecRad]*
00274 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3d].RadRecomb[ipRecEsc] ) *
00275 dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
00276 }
00277 else
00278 {
00279 flow = iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].RadRecomb[ipRecRad]*
00280 iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].RadRecomb[ipRecEsc]*
00281 dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
00282 }
00283 linadd(flow,0,"PA C",'i',
00284 "Paschen continuum emission ");
00285
00286
00287
00288
00289
00290 if( geometry.iEmissPower == 2 )
00291 {
00292 for( nBand=0; nBand < continuum.nContBand; ++nBand )
00293 {
00294 double EmergentContinuum = 0.;
00295 double DiffuseEmission = 0.;
00296 if( LineSave.ipass > 0 )
00297 {
00298
00299 for( i=continuum.ipContBandLow[nBand]; i<=continuum.ipContBandHi[nBand]; ++i )
00300 {
00301
00302
00303 double EdgeCorrection = 1.;
00304 if( i==continuum.ipContBandLow[nBand] )
00305 EdgeCorrection = continuum.BandEdgeCorrLow[nBand];
00306 else if( i==continuum.ipContBandHi[nBand])
00307 EdgeCorrection = continuum.BandEdgeCorrHi[nBand];
00308
00309 double xIntenOut =
00310
00311 rfield.flux[0][i-1] +
00312
00313
00314 (rfield.ConEmitOut[0][i-1] +
00315
00316
00317 rfield.outlin[0][i-1])*geometry.covgeo;
00318 xIntenOut *= EdgeCorrection;
00319
00320
00321
00322
00323 double xIntenIn = 0.;
00324 if( opac.E2TauAbsFace[i-1] > 0. )
00325 xIntenIn = (double)rfield.ConEmitReflec[0][i-1]/
00326 (double)opac.E2TauAbsFace[i-1]*geometry.covgeo;
00327
00328 xIntenIn += rfield.reflin[0][i-1]*geometry.covgeo;
00329 xIntenIn *= EdgeCorrection;
00330
00331
00332 EmergentContinuum += rfield.anu[i-1] *
00333 emergent_line( xIntenIn , xIntenOut , i )
00334 / SDIV(opac.tmn[i-1]);
00335
00336
00337 DiffuseEmission += (rfield.ConEmitLocal[nzone][i-1] +
00338 rfield.DiffuseLineEmission[i-1])*rfield.anu[i-1]*
00339 EdgeCorrection;
00340 }
00341 }
00342
00343
00344
00345
00346 double corr = emergent_line( 0.5 , 0.5 ,
00347 (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 );
00348 if( corr < SMALLFLOAT )
00349 EmergentContinuum = 0.;
00350 else
00351 EmergentContinuum /= corr;
00352
00353
00354 EmergentContinuum *= EN1RYD*radius.r1r0sq/radius.dVeffAper;
00355 DiffuseEmission *= EN1RYD;
00356
00357 if( LineSave.ipass > 0 )
00358 {
00359 LineSv[LineSave.nsum].SumLine[0] = 0.;
00360 LineSv[LineSave.nsum].SumLine[1] = 0.;
00361 }
00362
00363 lindst( EmergentContinuum ,
00364
00365
00366
00367
00368
00369
00370 -continuum.ContBandWavelength[nBand] ,
00371 continuum.chContBandLabels[nBand],
00372 (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 ,
00373 'i' , false,
00374 "continuum bands defined in continuum_bands.ini");
00375
00376
00377
00378 if( LineSave.ipass > 0 )
00379 {
00380 LineSv[LineSave.nsum-1].emslin[0] = DiffuseEmission;
00381 LineSv[LineSave.nsum-1].emslin[1] = DiffuseEmission;
00382 }
00383 }
00384 }
00385
00386 linadd(MAX2(0.,CoolHeavy.brems_cool_net),0,"HFFc",'c',
00387 "net free-free cooling, ALL species, free-free heating subtracted, so nearly cancels with cooling in LTE ");
00388
00389 linadd(MAX2(0.,-CoolHeavy.brems_cool_net),0,"HFFh",'h',
00390 "net free-free heating, nearly cancels with cooling in LTE ");
00391
00392 linadd(CoolHeavy.brems_cool_h,0,"H FF",'i',
00393 " H brems (free-free) cooling ");
00394
00395 linadd(CoolHeavy.brems_heat_total,0,"FF H",'i',
00396 "total free-free heating ");
00397
00398 linadd(CoolHeavy.brems_cool_he,0,"HeFF",'i',
00399 "He brems emission ");
00400
00401 linadd(CoolHeavy.heavfb,0,"MeFB",'c',
00402 "heavy element recombination cooling ");
00403
00404 linadd(CoolHeavy.brems_cool_metals,0,"MeFF",'i',
00405 "heavy elements (metals) brems cooling, heat not subtracted ");
00406
00407 linadd(CoolHeavy.brems_cool_h+CoolHeavy.brems_cool_he+CoolHeavy.brems_cool_metals,0,"ToFF",'i',
00408 "total brems emission - total cooling but not minus heating ");
00409
00410 linadd((CoolHeavy.brems_cool_h+CoolHeavy.brems_cool_he)*sexp(5.8e6/phycon.te),0,"FF X",'i',
00411 "part of H brems, in x-ray beyond 0.5KeV ");
00412
00413 linadd(CoolHeavy.eebrm,0,"eeff",'c',
00414 "electron - electron brems ");
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428 t_PredCont& PredCont = t_PredCont::Inst();
00429 if( LineSave.ipass == 0 )
00430 PredCont.set_offset(LineSave.nsum);
00431
00432
00433 if( geometry.iEmissPower == 2 )
00434 {
00435 for( i=0; i < long(PredCont.size()); i++ )
00436 {
00437 double SourceTransmitted , Cont_nInu;
00438 double SourceReflected, DiffuseOutward, DiffuseInward;
00439 double renorm;
00440
00441
00442
00443 (*TauDummy).WLAng() = (realnum)PredCont[i].Angstrom();
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454 renorm = rfield.anu2[PredCont[i].ip_C()]*EN1RYD/rfield.widflx[PredCont[i].ip_C()];
00455
00456
00457 if( prt.lgDiffuseInward )
00458 {
00459 DiffuseInward = rfield.ConEmitReflec[0][PredCont[i].ip_C()]*renorm;
00460 }
00461 else
00462 {
00463 DiffuseInward = 0.;
00464 }
00465
00466
00467 if( prt.lgDiffuseOutward )
00468 {
00469 DiffuseOutward = rfield.ConEmitOut[0][PredCont[i].ip_C()]*renorm*radius.r1r0sq;
00470 }
00471 else
00472 {
00473 DiffuseOutward = 0.;
00474 }
00475
00476
00477 if( prt.lgSourceReflected )
00478 {
00479 SourceReflected = rfield.ConRefIncid[0][PredCont[i].ip_C()]*renorm;
00480 }
00481 else
00482 {
00483 SourceReflected = 0.;
00484 }
00485
00486
00487 if( prt.lgSourceTransmitted )
00488 {
00489 SourceTransmitted = rfield.flux[0][PredCont[i].ip_C()]*renorm*radius.r1r0sq;
00490 }
00491 else
00492 {
00493 SourceTransmitted = 0.;
00494 }
00495
00496
00497
00498 if( LineSave.ipass > 0 )
00499 {
00500 LineSv[LineSave.nsum].SumLine[0] = 0.;
00501 LineSv[LineSave.nsum].SumLine[1] = 0.;
00502 }
00503
00504 linadd((DiffuseInward+SourceReflected+DiffuseOutward+SourceTransmitted)/radius.dVeffAper,
00505 (*TauDummy).WLAng(),"nFnu",'i',
00506 "total continuum at selected energy points " );
00507
00508
00509
00510
00511 if( KILL_CONT && LineSave.ipass > 0 )
00512 {
00513 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00514 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00515 }
00516
00517
00518 if( LineSave.ipass > 0 )
00519 {
00520 LineSv[LineSave.nsum].SumLine[0] = 0.;
00521 LineSv[LineSave.nsum].SumLine[1] = 0.;
00522 }
00523
00524
00525
00526 Cont_nInu = rfield.flux[0][PredCont[i].ip_C()]*renorm*radius.r1r0sq +
00527 rfield.ConRefIncid[0][PredCont[i].ip_C()]*renorm;
00528
00529 # if 0
00530
00531 if( !i )
00532 fprintf(ioQQQ,"\n");
00533 char chWL[1000];
00534 sprt_wl( chWL , (*TauDummy).WLAng() );
00535 fprintf( ioQQQ,"assert line luminosity \"nInu\" %s %.3f\n",
00536 chWL ,
00537 log10(SDIV(Cont_nInu/radius.dVeffAper)) + radius.Conv2PrtInten );
00538 # endif
00539
00540 linadd( Cont_nInu/radius.dVeffAper,(*TauDummy).WLAng(),"nInu",'i',
00541 "transmitted and reflected incident continuum at selected energy points " );
00542
00543
00544
00545 if( KILL_CONT && LineSave.ipass > 0 )
00546 {
00547 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00548 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00549 }
00550
00551
00552 if( LineSave.ipass > 0 )
00553 {
00554 LineSv[LineSave.nsum].SumLine[0] = 0.;
00555 LineSv[LineSave.nsum].SumLine[1] = 0.;
00556 }
00557
00558 linadd( (DiffuseInward+SourceReflected)/radius.dVeffAper,(*TauDummy).WLAng(),"InwT",'i',
00559 "total reflected continuum, total inward emission plus reflected (XXdiffuseXX) total continuum ");
00560
00561 if( KILL_CONT && LineSave.ipass > 0 )
00562 {
00563 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00564 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00565 }
00566
00567
00568 if( LineSave.ipass > 0 )
00569 {
00570 LineSv[LineSave.nsum].SumLine[0] = 0.;
00571 LineSv[LineSave.nsum].SumLine[1] = 0.;
00572 }
00573
00574 linadd(SourceReflected/radius.dVeffAper,(*TauDummy).WLAng(),"InwC",'i',
00575 "reflected incident continuum (only incident) ");
00576
00577 if( KILL_CONT && LineSave.ipass > 0 )
00578 {
00579 LineSv[LineSave.nsum-1].emslin[0] = 0.;
00580 LineSv[LineSave.nsum-1].emslin[1] = 0.;
00581 }
00582 }
00583 }
00584
00585 i = StuffComment( "RRC" );
00586 linadd( 0., (realnum)i , "####", 'i',"radiative recombination continua");
00587
00588
00589 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00590 {
00591 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00592 {
00593 if( nelem < 2 || dense.lgElmtOn[nelem] )
00594 {
00595 char chLabel[5]=" ";
00596 for( long n=0; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
00597 {
00598 if( LineSave.ipass < 0 )
00599
00600 linadd(0.,0.,"dumy",'i',"radiative recombination continuum");
00601 else if( LineSave.ipass == 0 )
00602 {
00603
00604
00605
00606 chIonLbl(chLabel, iso_sp[ipISO][nelem].trans(1,0));
00607 realnum wl = (realnum)(RYDLAM / iso_sp[ipISO][nelem].fb[n].xIsoLevNIonRyd);
00608 wl /= (realnum)RefIndex( 1e8/wl );
00609 linadd( 0. , wl ,chLabel,'i',
00610 "radiative recombination continuum");
00611 }
00612 else
00613 {
00614
00615 linadd(iso_sp[ipISO][nelem].fb[n].RadRecCon,0,"dumy",'i',
00616 "radiative recombination continuum");
00617 }
00618 }
00619 }
00620 }
00621 }
00622
00623
00624
00625
00626 for( long nelem=NISO; nelem < LIMELM; nelem++ )
00627 {
00628
00629
00630 for( long ion=0; ion < nelem-NISO+1; ion++ )
00631 {
00632 if( dense.lgElmtOn[nelem] )
00633 {
00634 if( LineSave.ipass < 0 )
00635
00636 linadd(0.,0.,"dumy",'i',"radiative recombination continuum");
00637 else if( LineSave.ipass == 0 )
00638 {
00639 char chLabel[5];
00640 strcpy( chLabel , elementnames.chElementSym[nelem] );
00641 strcat( chLabel , elementnames.chIonStage[ion]);
00642
00643 realnum wl = (realnum)(RYDLAM / Heavy.Valence_IP_Ryd[nelem][ion]);
00644 wl /= (realnum)RefIndex( 1e8/wl );
00645 linadd( 0. , wl ,chLabel,'i',
00646 "radiative recombination continuum");
00647 }
00648 else
00649 {
00650
00651 linadd(Heavy.RadRecCon[nelem][ion],0,"dumy",'i',
00652 "radiative recombination continuum");
00653 }
00654 }
00655 }
00656 }
00657
00658 return;
00659 }