00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "taulines.h"
00007 #include "hydrogenic.h"
00008 #include "elementnames.h"
00009 #include "phycon.h"
00010 #include "dense.h"
00011 #include "thermal.h"
00012 #include "cooling.h"
00013 #include "iso.h"
00014 #include "rt.h"
00015
00016 STATIC double iso_rad_rec_cooling_approx( long ipISO, long nelem );
00017 STATIC double iso_rad_rec_cooling_extra( long ipISO, long nelem, const double& ThinCoolingSum );
00018
00019
00020 #if defined(__HP_aCC)
00021 # pragma OPT_LEVEL 1
00022 #endif
00023
00024
00025 const bool lgPrintIonizCooling = false;
00026
00027 void iso_cool(
00028
00029 long int ipISO ,
00030
00031 long int nelem)
00032 {
00033 long int ipbig;
00034 double biggest = 0.,
00035 dCdT_all,
00036 edenIonAbund,
00037 CollisIonizatCoolingTotal,
00038 CollisIonizatHeatingTotal,
00039 dCollisIonizatCoolingTotalDT,
00040 HeatExcited,
00041 heat_max,
00042 CollisIonizatCooling,
00043 CollisIonizatHeating,
00044 CollisIonizatCoolingDT,
00045 hlone;
00046
00047 valarray<double> CollisIonizatCoolingArr,
00048 CollisIonizatCoolingDTArr,
00049 SavePhotoHeat,
00050 SaveInducCool;
00051
00052 long int nlo_heat_max , nhi_heat_max;
00053 int coolnum = thermal.ncltot;
00054
00055
00056 char chLabel[NCOLNT_LAB_LEN+1];
00057
00058 DEBUG_ENTRY( "iso_cool()" );
00059
00060 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00061
00062
00063 ASSERT( nelem >= ipISO );
00064 ASSERT( ipISO < NISO );
00065 ASSERT( nelem < LIMELM );
00066
00067
00068 ASSERT( sp->numLevels_local <= sp->numLevels_max );
00069
00070 if( dense.xIonDense[nelem][nelem-ipISO]<=0. ||
00071 !dense.lgElmtOn[nelem] )
00072 {
00073
00074 sp->coll_ion = 0.;
00075 sp->cLya_cool = 0.;
00076 sp->cLyrest_cool = 0.;
00077 sp->cBal_cool = 0.;
00078 sp->cRest_cool = 0.;
00079 sp->xLineTotCool = 0.;
00080 sp->RadRecCool = 0.;
00081 sp->FreeBnd_net_Cool_Rate = 0.;
00082 sp->dLTot = 0.;
00083 sp->RecomInducCool_Rate = 0.;
00084
00085 for( long ipHi=1; ipHi < sp->numLevels_max; ++ipHi )
00086 {
00087 for( long ipLo=0; ipLo < ipHi; ++ipLo )
00088 CollisionZero( sp->trans(ipHi,ipLo).Coll() );
00089 }
00090 return;
00091 }
00092
00093
00094
00095 if( lgPrintIonizCooling )
00096 {
00097 CollisIonizatCoolingArr.resize( sp->numLevels_local );
00098 CollisIonizatCoolingDTArr.resize( sp->numLevels_local );
00099 }
00100 SavePhotoHeat.resize( sp->numLevels_local );
00101 SaveInducCool.resize( sp->numLevels_local );
00102
00103
00104
00105
00106
00107
00108
00109
00110 CollisIonizatCoolingTotal = 0.;
00111 CollisIonizatHeatingTotal = 0.;
00112 dCollisIonizatCoolingTotalDT = 0.;
00113
00114
00115 for( long n=0; n < sp->numLevels_local; ++n )
00116 {
00117
00118 CollisIonizatCooling =
00119 EN1RYD*sp->fb[n].xIsoLevNIonRyd*sp->fb[n].ColIoniz*dense.EdenHCorr*
00120 sp->st[n].Pop();
00121 CollisIonizatCoolingTotal += CollisIonizatCooling;
00122 CollisIonizatHeating =
00123 EN1RYD*sp->fb[n].xIsoLevNIonRyd*sp->fb[n].ColIoniz*dense.EdenHCorr*
00124 sp->fb[n].PopLTE* dense.eden*
00125 dense.xIonDense[nelem][nelem+1-ipISO];
00126 CollisIonizatHeatingTotal += CollisIonizatHeating;
00127
00128 double CollisIonizatDiff = CollisIonizatCooling-CollisIonizatHeating;
00129
00130 CollisIonizatCoolingDT = CollisIonizatDiff*
00131 (sp->fb[n].xIsoLevNIonRyd*TE1RYD/POW2(phycon.te)- thermal.halfte);
00132
00133
00134 dCollisIonizatCoolingTotalDT += CollisIonizatCoolingDT;
00135
00136 if( lgPrintIonizCooling )
00137 {
00138 CollisIonizatCoolingArr[n] = CollisIonizatDiff;
00139 CollisIonizatCoolingDTArr[n] = CollisIonizatCoolingDT;
00140 }
00141 }
00142
00143 double CollisIonizatNetCooling =
00144 CollisIonizatCoolingTotal-CollisIonizatHeatingTotal;
00145
00146
00147 sp->coll_ion = CollisIonizatNetCooling;
00148
00149
00150 thermal.dCooldT += dCollisIonizatCoolingTotalDT;
00151
00152
00153 sprintf(chLabel , "IScion%2s%2s" , elementnames.chElementSym[ipISO] ,
00154 elementnames.chElementSym[nelem] );
00155
00156
00157
00158
00159 if (0)
00160 {
00161
00162
00163 CoolAdd(chLabel,(realnum)nelem,CollisIonizatCoolingTotal);
00164
00165
00166
00167
00168
00169 thermal.heating[0][3] += CollisIonizatHeatingTotal;
00170 }
00171 else
00172 {
00173
00174 CoolAdd(chLabel,(realnum)nelem,MAX2(0.,CollisIonizatNetCooling));
00175 if (CollisIonizatNetCooling < 0)
00176 thermal.heating[0][3] -= CollisIonizatNetCooling;
00177 }
00178
00179
00180 if( lgPrintIonizCooling && nelem==1 && ipISO==1 )
00181 {
00182 fprintf(ioQQQ,"DEBUG coll ioniz cool contributors:");
00183 for( long n=0; n < sp->numLevels_local; n++ )
00184 {
00185 if( CollisIonizatCoolingArr[n] / SDIV( CollisIonizatNetCooling ) > 0.1 )
00186 fprintf(ioQQQ," %2li %.1e",
00187 n,
00188 CollisIonizatCoolingArr[n]/ SDIV( CollisIonizatNetCooling ) );
00189 }
00190 fprintf(ioQQQ,"\n");
00191 fprintf(ioQQQ,"DEBUG coll ioniz derivcontributors:");
00192 for( long n=0; n < sp->numLevels_local; n++ )
00193 {
00194 if( CollisIonizatCoolingDTArr[n] / SDIV( dCollisIonizatCoolingTotalDT ) > 0.1 )
00195 fprintf(ioQQQ," %2li %.1e",
00196 n,
00197 CollisIonizatCoolingDTArr[n]/ SDIV( dCollisIonizatCoolingTotalDT ) );
00198 }
00199 fprintf(ioQQQ,"\n");
00200 }
00201
00202
00203
00204
00205
00206
00207
00208
00209 edenIonAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];
00210
00211 sp->RadRecCool = 0.;
00212 if( dense.gas_phase[nelem] > 1e-3 * dense.xNucleiTotal )
00213 {
00214 RT_iso_integrate_RRC( ipISO, nelem, false );
00215 }
00216 else
00217 {
00218 double ThinCoolingSum = iso_rad_rec_cooling_approx( ipISO, nelem );
00219
00220 sp->RadRecCool = iso_rad_rec_cooling_extra( ipISO, nelem, ThinCoolingSum );
00221 }
00222
00223
00224 for( long n=0; n < sp->numLevels_local; n++ )
00225 {
00226 sp->RadRecCool += sp->fb[n].RadRecCoolCoef;
00227 }
00228
00229
00230 sp->RadRecCool *= edenIonAbund;
00231
00232 {
00233
00234 enum {DEBUG_LOC=false};
00235 if( DEBUG_LOC )
00236 {
00237 if( nelem==ipISO )
00238 {
00239 for( long n=0; n < sp->numLevels_local; n++ )
00240 {
00241 fprintf(ioQQQ,"\t%.2f",sp->fb[n].RadRecCoolCoef/sp->RadRecCool);
00242 }
00243 fprintf(ioQQQ,"\n");
00244 }
00245 }
00246 }
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 HeatExcited = 0.;
00258 ipbig = -1000;
00259 for( long n=1; n < sp->numLevels_local; ++n )
00260 {
00261 ASSERT( sp->fb[n].PhotoHeat >= 0. );
00262 ASSERT( sp->st[n].Pop() >= 0. );
00263
00264 SavePhotoHeat[n] = sp->st[n].Pop()*
00265 sp->fb[n].PhotoHeat;
00266 HeatExcited += SavePhotoHeat[n];
00267 if( SavePhotoHeat[n] > biggest )
00268 {
00269 biggest = SavePhotoHeat[n];
00270 ipbig = (int)n;
00271 }
00272 }
00273 {
00274
00275 enum {DEBUG_LOC=false};
00276 if( DEBUG_LOC && ipISO==0 && nelem==0 && nzone > 700)
00277 {
00278
00279 SavePhotoHeat[ipH1s] = sp->st[ipH1s].Pop()*
00280 sp->fb[ipH1s].PhotoHeat;
00281 fprintf(ioQQQ,"ipISO = %li nelem=%li ipbig=%li biggest=%g HeatExcited=%.2e ctot=%.2e\n",
00282 ipISO , nelem,
00283 ipbig ,
00284 biggest,
00285 HeatExcited ,
00286 thermal.ctot);
00287
00288 for( long n=ipH1s; n< sp->numLevels_local; ++n )
00289 {
00290 fprintf(ioQQQ,"DEBUG phot heat%2li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00291 n,
00292 SavePhotoHeat[n]/HeatExcited,
00293 dense.xIonDense[nelem][nelem+1-ipISO],
00294 sp->st[n].Pop(),
00295 sp->fb[n].PhotoHeat,
00296 sp->fb[n].gamnc);
00297 }
00298 }
00299 }
00300
00301
00302
00303
00304 sp->FreeBnd_net_Cool_Rate = sp->RadRecCool - HeatExcited;
00305
00306
00307
00308
00309
00310
00311
00312 thermal.heating[0][1] += MAX2(0.,-sp->FreeBnd_net_Cool_Rate);
00313
00314
00315
00316 sprintf(chLabel , "ISrcol%2s%2s" , elementnames.chElementSym[ipISO] ,
00317 elementnames.chElementSym[nelem]);
00318 CoolAdd(chLabel, (realnum)nelem, MAX2(0.,sp->FreeBnd_net_Cool_Rate));
00319
00320
00321 thermal.dCooldT += 0.2*sp->FreeBnd_net_Cool_Rate*phycon.teinv;
00322
00323
00324
00325
00326
00327
00328
00329 sp->RecomInducCool_Rate = 0.;
00330
00331 for( long n=0; n < sp->numLevels_local; ++n )
00332 {
00333
00334 SaveInducCool[n] = sp->fb[n].RecomInducCool_Coef*sp->fb[n].PopLTE*edenIonAbund;
00335 sp->RecomInducCool_Rate += SaveInducCool[n];
00336 thermal.dCooldT += SaveInducCool[n]*
00337 (sp->fb[n].xIsoLevNIonRyd/phycon.te_ryd - 1.5)*phycon.teinv;
00338 }
00339
00340 {
00341
00342 enum {DEBUG_LOC=false};
00343 if( DEBUG_LOC && ipISO==0 && nelem==5 )
00344 {
00345 fprintf(ioQQQ," ipISO=%li nelem=%li ctot = %.2e\n",
00346 ipISO,
00347 nelem,
00348 thermal.ctot);
00349 fprintf(ioQQQ,"sum\t%.2e\t%.2e\t%.2e\n",
00350 HeatExcited,
00351 sp->RadRecCool,
00352 sp->RecomInducCool_Rate);
00353 fprintf(ioQQQ,"sum\tp ht\tr cl\ti cl\n");
00354
00355
00356 for( long n=0; n< sp->numLevels_local; ++n )
00357 {
00358 fprintf(ioQQQ,"%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e \n",
00359 n,
00360 SavePhotoHeat[n],
00361 sp->fb[n].RadRecCoolCoef,
00362 SaveInducCool[n] ,
00363 sp->fb[n].RecomInducCool_Coef,
00364 sp->fb[n].PopLTE,
00365 sp->fb[n].RecomInducRate);
00366 }
00367 fprintf(ioQQQ," \n");
00368 }
00369 }
00370
00371 sprintf(chLabel , "ISicol%2s%2s" , elementnames.chElementSym[ipISO] ,
00372 elementnames.chElementSym[nelem]);
00373
00374 CoolAdd(chLabel,(realnum)nelem,sp->RecomInducCool_Rate);
00375
00376
00377 sp->xLineTotCool = 0.;
00378 dCdT_all = 0.;
00379 heat_max = 0.;
00380 nlo_heat_max = -1;
00381 nhi_heat_max = -1;
00382
00383 if (0)
00384 {
00385 long ipHi = 0;
00386 fprintf(ioQQQ,"%ld %ld %ld %g\n",nelem,ipISO,ipHi,
00387 sp->st[ipHi].Pop()/sp->st[ipHi].g());
00388 for( ipHi=1; ipHi < sp->numLevels_local; ++ipHi )
00389 {
00390 fprintf(ioQQQ,"%ld %ld %ld %g\n",nelem,ipISO,ipHi,
00391 sp->st[ipHi].Pop()/(sp->st[ipHi].Boltzmann()*
00392 sp->st[ipHi].g()));
00393 }
00394
00395 }
00396
00397
00398
00399 vector<double> Boltzmann(sp->numLevels_local-1);
00400 for (unsigned lev = 0; lev < Boltzmann.size(); ++lev)
00401 Boltzmann[lev] = 1.0;
00402
00403 for( long ipHi=1; ipHi < sp->numLevels_local; ++ipHi )
00404 {
00405 double arg = (sp->st[ipHi].energy().K()-sp->st[ipHi-1].energy().K()) /
00406 phycon.te;
00407 arg = MAX2( -38., arg); fixit();
00408 double bstep = sexp( arg );
00409 for (long ipLo = 0; ipLo < ipHi; ++ipLo )
00410 Boltzmann[ipLo] *= bstep;
00411
00412 for( long ipLo=0; ipLo < ipHi; ++ipLo )
00413 {
00414
00415 hlone =
00416 sp->trans(ipHi,ipLo).Coll().ColUL( colliders ) *
00417 (sp->st[ipLo].Pop()*
00418 Boltzmann[ipLo]*
00419 sp->st[ipHi].g()/sp->st[ipLo].g() -
00420 sp->st[ipHi].Pop())*
00421 sp->trans(ipHi,ipLo).EnergyErg();
00422
00423 if( hlone > 0. )
00424 sp->trans(ipHi,ipLo).Coll().cool() = hlone;
00425 else
00426 sp->trans(ipHi,ipLo).Coll().heat() = -1.*hlone;
00427
00428 sp->xLineTotCool += hlone;
00429
00430
00431 if( hlone > 0. )
00432 {
00433
00434
00435
00436 dCdT_all += hlone*
00437 (sp->trans(ipHi,ipH1s).EnergyK()*thermal.tsq1 - thermal.halfte);
00438 }
00439 else
00440 {
00441
00442
00443
00444 if( hlone < heat_max )
00445 {
00446 heat_max = hlone;
00447 nlo_heat_max = ipLo;
00448 nhi_heat_max = ipHi;
00449 }
00450 dCdT_all -= hlone*thermal.halfte;
00451 }
00452 }
00453 }
00454 {
00455
00456 enum {DEBUG_LOC=false};
00457 if( DEBUG_LOC )
00458 {
00459 if( nelem==ipISO )
00460 fprintf(ioQQQ,"%li %li %.2e\n", nlo_heat_max, nhi_heat_max, heat_max );
00461 }
00462 }
00463
00464
00465
00466 sp->cLya_cool = 0.;
00467
00468 sp->cLyrest_cool = 0.;
00469
00470 for( long ipHi=1; ipHi < sp->numLevels_local; ipHi++ )
00471 {
00472 hlone = sp->trans(ipHi,ipH1s).Coll().ColUL( colliders ) *
00473 (sp->st[0].Pop()*
00474 sp->st[ipHi].Boltzmann()*
00475 sp->st[ipHi].g()/sp->st[0].g() -
00476 sp->st[ipHi].Pop())*
00477 sp->trans(ipHi,0).EnergyErg();
00478
00479 if( ipHi == iso_ctrl.nLyaLevel[ipISO] )
00480 {
00481 sp->cLya_cool = hlone;
00482
00483 }
00484 else
00485 {
00486
00487 sp->cLyrest_cool += hlone;
00488 }
00489 }
00490
00491
00492
00493 sp->cBal_cool = 0.;
00494 if (nelem == ipHYDROGEN)
00495 {
00496 double boltzH2s =
00497 sexp( (sp->st[2].energy().K()-sp->st[ipH2s].energy().K()) /
00498 phycon.te );
00499 double boltzH2p =
00500 sexp( (sp->st[2].energy().K()-sp->st[ipH2p].energy().K()) /
00501 phycon.te );
00502 for( long ipHi=3; ipHi < sp->numLevels_local; ipHi++ )
00503 {
00504 double bstep =
00505 sexp( (sp->st[ipHi].energy().K()-sp->st[ipHi-1].energy().K()) /
00506 phycon.te );
00507 boltzH2s *= bstep;
00508 boltzH2p *= bstep;
00509
00510
00511 long ipLo = ipH2s;
00512 hlone =
00513 sp->trans(ipHi,ipLo).Coll().ColUL( colliders ) *
00514 (sp->st[ipLo].Pop()*
00515 boltzH2s*
00516 sp->st[ipHi].g()/sp->st[ipLo].g() -
00517 sp->st[ipHi].Pop())*
00518 sp->trans(ipHi,ipLo).EnergyErg();
00519 ASSERT( sp->st[ipHi].energy().Erg() >
00520 sp->st[ipLo].energy().Erg() );
00521
00522 ipLo = ipH2p;
00523 hlone +=
00524 sp->trans(ipHi,ipLo).Coll().ColUL( colliders ) *
00525 (sp->st[ipLo].Pop()*
00526 boltzH2p*
00527 sp->st[ipHi].g()/sp->st[ipLo].g() -
00528 sp->st[ipHi].Pop())*
00529 sp->trans(ipHi,ipLo).EnergyErg();
00530 ASSERT( sp->st[ipHi].energy().Erg() >
00531 sp->st[ipLo].energy().Erg() );
00532
00533
00534 sp->cBal_cool += hlone;
00535 }
00536 }
00537
00538
00539
00540 sp->cRest_cool = 0.;
00541 if (nelem == ipHYDROGEN)
00542 {
00543 for (unsigned lev = 0; lev < Boltzmann.size(); ++lev)
00544 Boltzmann[lev] = 1.;
00545
00546 for( long ipHi=4; ipHi < sp->numLevels_local; ++ipHi )
00547 {
00548 double bstep =
00549 sexp( (sp->st[ipHi].energy().K()-sp->st[ipHi-1].energy().K()) /
00550 phycon.te );
00551 for ( long ipLo = 0; ipLo < ipHi; ++ipLo )
00552 Boltzmann[ipLo] *= bstep;
00553
00554 for( long ipLo=3; ipLo < ipHi; ++ipLo )
00555 {
00556 hlone =
00557 sp->trans(ipHi,ipLo).Coll().ColUL( colliders ) *
00558 (sp->st[ipLo].Pop()*
00559 Boltzmann[ipLo]*
00560 sp->st[ipHi].g()/sp->st[ipLo].g() -
00561 sp->st[ipHi].Pop())*
00562 sp->trans(ipHi,ipLo).EnergyErg();
00563
00564 sp->cRest_cool += hlone;
00565 }
00566 }
00567 }
00568
00569
00570
00571
00572
00573 sprintf(chLabel , "ISclin%2s%2s" , elementnames.chElementSym[ipISO] ,
00574 elementnames.chElementSym[nelem]);
00575 if( sp->xLineTotCool > 0. )
00576 {
00577
00578 CoolAdd(chLabel,(realnum)nelem,sp->xLineTotCool);
00579 thermal.dCooldT += dCdT_all;
00580 sp->dLTot = 0.;
00581 }
00582 else
00583 {
00584
00585 thermal.heating[0][23] -= sp->xLineTotCool;
00586 CoolAdd(chLabel,(realnum)nelem,0.);
00587 sp->dLTot = -dCdT_all;
00588 }
00589
00590 {
00591
00592 enum {DEBUG_LOC=false};
00593 if( DEBUG_LOC )
00594 {
00595 if( nelem == 0 )
00596 {
00597 fprintf(ioQQQ,"%.2e la %.2f restly %.2f barest %.2f hrest %.2f\n",
00598 sp->xLineTotCool ,
00599 sp->cLya_cool/sp->xLineTotCool ,
00600 sp->cLyrest_cool/sp->xLineTotCool ,
00601 sp->cBal_cool/sp->xLineTotCool ,
00602 sp->cRest_cool/sp->xLineTotCool );
00603 }
00604 }
00605 }
00606 {
00607
00608
00609 enum {DEBUG_LOC=false};
00610 enum {LTEPOP=true};
00611
00612 if( DEBUG_LOC && (nelem==1 || nelem==0) )
00613 {
00614
00615 if( LTEPOP )
00616 {
00617
00618 for( long n=ipH1s; n<sp->numLevels_local; ++n )
00619 {
00620 fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
00621 sp->fb[n].gamnc *sp->fb[n].PopLTE,
00622
00623
00624
00625 sp->fb[n].RadRecomb[ipRecRad]+
00626 sp->fb[n].RecomInducRate*sp->fb[n].PopLTE ,
00627
00628 sp->fb[n].RecomInducRate*sp->fb[n].PopLTE ,
00629
00630 sp->fb[n].RadRecomb[ipRecRad] ,
00631
00632 sp->fb[n].PhotoHeat*sp->fb[n].PopLTE,
00633 sp->fb[n].RecomInducCool_Coef*sp->fb[n].PopLTE+sp->fb[n].RadRecCoolCoef,
00634
00635 sp->fb[n].RecomInducCool_Coef*sp->fb[n].PopLTE ,
00636
00637 sp->fb[n].RadRecCoolCoef );
00638 }
00639 }
00640 else
00641 {
00642
00643 for( long n=ipH1s; n<sp->numLevels_local; ++n )
00644 {
00645 fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
00646 sp->fb[n].gamnc*sp->st[n].Pop(),
00647
00648 sp->fb[n].RadRecomb[ipRecRad]*edenIonAbund+
00649 sp->fb[n].RecomInducRate*sp->fb[n].PopLTE *edenIonAbund ,
00650 sp->fb[n].RadRecomb[ipRecRad]*edenIonAbund ,
00651 sp->fb[n].RecomInducRate*sp->fb[n].PopLTE *edenIonAbund ,
00652
00653 SavePhotoHeat[n],
00654 SaveInducCool[n]+sp->fb[n].RadRecCoolCoef*edenIonAbund ,
00655
00656 SaveInducCool[n] ,
00657
00658 sp->fb[n].RadRecCoolCoef*edenIonAbund );
00659 }
00660 }
00661 }
00662 }
00663
00664
00665 for( int i = coolnum; i < thermal.ncltot ; i++ )
00666 thermal.elementcool[nelem] += thermal.cooling[i];
00667
00668 return;
00669 }
00670 #if defined(__HP_aCC)
00671 #pragma OPTIMIZE OFF
00672 #pragma OPTIMIZE ON
00673 #endif
00674
00675 STATIC double iso_rad_rec_cooling_approx( long ipISO, long nelem )
00676 {
00677 DEBUG_ENTRY( "iso_rad_rec_cooling_approx()" );
00678
00679 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00680
00681
00682 double ThinCoolingSum = 0.;
00683
00684
00685 for( long n=0; n < sp->numLevels_local; n++ )
00686 {
00687 double thin = 0.;
00688
00689 if( n==0 && ipISO==ipH_LIKE )
00690 {
00691
00692 thin = HydroRecCool(
00693
00694 1 ,
00695
00696 nelem);
00697 }
00698 else
00699 {
00700
00701 thin = sp->fb[n].RadRecomb[ipRecRad]*
00702
00703
00704 HCoolRatio(
00705 phycon.te * POW2( (double)sp->st[n].n() / (double)(nelem+1-ipISO) ))*
00706
00707 phycon.te * BOLTZMANN;
00708 }
00709
00710
00711 sp->fb[n].RadRecCoolCoef = sp->fb[n].RadRecomb[ipRecNetEsc] * thin;
00712
00713
00714 if( n > 0 )
00715 ThinCoolingSum += thin;
00716 }
00717
00718 return ThinCoolingSum;
00719 }
00720
00721 STATIC double iso_rad_rec_cooling_extra( long ipISO, long nelem, const double& ThinCoolingSum )
00722 {
00723 DEBUG_ENTRY( "iso_rad_rec_cooling_extra()" );
00724
00725 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00726
00727 double RecCoolExtra = 0.;
00728 double ThinCoolingCaseB = 0.;
00729
00730
00731
00732
00733
00734 if( ipISO == ipH_LIKE )
00735 {
00736
00737 if( nelem == 0 )
00738 {
00739
00740 ThinCoolingCaseB = (-25.859117 +
00741 0.16229407*phycon.telogn[0] +
00742 0.34912863*phycon.telogn[1] -
00743 0.10615964*phycon.telogn[2])/(1. +
00744 0.050866793*phycon.telogn[0] -
00745 0.014118924*phycon.telogn[1] +
00746 0.0044980897*phycon.telogn[2] +
00747 6.0969594e-5*phycon.telogn[3]);
00748 }
00749 else
00750 {
00751
00752 ThinCoolingCaseB = (-25.859117 +
00753 0.16229407*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) +
00754 0.34912863*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) -
00755 0.10615964*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) )/(1. +
00756 0.050866793*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) -
00757 0.014118924*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) +
00758 0.0044980897*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) +
00759 6.0969594e-5*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),4) );
00760 }
00761
00762
00763 ThinCoolingCaseB = POW3(1.+nelem-ipISO)*pow(10.,ThinCoolingCaseB)/(phycon.te/POW2(1.+nelem-ipISO) );
00764
00765
00766 RecCoolExtra = ThinCoolingCaseB - ThinCoolingSum;
00767 }
00768 else
00769 {
00770 ThinCoolingCaseB = 0.;
00771 RecCoolExtra = 0.;
00772 }
00773
00774
00775
00776 RecCoolExtra = MAX2(0., RecCoolExtra );
00777
00778 return RecCoolExtra * sp->fb[sp->numLevels_local-1].RadRecomb[ipRecNetEsc];
00779 }
00780