00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "cddefines.h"
00012 #include "physconst.h"
00013 #include "iso.h"
00014 #include "thermal.h"
00015 #include "secondaries.h"
00016 #include "opacity.h"
00017 #include "rfield.h"
00018 #include "ionbal.h"
00019 #include "atmdat.h"
00020 #include "heavy.h"
00021 #include "gammas.h"
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 double GammaBn(
00033
00034 long int ipLoEnr,
00035
00036 long int ipHiEnr,
00037
00038 long int ipOpac,
00039
00040 double thresh,
00041
00042 double *ainduc,
00043
00044 double *rcool)
00045 {
00046 long int i,
00047 ilo,
00048 iup,
00049 limit;
00050 double GamHi,
00051 bnfun_v,
00052 emin,
00053 g,
00054 phisig,
00055 RateInducRec,
00056 RateInducRecCool,
00057 prod;
00058
00059 DEBUG_ENTRY( "GammaBn()" );
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
00071 {
00072 bnfun_v = 0.;
00073 thermal.HeatNet = 0.;
00074 thermal.HeatLowEnr = 0.;
00075 thermal.HeatHiEnr = 0.;
00076 *ainduc = 0.;
00077 *rcool = 0.;
00078 return( bnfun_v );
00079 }
00080
00081 ASSERT( ipLoEnr >= 0 && ipHiEnr >= 0 );
00082
00083
00084
00085
00086 emin = MAX2(thresh,rfield.anu[ipLoEnr-1]);
00087
00088
00089 emin = thresh;
00090
00091 thermal.HeatNet = 0.;
00092 g = 0.;
00093 RateInducRec = 0.;
00094 RateInducRecCool = 0.;
00095
00096
00097
00098 i = ipLoEnr;
00099
00100
00101 phisig = (rfield.flux[i-1] + rfield.otslin[i-1] +
00102 rfield.ConInterOut[i-1]*rfield.lgOutOnly)*
00103 opac.OpacStack[i-ipLoEnr+ipOpac-1];
00104
00105 g += phisig;
00106 thermal.HeatNet += phisig*rfield.anu[i-1];
00107
00108
00109 prod = phisig*rfield.ContBoltz[i-1];
00110 RateInducRec += prod;
00111
00112 RateInducRecCool += prod*(rfield.anu[i-1] - emin);
00113
00114 iup = MIN2(ipHiEnr,rfield.nflux);
00115 limit = MIN2(iup,secondaries.ipSecIon-1);
00116
00117
00118 for( i=ipLoEnr; i < limit; i++ )
00119 {
00120
00121
00122 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
00123
00124 g += phisig;
00125 thermal.HeatNet += phisig*rfield.anu[i];
00126
00127
00128 prod = phisig*rfield.ContBoltz[i];
00129
00130 RateInducRec += prod;
00131
00132 RateInducRecCool += prod*(rfield.anu[i] - emin);
00133 }
00134
00135
00136
00137
00138 thermal.HeatNet = MAX2(0.,thermal.HeatNet - g*thresh);
00139
00140
00141 thermal.HeatLowEnr = thermal.HeatNet;
00142
00143
00144 thermal.HeatHiEnr = 0.;
00145 GamHi = 0.;
00146
00147
00148 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
00149 for( i=ilo-1; i < iup; i++ )
00150 {
00151 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
00152
00153 GamHi += phisig;
00154 thermal.HeatHiEnr += phisig*rfield.anu[i];
00155
00156
00157 prod = phisig*rfield.ContBoltz[i];
00158 RateInducRec += prod;
00159
00160 RateInducRecCool += prod*(rfield.anu[i] - emin);
00161 }
00162
00163
00164 bnfun_v = g + GamHi;
00165
00166
00167 thermal.HeatHiEnr = thermal.HeatHiEnr - GamHi*thresh;
00168
00169
00170 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;
00171
00172
00173 thermal.HeatNet *= EN1RYD;
00174 thermal.HeatHiEnr *= EN1RYD;
00175 thermal.HeatLowEnr *= EN1RYD;
00176
00177
00178 if( rfield.lgInducProcess )
00179 {
00180 *rcool = RateInducRecCool*EN1RYD;
00181 *ainduc = RateInducRec;
00182 }
00183 else
00184 {
00185
00186 *rcool = 0.;
00187 *ainduc = 0.;
00188 }
00189
00190 ASSERT( bnfun_v >= 0. );
00191 ASSERT( thermal.HeatNet>= 0. );
00192 return( bnfun_v );
00193 }
00194
00195
00196
00197 void GammaPrtShells( long nelem , long ion )
00198 {
00199 double sum;
00200 long int ns;
00201
00202 DEBUG_ENTRY( "GammaPrtShells()" );
00203
00204 fprintf(ioQQQ," GammaPrtShells nz\t%.2f \t%.2li %.2li ",
00205 fnzone ,
00206 nelem,
00207 ion
00208 );
00209 sum = 0;
00210 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
00211 {
00212 sum += ionbal.PhotoRate_Shell[nelem][ion][ns][0];
00213 }
00214 fprintf(ioQQQ,"\ttot\t%.2e", sum);
00215
00216 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
00217 {
00218 fprintf(ioQQQ,"\t%.2e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
00219 # if 0
00220
00221 if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
00222 {
00223
00224 iplow = opac.ipElement[nelem][ion][ns][0];
00225 iphi = opac.ipElement[nelem][ion][ns][1];
00226 ipop = opac.ipElement[nelem][ion][ns][2];
00227
00228
00229
00230 ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
00231 GammaK(iplow,iphi,
00232 ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0))*ionbal.lgPhotoIoniz_On;
00233 }
00234 # endif
00235 }
00236 fprintf(ioQQQ,"\n");
00237 return;
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 void GammaPrt(
00249
00250 long int ipLoEnr,
00251 long int ipHiEnr,
00252 long int ipOpac,
00253
00254 FILE * ioFILE,
00255
00256 double total,
00257
00258 double threshold)
00259 {
00260 long int i,
00261 j,
00262 k;
00263 double flxcor,
00264 phisig;
00265
00266 DEBUG_ENTRY( "GammaPrt()" );
00267
00268
00269
00270
00271
00272 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
00273 {
00274 return;
00275 }
00276
00277 fprintf( ioFILE, " GammaPrt %.2f from ",fnzone);
00278 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipLoEnr-1]));
00279 fprintf( ioFILE, " to ");
00280 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipHiEnr-1]));
00281 fprintf( ioFILE, "R rates >");
00282 fprintf( ioFILE,PrintEfmt("%9.2e",threshold));
00283 fprintf( ioFILE, " of total=");
00284 fprintf( ioFILE,PrintEfmt("%9.2e",total));
00285 fprintf( ioFILE, " (frac inc, otslin, otscon, ConInterOut, outlin ConOTS_local_OTS_rate ) chL, C\n");
00286
00287 if( threshold <= 0. || total <= 0. )
00288 {
00289 return;
00290 }
00291
00292 k = ipLoEnr;
00293 j = MIN2(ipHiEnr,rfield.nflux);
00294
00295
00296 i = k-1;
00297
00298
00299 phisig = (rfield.flux[i] + rfield.otslin[i]+ rfield.ConInterOut[i]*rfield.lgOutOnly)*
00300 opac.OpacStack[i-k+ipOpac];
00301 if( phisig > threshold || phisig < 0.)
00302 {
00303 flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly;
00304
00305
00306 fprintf( ioFILE, "[%5ld]" , i );
00307 fprintf( ioFILE, PrintEfmt("%9.2e",rfield.anu[i]));
00308 fprintf( ioFILE, PrintEfmt("%9.2e",phisig/total));
00309 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
00310 rfield.flux[i]/SDIV(flxcor),
00311 rfield.otslin[i]/SDIV(flxcor),
00312
00313 0./SDIV(flxcor),
00314 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
00315 (rfield.outlin[i]+rfield.outlin_noplot[i])/SDIV(flxcor),
00316 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
00317 rfield.chLineLabel[i],
00318 rfield.chContLabel[i],
00319 opac.OpacStack[i-k+ipOpac]);
00320 }
00321
00322 for( i=k; i < j; i++ )
00323 {
00324 phisig = rfield.SummedCon[i]*opac.OpacStack[i-k+ipOpac];
00325 if( phisig > threshold || phisig < 0.)
00326 {
00327
00328 flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.otscon[i] +
00329 rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i]*rfield.lgOutOnly;
00330
00331 fprintf( ioFILE, "%5ld", i );
00332 fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i]));
00333 fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total));
00334 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
00335 rfield.flux[i]/SDIV(flxcor),
00336 rfield.otslin[i]/SDIV(flxcor),
00337 rfield.otscon[i]/SDIV(flxcor),
00338 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
00339 (rfield.outlin[i] + rfield.outlin_noplot[i])/SDIV(flxcor),
00340 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
00341 rfield.chLineLabel[i],
00342 rfield.chContLabel[i],
00343 opac.OpacStack[i-k+ipOpac]);
00344 }
00345 }
00346 return;
00347 }
00348
00349
00350
00351
00352
00353
00354 double GammaK(
00355
00356
00357 long int ipLoEnr,
00358 long int ipHiEnr,
00359
00360 long int ipOpac,
00361
00362
00363 double yield1)
00364 {
00365 long int i,
00366 ilo,
00367 ipOffset,
00368 iup,
00369 limit;
00370 double GamHi,
00371 eauger;
00372
00373 double gamk_v ,
00374 phisig;
00375
00376 DEBUG_ENTRY( "GammaK()" );
00377
00378
00379
00380
00381
00382
00383
00384 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr)
00385 {
00386 gamk_v = 0.;
00387 thermal.HeatNet = 0.;
00388 thermal.HeatHiEnr = 0.;
00389 thermal.HeatLowEnr = 0.;
00390 return( gamk_v );
00391 }
00392
00393 iup = MIN2(ipHiEnr,rfield.nflux);
00394
00395
00396
00397
00398 eauger = rfield.anu[ipLoEnr-1]*yield1;
00399
00400
00401
00402 gamk_v = 0.;
00403 thermal.HeatNet = 0.;
00404
00405
00406 ipOffset = ipOpac - ipLoEnr;
00407
00408
00409
00410
00411
00412
00413 phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+
00414 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly) * opac.OpacStack[ipOpac-1];
00415 gamk_v += phisig;
00416 thermal.HeatNet += phisig*rfield.anu[ipLoEnr-1];
00417
00418
00419 limit = MIN2(iup,secondaries.ipSecIon-1);
00420 for( i=ipLoEnr; i < limit; i++ )
00421 {
00422 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
00423 gamk_v += phisig;
00424 thermal.HeatNet += phisig*rfield.anu[i];
00425 }
00426
00427
00428
00429
00430 ASSERT( thermal.HeatNet >= 0. );
00431 thermal.HeatNet = MAX2(0.,thermal.HeatNet - gamk_v*eauger);
00432
00433
00434 thermal.HeatLowEnr = thermal.HeatNet;
00435
00436
00437 thermal.HeatHiEnr = 0.;
00438 GamHi = 0.;
00439
00440 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
00441 for( i=ilo-1; i < iup; i++ )
00442 {
00443
00444 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
00445 GamHi += phisig;
00446 thermal.HeatHiEnr += phisig*rfield.anu[i];
00447 }
00448
00449
00450 gamk_v += GamHi;
00451
00452 thermal.HeatHiEnr -= GamHi*eauger;
00453
00454
00455
00456 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;
00457
00458
00459 thermal.HeatNet *= EN1RYD;
00460 thermal.HeatLowEnr *= EN1RYD;
00461 thermal.HeatHiEnr *= EN1RYD;
00462
00463 ASSERT( gamk_v >= 0. );
00464 ASSERT( thermal.HeatNet>= 0. );
00465 return( gamk_v );
00466 }
00467
00468 #if 0
00469
00470
00471
00472 double GammaPL(
00473
00474 long int n,
00475
00476 long int nelem)
00477 {
00478 long int i,
00479 ilo,
00480 iup,
00481 limit,
00482 ipLoEnr,
00483 ipHiEnr;
00484 realnum GamHi,
00485 csThresh,
00486 gamPL_v,
00487 phisig;
00488
00489 DEBUG_ENTRY( "GammaPL()" );
00490
00491
00492
00493
00494
00495
00496
00497 ASSERT( nelem >= 0);
00498 ASSERT( nelem < LIMELM);
00499 ASSERT( n >= 0);
00500
00501
00502 ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n];
00503 ASSERT( ipLoEnr>0 );
00504
00505 if( ipLoEnr >= rfield.nflux )
00506 {
00507 gamPL_v = 0.;
00508 thermal.HeatNet = 0.;
00509 thermal.HeatHiEnr = 0.;
00510 thermal.HeatLowEnr = 0.;
00511 return( gamPL_v );
00512 }
00513
00514 ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
00515 iup = MIN2(ipHiEnr,rfield.nflux);
00516
00517 csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1]/ POW2(nelem+1.f);
00518
00519
00520
00521 gamPL_v = 0.;
00522 thermal.HeatNet = 0.;
00523
00524
00525
00526
00527 phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1] +
00528 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly ) /rfield.anu3[ipLoEnr-1];
00529 gamPL_v += phisig;
00530
00531
00532 limit = MIN2(iup,secondaries.ipSecIon-1);
00533 for( i=ipLoEnr; i < limit; i++ )
00534 {
00535 phisig = rfield.SummedCon[i]/rfield.anu3[i];
00536 gamPL_v += phisig;
00537 thermal.HeatNet += phisig*rfield.anu[i];
00538 }
00539
00540
00541 gamPL_v *= csThresh;
00542
00543
00544 thermal.HeatNet *= csThresh;
00545 thermal.HeatNet += -gamPL_v*rfield.anu[ipLoEnr-1];
00546 thermal.HeatLowEnr = thermal.HeatNet;
00547
00548
00549 thermal.HeatHiEnr = 0.;
00550 GamHi = 0.;
00551
00552
00553 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
00554 for( i=ilo-1; i < iup; i++ )
00555 {
00556
00557 phisig = rfield.SummedCon[i]/rfield.anu3[i];
00558 GamHi += phisig;
00559 thermal.HeatHiEnr += phisig*rfield.anu[i];
00560 }
00561
00562
00563 gamPL_v += GamHi*csThresh;
00564
00565
00566 thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*rfield.anu[ipLoEnr-1])*csThresh;
00567
00568
00569 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;
00570
00571
00572 thermal.HeatNet *= EN1RYD;
00573 thermal.HeatHiEnr *= EN1RYD;
00574 thermal.HeatLowEnr *= EN1RYD;
00575
00576 ASSERT( gamPL_v>= 0.);
00577 return( gamPL_v );
00578 }
00579 #endif
00580
00581
00582
00583 double GammaBnPL(long int n,
00584 long int nelem,
00585 double *ainduc,
00586 double *rcool)
00587 {
00588 long int i,
00589 ilo,
00590 iup,
00591 limit,
00592 ipLoEnr,
00593 ipHiEnr;
00594 double GamHi,
00595 bnfunPL_v,
00596 csThresh,
00597 emin,
00598 g,
00599 phisig,
00600 RateInducRec,
00601 RateInducRecCool,
00602 prod,
00603 thresh;
00604
00605 DEBUG_ENTRY( "GammaBnPL()" );
00606
00607
00608 ASSERT( nelem >= 0);
00609 ASSERT( nelem < LIMELM );
00610
00611
00612
00613
00614
00615 ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n];
00616 ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00617
00618 ilo = ipLoEnr;
00619 iup = MIN2(ipHiEnr,rfield.nflux);
00620 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
00621 {
00622 bnfunPL_v = 0.;
00623 thermal.HeatNet = 0.;
00624 thermal.HeatHiEnr = 0.;
00625 thermal.HeatLowEnr = 0.;
00626 *ainduc = 0.;
00627 *rcool = 0.;
00628 return( bnfunPL_v );
00629 }
00630
00631
00632
00633
00634 thresh = iso.xIsoLevNIonRyd[ipH_LIKE][nelem][n];
00635 emin = MAX2(thresh,rfield.anu[ipLoEnr-1]);
00636
00637
00638 emin = thresh;
00639
00640 thermal.HeatNet = 0.;
00641 g = 0.;
00642 RateInducRec = 0.;
00643 RateInducRecCool = 0.;
00644
00645
00646 csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1] / POW2(nelem+1.f);
00647
00648
00649
00650 phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+
00651 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly)/
00652 rfield.anu3[ipLoEnr-1];
00653
00654 g += phisig;
00655 thermal.HeatNet += phisig*rfield.anu[ipLoEnr-1];
00656
00657
00658 prod = phisig*rfield.ContBoltz[ipLoEnr-1];
00659 RateInducRec += prod;
00660
00661 RateInducRecCool += prod*(rfield.anu[ipLoEnr-1] - emin);
00662
00663 limit = MIN2(iup,secondaries.ipSecIon-1);
00664 for( i=ipLoEnr; i < limit; i++ )
00665 {
00666 phisig = rfield.SummedCon[i]/rfield.anu3[i];
00667
00668 g += phisig;
00669 thermal.HeatNet += phisig*rfield.anu[i];
00670
00671
00672 prod = phisig*rfield.ContBoltz[i];
00673 RateInducRec += prod;
00674
00675 RateInducRecCool += prod*(rfield.anu[i] - emin);
00676 }
00677
00678
00679 thermal.HeatNet = MAX2(0.,(thermal.HeatNet - g*thresh)*csThresh);
00680 thermal.HeatLowEnr = thermal.HeatNet;
00681
00682
00683 thermal.HeatHiEnr = 0.;
00684 GamHi = 0.;
00685
00686 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
00687 for( i=ilo-1; i < iup; i++ )
00688 {
00689 phisig = rfield.SummedCon[i]/rfield.anu3[i];
00690 GamHi += phisig;
00691 thermal.HeatHiEnr += phisig*rfield.anu[i];
00692
00693
00694 prod = phisig*rfield.ContBoltz[i];
00695 RateInducRec += prod;
00696
00697 RateInducRecCool += prod*(rfield.anu[i] - emin);
00698 }
00699 RateInducRec *= csThresh;
00700 RateInducRecCool *= csThresh;
00701
00702
00703 bnfunPL_v = (g + GamHi)*csThresh;
00704
00705
00706 thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*thresh)*csThresh;
00707
00708
00709 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;
00710
00711
00712 thermal.HeatNet *= EN1RYD;
00713 thermal.HeatHiEnr *= EN1RYD;
00714 thermal.HeatLowEnr *= EN1RYD;
00715
00716
00717 if( rfield.lgInducProcess )
00718 {
00719 *rcool = RateInducRecCool*EN1RYD;
00720 *ainduc = RateInducRec;
00721 }
00722 else
00723 {
00724
00725 *rcool = 0.;
00726 *ainduc = 0.;
00727 }
00728
00729 ASSERT( bnfunPL_v>= 0.);
00730 return( bnfunPL_v );
00731 }
00732
00733
00734 void GammaPrtRate(
00735
00736 FILE * ioFILE,
00737
00738 long int ion ,
00739
00740 long int nelem,
00741
00742 bool lgPRT )
00743 {
00744 long int nshell , ns;
00745
00746 DEBUG_ENTRY( "GammaPrtRate()" );
00747
00748
00749 nshell = Heavy.nsShells[nelem][ion];
00750
00751
00752 fprintf(ioFILE , "GammaPrtRate: %li %li",ion , nelem );
00753 for( ns=nshell-1; ns>=0; --ns )
00754 {
00755 fprintf(ioFILE , " %.2e" , ionbal.PhotoRate_Shell[nelem][ion][ns][0]);
00756
00757
00758 if( lgPRT )
00759 {
00760 fprintf(ioFILE , "\n");
00761 GammaPrt(
00762
00763 opac.ipElement[nelem][ion][ns][0],
00764 opac.ipElement[nelem][ion][ns][1],
00765 opac.ipElement[nelem][ion][ns][2],
00766
00767 ioFILE,
00768
00769 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
00770
00771 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.05);
00772 }
00773 }
00774 fprintf(ioFILE , "\n");
00775 return;
00776 }