00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #include "physconst.h"
00011 #include "iso.h"
00012 #include "thermal.h"
00013 #include "secondaries.h"
00014 #include "opacity.h"
00015 #include "rfield.h"
00016 #include "ionbal.h"
00017 #include "atmdat.h"
00018 #include "heavy.h"
00019 #include "gammas.h"
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 double GammaBn(
00036
00037 long int ipLoEnr,
00038
00039 long int ipHiEnr,
00040
00041 long int ipOpac,
00042
00043 double thresh,
00044
00045 double *ainduc,
00046
00047 double *rcool,
00048 t_phoHeat *photoHeat)
00049 {
00050 long int i,
00051 ilo,
00052 iup,
00053 limit;
00054 double GamHi,
00055 bnfun_v,
00056 emin,
00057 g,
00058 phisig,
00059 RateInducRec,
00060 RateInducRecCool,
00061 prod;
00062
00063 DEBUG_ENTRY( "GammaBn()" );
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
00075 {
00076 bnfun_v = 0.;
00077 photoHeat->HeatNet = 0.;
00078 photoHeat->HeatLowEnr = 0.;
00079 photoHeat->HeatHiEnr = 0.;
00080 *ainduc = 0.;
00081 *rcool = 0.;
00082 return( bnfun_v );
00083 }
00084
00085 ASSERT( ipLoEnr >= 0 && ipHiEnr >= 0 );
00086
00087
00088
00089
00090 emin = MAX2(thresh,rfield.anu[ipLoEnr-1]);
00091
00092
00093 emin = thresh;
00094
00095 photoHeat->HeatNet = 0.;
00096 g = 0.;
00097 RateInducRec = 0.;
00098 RateInducRecCool = 0.;
00099
00100
00101
00102 i = ipLoEnr;
00103
00104
00105 phisig = (rfield.flux[0][i-1] + rfield.otslin[i-1] +
00106 rfield.ConInterOut[i-1]*rfield.lgOutOnly)*
00107 opac.OpacStack[i-ipLoEnr+ipOpac-1];
00108
00109 g += phisig;
00110 photoHeat->HeatNet += phisig*rfield.anu[i-1];
00111
00112
00113 prod = phisig*rfield.ContBoltz[i-1];
00114 RateInducRec += prod;
00115
00116 RateInducRecCool += prod*(rfield.anu[i-1] - emin);
00117
00118 iup = MIN2(ipHiEnr,rfield.nflux);
00119 limit = MIN2(iup,secondaries.ipSecIon-1);
00120
00121
00122 for( i=ipLoEnr; i < limit; i++ )
00123 {
00124
00125
00126 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
00127
00128 g += phisig;
00129 photoHeat->HeatNet += phisig*rfield.anu[i];
00130
00131
00132 prod = phisig*rfield.ContBoltz[i];
00133
00134 RateInducRec += prod;
00135
00136 RateInducRecCool += prod*(rfield.anu[i] - emin);
00137 }
00138
00139
00140
00141
00142 photoHeat->HeatNet = MAX2(0.,photoHeat->HeatNet - g*thresh);
00143
00144
00145 photoHeat->HeatLowEnr = photoHeat->HeatNet;
00146
00147
00148 photoHeat->HeatHiEnr = 0.;
00149 GamHi = 0.;
00150
00151
00152 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
00153 for( i=ilo-1; i < iup; i++ )
00154 {
00155 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
00156
00157 GamHi += phisig;
00158 photoHeat->HeatHiEnr += phisig*rfield.anu[i];
00159
00160
00161 prod = phisig*rfield.ContBoltz[i];
00162 RateInducRec += prod;
00163
00164 RateInducRecCool += prod*(rfield.anu[i] - emin);
00165 }
00166
00167
00168 bnfun_v = g + GamHi;
00169
00170
00171 photoHeat->HeatHiEnr = photoHeat->HeatHiEnr - GamHi*thresh;
00172
00173
00174 photoHeat->HeatNet += photoHeat->HeatHiEnr*secondaries.HeatEfficPrimary;
00175
00176
00177 photoHeat->HeatNet *= EN1RYD;
00178 photoHeat->HeatHiEnr *= EN1RYD;
00179 photoHeat->HeatLowEnr *= EN1RYD;
00180
00181
00182 if( rfield.lgInducProcess )
00183 {
00184 *rcool = RateInducRecCool*EN1RYD;
00185 *ainduc = RateInducRec;
00186 }
00187 else
00188 {
00189
00190 *rcool = 0.;
00191 *ainduc = 0.;
00192 }
00193
00194 ASSERT( bnfun_v >= 0. );
00195 ASSERT( photoHeat->HeatNet>= 0. );
00196 return( bnfun_v );
00197 }
00198
00199
00200
00201 void GammaPrtShells( long nelem , long ion )
00202 {
00203 double sum;
00204 long int ns;
00205
00206 DEBUG_ENTRY( "GammaPrtShells()" );
00207
00208 fprintf(ioQQQ," GammaPrtShells nz\t%.2f \t%.2li %.2li ",
00209 fnzone ,
00210 nelem,
00211 ion
00212 );
00213 sum = 0;
00214 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
00215 {
00216 sum += ionbal.PhotoRate_Shell[nelem][ion][ns][0];
00217 }
00218 fprintf(ioQQQ,"\ttot\t%.2e", sum);
00219
00220 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
00221 {
00222 fprintf(ioQQQ,"\t%.2e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
00223 # if 0
00224
00225 if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
00226 {
00227
00228 iplow = opac.ipElement[nelem][ion][ns][0];
00229 iphi = opac.ipElement[nelem][ion][ns][1];
00230 ipop = opac.ipElement[nelem][ion][ns][2];
00231
00232
00233
00234 ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
00235 GammaK(iplow,iphi,
00236 ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0),
00237 &phoHeat)*ionbal.lgPhotoIoniz_On;
00238 }
00239 # endif
00240 }
00241 fprintf(ioQQQ,"\n");
00242 return;
00243 }
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 void GammaPrt(
00254
00255 long int ipLoEnr,
00256 long int ipHiEnr,
00257 long int ipOpac,
00258
00259 FILE * ioFILE,
00260
00261 double total,
00262
00263 double threshold)
00264 {
00265 long int i,
00266 j,
00267 k;
00268 double flxcor,
00269 phisig;
00270
00271 DEBUG_ENTRY( "GammaPrt()" );
00272
00273
00274
00275
00276
00277 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
00278 {
00279 return;
00280 }
00281
00282 fprintf( ioFILE, " GammaPrt %.2f from ",fnzone);
00283 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipLoEnr-1]));
00284 fprintf( ioFILE, " to ");
00285 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipHiEnr-1]));
00286 fprintf( ioFILE, "R rates >");
00287 fprintf( ioFILE,PrintEfmt("%9.2e",threshold));
00288 fprintf( ioFILE, " of total=");
00289 fprintf( ioFILE,PrintEfmt("%9.2e",total));
00290 fprintf( ioFILE, " (frac inc, otslin, otscon, ConInterOut, outlin ConOTS_local_OTS_rate ) chL, C\n");
00291
00292 if( threshold <= 0. || total <= 0. )
00293 {
00294 return;
00295 }
00296
00297 k = ipLoEnr;
00298 j = MIN2(ipHiEnr,rfield.nflux);
00299
00300
00301 i = k-1;
00302
00303
00304 phisig = (rfield.flux[0][i] + rfield.otslin[i]+ rfield.ConInterOut[i]*rfield.lgOutOnly)*
00305 opac.OpacStack[i-k+ipOpac];
00306 if( phisig > threshold || phisig < 0.)
00307 {
00308 flxcor = rfield.flux[0][i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly;
00309
00310
00311 fprintf( ioFILE, "[%5ld]" , i );
00312 fprintf( ioFILE, PrintEfmt("%9.2e",rfield.anu[i]));
00313 fprintf( ioFILE, PrintEfmt("%9.2e",phisig/total));
00314 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
00315 rfield.flux[0][i]/SDIV(flxcor),
00316 rfield.otslin[i]/SDIV(flxcor),
00317
00318 0./SDIV(flxcor),
00319 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
00320 (rfield.outlin[0][i]+rfield.outlin_noplot[i])/SDIV(flxcor),
00321 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
00322 rfield.chLineLabel[i],
00323 rfield.chContLabel[i],
00324 opac.OpacStack[i-k+ipOpac]);
00325 }
00326
00327 for( i=k; i < j; i++ )
00328 {
00329 phisig = rfield.SummedCon[i]*opac.OpacStack[i-k+ipOpac];
00330 if( phisig > threshold || phisig < 0.)
00331 {
00332
00333 flxcor = rfield.flux[0][i] + rfield.otslin[i] + rfield.otscon[i] +
00334 rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i]*rfield.lgOutOnly;
00335
00336 fprintf( ioFILE, "%5ld", i );
00337 fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i]));
00338 fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total));
00339 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
00340 rfield.flux[0][i]/SDIV(flxcor),
00341 rfield.otslin[i]/SDIV(flxcor),
00342 rfield.otscon[i]/SDIV(flxcor),
00343 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
00344 (rfield.outlin[0][i] + rfield.outlin_noplot[i])/SDIV(flxcor),
00345 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
00346 rfield.chLineLabel[i],
00347 rfield.chContLabel[i],
00348 opac.OpacStack[i-k+ipOpac]);
00349 }
00350 }
00351 return;
00352 }
00353
00354
00355
00356
00357
00358
00359 double GammaK(
00360
00361
00362 long int ipLoEnr,
00363 long int ipHiEnr,
00364
00365 long int ipOpac,
00366
00367
00368 double yield1,
00369 t_phoHeat *photoHeat)
00370 {
00371 long int i,
00372 ilo,
00373 ipOffset,
00374 iup,
00375 limit;
00376 double GamHi,
00377 eauger;
00378
00379 double gamk_v ,
00380 phisig;
00381
00382 DEBUG_ENTRY( "GammaK()" );
00383
00384
00385
00386
00387
00388
00389
00390 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr)
00391 {
00392 gamk_v = 0.;
00393 photoHeat->HeatNet = 0.;
00394 photoHeat->HeatHiEnr = 0.;
00395 photoHeat->HeatLowEnr = 0.;
00396 return( gamk_v );
00397 }
00398
00399 iup = MIN2(ipHiEnr,rfield.nflux);
00400
00401
00402
00403
00404 eauger = rfield.anu[ipLoEnr-1]*yield1;
00405
00406
00407
00408 gamk_v = 0.;
00409 photoHeat->HeatNet = 0.;
00410
00411
00412 ipOffset = ipOpac - ipLoEnr;
00413
00414
00415
00416
00417
00418
00419 phisig = (rfield.flux[0][ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+
00420 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly) * opac.OpacStack[ipOpac-1];
00421 gamk_v += phisig;
00422 photoHeat->HeatNet += phisig*rfield.anu[ipLoEnr-1];
00423
00424
00425 limit = MIN2(iup,secondaries.ipSecIon-1);
00426 for( i=ipLoEnr; i < limit; i++ )
00427 {
00428 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
00429 gamk_v += phisig;
00430 photoHeat->HeatNet += phisig*rfield.anu[i];
00431 }
00432
00433
00434
00435
00436 ASSERT( photoHeat->HeatNet >= 0. );
00437 photoHeat->HeatNet = MAX2(0.,photoHeat->HeatNet - gamk_v*eauger);
00438
00439
00440 photoHeat->HeatLowEnr = photoHeat->HeatNet;
00441
00442
00443 photoHeat->HeatHiEnr = 0.;
00444 GamHi = 0.;
00445
00446 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
00447 for( i=ilo-1; i < iup; i++ )
00448 {
00449
00450 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
00451 GamHi += phisig;
00452 photoHeat->HeatHiEnr += phisig*rfield.anu[i];
00453 }
00454
00455
00456 gamk_v += GamHi;
00457
00458 photoHeat->HeatHiEnr -= GamHi*eauger;
00459
00460
00461
00462 photoHeat->HeatNet += photoHeat->HeatHiEnr*secondaries.HeatEfficPrimary;
00463
00464
00465 photoHeat->HeatNet *= EN1RYD;
00466 photoHeat->HeatLowEnr *= EN1RYD;
00467 photoHeat->HeatHiEnr *= EN1RYD;
00468
00469 ASSERT( gamk_v >= 0. );
00470 ASSERT( photoHeat->HeatNet>= 0. );
00471 return( gamk_v );
00472 }
00473
00474
00475 void GammaPrtRate(
00476
00477 FILE * ioFILE,
00478
00479 long int ion ,
00480
00481 long int nelem,
00482
00483 bool lgPRT )
00484 {
00485 long int nshell , ns;
00486
00487 DEBUG_ENTRY( "GammaPrtRate()" );
00488
00489
00490 nshell = Heavy.nsShells[nelem][ion];
00491
00492
00493 fprintf(ioFILE , "GammaPrtRate: %li %li",ion , nelem );
00494 for( ns=nshell-1; ns>=0; --ns )
00495 {
00496 fprintf(ioFILE , " %.2e" , ionbal.PhotoRate_Shell[nelem][ion][ns][0]);
00497
00498
00499 if( lgPRT )
00500 {
00501 fprintf(ioFILE , "\n");
00502 GammaPrt(
00503
00504 opac.ipElement[nelem][ion][ns][0],
00505 opac.ipElement[nelem][ion][ns][1],
00506 opac.ipElement[nelem][ion][ns][2],
00507
00508 ioFILE,
00509
00510 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
00511
00512 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.05);
00513 }
00514 }
00515 fprintf(ioFILE , "\n");
00516 return;
00517 }