00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "atmdat.h"
00007 #include "rfield.h"
00008 #include "hmi.h"
00009 #include "trace.h"
00010 #include "conv.h"
00011 #include "ionbal.h"
00012 #include "thermal.h"
00013 #include "phycon.h"
00014 #include "doppvel.h"
00015 #include "taulines.h"
00016 #include "mole.h"
00017 #include "heavy.h"
00018 #include "thirdparty.h"
00019 #include "dense.h"
00020 #include "ipoint.h"
00021 #include "elementnames.h"
00022 #include "grainvar.h"
00023 #include "grains.h"
00024 #include "iso.h"
00025
00026
00027
00028
00029
00030
00034 inline double ASINH(double x)
00035 {
00036 if( abs(x) <= 8.e-3 )
00037 return ((0.075*pow2(x) - 1./6.)*pow2(x) + 1.)*x;
00038 else if( abs(x) <= 1./sqrt(DBL_EPSILON) )
00039 {
00040 if( x < 0. )
00041 return -log(sqrt(1. + pow2(x)) - x);
00042 else
00043 return log(sqrt(1. + pow2(x)) + x);
00044 }
00045 else
00046 {
00047 if( x < 0. )
00048 return -(log(-x)+LN_TWO);
00049 else
00050 return log(x)+LN_TWO;
00051 }
00052 }
00053
00054
00055 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
00056 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
00057
00058 static const long MAGIC_AUGER_DATA = 20060126L;
00059
00060 static const bool INCL_TUNNEL = true;
00061 static const bool NO_TUNNEL = false;
00062
00063 static const bool ALL_STAGES = true;
00064
00065
00066
00067 static long int nCalledGrainDrive;
00068
00069
00070
00071
00072
00073 static const long NTOP = NDEMS/5;
00074
00075
00076
00077 static const double TOLER = CONSERV_TOL/10.;
00078 static const long BRACKET_MAX = 50L;
00079
00080
00081
00082 static const int NCHU = NCHS/3;
00083
00084
00085
00086
00087 static const long CT_LOOP_MAX = 25L;
00088
00089
00090 static const long T_LOOP_MAX = 50L;
00091
00092
00093 static double HEAT_TOLER = DBL_MAX;
00094 static double HEAT_TOLER_BIN = DBL_MAX;
00095 static double CHRG_TOLER = DBL_MAX;
00096
00097
00098
00099
00100
00101
00102 static const double AC0 = 3.e-9;
00103 static const double AC1G = 4.e-8;
00104 static const double AC2G = 7.e-8;
00105
00106
00107 static const double ETILDE = 2.*SQRT2/EVRYD;
00108
00109
00110 static const double THERMCONST = PI4*ELECTRON_MASS*POW2(BOLTZMANN)/POW3(HPLANCK);
00111
00112
00113 static const double STICK_ELEC = 0.5;
00114 static const double STICK_ION = 1.0;
00115
00117 inline double one_elec(long nd)
00118 {
00119 return ELEM_CHARGE/EVRYD/gv.bin[nd]->Capacity;
00120 }
00121
00123 inline double pot2chrg(double x,
00124 long nd)
00125 {
00126 return x/one_elec(nd) - 1.;
00127 }
00128
00130 inline double chrg2pot(double x,
00131 long nd)
00132 {
00133 return (x+1.)*one_elec(nd);
00134 }
00135
00137 inline double elec_esc_length(double e,
00138 long nd)
00139 {
00140
00141 if( e <= gv.bin[nd]->le_thres )
00142 return 1.e-7;
00143 else
00144 return 3.e-6*gv.bin[nd]->eec*sqrt(pow3(e*EVRYD*1.e-3));
00145 }
00146
00147
00148 STATIC void ReadAugerData();
00149
00150 STATIC void InitBinAugerData(size_t,long,long);
00151
00152 STATIC void GetNextLine(const char*, FILE*, char[]);
00153
00154 STATIC void InitEmissivities(void);
00155
00156 STATIC double PlanckIntegral(double,size_t,long);
00157
00158 STATIC void NewChargeData(long);
00159
00160 STATIC double GrnStdDpth(long);
00161
00162 STATIC void GrainChargeTemp(void);
00163
00164 STATIC void GrainCharge(size_t,double*);
00165
00166 STATIC double GrainElecRecomb1(size_t,long,double*,double*);
00167
00168 STATIC double GrainElecEmis1(size_t,long,double*,double*,double*,double*);
00169
00170
00171 STATIC void GrainScreen(long,size_t,long,double*,double*);
00172
00173 STATIC double ThetaNu(double);
00174
00175 STATIC void UpdatePot(size_t,long,long,double[],double[]);
00176
00177 STATIC void GetFracPop(size_t,long,double[],double[],long*);
00178
00179 STATIC void UpdatePot1(size_t,long,long,long);
00180
00181 STATIC void UpdatePot2(size_t,long);
00182
00183 inline void Yfunc(long,long,double,double,double,double,double,double*,double*,
00184 double*,double*);
00185
00186 STATIC double y0b(size_t,long,long);
00187
00188 STATIC double y0b01(size_t,long,long);
00189
00190 STATIC double y0psa(size_t,long,long,double);
00191
00192 STATIC double y1psa(size_t,long,double);
00193
00194 inline double y2pa(double,double,long,double*);
00195
00196 inline double y2s(double,double,long,double*);
00197
00198 STATIC long HighestIonStage(void);
00199
00200 STATIC void UpdateRecomZ0(size_t,long,bool);
00201
00202 STATIC void GetPotValues(size_t,long,double*,double*,double*,
00203 double*,double*,double*,bool);
00204
00205
00206
00207 STATIC void GrainIonColl(size_t,long,long,long,const double[],const double[],long*,
00208 realnum*,realnum*);
00209
00210 STATIC void GrainChrgTransferRates(long);
00211
00212 STATIC void GrainUpdateRadius1(void);
00213
00214 STATIC void GrainUpdateRadius2();
00215
00216 STATIC void GrainTemperature(size_t,realnum*,double*,double*,
00217 double*);
00218
00219 STATIC void PE_init(size_t,long,long,double*,double*,double*,
00220 double*,double*,double*,double*);
00221
00222 STATIC void GrainCollHeating(size_t,realnum*,realnum*);
00223
00224 STATIC double GrnVryDpth(size_t);
00225
00226
00227 void AEInfo::p_clear0()
00228 {
00229 nData.clear();
00230 IonThres.clear();
00231 AvNumber.clear();
00232 Energy.clear();
00233 }
00234
00235 void AEInfo::p_clear1()
00236 {
00237 nSubShell = 0;
00238 }
00239
00240 void ShellData::p_clear0()
00241 {
00242 p.clear();
00243 y01.clear();
00244 AvNr.clear();
00245 Ener.clear();
00246 y01A.clear();
00247 }
00248
00249 void ShellData::p_clear1()
00250 {
00251 nelem = LONG_MIN;
00252 ns = LONG_MIN;
00253 ionPot = -DBL_MAX;
00254 ipLo = LONG_MIN;
00255 nData = 0;
00256 }
00257
00258 void ChargeBin::p_clear0()
00259 {
00260 yhat.clear();
00261 yhat_primary.clear();
00262 ehat.clear();
00263 cs_pdt.clear();
00264 fac1.clear();
00265 fac2.clear();
00266 }
00267
00268 void ChargeBin::p_clear1()
00269 {
00270 DustZ = LONG_MIN;
00271 nfill = 0;
00272 FracPop = -DBL_MAX;
00273 tedust = 1.f;
00274 }
00275
00276 void GrainBin::p_clear0()
00277 {
00278 dstab1.clear();
00279 pure_sc1.clear();
00280 asym.clear();
00281 y0b06.clear();
00282 inv_att_len.clear();
00283
00284 for( unsigned int ns=0; ns < sd.size(); ns++ )
00285 delete sd[ns];
00286 sd.clear();
00287
00288 for( int nz=0; nz < NCHS; nz++ )
00289 {
00290 delete chrg[nz];
00291 chrg[nz] = NULL;
00292 }
00293 }
00294
00295 void GrainBin::p_clear1()
00296 {
00297 nDustFunc = DF_STANDARD;
00298 lgPAHsInIonizedRegion = false;
00299 avDGRatio = 0.;
00300 dstfactor = 1.f;
00301 dstAbund = -FLT_MAX;
00302 GrnDpth = 1.f;
00303 cnv_H_pGR = -DBL_MAX;
00304 cnv_H_pCM3 = -DBL_MAX;
00305 cnv_CM3_pGR = -DBL_MAX;
00306 cnv_CM3_pH = -DBL_MAX;
00307 cnv_GR_pH = -DBL_MAX;
00308 cnv_GR_pCM3 = -DBL_MAX;
00309
00310
00311 RSFCheck = 0.;
00312 memset( dstems, 0, NDEMS*sizeof(double) );
00313 memset( dstslp, 0, NDEMS*sizeof(double) );
00314 memset( dstslp2, 0, NDEMS*sizeof(double) );
00315 lgTdustConverged = false;
00316
00317
00318 tedust = 1.f;
00319 TeGrainMax = FLT_MAX;
00320 avdust = 0.;
00321 lgChrgConverged = false;
00322 LowestZg = LONG_MIN;
00323 nfill = 0;
00324 sd.reserve(15);
00325 AveDustZ = -DBL_MAX;
00326 dstpot = -DBL_MAX;
00327 dstpotsav = -DBL_MAX;
00328 LowestPot = -DBL_MAX;
00329 RateUp = -DBL_MAX;
00330 RateDn = -DBL_MAX;
00331 StickElecNeg = -DBL_MAX;
00332 StickElecPos = -DBL_MAX;
00333 avdpot = 0.;
00334 le_thres = FLT_MAX;
00335 BolFlux = -DBL_MAX;
00336 GrainCoolTherm = -DBL_MAX;
00337 GasHeatPhotoEl = -DBL_MAX;
00338 GrainHeat = DBL_MAX/10.;
00339 GrainHeatColl = -DBL_MAX;
00340 GrainGasCool = DBL_MAX/10.;
00341 ChemEn = -DBL_MAX;
00342 ChemEnH2 = -DBL_MAX;
00343 thermionic = -DBL_MAX;
00344 lgQHeat = false;
00345 lgUseQHeat = false;
00346 lgEverQHeat = false;
00347 lgQHTooWide = false;
00348 QHeatFailures = 0;
00349 qnflux = LONG_MAX;
00350 qnflux2 = LONG_MAX;
00351 qtmin = -DBL_MAX;
00352 qtmin_zone1 = -DBL_MAX;
00353 HeatingRate1 = -DBL_MAX;
00354 memset( DustEnth, 0, NDEMS*sizeof(double) );
00355 memset( EnthSlp, 0, NDEMS*sizeof(double) );
00356 memset( EnthSlp2, 0, NDEMS*sizeof(double) );
00357 rate_h2_form_grains_HM79 = 0.;
00358 rate_h2_form_grains_CT02 = 0.;
00359
00360 rate_h2_form_grains_used = 0.;
00361 DustDftVel = 1.e3f;
00362 avdft = 0.;
00363
00364 nChrgOrg = gv.nChrgRequested;
00365 nChrg = nChrgOrg;
00366 for( int nz=0; nz < NCHS; nz++ )
00367 chrg[nz] = NULL;
00368 }
00369
00370 void GrainVar::p_clear0()
00371 {
00372 for( size_t nd=0; nd < bin.size(); nd++ )
00373 delete bin[nd];
00374 bin.clear();
00375
00376 for( int nelem=0; nelem < LIMELM; nelem++ )
00377 {
00378 delete AugerData[nelem];
00379 AugerData[nelem] = NULL;
00380 }
00381
00382 ReadRecord.clear();
00383 anumin.clear();
00384 anumax.clear();
00385 dstab.clear();
00386 dstsc.clear();
00387 GrainEmission.clear();
00388 GraphiteEmission.clear();
00389 SilicateEmission.clear();
00390 }
00391
00392 void GrainVar::p_clear1()
00393 {
00394 bin.reserve(50);
00395
00396 for( int nelem=0; nelem < LIMELM; nelem++ )
00397 AugerData[nelem] = NULL;
00398
00399 lgAnyDustVary = false;
00400 TotalEden = 0.;
00401 dHeatdT = 0.;
00402 lgQHeatAll = false;
00403
00404
00405 lgGrainElectrons = true;
00406 lgQHeatOn = true;
00407 lgDHetOn = true;
00408 lgDColOn = true;
00409 GrainMetal = 1.;
00410 dstAbundThresholdNear = 1.e-6f;
00411 dstAbundThresholdFar = 1.e-3f;
00412 lgWD01 = false;
00413 nChrgRequested = NCHRG_DEFAULT;
00414
00415 lgReevaluate = true;
00416
00417 lgNegGrnDrg = false;
00418
00419
00420 nCalledGrainDrive = 0;
00421
00422
00423
00424 lgBakesPAH_heat = false;
00425
00426
00427
00428 lgGrainPhysicsOn = true;
00429
00430
00431
00432 GrainHeatScaleFactor = 1.f;
00433
00434
00435
00436 which_enth[MAT_CAR] = ENTH_CAR;
00437 which_zmin[MAT_CAR] = ZMIN_CAR;
00438 which_pot[MAT_CAR] = POT_CAR;
00439 which_ial[MAT_CAR] = IAL_CAR;
00440 which_pe[MAT_CAR] = PE_CAR;
00441 which_strg[MAT_CAR] = STRG_CAR;
00442 which_H2distr[MAT_CAR] = H2_CAR;
00443
00444 which_enth[MAT_SIL] = ENTH_SIL;
00445 which_zmin[MAT_SIL] = ZMIN_SIL;
00446 which_pot[MAT_SIL] = POT_SIL;
00447 which_ial[MAT_SIL] = IAL_SIL;
00448 which_pe[MAT_SIL] = PE_SIL;
00449 which_strg[MAT_SIL] = STRG_SIL;
00450 which_H2distr[MAT_SIL] = H2_SIL;
00451
00452 which_enth[MAT_PAH] = ENTH_PAH;
00453 which_zmin[MAT_PAH] = ZMIN_CAR;
00454 which_pot[MAT_PAH] = POT_CAR;
00455 which_ial[MAT_PAH] = IAL_CAR;
00456 which_pe[MAT_PAH] = PE_CAR;
00457 which_strg[MAT_PAH] = STRG_CAR;
00458 which_H2distr[MAT_PAH] = H2_CAR;
00459
00460 which_enth[MAT_CAR2] = ENTH_CAR2;
00461 which_zmin[MAT_CAR2] = ZMIN_CAR;
00462 which_pot[MAT_CAR2] = POT_CAR;
00463 which_ial[MAT_CAR2] = IAL_CAR;
00464 which_pe[MAT_CAR2] = PE_CAR;
00465 which_strg[MAT_CAR2] = STRG_CAR;
00466 which_H2distr[MAT_CAR2] = H2_CAR;
00467
00468 which_enth[MAT_SIL2] = ENTH_SIL2;
00469 which_zmin[MAT_SIL2] = ZMIN_SIL;
00470 which_pot[MAT_SIL2] = POT_SIL;
00471 which_ial[MAT_SIL2] = IAL_SIL;
00472 which_pe[MAT_SIL2] = PE_SIL;
00473 which_strg[MAT_SIL2] = STRG_SIL;
00474 which_H2distr[MAT_SIL2] = H2_SIL;
00475
00476 which_enth[MAT_PAH2] = ENTH_PAH2;
00477 which_zmin[MAT_PAH2] = ZMIN_CAR;
00478 which_pot[MAT_PAH2] = POT_CAR;
00479 which_ial[MAT_PAH2] = IAL_CAR;
00480 which_pe[MAT_PAH2] = PE_CAR;
00481 which_strg[MAT_PAH2] = STRG_CAR;
00482 which_H2distr[MAT_PAH2] = H2_CAR;
00483
00484 for( int nelem=0; nelem < LIMELM; nelem++ )
00485 {
00486 for( int ion=0; ion <= nelem+1; ion++ )
00487 {
00488 for( int ion_to=0; ion_to <= nelem+1; ion_to++ )
00489 {
00490 GrainChTrRate[nelem][ion][ion_to] = 0.f;
00491 }
00492 }
00493 }
00494
00495
00496
00497
00498 chPAH_abundance = "H";
00499 }
00500
00501
00502
00503
00504 void GrainZero(void)
00505 {
00506 DEBUG_ENTRY( "GrainZero()" );
00507
00508
00509
00510 gv.clear();
00511 return;
00512 }
00513
00514
00515
00516
00517 void GrainStartIter(void)
00518 {
00519 DEBUG_ENTRY( "GrainStartIter()" );
00520
00521 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
00522 {
00523 gv.lgNegGrnDrg = false;
00524 gv.TotalDustHeat = 0.;
00525 gv.GrnElecDonateMax = 0.;
00526 gv.GrnElecHoldMax = 0.;
00527 gv.dphmax = 0.f;
00528 gv.dclmax = 0.f;
00529
00530 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00531 {
00532
00533
00534 gv.bin[nd]->dstpotsav = gv.bin[nd]->dstpot;
00535 gv.bin[nd]->qtmin = ( gv.bin[nd]->qtmin_zone1 > 0. ) ?
00536 gv.bin[nd]->qtmin_zone1 : DBL_MAX;
00537 gv.bin[nd]->avdust = 0.;
00538 gv.bin[nd]->avdpot = 0.;
00539 gv.bin[nd]->avdft = 0.;
00540 gv.bin[nd]->avDGRatio = 0.;
00541 gv.bin[nd]->TeGrainMax = -1.f;
00542 gv.bin[nd]->lgEverQHeat = false;
00543 gv.bin[nd]->QHeatFailures = 0L;
00544 gv.bin[nd]->lgQHTooWide = false;
00545 gv.bin[nd]->lgPAHsInIonizedRegion = false;
00546 gv.bin[nd]->nChrgOrg = gv.bin[nd]->nChrg;
00547 }
00548 }
00549 return;
00550 }
00551
00552
00553
00554
00555 void GrainRestartIter(void)
00556 {
00557 DEBUG_ENTRY( "GrainRestartIter()" );
00558
00559 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
00560 {
00561 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00562 {
00563
00564
00565 gv.bin[nd]->dstpot = gv.bin[nd]->dstpotsav;
00566 gv.bin[nd]->nChrg = gv.bin[nd]->nChrgOrg;
00567 }
00568 }
00569 return;
00570 }
00571
00572
00573
00574 void SetNChrgStates(long nChrg)
00575 {
00576 DEBUG_ENTRY( "SetNChrgStates()" );
00577
00578 ASSERT( nChrg >= 2 && nChrg <= NCHU );
00579 gv.nChrgRequested = nChrg;
00580 return;
00581 }
00582
00583
00584
00585
00586
00587 void GrainsInit(void)
00588 {
00589 long int i,
00590 nelem;
00591 unsigned int ns;
00592
00593 DEBUG_ENTRY( "GrainsInit()" );
00594
00595 if( trace.lgTrace && trace.lgDustBug )
00596 {
00597 fprintf( ioQQQ, " GrainsInit called.\n" );
00598 }
00599
00600 gv.anumin.resize( rfield.nupper );
00601 gv.anumax.resize( rfield.nupper );
00602 gv.dstab.resize( rfield.nupper );
00603 gv.dstsc.resize( rfield.nupper );
00604 gv.GrainEmission.resize( rfield.nupper );
00605 gv.GraphiteEmission.resize( rfield.nupper );
00606 gv.SilicateEmission.resize( rfield.nupper );
00607
00608
00609 for( nelem=0; nelem < LIMELM; nelem++ )
00610 {
00611 gv.elmSumAbund[nelem] = 0.f;
00612 }
00613
00614 for( i=0; i < rfield.nupper; i++ )
00615 {
00616 gv.dstab[i] = 0.;
00617 gv.dstsc[i] = 0.;
00618
00619 gv.GrainEmission[i] = 0.;
00620 gv.SilicateEmission[i] = 0.;
00621 gv.GraphiteEmission[i] = 0.;
00622 }
00623
00624 if( !gv.lgDustOn() )
00625 {
00626
00627 gv.GrainHeatInc = 0.;
00628 gv.GrainHeatDif = 0.;
00629 gv.GrainHeatLya = 0.;
00630 gv.GrainHeatCollSum = 0.;
00631 gv.GrainHeatSum = 0.;
00632 gv.GasCoolColl = 0.;
00633 thermal.heating[0][13] = 0.;
00634 thermal.heating[0][14] = 0.;
00635 thermal.heating[0][25] = 0.;
00636
00637 if( trace.lgTrace && trace.lgDustBug )
00638 {
00639 fprintf( ioQQQ, " GrainsInit exits.\n" );
00640 }
00641 return;
00642 }
00643
00644 #ifdef WD_TEST2
00645 gv.lgWD01 = true;
00646 #endif
00647
00648 HEAT_TOLER = conv.HeatCoolRelErrorAllowed / 3.;
00649 HEAT_TOLER_BIN = HEAT_TOLER / sqrt((double)gv.bin.size());
00650 CHRG_TOLER = conv.EdenErrorAllowed / 3.;
00651
00652
00653 gv.anumin[0] = 0.f;
00654 for( i=1; i < rfield.nupper; i++ )
00655 gv.anumax[i-1] = gv.anumin[i] = (realnum)sqrt(rfield.anu[i-1]*rfield.anu[i]);
00656 gv.anumax[rfield.nupper-1] = FLT_MAX;
00657
00658 ReadAugerData();
00659
00660 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00661 {
00662 double help,atoms,p_rad,ThresInf,ThresInfVal,Emin,d[5];
00663 long low1,low2,low3,lowm;
00664
00665
00666 ASSERT( gv.bin[nd] != NULL );
00667 ASSERT( gv.bin[nd]->nChrg >= 2 && gv.bin[nd]->nChrg <= NCHU );
00668
00669 if( gv.bin[nd]->DustWorkFcn < rfield.anu[0] || gv.bin[nd]->DustWorkFcn > rfield.anu[rfield.nupper] )
00670 {
00671 fprintf( ioQQQ, " Grain work function for %s has insane value: %.4e\n",
00672 gv.bin[nd]->chDstLab,gv.bin[nd]->DustWorkFcn );
00673 cdEXIT(EXIT_FAILURE);
00674 }
00675
00676
00677 if( gv.lgQHeatAll )
00678 {
00679 gv.bin[nd]->lgQHeat = true;
00680 }
00681
00682
00683 if( !gv.lgQHeatOn )
00684 {
00685 gv.bin[nd]->lgQHeat = false;
00686 }
00687
00688
00689 if( thermal.ConstGrainTemp > 0. )
00690 {
00691 gv.bin[nd]->lgQHeat = false;
00692 }
00693
00694 #ifndef IGNORE_QUANTUM_HEATING
00695 gv.bin[nd]->lgQHTooWide = false;
00696 gv.bin[nd]->qtmin = DBL_MAX;
00697 #endif
00698
00699 if( gv.bin[nd]->nDustFunc>DF_STANDARD || gv.bin[nd]->matType == MAT_PAH ||
00700 gv.bin[nd]->matType == MAT_PAH2 )
00701 gv.lgAnyDustVary = true;
00702
00703
00704
00705 gv.bin[nd]->dstAbund = -FLT_MAX;
00706
00707 gv.bin[nd]->GrnDpth = 1.f;
00708
00709 gv.bin[nd]->qtmin_zone1 = -1.;
00710
00711
00712 gv.bin[nd]->le_thres = gv.lgWD01 ? FLT_MAX :
00713 (realnum)(pow(pow((double)gv.bin[nd]->dustp[0],0.85)/30.,2./3.)*1.e3/EVRYD);
00714
00715 for( long nz=0; nz < NCHS; nz++ )
00716 {
00717 ASSERT( gv.bin[nd]->chrg[nz] == NULL );
00718 gv.bin[nd]->chrg[nz] = new ChargeBin;
00719 }
00720
00721
00722
00723 zmin_type zcase = gv.which_zmin[gv.bin[nd]->matType];
00724 switch( zcase )
00725 {
00726 case ZMIN_CAR:
00727
00728 help = gv.bin[nd]->AvRadius*1.e7;
00729 help = ceil(-(1.2*POW2(help)+3.9*help+0.2)/1.44);
00730 break;
00731 case ZMIN_SIL:
00732
00733 help = gv.bin[nd]->AvRadius*1.e7;
00734 help = ceil(-(0.7*POW2(help)+2.5*help+0.8)/1.44);
00735 break;
00736 default:
00737 fprintf( ioQQQ, " GrainsInit detected unknown Zmin type: %d\n" , zcase );
00738 cdEXIT(EXIT_FAILURE);
00739 }
00740
00741
00742 ASSERT( help > (double)(LONG_MIN+1) );
00743 low1 = nint(help);
00744
00745
00746
00747
00748
00749
00750 low2 = low1;
00751 GetPotValues(nd,low2,&ThresInf,&d[0],&d[1],&d[2],&d[3],&d[4],INCL_TUNNEL);
00752 if( ThresInf < 0. )
00753 {
00754 low3 = 0;
00755
00756
00757 while( low3-low2 > 1 )
00758 {
00759 lowm = (low2+low3)/2;
00760 GetPotValues(nd,lowm,&ThresInf,&d[0],&d[1],&d[2],&d[3],&d[4],INCL_TUNNEL);
00761 if( ThresInf < 0. )
00762 low2 = lowm;
00763 else
00764 low3 = lowm;
00765 }
00766 low2 = low3;
00767 }
00768
00769
00770
00771
00772 gv.bin[nd]->LowestZg = MAX2(low1,low2);
00773 gv.bin[nd]->LowestPot = chrg2pot(gv.bin[nd]->LowestZg,nd);
00774
00775 ns = 0;
00776
00777 ASSERT( gv.bin[nd]->sd.size() == 0 );
00778 gv.bin[nd]->sd.push_back( new ShellData );
00779
00780
00781 gv.bin[nd]->sd[ns]->nelem = -1;
00782 gv.bin[nd]->sd[ns]->ns = -1;
00783 gv.bin[nd]->sd[ns]->ionPot = gv.bin[nd]->DustWorkFcn;
00784
00785
00786 for( nelem=ipLITHIUM; nelem < LIMELM && !gv.lgWD01; nelem++ )
00787 {
00788 if( gv.bin[nd]->elmAbund[nelem] > 0. )
00789 {
00790 if( gv.AugerData[nelem] == NULL )
00791 {
00792 fprintf( ioQQQ, " Grain Auger data are missing for element %s\n",
00793 elementnames.chElementName[nelem] );
00794 fprintf( ioQQQ, " Please include the NO GRAIN X-RAY TREATMENT command "
00795 "to disable the Auger treatment in grains.\n" );
00796 cdEXIT(EXIT_FAILURE);
00797 }
00798
00799 for( unsigned int j=0; j < gv.AugerData[nelem]->nSubShell; j++ )
00800 {
00801 ++ns;
00802
00803 gv.bin[nd]->sd.push_back( new ShellData );
00804
00805 gv.bin[nd]->sd[ns]->nelem = nelem;
00806 gv.bin[nd]->sd[ns]->ns = j;
00807 gv.bin[nd]->sd[ns]->ionPot = gv.AugerData[nelem]->IonThres[j];
00808 }
00809 }
00810 }
00811
00812 GetPotValues(nd,gv.bin[nd]->LowestZg,&d[0],&ThresInfVal,&d[1],&d[2],&d[3],&Emin,INCL_TUNNEL);
00813
00814 for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
00815 {
00816 long ipLo;
00817 double Ethres = ( ns == 0 ) ? ThresInfVal : gv.bin[nd]->sd[ns]->ionPot;
00818 ShellData *sptr = gv.bin[nd]->sd[ns];
00819
00820 sptr->ipLo = hunt_bisect( rfield.anu, rfield.nupper, Ethres ) + 1;
00821
00822 ipLo = sptr->ipLo;
00823
00824
00825 long len = rfield.nflux + 10 - ipLo;
00826
00827 sptr->p.reserve( len );
00828 sptr->p.alloc( ipLo, rfield.nflux );
00829
00830 sptr->y01.reserve( len );
00831 sptr->y01.alloc( ipLo, rfield.nflux );
00832
00833
00834 if( ns > 0 )
00835 {
00836 sptr->nData = gv.AugerData[sptr->nelem]->nData[sptr->ns];
00837 sptr->AvNr.resize( sptr->nData );
00838 sptr->Ener.resize( sptr->nData );
00839 sptr->y01A.resize( sptr->nData );
00840
00841 for( long n=0; n < sptr->nData; n++ )
00842 {
00843 sptr->AvNr[n] = gv.AugerData[sptr->nelem]->AvNumber[sptr->ns][n];
00844 sptr->Ener[n] = gv.AugerData[sptr->nelem]->Energy[sptr->ns][n];
00845
00846 sptr->y01A[n].reserve( len );
00847 sptr->y01A[n].alloc( ipLo, rfield.nflux );
00848 }
00849 }
00850 }
00851
00852 gv.bin[nd]->y0b06.resize( rfield.nupper );
00853
00854 InitBinAugerData( nd, 0, rfield.nflux );
00855
00856 gv.bin[nd]->nfill = rfield.nflux;
00857
00858
00859
00860
00861
00862
00864 gv.bin[nd]->StickElecPos = STICK_ELEC*(1. - exp(-gv.bin[nd]->AvRadius/elec_esc_length(0.,nd)));
00865 atoms = gv.bin[nd]->AvVol*gv.bin[nd]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[nd]->atomWeight;
00866 p_rad = 1./(1.+exp(20.-atoms));
00867 gv.bin[nd]->StickElecNeg = gv.bin[nd]->StickElecPos*p_rad;
00868
00869
00870
00871
00872 gv.bin[nd]->GrnDpth = (realnum)GrnStdDpth(nd);
00873 gv.bin[nd]->dstAbund = (realnum)(gv.bin[nd]->dstfactor*gv.GrainMetal*gv.bin[nd]->GrnDpth);
00874 ASSERT( gv.bin[nd]->dstAbund > 0.f );
00875
00876 gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
00877 gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
00878
00879 gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
00880 gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;
00881 }
00882
00883
00884
00885
00886
00887
00888 for( nelem=0; nelem < LIMELM; nelem++ )
00889 {
00890 gv.elmSumAbund[nelem] = 0.f;
00891 }
00892
00893 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00894 {
00895 for( nelem=0; nelem < LIMELM; nelem++ )
00896 {
00897 gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(realnum)gv.bin[nd]->cnv_H_pCM3;
00898 }
00899 }
00900
00901 gv.nzone = -1;
00902 gv.GrnRecomTe = -1.;
00903
00904
00905
00906 for( i=0; i < rfield.nupper; i++ )
00907 {
00908
00909
00910 gv.dstab[i] = -DBL_MAX;
00911 gv.dstsc[i] = -DBL_MAX;
00912 }
00913
00914 InitEmissivities();
00915 InitEnthalpy();
00916
00917 if( trace.lgDustBug && trace.lgTrace )
00918 {
00919 fprintf( ioQQQ, " There are %ld grain types turned on.\n", (unsigned long)gv.bin.size() );
00920
00921 fprintf( ioQQQ, " grain depletion factors, dstfactor*GrainMetal=" );
00922 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00923 {
00924 fprintf( ioQQQ, "%10.2e", gv.bin[nd]->dstfactor*gv.GrainMetal );
00925 }
00926 fprintf( ioQQQ, "\n" );
00927
00928 fprintf( ioQQQ, " nChrg =" );
00929 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00930 {
00931 fprintf( ioQQQ, " %ld", gv.bin[nd]->nChrg );
00932 }
00933 fprintf( ioQQQ, "\n" );
00934
00935 fprintf( ioQQQ, " lowest charge (e) =" );
00936 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00937 {
00938 fprintf( ioQQQ, "%10.2e", pot2chrg(gv.bin[nd]->LowestPot,nd) );
00939 }
00940 fprintf( ioQQQ, "\n" );
00941
00942 fprintf( ioQQQ, " nDustFunc flag for user requested custom depth dependence:" );
00943 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00944 {
00945 fprintf( ioQQQ, "%2i", gv.bin[nd]->nDustFunc );
00946 }
00947 fprintf( ioQQQ, "\n" );
00948
00949 fprintf( ioQQQ, " Quantum heating flag:" );
00950 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00951 {
00952 fprintf( ioQQQ, "%2c", TorF(gv.bin[nd]->lgQHeat) );
00953 }
00954 fprintf( ioQQQ, "\n" );
00955
00956
00957 fprintf( ioQQQ, " NU(Ryd), Abs cross sec per proton\n" );
00958
00959 fprintf( ioQQQ, " Ryd " );
00960 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00961 {
00962 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00963 }
00964 fprintf( ioQQQ, "\n" );
00965
00966 for( i=0; i < rfield.nupper; i += 40 )
00967 {
00968 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00969 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00970 {
00971 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->dstab1[i] );
00972 }
00973 fprintf( ioQQQ, "\n" );
00974 }
00975
00976 fprintf( ioQQQ, " NU(Ryd), Sct cross sec per proton\n" );
00977
00978 fprintf( ioQQQ, " Ryd " );
00979 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00980 {
00981 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00982 }
00983 fprintf( ioQQQ, "\n" );
00984
00985 for( i=0; i < rfield.nupper; i += 40 )
00986 {
00987 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00988 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00989 {
00990 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->pure_sc1[i] );
00991 }
00992 fprintf( ioQQQ, "\n" );
00993 }
00994
00995 fprintf( ioQQQ, " NU(Ryd), Q abs\n" );
00996
00997 fprintf( ioQQQ, " Ryd " );
00998 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00999 {
01000 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
01001 }
01002 fprintf( ioQQQ, "\n" );
01003
01004 for( i=0; i < rfield.nupper; i += 40 )
01005 {
01006 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
01007 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01008 {
01009 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->dstab1[i]*4./gv.bin[nd]->IntArea );
01010 }
01011 fprintf( ioQQQ, "\n" );
01012 }
01013
01014 fprintf( ioQQQ, " NU(Ryd), Q sct\n" );
01015
01016 fprintf( ioQQQ, " Ryd " );
01017 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01018 {
01019 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
01020 }
01021 fprintf( ioQQQ, "\n" );
01022
01023 for( i=0; i < rfield.nupper; i += 40 )
01024 {
01025 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
01026 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01027 {
01028 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->pure_sc1[i]*4./gv.bin[nd]->IntArea );
01029 }
01030 fprintf( ioQQQ, "\n" );
01031 }
01032
01033 fprintf( ioQQQ, " NU(Ryd), asymmetry factor\n" );
01034
01035 fprintf( ioQQQ, " Ryd " );
01036 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01037 {
01038 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
01039 }
01040 fprintf( ioQQQ, "\n" );
01041
01042 for( i=0; i < rfield.nupper; i += 40 )
01043 {
01044 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
01045 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01046 {
01047 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->asym[i] );
01048 }
01049 fprintf( ioQQQ, "\n" );
01050 }
01051
01052 fprintf( ioQQQ, " GrainsInit exits.\n" );
01053 }
01054 return;
01055 }
01056
01057
01058 STATIC void ReadAugerData()
01059 {
01060 char chString[FILENAME_PATH_LENGTH_2];
01061 long version;
01062 FILE *fdes;
01063
01064 DEBUG_ENTRY( "ReadAugerData()" );
01065
01066 static const char chFile[] = "auger_spectrum.dat";
01067 fdes = open_data( chFile, "r" );
01068
01069 GetNextLine( chFile, fdes, chString );
01070 sscanf( chString, "%ld", &version );
01071 if( version != MAGIC_AUGER_DATA )
01072 {
01073 fprintf( ioQQQ, " File %s has wrong version number\n", chFile );
01074 fprintf( ioQQQ, " please update your installation...\n" );
01075 cdEXIT(EXIT_FAILURE);
01076 }
01077
01078 while( true )
01079 {
01080 int res;
01081 long nelem;
01082 unsigned int ns;
01083 AEInfo *ptr;
01084
01085 GetNextLine( chFile, fdes, chString );
01086 res = sscanf( chString, "%ld", &nelem );
01087 ASSERT( res == 1 );
01088
01089 if( nelem < 0 )
01090 break;
01091
01092 ASSERT( nelem < LIMELM );
01093
01094 ptr = new AEInfo;
01095
01096 GetNextLine( chFile, fdes, chString );
01097 res = sscanf( chString, "%u", &ptr->nSubShell );
01098 ASSERT( res == 1 && ptr->nSubShell > 0 );
01099
01100 ptr->nData.resize( ptr->nSubShell );
01101 ptr->IonThres.resize( ptr->nSubShell );
01102 ptr->Energy.resize( ptr->nSubShell );
01103 ptr->AvNumber.resize( ptr->nSubShell );
01104
01105 for( ns=0; ns < ptr->nSubShell; ns++ )
01106 {
01107 unsigned int ss;
01108
01109 GetNextLine( chFile, fdes, chString );
01110 res = sscanf( chString, "%u", &ss );
01111 ASSERT( res == 1 && ns == ss );
01112
01113 GetNextLine( chFile, fdes, chString );
01114 res = sscanf( chString, "%le", &ptr->IonThres[ns] );
01115 ASSERT( res == 1 );
01116 ptr->IonThres[ns] /= EVRYD;
01117
01118 GetNextLine( chFile, fdes, chString );
01119 res = sscanf( chString, "%u", &ptr->nData[ns] );
01120 ASSERT( res == 1 && ptr->nData[ns] > 0 );
01121
01122 ptr->Energy[ns].resize( ptr->nData[ns] );
01123 ptr->AvNumber[ns].resize( ptr->nData[ns] );
01124
01125 for( unsigned int n=0; n < ptr->nData[ns]; n++ )
01126 {
01127 GetNextLine( chFile, fdes, chString );
01128 res = sscanf(chString,"%le %le",&ptr->Energy[ns][n],&ptr->AvNumber[ns][n]);
01129 ASSERT( res == 2 );
01130 ptr->Energy[ns][n] /= EVRYD;
01131 ASSERT( ptr->Energy[ns][n] < ptr->IonThres[ns] );
01132 }
01133 }
01134
01135 ASSERT( gv.AugerData[nelem] == NULL );
01136 gv.AugerData[nelem] = ptr;
01137 }
01138
01139 fclose( fdes );
01140 }
01141
01143 STATIC void InitBinAugerData(size_t nd,
01144 long ipBegin,
01145 long ipEnd)
01146 {
01147 DEBUG_ENTRY( "InitBinAugerData()" );
01148
01149 long i, ipLo, nelem;
01150 unsigned int ns;
01151
01152 flex_arr<realnum> temp( ipBegin, ipEnd );
01153 temp.zero();
01154
01155
01156 double norm = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->AvVol;
01157
01158
01159 for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
01160 {
01161 ipLo = max( gv.bin[nd]->sd[ns]->ipLo, ipBegin );
01162
01163 gv.bin[nd]->sd[ns]->p.realloc( ipEnd );
01164
01167 for( i=ipLo; i < ipEnd; i++ )
01168 {
01169 long nel,nsh;
01170 double phot_ev,cs;
01171
01172 phot_ev = rfield.anu[i]*EVRYD;
01173
01174 if( ns == 0 )
01175 {
01176
01177
01178 gv.bin[nd]->sd[ns]->p[i] = 0.;
01179
01180 for( nelem=ipHYDROGEN; nelem < LIMELM && !gv.lgWD01; nelem++ )
01181 {
01182 if( gv.bin[nd]->elmAbund[nelem] == 0. )
01183 continue;
01184
01185 long nshmax = Heavy.nsShells[nelem][0];
01186
01187 for( nsh = gv.AugerData[nelem]->nSubShell; nsh < nshmax; nsh++ )
01188 {
01189 nel = nelem+1;
01190 cs = t_ADfA::Inst().phfit(nelem+1,nel,nsh+1,phot_ev);
01191 gv.bin[nd]->sd[ns]->p[i] +=
01192 (realnum)(norm*gv.bin[nd]->elmAbund[nelem]*cs*1e-18);
01193 }
01194 }
01195
01196 temp[i] += gv.bin[nd]->sd[ns]->p[i];
01197 }
01198 else
01199 {
01200
01201 nelem = gv.bin[nd]->sd[ns]->nelem;
01202 nel = nelem+1;
01203 nsh = gv.bin[nd]->sd[ns]->ns;
01204 cs = t_ADfA::Inst().phfit(nelem+1,nel,nsh+1,phot_ev);
01205 gv.bin[nd]->sd[ns]->p[i] =
01206 (realnum)(norm*gv.bin[nd]->elmAbund[nelem]*cs*1e-18);
01207 temp[i] += gv.bin[nd]->sd[ns]->p[i];
01208 }
01209 }
01210 }
01211
01212 for( i=ipBegin; i < ipEnd && !gv.lgWD01; i++ )
01213 {
01214
01215 if( rfield.anu[i] > 20./EVRYD )
01216 gv.bin[nd]->inv_att_len[i] = temp[i];
01217 }
01218
01219 for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
01220 {
01221 ipLo = max( gv.bin[nd]->sd[ns]->ipLo, ipBegin );
01222
01223 for( i=ipLo; i < ipEnd; i++ )
01224 {
01225 if( temp[i] > 0. )
01226 gv.bin[nd]->sd[ns]->p[i] /= temp[i];
01227 else
01228 gv.bin[nd]->sd[ns]->p[i] = ( ns == 0 ) ? 1.f : 0.f;
01229 }
01230 }
01231
01232 temp.clear();
01233
01234 for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
01235 {
01236 long n;
01237 ShellData *sptr = gv.bin[nd]->sd[ns];
01238
01239 ipLo = max( sptr->ipLo, ipBegin );
01240
01241
01242 sptr->y01.realloc( ipEnd );
01243
01244 for( i=ipLo; i < ipEnd; i++ )
01245 {
01246 double elec_en,yzero,yone;
01247
01248 elec_en = MAX2(rfield.anu[i] - sptr->ionPot,0.);
01249 yzero = y0psa( nd, ns, i, elec_en );
01250
01251
01252
01253 yone = y1psa( nd, i, elec_en );
01254
01255 if( ns == 0 )
01256 {
01257 gv.bin[nd]->y0b06[i] = (realnum)yzero;
01258 sptr->y01[i] = (realnum)yone;
01259 }
01260 else
01261 {
01262 sptr->y01[i] = (realnum)(yzero*yone);
01263 }
01264 }
01265
01266
01267 if( ns > 0 )
01268 {
01269
01270 for( n=0; n < sptr->nData; n++ )
01271 {
01272 sptr->y01A[n].realloc( ipEnd );
01273
01274 for( i=ipLo; i < ipEnd; i++ )
01275 {
01276 double yzero = sptr->AvNr[n] * y0psa( nd, ns, i, sptr->Ener[n] );
01277
01278
01279
01280 double yone = y1psa( nd, i, sptr->Ener[n] );
01281
01282 sptr->y01A[n][i] = (realnum)(yzero*yone);
01283 }
01284 }
01285 }
01286 }
01287 }
01288
01289
01290 STATIC void GetNextLine(const char *chFile,
01291 FILE *io,
01292 char chLine[])
01293 {
01294 char *str;
01295
01296 DEBUG_ENTRY( "GetNextLine()" );
01297
01298 do
01299 {
01300 if( read_whole_line( chLine, FILENAME_PATH_LENGTH_2, io ) == NULL )
01301 {
01302 fprintf( ioQQQ, " Could not read from %s\n", chFile );
01303 if( feof(io) )
01304 fprintf( ioQQQ, " EOF reached\n");
01305 cdEXIT(EXIT_FAILURE);
01306 }
01307 }
01308 while( chLine[0] == '#' );
01309
01310
01311 str = strstr_s(chLine,"#");
01312 if( str != NULL )
01313 *str = '\0';
01314 return;
01315 }
01316
01317 STATIC void InitEmissivities(void)
01318 {
01319 double fac,
01320 fac2,
01321 mul,
01322 tdust;
01323 long int i;
01324
01325 DEBUG_ENTRY( "InitEmissivities()" );
01326
01327 if( trace.lgTrace && trace.lgDustBug )
01328 {
01329 fprintf( ioQQQ, " InitEmissivities starts\n" );
01330 fprintf( ioQQQ, " ND Tdust Emis BB Check 4pi*a^2*<Q>\n" );
01331 }
01332
01333 ASSERT( NTOP >= 2 && NDEMS > 2*NTOP );
01334 fac = exp(log(GRAIN_TMID/GRAIN_TMIN)/(double)(NDEMS-NTOP));
01335 tdust = GRAIN_TMIN;
01336 for( i=0; i < NDEMS-NTOP; i++ )
01337 {
01338 gv.dsttmp[i] = log(tdust);
01339 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01340 {
01341 gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
01342 }
01343 tdust *= fac;
01344 }
01345
01346
01347 fac2 = exp(log(GRAIN_TMAX/GRAIN_TMID/powi(fac,NTOP-1))/(double)((NTOP-1)*NTOP/2));
01348 for( i=NDEMS-NTOP; i < NDEMS; i++ )
01349 {
01350 gv.dsttmp[i] = log(tdust);
01351 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01352 {
01353 gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
01354 }
01355 fac *= fac2;
01356 tdust *= fac;
01357 }
01358
01359
01360 mul = 1.;
01361 ASSERT( fabs(gv.dsttmp[0] - log(GRAIN_TMIN)) < 10.*mul*DBL_EPSILON );
01362 mul = sqrt((double)(NDEMS-NTOP));
01363 ASSERT( fabs(gv.dsttmp[NDEMS-NTOP] - log(GRAIN_TMID)) < 10.*mul*DBL_EPSILON );
01364 mul = (double)NTOP + sqrt((double)NDEMS);
01365 ASSERT( fabs(gv.dsttmp[NDEMS-1] - log(GRAIN_TMAX)) < 10.*mul*DBL_EPSILON );
01366
01367
01368 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01369 {
01370
01371 spline(gv.bin[nd]->dstems,gv.dsttmp,NDEMS,2e31,2e31,gv.bin[nd]->dstslp);
01372 spline(gv.dsttmp,gv.bin[nd]->dstems,NDEMS,2e31,2e31,gv.bin[nd]->dstslp2);
01373 }
01374
01375 # if 0
01376
01377 nd = nint(fudge(0));
01378 ASSERT( nd >= 0 && nd < gv.bin.size() );
01379 for( i=0; i < 2000; i++ )
01380 {
01381 double x,y,z;
01382 z = pow(10.,-40. + (double)i/50.);
01383 splint(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,log(z),&y);
01384 if( exp(y) > GRAIN_TMIN && exp(y) < GRAIN_TMAX )
01385 {
01386 x = PlanckIntegral(exp(y),nd,3);
01387 printf(" input %.6e temp %.6e output %.6e rel. diff. %.6e\n",z,exp(y),x,(x-z)/z);
01388 }
01389 }
01390 cdEXIT(EXIT_SUCCESS);
01391 # endif
01392 return;
01393 }
01394
01395
01396
01397 STATIC double PlanckIntegral(double tdust,
01398 size_t nd,
01399 long int ip)
01400 {
01401 long int i;
01402
01403 double arg,
01404 ExpM1,
01405 integral1 = 0.,
01406 integral2 = 0.,
01407 Planck1,
01408 Planck2,
01409 TDustRyg,
01410 x;
01411
01412 DEBUG_ENTRY( "PlanckIntegral()" );
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424 TDustRyg = TE1RYD/tdust;
01425
01426 x = 0.999*log(DBL_MAX);
01427
01428 for( i=0; i < rfield.nupper; i++ )
01429 {
01430
01431 arg = TDustRyg*rfield.anu[i];
01432
01433
01434 if( arg < 1.e-5 )
01435 {
01436
01437 ExpM1 = arg*(1. + arg/2.);
01438 }
01439 else
01440 {
01441
01442 ExpM1 = exp(MIN2(x,arg)) - 1.;
01443 }
01444
01445 Planck1 = PI4*2.*HPLANCK/POW2(SPEEDLIGHT)*POW2(FR1RYD)*POW2(FR1RYD)*
01446 rfield.anu3[i]/ExpM1*rfield.widflx[i];
01447 Planck2 = Planck1*gv.bin[nd]->dstab1[i];
01448
01449
01450 if( i == 0 )
01451 {
01452 integral1 = Planck1/rfield.widflx[0]*rfield.anu[0]/3.;
01453 integral2 = Planck2/rfield.widflx[0]*rfield.anu[0]/5.;
01454 }
01455
01456 if( Planck1/integral1 < DBL_EPSILON && Planck2/integral2 < DBL_EPSILON )
01457 break;
01458
01459 integral1 += Planck1;
01460 integral2 += Planck2;
01461 }
01462
01463
01464 if( trace.lgDustBug && trace.lgTrace && ip%10 == 0 )
01465 {
01466 fprintf( ioQQQ, " %4ld %11.4e %11.4e %11.4e %11.4e\n", (unsigned long)nd, tdust,
01467 integral2, integral1/4./5.67051e-5/powi(tdust,4), integral2*4./integral1 );
01468 }
01469
01470 ASSERT( integral2 > 0. );
01471 return integral2;
01472 }
01473
01474
01475
01476 STATIC void NewChargeData(long nd)
01477 {
01478 long nz;
01479
01480 DEBUG_ENTRY( "NewChargeData()" );
01481
01482 for( nz=0; nz < NCHS; nz++ )
01483 {
01484 gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
01485 gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
01486 gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
01487 gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
01488 gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
01489
01491 gv.bin[nd]->chrg[nz]->ThermRate = -DBL_MAX;
01492 gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
01493 gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
01494
01495 gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
01496 gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
01497 gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
01498 }
01499
01500 if( !fp_equal(phycon.te,gv.GrnRecomTe) )
01501 {
01502 for( nz=0; nz < NCHS; nz++ )
01503 {
01504 memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
01505 memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );
01506 }
01507 }
01508
01509 if( nzone != gv.nzone )
01510 {
01511 for( nz=0; nz < NCHS; nz++ )
01512 {
01513 gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
01514 }
01515 }
01516 return;
01517 }
01518
01519
01520
01521
01522 STATIC double GrnStdDpth(long int nd)
01523 {
01524 double GrnStdDpth_v;
01525
01526 DEBUG_ENTRY( "GrnStdDpth()" );
01527
01528
01529
01530
01531
01532
01533 if( gv.bin[nd]->nDustFunc == DF_STANDARD )
01534 {
01535 if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
01536 {
01537
01538 if( gv.chPAH_abundance == "H" )
01539 {
01540
01541
01542
01543
01544
01545 GrnStdDpth_v = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
01546 }
01547 else if( gv.chPAH_abundance == "H,H2" )
01548 {
01549
01550
01551
01552
01553
01554 GrnStdDpth_v = (dense.xIonDense[ipHYDROGEN][0]+2*hmi.H2_total)/dense.gas_phase[ipHYDROGEN];
01555 }
01556 else if( gv.chPAH_abundance == "CON" )
01557 {
01558
01559 GrnStdDpth_v = 1.;
01560 }
01561 else
01562 {
01563 fprintf(ioQQQ,"Invalid argument to SET PAH: %s\n",gv.chPAH_abundance.c_str());
01564 TotalInsanity();
01565 }
01566 }
01567 else
01568 {
01569
01570 GrnStdDpth_v = 1.;
01571 }
01572 }
01573 else if( gv.bin[nd]->nDustFunc == DF_USER_FUNCTION )
01574 {
01575 GrnStdDpth_v = GrnVryDpth(nd);
01576 }
01577 else if( gv.bin[nd]->nDustFunc == DF_SUBLIMATION )
01578 {
01579
01580
01581 GrnStdDpth_v = sexp( pow3( gv.bin[nd]->tedust / gv.bin[nd]->Tsublimat ) );
01582 }
01583 else
01584 {
01585 TotalInsanity();
01586 }
01587
01588 GrnStdDpth_v = min(max(1.e-10,GrnStdDpth_v),1.);
01589
01590 return GrnStdDpth_v;
01591 }
01592
01593
01594
01595 void GrainDrive(void)
01596 {
01597 DEBUG_ENTRY( "GrainDrive()" );
01598
01599
01600 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
01601 {
01602 static double tesave = -1.;
01603 static long int nzonesave = -1;
01604
01605
01606
01607
01608
01609
01610
01611 if( gv.lgReevaluate || conv.lgSearch || nzonesave != nzone ||
01612
01613 ! fp_equal( phycon.te, tesave ) ||
01614
01615
01616 fabs(gv.TotalEden)/dense.eden > conv.EdenErrorAllowed/5. ||
01617
01618
01619 (fabs( gv.GasCoolColl ) + fabs( thermal.heating[0][13] ))/SDIV(thermal.ctot)>0.1 )
01620 {
01621 nzonesave = nzone;
01622 tesave = phycon.te;
01623
01624 if( trace.nTrConvg >= 5 )
01625 {
01626 fprintf( ioQQQ, " grain5 calling GrainChargeTemp\n");
01627 }
01628
01629
01630 GrainChargeTemp();
01631
01632
01633 }
01634 }
01635 else if( gv.lgDustOn() && !gv.lgGrainPhysicsOn )
01636 {
01637
01638
01639 if( nCalledGrainDrive == 0 )
01640 {
01641 long nelem, ion, ion_to;
01642
01643
01644
01645 gv.GasHeatPhotoEl = 0.;
01646 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01647 {
01648 long nz;
01649
01650
01651 gv.bin[nd]->lgPAHsInIonizedRegion = false;
01652
01653 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
01654 {
01655 gv.bin[nd]->chrg[nz]->DustZ = nz;
01656 gv.bin[nd]->chrg[nz]->FracPop = ( nz == 0 ) ? 1. : 0.;
01657 gv.bin[nd]->chrg[nz]->nfill = 0;
01658 gv.bin[nd]->chrg[nz]->tedust = 100.f;
01659 }
01660
01661 gv.bin[nd]->AveDustZ = 0.;
01662 gv.bin[nd]->dstpot = chrg2pot(0.,nd);
01663
01664 gv.bin[nd]->tedust = 100.f;
01665 gv.bin[nd]->TeGrainMax = 100.;
01666
01667
01668 gv.bin[nd]->BolFlux = 0.;
01669 gv.bin[nd]->GrainCoolTherm = 0.;
01670 gv.bin[nd]->GasHeatPhotoEl = 0.;
01671 gv.bin[nd]->GrainHeat = 0.;
01672 gv.bin[nd]->GrainHeatColl = 0.;
01673 gv.bin[nd]->ChemEn = 0.;
01674 gv.bin[nd]->ChemEnH2 = 0.;
01675 gv.bin[nd]->thermionic = 0.;
01676
01677 gv.bin[nd]->lgUseQHeat = false;
01678 gv.bin[nd]->lgEverQHeat = false;
01679 gv.bin[nd]->QHeatFailures = 0;
01680
01681 gv.bin[nd]->DustDftVel = 0.;
01682
01683 gv.bin[nd]->avdust = gv.bin[nd]->tedust;
01684 gv.bin[nd]->avdft = 0.f;
01685 gv.bin[nd]->avdpot = (realnum)(gv.bin[nd]->dstpot*EVRYD);
01686 gv.bin[nd]->avDGRatio = -1.f;
01687
01688
01689
01690 if( 0 && gv.lgBakesPAH_heat )
01691 {
01692
01693
01694
01695
01696
01697 gv.bin[nd]->GasHeatPhotoEl = 1.e-24*hmi.UV_Cont_rel2_Habing_TH85_depth*
01698 dense.gas_phase[ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((hmi.UV_Cont_rel2_Habing_TH85_depth*
01699 sqrt(phycon.te)/dense.eden),0.73)) + 3.65e-2*pow(phycon.te/1.e4,0.7)/
01700 (1.+2.e-4*(hmi.UV_Cont_rel2_Habing_TH85_depth*sqrt(phycon.te)/dense.eden)))/gv.bin.size() *
01701 gv.GrainHeatScaleFactor;
01702 gv.GasHeatPhotoEl += gv.bin[nd]->GasHeatPhotoEl;
01703 }
01704 }
01705
01706 gv.TotalEden = 0.;
01707 gv.GrnElecDonateMax = 0.f;
01708 gv.GrnElecHoldMax = 0.f;
01709
01710 for( nelem=0; nelem < LIMELM; nelem++ )
01711 {
01712 for( ion=0; ion <= nelem+1; ion++ )
01713 {
01714 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
01715 {
01716 gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
01717 }
01718 }
01719 }
01720
01721
01722 gv.GrainHeatInc = 0.;
01723 gv.GrainHeatDif = 0.;
01724 gv.GrainHeatLya = 0.;
01725 gv.GrainHeatCollSum = 0.;
01726 gv.GrainHeatSum = 0.;
01727 gv.GrainHeatChem = 0.;
01728 gv.GasCoolColl = 0.;
01729 gv.TotalDustHeat = 0.f;
01730 gv.dphmax = 0.f;
01731 gv.dclmax = 0.f;
01732
01733 thermal.heating[0][13] = 0.;
01734 thermal.heating[0][14] = 0.;
01735 thermal.heating[0][25] = 0.;
01736 }
01737
01738 if( nCalledGrainDrive == 0 || gv.lgAnyDustVary )
01739 {
01740 GrainUpdateRadius1();
01741 GrainUpdateRadius2();
01742 }
01743 }
01744
01745 ++nCalledGrainDrive;
01746 return;
01747 }
01748
01749
01750 STATIC void GrainChargeTemp(void)
01751 {
01752 long int i,
01753 ion,
01754 ion_to,
01755 nelem,
01756 nz;
01757 realnum dccool = FLT_MAX;
01758 double delta,
01759 GasHeatNet,
01760 hcon = DBL_MAX,
01761 hla = DBL_MAX,
01762 hots = DBL_MAX,
01763 oldtemp,
01764 oldTotalEden,
01765 ratio,
01766 ThermRatio;
01767
01768 static long int oldZone = -1;
01769 static double oldTe = -DBL_MAX,
01770 oldHeat = -DBL_MAX;
01771
01772 DEBUG_ENTRY( "GrainChargeTemp()" );
01773
01774 if( trace.lgTrace && trace.lgDustBug )
01775 {
01776 fprintf( ioQQQ, "\n GrainChargeTemp called lgSearch%2c\n\n", TorF(conv.lgSearch) );
01777 }
01778
01779 oldTotalEden = gv.TotalEden;
01780
01781
01782 gv.GrainHeatInc = 0.;
01783 gv.GrainHeatDif = 0.;
01784 gv.GrainHeatLya = 0.;
01785 gv.GrainHeatCollSum = 0.;
01786 gv.GrainHeatSum = 0.;
01787 gv.GrainHeatChem = 0.;
01788
01789 gv.GasCoolColl = 0.;
01790 gv.GasHeatPhotoEl = 0.;
01791 gv.GasHeatTherm = 0.;
01792
01793 gv.TotalEden = 0.;
01794
01795 for( nelem=0; nelem < LIMELM; nelem++ )
01796 {
01797 for( ion=0; ion <= nelem+1; ion++ )
01798 {
01799 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
01800 {
01801 gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
01802 }
01803 }
01804 }
01805
01806 gv.HighestIon = HighestIonStage();
01807
01808
01809 GrainUpdateRadius1();
01810
01811 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01812 {
01813 double one;
01814 double ChTdBracketLo = 0., ChTdBracketHi = -DBL_MAX;
01815 long relax = ( conv.lgSearch ) ? 3 : 1;
01816
01817
01818 if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
01819 {
01820 if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.50 )
01821 {
01822 gv.bin[nd]->lgPAHsInIonizedRegion = true;
01823 }
01824 }
01825
01826
01827
01828
01829 NewChargeData(nd);
01830
01831 if( trace.lgTrace && trace.lgDustBug )
01832 {
01833 fprintf( ioQQQ, " >>GrainChargeTemp starting grain %s\n",
01834 gv.bin[nd]->chDstLab );
01835 }
01836
01837 delta = 2.*TOLER;
01838
01839 for( i=0; i < relax*CT_LOOP_MAX && delta > TOLER; ++i )
01840 {
01841 string which;
01842 long j;
01843 double TdBracketLo = 0., TdBracketHi = -DBL_MAX;
01844 double ThresEst = 0.;
01845 oldtemp = gv.bin[nd]->tedust;
01846
01847
01848
01849
01850
01851 GrainCharge(nd,&ThermRatio);
01852
01853 ASSERT( gv.bin[nd]->GrainHeat > 0. );
01854 ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
01855
01856
01857
01858
01859
01860 gv.bin[nd]->lgTdustConverged = false;
01861 for( j=0; j < relax*T_LOOP_MAX; ++j )
01862 {
01863 double oldTemp2 = gv.bin[nd]->tedust;
01864 double oldHeat2 = gv.bin[nd]->GrainHeat;
01865 double oldCool = gv.bin[nd]->GrainGasCool;
01866
01867
01868 GrainTemperature(nd,&dccool,&hcon,&hots,&hla);
01869
01870 gv.bin[nd]->GrainGasCool = dccool;
01871
01872 if( trace.lgTrace && trace.lgDustBug )
01873 {
01874 fprintf( ioQQQ, " >>loop %ld BracketLo %.6e BracketHi %.6e",
01875 j, TdBracketLo, TdBracketHi );
01876 }
01877
01878
01879
01880
01881
01882 if( fabs(gv.bin[nd]->GrainHeat-oldHeat2) < HEAT_TOLER*gv.bin[nd]->GrainHeat &&
01883 fabs(gv.bin[nd]->GrainGasCool-oldCool) < HEAT_TOLER_BIN*thermal.ctot )
01884 {
01885 gv.bin[nd]->lgTdustConverged = true;
01886 if( trace.lgTrace && trace.lgDustBug )
01887 fprintf( ioQQQ, " converged\n" );
01888 break;
01889 }
01890
01891
01892 if( gv.bin[nd]->tedust < oldTemp2 )
01893 TdBracketHi = oldTemp2;
01894 else
01895 TdBracketLo = oldTemp2;
01896
01897
01898
01899
01900
01901
01904
01905 if( TdBracketHi > TdBracketLo )
01906 {
01907
01908
01909 if( ( j >= 2 && TdBracketLo > 0. ) ||
01910 gv.bin[nd]->tedust <= TdBracketLo ||
01911 gv.bin[nd]->tedust >= TdBracketHi )
01912 {
01913 gv.bin[nd]->tedust = (realnum)(0.5*(TdBracketLo + TdBracketHi));
01914 if( trace.lgTrace && trace.lgDustBug )
01915 fprintf( ioQQQ, " bisection\n" );
01916 }
01917 else
01918 {
01919 if( trace.lgTrace && trace.lgDustBug )
01920 fprintf( ioQQQ, " iteration\n" );
01921 }
01922 }
01923 else
01924 {
01925 if( trace.lgTrace && trace.lgDustBug )
01926 fprintf( ioQQQ, " iteration\n" );
01927 }
01928
01929 ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
01930 }
01931
01932 if( gv.bin[nd]->lgTdustConverged )
01933 {
01934
01935 if( gv.bin[nd]->tedust < oldtemp )
01936 ChTdBracketHi = oldtemp;
01937 else
01938 ChTdBracketLo = oldtemp;
01939 }
01940 else
01941 {
01942 bool lgBoundErr;
01943 double y, x = log(gv.bin[nd]->tedust);
01944
01945 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgBoundErr);
01946 gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;
01947
01948 fprintf( ioQQQ," PROBLEM temperature of grain species %s (Tg=%.3eK) not converged\n",
01949 gv.bin[nd]->chDstLab , gv.bin[nd]->tedust );
01950 ConvFail("grai","");
01951 if( lgAbort )
01952 {
01953 return;
01954 }
01955 }
01956
01957 ASSERT( gv.bin[nd]->GrainHeat > 0. );
01958 ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
01959
01960
01961
01962
01963
01964
01966 ratio = gv.bin[nd]->tedust/oldtemp;
01967 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
01968 {
01969 ThresEst += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->ThresInf;
01970 }
01971 delta = ThresEst*TE1RYD/gv.bin[nd]->tedust*(ratio - 1.);
01973 delta = ( delta < 0.9*log(DBL_MAX) ) ?
01974 ThermRatio*fabs(POW2(ratio)*exp(delta)-1.) : DBL_MAX;
01975
01976
01977 if( delta > TOLER )
01978 {
01979 if( trace.lgTrace && trace.lgDustBug )
01980 which = "iteration";
01981
01982
01983
01984
01985
01986
01989
01990 if( ChTdBracketHi > ChTdBracketLo )
01991 {
01992
01993
01994 if( ( i >= 2 && ChTdBracketLo > 0. ) ||
01995 gv.bin[nd]->tedust <= ChTdBracketLo ||
01996 gv.bin[nd]->tedust >= ChTdBracketHi )
01997 {
01998 gv.bin[nd]->tedust = (realnum)(0.5*(ChTdBracketLo + ChTdBracketHi));
01999 if( trace.lgTrace && trace.lgDustBug )
02000 which = "bisection";
02001 }
02002 }
02003 }
02004
02005 if( trace.lgTrace && trace.lgDustBug )
02006 {
02007 fprintf( ioQQQ, " >>GrainChargeTemp finds delta=%.4e, ", delta );
02008 fprintf( ioQQQ, " old/new temp=%.5e %.5e, ", oldtemp, gv.bin[nd]->tedust );
02009 if( delta > TOLER )
02010 fprintf( ioQQQ, "doing another %s\n", which.c_str() );
02011 else
02012 fprintf( ioQQQ, "converged\n" );
02013 }
02014 }
02015 if( delta > TOLER )
02016 {
02017 fprintf( ioQQQ, " PROBLEM charge/temperature not converged for %s zone %.2f\n",
02018 gv.bin[nd]->chDstLab , fnzone );
02019 ConvFail("grai","");
02020 }
02021
02022
02023
02024
02025 if( ionbal.lgGrainIonRecom )
02026 GrainChrgTransferRates(nd);
02027
02028
02029
02030
02031
02032
02033 gv.GrainHeatInc += hcon;
02034
02035 gv.GrainHeatDif += hots;
02036
02037
02038 gv.GrainHeatLya += hla;
02039
02040
02041 gv.GrainHeatCollSum += gv.bin[nd]->GrainHeatColl;
02042
02043
02044
02045
02046 gv.GrainHeatSum += gv.bin[nd]->GrainHeat;
02047
02048
02049 gv.GrainHeatChem += gv.bin[nd]->ChemEn + gv.bin[nd]->ChemEnH2;
02050
02051
02052
02053 gv.GasCoolColl += dccool;
02054 gv.GasHeatPhotoEl += gv.bin[nd]->GasHeatPhotoEl;
02055 gv.GasHeatTherm += gv.bin[nd]->thermionic;
02056
02057
02058
02059 one = 0.;
02060 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02061 {
02062 one += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
02063 gv.bin[nd]->cnv_GR_pCM3;
02064 }
02065
02066 gv.TotalEden += one;
02067 {
02068
02069 enum {DEBUG_LOC=false};
02070
02071 if( DEBUG_LOC )
02072 {
02073 fprintf(ioQQQ," DEBUG grn chr nz\t%.2f\teden\t%.3e\tnd\t%li",
02074 fnzone,
02075 dense.eden,
02076 (unsigned long)nd);
02077 fprintf(ioQQQ,"\tne\t%.2e\tAveDustZ\t%.2e\t%.2e\t%.2e\t%.2e",
02078 one,
02079 gv.bin[nd]->AveDustZ,
02080 gv.bin[nd]->chrg[0]->FracPop,(double)gv.bin[nd]->chrg[0]->DustZ,
02081 gv.bin[nd]->cnv_GR_pCM3);
02082 fprintf(ioQQQ,"\n");
02083 }
02084 }
02085
02086 if( trace.lgTrace && trace.lgDustBug )
02087 {
02088 fprintf(ioQQQ," %s Pot %.5e Thermal %.5e GasCoolColl %.5e" ,
02089 gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot, gv.bin[nd]->GrainHeat, dccool );
02090 fprintf(ioQQQ," GasPEHeat %.5e GasThermHeat %.5e ChemHeat %.5e\n\n" ,
02091 gv.bin[nd]->GasHeatPhotoEl, gv.bin[nd]->thermionic, gv.bin[nd]->ChemEn );
02092 }
02093 }
02094
02095
02096 GasHeatNet = gv.GasHeatPhotoEl + gv.GasHeatTherm - gv.GasCoolColl;
02097
02098 if( !fp_equal(phycon.te,gv.GrnRecomTe) )
02099 {
02100 oldZone = gv.nzone;
02101 oldTe = gv.GrnRecomTe;
02102 oldHeat = gv.GasHeatNet;
02103 }
02104
02105
02106 if( nzone == oldZone && !fp_equal(phycon.te,oldTe) )
02107 {
02108 gv.dHeatdT = (GasHeatNet-oldHeat)/(phycon.te-oldTe);
02109 }
02110
02111
02112 if( nzone != gv.nzone || !fp_equal(phycon.te,gv.GrnRecomTe) ||
02113 fabs(gv.GasHeatNet-GasHeatNet) > HEAT_TOLER*thermal.ctot ||
02114 fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
02115 {
02116
02117
02118
02119 if( fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
02120 {
02121 conv.setConvIonizFail( "grn eden chg" , oldTotalEden, gv.TotalEden);
02122 }
02123 else if( fabs(gv.GasHeatNet-GasHeatNet) > HEAT_TOLER*thermal.ctot )
02124 {
02125 conv.setConvIonizFail( "grn het chg" , gv.GasHeatNet, GasHeatNet);
02126 }
02127 else if( !fp_equal(phycon.te,gv.GrnRecomTe) )
02128 {
02129 conv.setConvIonizFail( "grn ter chg" , gv.GrnRecomTe, phycon.te);
02130 }
02131 else if( nzone != gv.nzone )
02132 {
02133 conv.setConvIonizFail( "grn zon chg" , gv.nzone, nzone );
02134 }
02135 else
02136 TotalInsanity();
02137 }
02138
02139
02140
02141
02142
02143
02144
02145
02146 GrainUpdateRadius2();
02147
02148 gv.nzone = nzone;
02149 gv.GrnRecomTe = phycon.te;
02150 gv.GasHeatNet = GasHeatNet;
02151
02152 #ifdef WD_TEST2
02153 printf("wd test: proton fraction %.5e Total DustZ %.6f heating/cooling rate %.5e %.5e\n",
02154 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
02155 gv.bin[0]->AveDustZ,gv.GasHeatPhotoEl/dense.gas_phase[ipHYDROGEN]/fudge(0),
02156 gv.GasCoolColl/dense.gas_phase[ipHYDROGEN]/fudge(0));
02157 #endif
02158
02159 if( trace.lgTrace )
02160 {
02161
02162 enum {DEBUG_LOC=false};
02163
02164 if( DEBUG_LOC )
02165 {
02166 fprintf( ioQQQ, " %.2f Grain surface charge transfer rates\n", fnzone );
02167
02168 for( nelem=0; nelem < LIMELM; nelem++ )
02169 {
02170 if( dense.lgElmtOn[nelem] )
02171 {
02172 printf( " %s:", elementnames.chElementSym[nelem] );
02173 for( ion=dense.IonLow[nelem]; ion <= dense.IonHigh[nelem]; ++ion )
02174 {
02175 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
02176 {
02177 if( gv.GrainChTrRate[nelem][ion][ion_to] > 0.f )
02178 {
02179 printf( " %ld->%ld %.2e", ion, ion_to,
02180 gv.GrainChTrRate[nelem][ion][ion_to] );
02181 }
02182 }
02183 }
02184 printf( "\n" );
02185 }
02186 }
02187 }
02188
02189 fprintf( ioQQQ, " %.2f Grain contribution to electron density %.2e\n",
02190 fnzone , gv.TotalEden );
02191
02192 fprintf( ioQQQ, " Grain electons: " );
02193 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02194 {
02195 double sum = 0.;
02196 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02197 {
02198 sum += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
02199 gv.bin[nd]->cnv_GR_pCM3;
02200 }
02201 fprintf( ioQQQ, " %.2e", sum );
02202 }
02203 fprintf( ioQQQ, "\n" );
02204
02205 fprintf( ioQQQ, " Grain potentials:" );
02206 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02207 {
02208 fprintf( ioQQQ, " %.2e", gv.bin[nd]->dstpot );
02209 }
02210 fprintf( ioQQQ, "\n" );
02211
02212 fprintf( ioQQQ, " Grain temperatures:" );
02213 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02214 {
02215 fprintf( ioQQQ, " %.2e", gv.bin[nd]->tedust );
02216 }
02217 fprintf( ioQQQ, "\n" );
02218
02219 fprintf( ioQQQ, " GrainCollCool: %.6e\n", gv.GasCoolColl );
02220 }
02221
02222
02223
02224
02225
02226
02227
02228 return;
02229 }
02230
02231
02232 STATIC void GrainCharge(size_t nd,
02233 double *ThermRatio)
02234 {
02235 bool lgBigError,
02236 lgInitial;
02237 long backup,
02238 i,
02239 loopMax,
02240 newZlo,
02241 nz,
02242 power,
02243 stride,
02244 stride0,
02245 Zlo;
02246 double crate,
02247 csum1,
02248 csum1a,
02249 csum1b,
02250 csum2,
02251 csum3,
02252 netloss0 = -DBL_MAX,
02253 netloss1 = -DBL_MAX,
02254 rate_up[NCHU],
02255 rate_dn[NCHU],
02256 step;
02257
02258 DEBUG_ENTRY( "GrainCharge()" );
02259
02260
02261 if( trace.lgTrace && trace.lgDustBug )
02262 {
02263 fprintf( ioQQQ, " Charge loop, search %c,", TorF(conv.lgSearch) );
02264 }
02265
02266 ASSERT( nd < gv.bin.size() );
02267
02268 for( nz=0; nz < NCHS; nz++ )
02269 {
02270 gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
02271 }
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321 lgInitial = ( gv.bin[nd]->chrg[0]->DustZ == LONG_MIN );
02322
02323 backup = gv.bin[nd]->nChrg;
02324 gv.bin[nd]->nChrg = MAX2(gv.bin[nd]->nChrg,3);
02325
02326 stride0 = gv.bin[nd]->nChrg-1;
02327
02328
02329 if( lgInitial )
02330 {
02331 double xxx;
02332 step = MAX2((double)(-gv.bin[nd]->LowestZg),1.);
02333 power = (int)(log(step)/log((double)stride0));
02334 power = MAX2(power,0);
02335 xxx = powi((double)stride0,power);
02336 stride = nint(xxx);
02337 Zlo = gv.bin[nd]->LowestZg;
02338 }
02339 else
02340 {
02341
02342 stride = 1;
02343 Zlo = gv.bin[nd]->chrg[0]->DustZ;
02344 }
02345 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
02346
02347
02348 for( i=0; i < BRACKET_MAX; i++ )
02349 {
02350 netloss0 = rate_up[0] - rate_dn[0];
02351 netloss1 = rate_up[gv.bin[nd]->nChrg-1] - rate_dn[gv.bin[nd]->nChrg-1];
02352
02353 if( netloss0*netloss1 <= 0. )
02354 break;
02355
02356 if( netloss1 > 0. )
02357 Zlo += (gv.bin[nd]->nChrg-1)*stride;
02358
02359 if( i > 0 )
02360 stride *= stride0;
02361
02362 if( netloss1 < 0. )
02363 Zlo -= (gv.bin[nd]->nChrg-1)*stride;
02364
02365 Zlo = MAX2(Zlo,gv.bin[nd]->LowestZg);
02366 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
02367 }
02368
02369 if( netloss0*netloss1 > 0. ) {
02370 fprintf( ioQQQ, " insanity: could not bracket grain charge for %s\n", gv.bin[nd]->chDstLab );
02371 ShowMe();
02372 cdEXIT(EXIT_FAILURE);
02373 }
02374
02375
02376 while( stride > 1 )
02377 {
02378 stride /= stride0;
02379
02380 netloss0 = rate_up[0] - rate_dn[0];
02381 for( nz=0; nz < gv.bin[nd]->nChrg-1; nz++ )
02382 {
02383 netloss1 = rate_up[nz+1] - rate_dn[nz+1];
02384
02385 if( netloss0*netloss1 <= 0. )
02386 {
02387 Zlo = gv.bin[nd]->chrg[nz]->DustZ;
02388 break;
02389 }
02390
02391 netloss0 = netloss1;
02392 }
02393 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
02394 }
02395
02396 ASSERT( netloss0*netloss1 <= 0. );
02397
02398 gv.bin[nd]->nChrg = backup;
02399
02400
02401 loopMax = ( lgInitial ) ? 4*gv.bin[nd]->nChrg : 2*gv.bin[nd]->nChrg;
02402
02403 lgBigError = true;
02404 for( i=0; i < loopMax; i++ )
02405 {
02406 GetFracPop( nd, Zlo, rate_up, rate_dn, &newZlo );
02407
02408 if( newZlo == Zlo )
02409 {
02410 lgBigError = false;
02411 break;
02412 }
02413
02414 Zlo = newZlo;
02415 UpdatePot( nd, Zlo, 1, rate_up, rate_dn );
02416 }
02417
02418 if( lgBigError ) {
02419 fprintf( ioQQQ, " insanity: could not converge grain charge for %s\n", gv.bin[nd]->chDstLab );
02420 ShowMe();
02421 cdEXIT(EXIT_FAILURE);
02422 }
02423
02426 gv.bin[nd]->lgChrgConverged = true;
02427
02428 gv.bin[nd]->AveDustZ = 0.;
02429 crate = csum3 = 0.;
02430 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02431 {
02432 double d[4];
02433
02434 gv.bin[nd]->AveDustZ += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->DustZ;
02435
02436 crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
02437 csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[3];
02438 }
02439 gv.bin[nd]->dstpot = chrg2pot(gv.bin[nd]->AveDustZ,nd);
02440 *ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;
02441
02442 ASSERT( *ThermRatio >= 0. );
02443
02444 if( trace.lgTrace && trace.lgDustBug )
02445 {
02446 double d[4];
02447
02448 fprintf( ioQQQ, "\n" );
02449
02450 crate = csum1a = csum1b = csum2 = csum3 = 0.;
02451 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02452 {
02453 crate += gv.bin[nd]->chrg[nz]->FracPop*
02454 GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
02455 csum1a += gv.bin[nd]->chrg[nz]->FracPop*d[0];
02456 csum1b += gv.bin[nd]->chrg[nz]->FracPop*d[1];
02457 csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[2];
02458 csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[3];
02459 }
02460
02461 fprintf( ioQQQ, " ElecEm rate1a=%.4e, rate1b=%.4e, ", csum1a, csum1b );
02462 fprintf( ioQQQ, "rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum2, csum3, crate );
02463 if( crate > 0. )
02464 {
02465 fprintf( ioQQQ, " rate1a/sum=%.4e, rate1b/sum=%.4e, ", csum1a/crate, csum1b/crate );
02466 fprintf( ioQQQ, "rate2/sum=%.4e, rate3/sum=%.4e\n", csum2/crate, csum3/crate );
02467 }
02468
02469 crate = csum1 = csum2 = 0.;
02470 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02471 {
02472 crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecRecomb1(nd,nz,&d[0],&d[1]);
02473 csum1 += gv.bin[nd]->chrg[nz]->FracPop*d[0];
02474 csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[1];
02475 }
02476
02477 fprintf( ioQQQ, " ElecRc rate1=%.4e, rate2=%.4e, sum=%.4e\n", csum1, csum2, crate );
02478 if( crate > 0. )
02479 {
02480 fprintf( ioQQQ, " rate1/sum=%.4e, rate2/sum=%.4e\n", csum1/crate, csum2/crate );
02481 }
02482
02483 fprintf( ioQQQ, " Charging rates:" );
02484 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02485 {
02486 fprintf( ioQQQ, " Zg %ld up %.4e dn %.4e",
02487 gv.bin[nd]->chrg[nz]->DustZ, rate_up[nz], rate_dn[nz] );
02488 }
02489 fprintf( ioQQQ, "\n" );
02490
02491 fprintf( ioQQQ, " FracPop: fnzone %.2f nd %ld AveZg %.5e",
02492 fnzone, (unsigned long)nd, gv.bin[nd]->AveDustZ );
02493 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02494 {
02495 fprintf( ioQQQ, " Zg %ld %.5f", Zlo+nz, gv.bin[nd]->chrg[nz]->FracPop );
02496 }
02497 fprintf( ioQQQ, "\n" );
02498
02499 fprintf( ioQQQ, " >Grain potential:%12.12s %.6fV",
02500 gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot*EVRYD );
02501 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02502 {
02503 fprintf( ioQQQ, " Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V",
02504 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInf*EVRYD,
02505 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInfVal*EVRYD );
02506 }
02507 fprintf( ioQQQ, "\n" );
02508 }
02509 return;
02510 }
02511
02512
02513
02514 STATIC double GrainElecRecomb1(size_t nd,
02515 long nz,
02516 double *sum1,
02517 double *sum2)
02518 {
02519 long ion,
02520 nelem;
02521 double eta,
02522 rate,
02523 Stick,
02524 ve,
02525 xi;
02526
02527 DEBUG_ENTRY( "GrainElecRecomb1()" );
02528
02529 ASSERT( nd < gv.bin.size() );
02530 ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
02531
02532
02533
02534 if( gv.bin[nd]->chrg[nz]->RSum1 >= 0. )
02535 {
02536 *sum1 = gv.bin[nd]->chrg[nz]->RSum1;
02537 *sum2 = gv.bin[nd]->chrg[nz]->RSum2;
02538 rate = *sum1 + *sum2;
02539 return rate;
02540 }
02541
02542
02543 ion = -1;
02544
02545
02546 ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);
02547
02548 Stick = ( gv.bin[nd]->chrg[nz]->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
02549
02550
02551 GrainScreen(ion,nd,nz,&eta,&xi);
02552
02553 *sum1 = ( gv.bin[nd]->chrg[nz]->DustZ > gv.bin[nd]->LowestZg ) ? Stick*dense.eden*ve*eta : 0.;
02554
02555
02556 *sum2 = 0.;
02557
02558 #ifndef IGNORE_GRAIN_ION_COLLISIONS
02559 for( ion=0; ion <= LIMELM; ion++ )
02560 {
02561 double CollisionRateAll = 0.;
02562
02563 for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
02564 {
02565 if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
02566 gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] > ion )
02567 {
02568
02569 CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*GetAveVelocity( dense.AtomicWeight[nelem] )*
02570 (double)(gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]-ion);
02571 }
02572 }
02573
02574 if( CollisionRateAll > 0. )
02575 {
02576
02577
02578
02579 GrainScreen(ion,nd,nz,&eta,&xi);
02580 *sum2 += CollisionRateAll*eta;
02581 }
02582 }
02583 #endif
02584
02585 rate = *sum1 + *sum2;
02586
02587
02588 gv.bin[nd]->chrg[nz]->RSum1 = *sum1;
02589 gv.bin[nd]->chrg[nz]->RSum2 = *sum2;
02590
02591 ASSERT( *sum1 >= 0. && *sum2 >= 0. );
02592 return rate;
02593 }
02594
02595
02596
02597 STATIC double GrainElecEmis1(size_t nd,
02598 long nz,
02599 double *sum1a,
02600 double *sum1b,
02601 double *sum2,
02602 double *sum3)
02603 {
02604 long int i,
02605 ion,
02606 ipLo,
02607 nelem;
02608 double eta,
02609 rate,
02610 xi;
02611
02612 DEBUG_ENTRY( "GrainElecEmis1()" );
02613
02614 ASSERT( nd < gv.bin.size() );
02615 ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
02616
02617
02618
02619 if( gv.bin[nd]->chrg[nz]->ESum1a >= 0. )
02620 {
02621 *sum1a = gv.bin[nd]->chrg[nz]->ESum1a;
02622 *sum1b = gv.bin[nd]->chrg[nz]->ESum1b;
02623 *sum2 = gv.bin[nd]->chrg[nz]->ESum2;
02624
02625 *sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;
02626 rate = *sum1a + *sum1b + *sum2 + *sum3;
02627 return rate;
02628 }
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640
02641 *sum1a = 0.;
02642 ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
02643 for( i=ipLo; i < rfield.nflux; i++ )
02644 {
02645 # ifdef WD_TEST2
02646 *sum1a += rfield.flux[0][i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i];
02647 # else
02648 *sum1a += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i];
02649 # endif
02650 }
02651
02652 *sum1a /= gv.bin[nd]->IntArea/4.;
02653
02654 *sum1b = 0.;
02655 if( gv.bin[nd]->chrg[nz]->DustZ <= -1 )
02656 {
02657 ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
02658 for( i=ipLo; i < rfield.nflux; i++ )
02659 {
02660
02661 # ifdef WD_TEST2
02662 *sum1b += rfield.flux[0][i]*gv.bin[nd]->chrg[nz]->cs_pdt[i];
02663 # else
02664 *sum1b += rfield.SummedCon[i]*gv.bin[nd]->chrg[nz]->cs_pdt[i];
02665 # endif
02666 }
02667 *sum1b /= gv.bin[nd]->IntArea/4.;
02668 }
02669
02670
02671 *sum2 = 0.;
02672 # ifndef IGNORE_GRAIN_ION_COLLISIONS
02673 for( ion=0; ion <= LIMELM; ion++ )
02674 {
02675 double CollisionRateAll = 0.;
02676
02677 for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
02678 {
02679 if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
02680 ion > gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] )
02681 {
02682
02683 CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*GetAveVelocity( dense.AtomicWeight[nelem] )*
02684 (double)(ion-gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]);
02685 }
02686 }
02687
02688 if( CollisionRateAll > 0. )
02689 {
02690
02691
02692
02693 GrainScreen(ion,nd,nz,&eta,&xi);
02694 *sum2 += CollisionRateAll*eta;
02695 }
02696 }
02697 # endif
02698
02699
02700
02701
02702 *sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;
02703
02706 rate = *sum1a + *sum1b + *sum2 + *sum3;
02707
02708
02709 gv.bin[nd]->chrg[nz]->ESum1a = *sum1a;
02710 gv.bin[nd]->chrg[nz]->ESum1b = *sum1b;
02711 gv.bin[nd]->chrg[nz]->ESum2 = *sum2;
02712
02713 ASSERT( *sum1a >= 0. && *sum1b >= 0. && *sum2 >= 0. && *sum3 >= 0. );
02714 return rate;
02715 }
02716
02717
02718
02719
02720 STATIC void GrainScreen(long ion,
02721 size_t nd,
02722 long nz,
02723 double *eta,
02724 double *xi)
02725 {
02726
02727
02728
02729
02730 long ind = ion+1;
02731
02732 DEBUG_ENTRY( "GrainScreen()" );
02733
02734 ASSERT( ind >= 0 && ind < LIMELM+2 );
02735
02736 if( gv.bin[nd]->chrg[nz]->eta[ind] > 0. )
02737 {
02738 *eta = gv.bin[nd]->chrg[nz]->eta[ind];
02739 *xi = gv.bin[nd]->chrg[nz]->xi[ind];
02740 return;
02741 }
02742
02743
02744
02745 if( ion == 0 )
02746 {
02747 *eta = 1.;
02748 *xi = 1.;
02749 }
02750 else
02751 {
02752
02753
02754
02755
02756
02757
02758 double nu = (double)gv.bin[nd]->chrg[nz]->DustZ/(double)ion;
02759 double tau = gv.bin[nd]->Capacity*BOLTZMANN*phycon.te*1.e-7/POW2((double)ion*ELEM_CHARGE);
02760 if( nu < 0. )
02761 {
02762 *eta = (1. - nu/tau)*(1. + sqrt(2./(tau - 2.*nu)));
02763 *xi = (1. - nu/(2.*tau))*(1. + 1./sqrt(tau - nu));
02764 }
02765 else if( nu == 0. )
02766 {
02767 *eta = 1. + sqrt(PI/(2.*tau));
02768 *xi = 1. + 0.75*sqrt(PI/(2.*tau));
02769 }
02770 else
02771 {
02772 double theta_nu = ThetaNu(nu);
02773
02774 double xxx = 1. + 1./sqrt(4.*tau+3.*nu);
02775 *eta = POW2(xxx)*exp(-theta_nu/tau);
02776 # ifdef WD_TEST2
02777 *xi = (1. + nu/(2.*tau))*(1. + 1./sqrt(3./(2.*tau)+3.*nu))*exp(-theta_nu/tau);
02778 # else
02779
02780
02781
02782 xxx = 0.25*pow(nu/tau,0.75)/(pow(nu/tau,0.75) + pow((25.+3.*nu)/5.,0.75)) +
02783 (1. + 0.75*sqrt(PI/(2.*tau)))/(1. + sqrt(PI/(2.*tau)));
02784 *xi = (MIN2(xxx,1.) + theta_nu/(2.*tau))*(*eta);
02785 # endif
02786 }
02787
02788 ASSERT( *eta >= 0. && *xi >= 0. );
02789 }
02790
02791 gv.bin[nd]->chrg[nz]->eta[ind] = *eta;
02792 gv.bin[nd]->chrg[nz]->xi[ind] = *xi;
02793
02794 return;
02795 }
02796
02797
02798 STATIC double ThetaNu(double nu)
02799 {
02800 double theta_nu;
02801
02802 DEBUG_ENTRY( "ThetaNu()" );
02803
02804 if( nu > 0. )
02805 {
02806 double xxx;
02807 const double REL_TOLER = 10.*DBL_EPSILON;
02808
02809
02810 double xi_nu = 1. + 1./sqrt(3.*nu);
02811 double xi_nu2 = POW2(xi_nu);
02812 do
02813 {
02814 double old = xi_nu;
02815
02816
02817 double fnu = 2.*xi_nu2 - 1. - nu*xi_nu*POW2(xi_nu2 - 1.);
02818 double dfdxi = 4.*xi_nu - nu*((5.*xi_nu2 - 6.)*xi_nu2 + 1.);
02819 xi_nu -= fnu/dfdxi;
02820 xi_nu2 = POW2(xi_nu);
02821 xxx = fabs(old-xi_nu);
02822 } while( xxx > REL_TOLER*xi_nu );
02823
02824 theta_nu = nu/xi_nu - 1./(2.*xi_nu2*(xi_nu2-1.));
02825 }
02826 else
02827 {
02828 theta_nu = 0.;
02829 }
02830 return theta_nu;
02831 }
02832
02833
02834
02835 STATIC void UpdatePot(size_t nd,
02836 long Zlo,
02837 long stride,
02838 double rate_up[],
02839 double rate_dn[])
02840 {
02841 long nz,
02842 Zg;
02843 double BoltzFac,
02844 HighEnergy;
02845
02846 DEBUG_ENTRY( "UpdatePot()" );
02847
02848 ASSERT( nd < gv.bin.size() );
02849 ASSERT( Zlo >= gv.bin[nd]->LowestZg );
02850 ASSERT( stride >= 1 );
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861 if( trace.lgTrace && trace.lgDustBug )
02862 {
02863 fprintf( ioQQQ, " %ld/%ld", Zlo, stride );
02864 }
02865
02866 if( gv.bin[nd]->nfill < rfield.nflux )
02867 {
02868 InitBinAugerData( nd, gv.bin[nd]->nfill, rfield.nflux );
02869 gv.bin[nd]->nfill = rfield.nflux;
02870 }
02871
02872 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02873 {
02874 long ind, zz;
02875 double d[4];
02876 ChargeBin *ptr;
02877
02878 Zg = Zlo + nz*stride;
02879
02880
02881 ind = NCHS-1;
02882 for( zz=0; zz < NCHS-1; zz++ )
02883 {
02884 if( gv.bin[nd]->chrg[zz]->DustZ == Zg )
02885 {
02886 ind = zz;
02887 break;
02888 }
02889 }
02890
02891
02892
02893 ptr = gv.bin[nd]->chrg[ind];
02894
02895
02896 for( zz=ind-1; zz >= nz; zz-- )
02897 gv.bin[nd]->chrg[zz+1] = gv.bin[nd]->chrg[zz];
02898
02899 gv.bin[nd]->chrg[nz] = ptr;
02900
02901 if( gv.bin[nd]->chrg[nz]->DustZ != Zg )
02902 UpdatePot1(nd,nz,Zg,0);
02903 else if( gv.bin[nd]->chrg[nz]->nfill < rfield.nflux )
02904 UpdatePot1(nd,nz,Zg,gv.bin[nd]->chrg[nz]->nfill);
02905
02906 UpdatePot2(nd,nz);
02907
02908 rate_up[nz] = GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
02909 rate_dn[nz] = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
02910
02911
02912 ASSERT( gv.bin[nd]->chrg[nz]->DustZ == Zg );
02913 ASSERT( gv.bin[nd]->chrg[nz]->nfill >= rfield.nflux );
02914 ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
02915 }
02916
02917
02918
02919
02920
02921
02922 BoltzFac = (-log(CONSERV_TOL) + 8.)*BOLTZMANN/EN1RYD;
02923 HighEnergy = 0.;
02924 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02925 {
02926
02927 HighEnergy = MAX2(HighEnergy,
02928 MAX2(gv.bin[nd]->chrg[nz]->ThresInfInc,0.) + BoltzFac*MAX2(phycon.te,gv.bin[nd]->tedust));
02929 }
02930 HighEnergy = MIN2(HighEnergy,rfield.anu[rfield.nupper-1]);
02931 gv.bin[nd]->qnflux2 = ipoint(HighEnergy);
02932 gv.bin[nd]->qnflux = MAX2(rfield.nflux,gv.bin[nd]->qnflux2);
02933
02934 ASSERT( gv.bin[nd]->qnflux <= rfield.nupper-1 );
02935 return;
02936 }
02937
02938
02939
02940 STATIC void GetFracPop(size_t nd,
02941 long Zlo,
02942 double rate_up[],
02943 double rate_dn[],
02944 long *newZlo)
02945 {
02946 bool lgRedo;
02947 long i,
02948 nz;
02949 double netloss[2],
02950 pop[2][NCHU-1];
02951
02952
02953 DEBUG_ENTRY( "GetFracPop()" );
02954
02955 ASSERT( nd < gv.bin.size() );
02956 ASSERT( Zlo >= gv.bin[nd]->LowestZg );
02957
02958
02959
02960
02961
02962
02963
02964
02965 do
02966 {
02967 for( i=0; i < 2; i++ )
02968 {
02969 long j, k;
02970 double sum;
02971
02972 sum = pop[i][0] = 1.;
02973 for( j=1; j < gv.bin[nd]->nChrg-1; j++ )
02974 {
02975 nz = i + j;
02976 if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
02977 {
02978 pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
02979 sum += pop[i][j];
02980 }
02981 else
02982 {
02983 for( k=0; k < j; k++ )
02984 {
02985 pop[i][k] = 0.;
02986 }
02987 pop[i][j] = 1.;
02988 sum = pop[i][j];
02989 }
02990
02991 if( pop[i][j] > sqrt(DBL_MAX) )
02992 {
02993 for( k=0; k <= j; k++ )
02994 {
02995 pop[i][k] /= DBL_MAX/10.;
02996 }
02997 sum /= DBL_MAX/10.;
02998 }
02999 }
03000 netloss[i] = 0.;
03001 for( j=0; j < gv.bin[nd]->nChrg-1; j++ )
03002 {
03003 nz = i + j;
03004 pop[i][j] /= sum;
03005 netloss[i] += pop[i][j]*(rate_up[nz] - rate_dn[nz]);
03006 }
03007 }
03008
03009
03010
03011 if( netloss[0]*netloss[1] > 0. )
03012 *newZlo = ( netloss[1] > 0. ) ? Zlo + 1 : Zlo - 1;
03013 else
03014 *newZlo = Zlo;
03015
03016
03017
03018
03019
03020
03021 if( gv.bin[nd]->nChrg > 2 &&
03022 ( *newZlo < gv.bin[nd]->LowestZg ||
03023 ( *newZlo == Zlo && pop[1][gv.bin[nd]->nChrg-2] < DBL_EPSILON ) ) )
03024 {
03025 gv.bin[nd]->nChrg--;
03026 lgRedo = true;
03027 }
03028 else
03029 {
03030 lgRedo = false;
03031 }
03032
03033 # if 0
03034 printf( " fnzone %.2f nd %ld Zlo %ld newZlo %ld netloss %.4e %.4e nChrg %ld lgRedo %d\n",
03035 fnzone, nd, Zlo, *newZlo, netloss[0], netloss[1], gv.bin[nd]->nChrg, lgRedo );
03036 # endif
03037 }
03038 while( lgRedo );
03039
03040 if( *newZlo < gv.bin[nd]->LowestZg )
03041 {
03042 fprintf( ioQQQ, " could not converge charge state populations for %s\n", gv.bin[nd]->chDstLab );
03043 ShowMe();
03044 cdEXIT(EXIT_FAILURE);
03045 }
03046
03047 if( *newZlo == Zlo )
03048 {
03049 double frac0, frac1;
03050 # ifndef NDEBUG
03051 double test1, test2, test3, x1, x2;
03052 # endif
03053
03054
03055 frac0 = netloss[1]/(netloss[1]-netloss[0]);
03056 frac1 = -netloss[0]/(netloss[1]-netloss[0]);
03057
03058 gv.bin[nd]->chrg[0]->FracPop = frac0*pop[0][0];
03059 gv.bin[nd]->chrg[gv.bin[nd]->nChrg-1]->FracPop = frac1*pop[1][gv.bin[nd]->nChrg-2];
03060 for( nz=1; nz < gv.bin[nd]->nChrg-1; nz++ )
03061 {
03062 gv.bin[nd]->chrg[nz]->FracPop = frac0*pop[0][nz] + frac1*pop[1][nz-1];
03063 }
03064
03065 # ifndef NDEBUG
03066 test1 = test2 = test3 = 0.;
03067 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
03068 {
03069 ASSERT( gv.bin[nd]->chrg[nz]->FracPop >= 0. );
03070 test1 += gv.bin[nd]->chrg[nz]->FracPop;
03071 test2 += gv.bin[nd]->chrg[nz]->FracPop*rate_up[nz];
03072 test3 += gv.bin[nd]->chrg[nz]->FracPop*(rate_up[nz]-rate_dn[nz]);
03073 }
03074 x1 = fabs(test1-1.);
03075 x2 = fabs( safe_div(test3, test2, 0.) );
03076 ASSERT( MAX2(x1,x2) < 10.*sqrt((double)gv.bin[nd]->nChrg)*DBL_EPSILON );
03077 # endif
03078 }
03079 return;
03080 }
03081
03082
03083
03084
03085
03086
03087
03088
03089
03090
03091
03092 STATIC void UpdatePot1(size_t nd,
03093 long nz,
03094 long Zg,
03095 long ipStart)
03096 {
03097 long i,
03098 ipLo,
03099 nfill;
03100 double c1,
03101 cnv_GR_pH,
03102 d[2],
03103 DustWorkFcn,
03104 Elo,
03105 Emin,
03106 ThresInf,
03107 ThresInfVal;
03108
03109 double *anu = rfield.anu;
03110
03111
03112
03113
03114 const double INV_DELTA_E = EVRYD/3.;
03115 const double CS_PDT = 1.2e-17;
03116
03117 DEBUG_ENTRY( "UpdatePot1()" );
03118
03119
03120
03121
03122
03123
03124
03125 if( ipStart == 0 )
03126 {
03127 gv.bin[nd]->chrg[nz]->DustZ = Zg;
03128
03129
03130 memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
03131 memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );
03132
03133 GetPotValues(nd,Zg,&gv.bin[nd]->chrg[nz]->ThresInf,&gv.bin[nd]->chrg[nz]->ThresInfVal,
03134 &gv.bin[nd]->chrg[nz]->ThresSurf,&gv.bin[nd]->chrg[nz]->ThresSurfVal,
03135 &gv.bin[nd]->chrg[nz]->PotSurf,&gv.bin[nd]->chrg[nz]->Emin,INCL_TUNNEL);
03136
03137
03138
03139 GetPotValues(nd,Zg-1,&gv.bin[nd]->chrg[nz]->ThresInfInc,&d[0],&gv.bin[nd]->chrg[nz]->ThresSurfInc,
03140 &d[1],&gv.bin[nd]->chrg[nz]->PotSurfInc,&gv.bin[nd]->chrg[nz]->EminInc,NO_TUNNEL);
03141
03142 gv.bin[nd]->chrg[nz]->ipThresInfVal =
03143 hunt_bisect( rfield.anu, rfield.nupper, gv.bin[nd]->chrg[nz]->ThresInfVal ) + 1;
03144 }
03145
03146 ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
03147
03148
03149 gv.bin[nd]->chrg[nz]->nfill = rfield.nflux;
03150 nfill = gv.bin[nd]->chrg[nz]->nfill;
03151
03152
03153 long len = nfill + 10 - ipLo;
03154 if( ipStart == 0 )
03155 {
03156 gv.bin[nd]->chrg[nz]->yhat.reserve( len );
03157 gv.bin[nd]->chrg[nz]->yhat_primary.reserve( len );
03158 gv.bin[nd]->chrg[nz]->ehat.reserve( len );
03159 gv.bin[nd]->chrg[nz]->yhat.alloc( ipLo, nfill );
03160 gv.bin[nd]->chrg[nz]->yhat_primary.alloc( ipLo, nfill );
03161 gv.bin[nd]->chrg[nz]->ehat.alloc( ipLo, nfill );
03162 }
03163 else
03164 {
03165 gv.bin[nd]->chrg[nz]->yhat.realloc( nfill );
03166 gv.bin[nd]->chrg[nz]->yhat_primary.realloc( nfill );
03167 gv.bin[nd]->chrg[nz]->ehat.realloc( nfill );
03168 }
03169
03170 double GrainPot = chrg2pot(Zg,nd);
03171
03172 if( nfill > ipLo )
03173 {
03174 DustWorkFcn = gv.bin[nd]->DustWorkFcn;
03175 Elo = -gv.bin[nd]->chrg[nz]->PotSurf;
03176 ThresInfVal = gv.bin[nd]->chrg[nz]->ThresInfVal;
03177 Emin = gv.bin[nd]->chrg[nz]->Emin;
03178
03182 ASSERT( gv.bin[nd]->sd[0]->ipLo <= ipLo );
03183
03184 for( i=max(ipLo,ipStart); i < nfill; i++ )
03185 {
03186 long n;
03187 unsigned int ns=0;
03188 double Yp1,Ys1,Ehp,Ehs,yp,ya,ys,eyp,eya,eys;
03189 double yzero,yone,Ehi,Ehi_band,Wcorr,Eel;
03190 ShellData *sptr = gv.bin[nd]->sd[ns];
03191
03192 yp = ya = ys = 0.;
03193 eyp = eya = eys = 0.;
03194
03195
03196 Ehi = Ehi_band = anu[i] - ThresInfVal - Emin;
03197 Wcorr = ThresInfVal + Emin - GrainPot;
03198 Eel = anu[i] - DustWorkFcn;
03199 yzero = y0b( nd, nz, i );
03200 yone = sptr->y01[i];
03201 Yfunc(nd,nz,yzero*yone,sptr->p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
03202 yp += Yp1;
03203 ys += Ys1;
03204 eyp += Yp1*Ehp;
03205 eys += Ys1*Ehs;
03206
03207
03208 unsigned int nsmax = gv.bin[nd]->sd.size();
03209 for( ns=1; ns < nsmax; ns++ )
03210 {
03211 sptr = gv.bin[nd]->sd[ns];
03212
03213 if( i < sptr->ipLo )
03214 continue;
03215
03216 Ehi = Ehi_band + Wcorr - sptr->ionPot;
03217 Eel = anu[i] - sptr->ionPot;
03218 Yfunc(nd,nz,sptr->y01[i],sptr->p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
03219 yp += Yp1;
03220 ys += Ys1;
03221 eyp += Yp1*Ehp;
03222 eys += Ys1*Ehs;
03223
03224
03225 long nmax = sptr->nData;
03226 for( n=0; n < nmax; n++ )
03227 {
03228 double max = sptr->AvNr[n]*sptr->p[i];
03229 Ehi = sptr->Ener[n] - GrainPot;
03230 Eel = sptr->Ener[n];
03231 Yfunc(nd,nz,sptr->y01A[n][i],max,Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
03232 ya += Yp1;
03233 ys += Ys1;
03234 eya += Yp1*Ehp;
03235 eys += Ys1*Ehs;
03236 }
03237 }
03238
03239 gv.bin[nd]->chrg[nz]->yhat[i] = (realnum)(yp + ya + ys);
03240 gv.bin[nd]->chrg[nz]->yhat_primary[i] = min((realnum)yp,1.f);
03241 gv.bin[nd]->chrg[nz]->ehat[i] = ( gv.bin[nd]->chrg[nz]->yhat[i] > 0. ) ?
03242 (realnum)((eyp + eya + eys)/gv.bin[nd]->chrg[nz]->yhat[i]) : 0.f;
03243
03244 ASSERT( yp <= 1.00001 );
03245 ASSERT( gv.bin[nd]->chrg[nz]->ehat[i] >= 0.f );
03246 }
03247 }
03248
03249 if( ipStart == 0 )
03250 {
03251
03252
03253
03254 gv.bin[nd]->chrg[nz]->ipThresInf =
03255 hunt_bisect( rfield.anu, rfield.nupper, gv.bin[nd]->chrg[nz]->ThresInf ) + 1;
03256 }
03257
03258 ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
03259
03260 len = nfill + 10 - ipLo;
03261
03262 if( Zg <= -1 )
03263 {
03264
03265 if( ipStart == 0 )
03266 {
03267 gv.bin[nd]->chrg[nz]->cs_pdt.reserve( len );
03268 gv.bin[nd]->chrg[nz]->cs_pdt.alloc( ipLo, nfill );
03269 }
03270 else
03271 {
03272 gv.bin[nd]->chrg[nz]->cs_pdt.realloc( nfill );
03273 }
03274
03275 if( nfill > ipLo )
03276 {
03277 c1 = -CS_PDT*(double)Zg;
03278 ThresInf = gv.bin[nd]->chrg[nz]->ThresInf;
03279 cnv_GR_pH = gv.bin[nd]->cnv_GR_pH;
03280
03281 for( i=max(ipLo,ipStart); i < nfill; i++ )
03282 {
03283 double x = (anu[i] - ThresInf)*INV_DELTA_E;
03284 double cs = c1*x/POW2(1.+(1./3.)*POW2(x));
03285
03286
03287 gv.bin[nd]->chrg[nz]->cs_pdt[i] = MAX2(cs,0.)*cnv_GR_pH;
03288 }
03289 }
03290 }
03291
03292
03293 if( ipStart == 0 )
03294 {
03295 gv.bin[nd]->chrg[nz]->fac1.reserve( len );
03296 gv.bin[nd]->chrg[nz]->fac2.reserve( len );
03297 gv.bin[nd]->chrg[nz]->fac1.alloc( ipLo, nfill );
03298 gv.bin[nd]->chrg[nz]->fac2.alloc( ipLo, nfill );
03299 }
03300 else
03301 {
03302 gv.bin[nd]->chrg[nz]->fac1.realloc( nfill );
03303 gv.bin[nd]->chrg[nz]->fac2.realloc( nfill );
03304 }
03305
03306 if( nfill > ipLo )
03307 {
03308 for( i=max(ipLo,ipStart); i < nfill; i++ )
03309 {
03310 double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;
03311
03312
03313 PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);
03314
03315 gv.bin[nd]->chrg[nz]->fac1[i] = cs_tot*anu[i]-cs1*cool1-cs2*cool2;
03316 gv.bin[nd]->chrg[nz]->fac2[i] = cs1*ehat1 + cs2*ehat2;
03317
03318 ASSERT( gv.bin[nd]->chrg[nz]->fac1[i] >= 0. && gv.bin[nd]->chrg[nz]->fac2[i] >= 0. );
03319 }
03320 }
03321
03322 if( ipStart == 0 )
03323 {
03324
03325
03326 UpdateRecomZ0(nd,nz,ALL_STAGES);
03327 }
03328
03329
03330 gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
03331
03332 gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
03333 gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
03334 gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
03335 gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
03336 gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
03337
03338 gv.bin[nd]->chrg[nz]->tedust = 1.f;
03339
03340 gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
03341 gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
03342 gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
03343 gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
03344
03345 gv.bin[nd]->chrg[nz]->BolFlux = -DBL_MAX;
03346 gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
03347 gv.bin[nd]->chrg[nz]->GrainHeatColl = -DBL_MAX;
03348 gv.bin[nd]->chrg[nz]->GasHeatPhotoEl = -DBL_MAX;
03349 gv.bin[nd]->chrg[nz]->GasHeatTherm = -DBL_MAX;
03350 gv.bin[nd]->chrg[nz]->GrainCoolTherm = -DBL_MAX;
03351 gv.bin[nd]->chrg[nz]->ChemEnIon = -DBL_MAX;
03352 gv.bin[nd]->chrg[nz]->ChemEnH2 = -DBL_MAX;
03353
03354 gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
03355
03356
03357 ASSERT( gv.bin[nd]->chrg[nz]->ipThresInf <= gv.bin[nd]->chrg[nz]->ipThresInfVal );
03358 return;
03359 }
03360
03361
03362
03363
03364
03365
03366
03367
03368
03369
03370
03371 STATIC void UpdatePot2(size_t nd,
03372 long nz)
03373 {
03374 double ThermExp;
03375
03376 DEBUG_ENTRY( "UpdatePot2()" );
03377
03378
03379 ThermExp = gv.bin[nd]->chrg[nz]->ThresInf*TE1RYD/gv.bin[nd]->tedust;
03380
03381 gv.bin[nd]->chrg[nz]->ThermRate = THERMCONST*gv.bin[nd]->ThermEff*POW2(gv.bin[nd]->tedust)*exp(-ThermExp);
03382 # if defined( WD_TEST2 ) || defined( IGNORE_THERMIONIC )
03383 gv.bin[nd]->chrg[nz]->ThermRate = 0.;
03384 # endif
03385 return;
03386 }
03387
03388
03389
03390 inline void Yfunc(long nd,
03391 long nz,
03392 double y01,
03393 double maxval,
03394 double Elo,
03395 double Ehi,
03396 double Eel,
03397 double *Yp,
03398 double *Ys,
03399 double *Ehp,
03400 double *Ehs)
03401 {
03402 DEBUG_ENTRY( "Yfunc()" );
03403
03404 long Zg = gv.bin[nd]->chrg[nz]->DustZ;
03405 double y2pr, y2sec;
03406
03407 ASSERT( Ehi >= Elo );
03408
03409 y2pr = y2pa( Elo, Ehi, Zg, Ehp );
03410
03411 if( y2pr > 0. )
03412 {
03413 pe_type pcase = gv.which_pe[gv.bin[nd]->matType];
03414 double eps, f3;
03415
03416 *Yp = y2pr*min(y01,maxval);
03417
03418 y2sec = y2s( Elo, Ehi, Zg, Ehs );
03419 if( pcase == PE_CAR )
03420 eps = 117./EVRYD;
03421 else if( pcase == PE_SIL )
03422 eps = 155./EVRYD;
03423 else
03424 {
03425 fprintf( ioQQQ, " Yfunc: unknown type for PE effect: %d\n" , pcase );
03426 cdEXIT(EXIT_FAILURE);
03427 }
03428
03429
03430 f3 = max(Eel,0.)/(eps*elec_esc_length(Eel,nd)*gv.bin[nd]->eyc);
03431 *Ys = y2sec*f3*min(y01,maxval);
03432 }
03433 else
03434 {
03435 *Yp = 0.;
03436 *Ys = 0.;
03437 *Ehp = 0.;
03438 *Ehs = 0.;
03439 }
03440 return;
03441 }
03442
03443
03444
03445 STATIC double y0b(size_t nd,
03446 long nz,
03447 long i)
03448 {
03449 double yzero;
03450
03451 DEBUG_ENTRY( "y0b()" );
03452
03453 if( gv.lgWD01 )
03454 yzero = y0b01( nd, nz, i );
03455 else
03456 {
03457 double Eph = rfield.anu[i];
03458
03459 if( Eph <= 20./EVRYD )
03460 yzero = y0b01( nd, nz, i );
03461 else if( Eph < 50./EVRYD )
03462 {
03463 double y0a = y0b01( nd, nz, i );
03464 double y0b = gv.bin[nd]->y0b06[i];
03465
03466 double frac = log(Eph*(EVRYD/20.))*1.0913566679372915;
03467
03468 yzero = y0a * exp(log(y0b/y0a)*frac);
03469 }
03470 else
03471 yzero = gv.bin[nd]->y0b06[i];
03472 }
03473
03474 ASSERT( yzero > 0. );
03475 return yzero;
03476 }
03477
03478
03479
03480 STATIC double y0b01(size_t nd,
03481 long nz,
03482 long i)
03483 {
03484 pe_type pcase = gv.which_pe[gv.bin[nd]->matType];
03485 double xv, yzero;
03486
03487 DEBUG_ENTRY( "y0b01()" );
03488
03489 xv = MAX2((rfield.anu[i] - gv.bin[nd]->chrg[nz]->ThresSurfVal)/gv.bin[nd]->DustWorkFcn,0.);
03490
03491 switch( pcase )
03492 {
03493 case PE_CAR:
03494
03495 xv = POW2(xv)*POW3(xv);
03496 yzero = xv/((1./9.e-3) + (3.7e-2/9.e-3)*xv);
03497 break;
03498 case PE_SIL:
03499
03500 yzero = xv/(2.+10.*xv);
03501 break;
03502 default:
03503 fprintf( ioQQQ, " y0b01: unknown type for PE effect: %d\n" , pcase );
03504 cdEXIT(EXIT_FAILURE);
03505 }
03506
03507 ASSERT( yzero > 0. );
03508 return yzero;
03509 }
03510
03511
03512
03513 STATIC double y0psa(size_t nd,
03514 long ns,
03515 long i,
03516 double Eel)
03517 {
03518 double yzero, leola;
03519
03520 DEBUG_ENTRY( "y0psa()" );
03521
03522 ASSERT( i >= gv.bin[nd]->sd[ns]->ipLo );
03523
03524
03525 leola = elec_esc_length(Eel,nd)*gv.bin[nd]->inv_att_len[i];
03526
03527 ASSERT( leola > 0. );
03528
03529
03530 if( leola < 1.e4 )
03531 yzero = gv.bin[nd]->sd[ns]->p[i]*leola*(1. - leola*log(1.+1./leola));
03532 else
03533 {
03534 double x = 1./leola;
03535 yzero = gv.bin[nd]->sd[ns]->p[i]*(((-1./5.*x+1./4.)*x-1./3.)*x+1./2.);
03536 }
03537
03538 ASSERT( yzero > 0. );
03539 return yzero;
03540 }
03541
03542
03543
03544 STATIC double y1psa(size_t nd,
03545 long i,
03546 double Eel)
03547 {
03548 double alpha, beta, af, bf, yone;
03549
03550 DEBUG_ENTRY( "y1psa()" );
03551
03552 beta = gv.bin[nd]->AvRadius*gv.bin[nd]->inv_att_len[i];
03553 if( beta > 1.e-4 )
03554 bf = pow2(beta) - 2.*beta + 2. - 2.*exp(-beta);
03555 else
03556 bf = ((1./60.*beta - 1./12.)*beta + 1./3.)*pow3(beta);
03557
03558 alpha = beta + gv.bin[nd]->AvRadius/elec_esc_length(Eel,nd);
03559 if( alpha > 1.e-4 )
03560 af = pow2(alpha) - 2.*alpha + 2. - 2.*exp(-alpha);
03561 else
03562 af = ((1./60.*alpha - 1./12.)*alpha + 1./3.)*pow3(alpha);
03563
03564 yone = pow2(beta/alpha)*af/bf;
03565
03566 ASSERT( yone > 0. );
03567 return yone;
03568 }
03569
03570
03571
03572 inline double y2pa(double Elo,
03573 double Ehi,
03574 long Zg,
03575 double *Ehp)
03576 {
03577 DEBUG_ENTRY( "y2pa()" );
03578
03579 double ytwo;
03580
03581 if( Zg > -1 )
03582 {
03583 if( Ehi > 0. )
03584 {
03585 double x = Elo/Ehi;
03586 *Ehp = 0.5*Ehi*(1.-2.*x)/(1.-3.*x);
03587
03588 ytwo = ( abs(x) > 1e-4 ) ? (1.-3.*x)/pow3(1.-x) : 1. - (3. + 8.*x)*x*x;
03589 ASSERT( *Ehp > 0. && *Ehp <= Ehi && ytwo > 0. && ytwo <= 1. );
03590 }
03591 else
03592 {
03593 *Ehp = 0.;
03594 ytwo = 0.;
03595 }
03596 }
03597 else
03598 {
03599 if( Ehi > Elo )
03600 {
03601 *Ehp = 0.5*(Elo+Ehi);
03602 ytwo = 1.;
03603 ASSERT( *Ehp >= Elo && *Ehp <= Ehi );
03604 }
03605 else
03606 {
03607 *Ehp = 0.;
03608 ytwo = 0.;
03609 }
03610 }
03611 return ytwo;
03612 }
03613
03614
03615
03616 inline double y2s(double Elo,
03617 double Ehi,
03618 long Zg,
03619 double *Ehs)
03620 {
03621 DEBUG_ENTRY( "y2s()" );
03622
03623 double ytwo;
03624
03625 if( Zg > -1 )
03626 {
03627 if( !gv.lgWD01 && Ehi > 0. )
03628 {
03629 double yl = Elo/ETILDE;
03630 double yh = Ehi/ETILDE;
03631 double x = yh - yl;
03632 double E0, N0;
03633 if( x < 0.01 )
03634 {
03635
03636 double x2 = x*x, x3 = x2*x, x4 = x3*x, x5 = x4*x;
03637 double yh2 = yh*yh, yh3 = yh2*yh, yh4 = yh3*yh, yh5 = yh4*yh;
03638 double help1 = 2.*x-yh;
03639 double help2 = (6.*x3-15.*yh*x2+12.*yh2*x-3.*yh3)/4.;
03640 double help3 = (22.*x5-95.*yh*x4+164.*yh2*x3-141.*yh3*x2+60.*yh4*x-10.*yh5)/16.;
03641 N0 = yh*(help1 - help2 + help3)/x2;
03642
03643 help1 = (3.*x-yh)/3.;
03644 help2 = (15.*x3-25.*yh*x2+15.*yh2*x-3.*yh3)/20.;
03645 help3 = (1155.*x5-3325.*yh*x4+4305.*yh2*x3-2961.*yh3*x2+1050.*yh4*x-150.*yh5)/1680.;
03646 E0 = ETILDE*yh2*(help1 - help2 + help3)/x2;
03647 }
03648 else
03649 {
03650 double sR0 = (1. + yl*yl);
03651 double sqR0 = sqrt(sR0);
03652 double sqRh = sqrt(1. + x*x);
03653 double alpha = sqRh/(sqRh - 1.);
03654 if( yh/sqR0 < 0.01 )
03655 {
03656
03657 double z = yh*(yh - 2.*yl)/sR0;
03658 N0 = ((((7./256.*z-5./128.)*z+1./16.)*z-1./8.)*z+1./2.)*z/(sqRh-1.);
03659
03660 double yl2 = yl*yl, yl3 = yl2*yl, yl4 = yl3*yl;
03661 double help1 = yl/2.;
03662 double help2 = (2.*yl2-1.)/3.;
03663 double help3 = (6.*yl3-9.*yl)/8.;
03664 double help4 = (8.*yl4-24.*yl2+3.)/10.;
03665 double h = yh/sR0;
03666 E0 = -alpha*Ehi*(((help4*h + help3)*h + help2)*h + help1)*h/sqR0;
03667 }
03668 else
03669 {
03670 N0 = alpha*(1./sqR0 - 1./sqRh);
03671 E0 = alpha*ETILDE*(ASINH(x*sqR0 + yl*sqRh) - yh/sqRh);
03672 }
03673 }
03674 ASSERT( N0 > 0. && N0 <= 1. );
03675
03676 *Ehs = E0/N0;
03677
03678 ASSERT( *Ehs > 0. && *Ehs <= Ehi );
03679
03680 ytwo = N0;
03681 }
03682 else
03683 {
03684 *Ehs = 0.;
03685 ytwo = 0.;
03686 }
03687 }
03688 else
03689 {
03690 if( !gv.lgWD01 && Ehi > Elo )
03691 {
03692 double yl = Elo/ETILDE;
03693 double yh = Ehi/ETILDE;
03694 double x = yh - yl;
03695 double x2 = x*x;
03696 if( x > 0.025 )
03697 {
03698 double sqRh = sqrt(1. + x2);
03699 double alpha = sqRh/(sqRh - 1.);
03700 *Ehs = alpha*ETILDE*(ASINH(x) - yh/sqRh + yl);
03701 }
03702 else
03703 {
03704
03705 *Ehs = Ehi - (Ehi-Elo)*((-37./840.*x2 + 1./10.)*x2 + 1./3.);
03706 }
03707
03708 ASSERT( *Ehs >= Elo && *Ehs <= Ehi );
03709
03710 ytwo = 1.;
03711 }
03712 else
03713 {
03714 *Ehs = 0.;
03715 ytwo = 0.;
03716 }
03717 }
03718 return ytwo;
03719 }
03720
03721
03722
03723 STATIC long HighestIonStage(void)
03724 {
03725 long high,
03726 ion,
03727 nelem;
03728
03729 DEBUG_ENTRY( "HighestIonStage()" );
03730
03731 high = 0;
03732 for( nelem=LIMELM-1; nelem >= 0; nelem-- )
03733 {
03734 if( dense.lgElmtOn[nelem] )
03735 {
03736 for( ion=nelem+1; ion >= 0; ion-- )
03737 {
03738 if( ion == high || dense.xIonDense[nelem][ion] > 0. )
03739 break;
03740 }
03741 high = MAX2(high,ion);
03742 }
03743 if( nelem <= high )
03744 break;
03745 }
03746 return high;
03747 }
03748
03749
03750 STATIC void UpdateRecomZ0(size_t nd,
03751 long nz,
03752 bool lgAllIonStages)
03753 {
03754 long hi_ion,
03755 i,
03756 ion,
03757 nelem,
03758 Zg;
03759 double d[5],
03760 phi_s_up[LIMELM+1],
03761 phi_s_dn[2];
03762
03763 DEBUG_ENTRY( "UpdateRecomZ0()" );
03764
03765 Zg = gv.bin[nd]->chrg[nz]->DustZ;
03766
03767 hi_ion = ( lgAllIonStages ) ? LIMELM : gv.HighestIon;
03768
03769 phi_s_up[0] = gv.bin[nd]->chrg[nz]->ThresSurf;
03770 for( i=1; i <= LIMELM; i++ )
03771 {
03772 if( i <= hi_ion )
03773 GetPotValues(nd,Zg+i,&d[0],&d[1],&phi_s_up[i],&d[2],&d[3],&d[4],INCL_TUNNEL);
03774 else
03775 phi_s_up[i] = -DBL_MAX;
03776 }
03777 phi_s_dn[0] = gv.bin[nd]->chrg[nz]->ThresSurfInc;
03778 GetPotValues(nd,Zg-2,&d[0],&d[1],&phi_s_dn[1],&d[2],&d[3],&d[4],NO_TUNNEL);
03779
03780
03781 for( nelem=0; nelem < LIMELM; nelem++ )
03782 {
03783 if( dense.lgElmtOn[nelem] )
03784 {
03785 for( ion=0; ion <= nelem+1; ion++ )
03786 {
03787 if( lgAllIonStages || dense.xIonDense[nelem][ion] > 0. )
03788 {
03789 GrainIonColl(nd,nz,nelem,ion,phi_s_up,phi_s_dn,
03790 &gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion],
03791 &gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion],
03792 &gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion]);
03793 }
03794 else
03795 {
03796 gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] = ion;
03797 gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion] = 0.f;
03798 gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion] = 0.f;
03799 }
03800 }
03801 }
03802 }
03803 return;
03804 }
03805
03806 STATIC void GetPotValues(size_t nd,
03807 long Zg,
03808 double *ThresInf,
03809 double *ThresInfVal,
03810 double *ThresSurf,
03811 double *ThresSurfVal,
03812 double *PotSurf,
03813 double *Emin,
03814 bool lgUseTunnelCorr)
03815 {
03816 double dstpot,
03817 dZg = (double)Zg,
03818 IP_v;
03819
03820 DEBUG_ENTRY( "GetPotValues()" );
03821
03822
03823
03824
03825
03826
03827 dstpot = chrg2pot(dZg,nd);
03828
03829
03830
03831
03832 IP_v = gv.bin[nd]->DustWorkFcn + dstpot - 0.5*one_elec(nd) + (dZg+2.)*AC0/gv.bin[nd]->AvRadius*one_elec(nd);
03833
03834
03835
03836
03837
03838 if( Zg <= -1 )
03839 {
03840 pot_type pcase = gv.which_pot[gv.bin[nd]->matType];
03841 double IP;
03842
03843 IP = gv.bin[nd]->DustWorkFcn - gv.bin[nd]->BandGap + dstpot - 0.5*one_elec(nd);
03844 switch( pcase )
03845 {
03846 case POT_CAR:
03847 IP -= AC1G/(gv.bin[nd]->AvRadius+AC2G)*one_elec(nd);
03848 break;
03849 case POT_SIL:
03850
03851 break;
03852 default:
03853 fprintf( ioQQQ, " GetPotValues detected unknown type for ionization pot: %d\n" , pcase );
03854 cdEXIT(EXIT_FAILURE);
03855 }
03856
03857
03858
03859 IP_v = MAX2(IP,IP_v);
03860
03861 if( Zg < -1 )
03862 {
03863
03864 double help = fabs(dZg+1);
03865
03866 *Emin = -ThetaNu(help)*one_elec(nd);
03867 if( lgUseTunnelCorr )
03868 {
03869
03870 *Emin *= 1. - 2.124e-4/(pow(gv.bin[nd]->AvRadius,(realnum)0.45)*pow(help,0.26));
03871 }
03872 }
03873 else
03874 {
03875 *Emin = 0.;
03876 }
03877
03878 *ThresInf = IP - *Emin;
03879 *ThresInfVal = IP_v - *Emin;
03880 *ThresSurf = *ThresInf;
03881 *ThresSurfVal = *ThresInfVal;
03882 *PotSurf = *Emin;
03883 }
03884 else
03885 {
03886 *ThresInf = IP_v;
03887 *ThresInfVal = IP_v;
03888 *ThresSurf = *ThresInf - dstpot;
03889 *ThresSurfVal = *ThresInfVal - dstpot;
03890 *PotSurf = dstpot;
03891 *Emin = 0.;
03892 }
03893 return;
03894 }
03895
03896
03897
03898
03899
03900 STATIC void GrainIonColl(size_t nd,
03901 long int nz,
03902 long int nelem,
03903 long int ion,
03904 const double phi_s_up[],
03905 const double phi_s_dn[],
03906 long *Z0,
03907 realnum *ChEn,
03908 realnum *ChemEn)
03909 {
03910 long Zg;
03911 double d[5];
03912 double phi_s;
03913
03914 long save = ion;
03915
03916 DEBUG_ENTRY( "GrainIonColl()" );
03917 if( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (realnum)phi_s_up[0] )
03918 {
03919
03920 *ChEn = 0.f;
03921 *ChemEn = 0.f;
03922 Zg = gv.bin[nd]->chrg[nz]->DustZ;
03923 phi_s = phi_s_up[0];
03924 do
03925 {
03926 *ChEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] - (realnum)phi_s;
03927 *ChemEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1];
03928
03929
03930
03931 *ChemEn -= (realnum)(phi_s - phi_s_up[0]);
03932 --ion;
03933 ++Zg;
03934 phi_s = phi_s_up[save-ion];
03935 } while( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (realnum)phi_s );
03936
03937 *Z0 = ion;
03938 }
03939 else if( ion <= nelem && gv.bin[nd]->chrg[nz]->DustZ > gv.bin[nd]->LowestZg &&
03940 rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (realnum)phi_s_dn[0] )
03941 {
03942
03943 *ChEn = 0.f;
03944 *ChemEn = 0.f;
03945 Zg = gv.bin[nd]->chrg[nz]->DustZ;
03946 phi_s = phi_s_dn[0];
03947 do
03948 {
03949 *ChEn += (realnum)phi_s - rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
03950 *ChemEn -= rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
03951
03952
03953
03954 *ChemEn += (realnum)(phi_s - phi_s_dn[0]);
03955 ++ion;
03956 --Zg;
03957
03958 if( ion-save < 2 )
03959 phi_s = phi_s_dn[ion-save];
03960 else
03961 GetPotValues(nd,Zg-1,&d[0],&d[1],&phi_s,&d[2],&d[3],&d[4],NO_TUNNEL);
03962
03963 } while( ion <= nelem && Zg > gv.bin[nd]->LowestZg &&
03964 rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (realnum)phi_s );
03965 *Z0 = ion;
03966 }
03967 else
03968 {
03969
03970 *ChEn = 0.f;
03971 *ChemEn = 0.f;
03972 *Z0 = ion;
03973 }
03974
03975 return;
03976 }
03977
03978
03979
03980 STATIC void GrainChrgTransferRates(long nd)
03981 {
03982 long nz;
03983 double fac0 = STICK_ION*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
03984
03985 DEBUG_ENTRY( "GrainChrgTransferRates()" );
03986
03987 # ifndef IGNORE_GRAIN_ION_COLLISIONS
03988
03989 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
03990 {
03991 long ion;
03992 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
03993 double fac1 = gptr->FracPop*fac0;
03994
03995 if( fac1 == 0. )
03996 continue;
03997
03998 for( ion=0; ion <= LIMELM; ion++ )
03999 {
04000 long nelem;
04001 double eta, fac2, xi;
04002
04003
04004
04005 GrainScreen(ion,nd,nz,&eta,&xi);
04006
04007 fac2 = eta*fac1;
04008
04009 if( fac2 == 0. )
04010 continue;
04011
04012 for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
04013 {
04014 if( dense.lgElmtOn[nelem] && ion != gptr->RecomZ0[nelem][ion] )
04015 {
04016 gv.GrainChTrRate[nelem][ion][gptr->RecomZ0[nelem][ion]] +=
04017 (realnum)(fac2*GetAveVelocity( dense.AtomicWeight[nelem] )*atmdat.lgCTOn);
04018 }
04019 }
04020 }
04021 }
04022 # endif
04023 return;
04024 }
04025
04026
04027
04028
04029 STATIC void GrainUpdateRadius1(void)
04030 {
04031 long nelem;
04032
04033 DEBUG_ENTRY( "GrainUpdateRadius1()" );
04034
04035 for( nelem=0; nelem < LIMELM; nelem++ )
04036 {
04037 gv.elmSumAbund[nelem] = 0.f;
04038 }
04039
04040
04041 for( size_t nd=0; nd < gv.bin.size(); nd++ )
04042 {
04043 gv.bin[nd]->GrnDpth = (realnum)GrnStdDpth(nd);
04044 gv.bin[nd]->dstAbund = (realnum)(gv.bin[nd]->dstfactor*gv.GrainMetal*gv.bin[nd]->GrnDpth);
04045 ASSERT( gv.bin[nd]->dstAbund > 0.f );
04046
04047
04048 gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
04049 gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
04050
04051 gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
04052 gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;
04053
04054
04055
04056
04057 for( nelem=0; nelem < LIMELM; nelem++ )
04058 {
04059 gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(realnum)gv.bin[nd]->cnv_H_pCM3;
04060 }
04061 }
04062 return;
04063 }
04064
04065
04066
04067
04068 STATIC void GrainUpdateRadius2()
04069 {
04070 DEBUG_ENTRY( "GrainUpdateRadius2()" );
04071
04072 for( long i=0; i < rfield.nupper; i++ )
04073 {
04074 gv.dstab[i] = 0.;
04075 gv.dstsc[i] = 0.;
04076 }
04077
04078
04079
04080 for( size_t nd=0; nd < gv.bin.size(); nd++ )
04081 {
04082 realnum dstAbund = gv.bin[nd]->dstAbund;
04083
04084
04085 for( long i=0; i < rfield.nflux; i++ )
04086 {
04087
04088
04089
04090
04091
04092
04093 gv.dstab[i] += gv.bin[nd]->dstab1[i]*dstAbund;
04094 gv.dstsc[i] += gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]*dstAbund;
04095 }
04096
04097 for( long nz=0; nz < gv.bin[nd]->nChrg; nz++ )
04098 {
04099 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
04100 if( gptr->DustZ <= -1 )
04101 {
04102 double FracPop = gptr->FracPop;
04103
04104 for( long i=gptr->ipThresInf; i < rfield.nflux; i++ )
04105 gv.dstab[i] += FracPop*gptr->cs_pdt[i]*dstAbund;
04106 }
04107 }
04108 }
04109
04110 for( long i=0; i < rfield.nflux; i++ )
04111 {
04112
04113 ASSERT( gv.dstab[i] > 0. && gv.dstsc[i] > 0. );
04114 }
04115 return;
04116 }
04117
04118
04119
04120 STATIC void GrainTemperature(size_t nd,
04121 realnum *dccool,
04122 double *hcon,
04123 double *hots,
04124 double *hla)
04125 {
04126 long int i,
04127 ipLya,
04128 nz;
04129 double EhatThermionic,
04130 norm,
04131 rate,
04132 x,
04133 y;
04134 realnum dcheat;
04135
04136 DEBUG_ENTRY( "GrainTemperature()" );
04137
04138
04139 ASSERT( nd < gv.bin.size() );
04140
04141 if( trace.lgTrace && trace.lgDustBug )
04142 {
04143 fprintf( ioQQQ, " GrainTemperature starts for grain %s\n", gv.bin[nd]->chDstLab );
04144 }
04145
04146
04147
04148
04149
04150 *hcon = 0.;
04151
04152 *hots = 0.;
04153
04154 *hla = 0.;
04155
04156 ipLya = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont() - 1;
04157
04158
04159
04160 gv.bin[nd]->GasHeatPhotoEl = 0.;
04161
04162 gv.bin[nd]->GrainCoolTherm = 0.;
04163 gv.bin[nd]->thermionic = 0.;
04164
04165 dcheat = 0.f;
04166 *dccool = 0.f;
04167
04168 gv.bin[nd]->BolFlux = 0.;
04169
04170
04171
04172 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
04173 {
04174 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
04175
04176 double hcon1 = 0.;
04177 double hots1 = 0.;
04178 double hla1 = 0.;
04179 double bolflux1 = 0.;
04180 double pe1 = 0.;
04181
04182
04183 bool lgReEvaluate1 = gptr->hcon1 < 0.;
04184 bool lgReEvaluate2 = gptr->hots1 < 0.;
04185 bool lgReEvaluate = lgReEvaluate1 || lgReEvaluate2;
04186
04187
04188 if( lgReEvaluate )
04189 {
04190 long loopmax = MIN2(gptr->ipThresInf,rfield.nflux);
04191 for( i=0; i < loopmax; i++ )
04192 {
04193 double fac = gv.bin[nd]->dstab1[i]*rfield.anu[i];
04194
04195 if( lgReEvaluate1 )
04196 hcon1 += rfield.flux[0][i]*fac;
04197
04198 if( lgReEvaluate2 )
04199 {
04200 hots1 += rfield.SummedDif[i]*fac;
04201 # ifndef NDEBUG
04202 bolflux1 += rfield.SummedCon[i]*fac;
04203 # endif
04204 }
04205 }
04206 }
04207
04208
04209
04210
04211
04212
04213
04214 if( lgReEvaluate1 )
04215 {
04216 for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
04217 {
04218 hcon1 += rfield.flux[0][i]*gptr->fac1[i];
04219 }
04220
04221 gptr->hcon1 = hcon1;
04222 }
04223 else
04224 {
04225 hcon1 = gptr->hcon1;
04226 }
04227
04228 if( lgReEvaluate2 )
04229 {
04230 for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
04231 {
04232
04233
04234 hots1 += rfield.SummedDif[i]*gptr->fac1[i];
04235
04236 #ifdef WD_TEST2
04237 pe1 += rfield.flux[0][i]*gptr->fac2[i];
04238 #else
04239 pe1 += rfield.SummedCon[i]*gptr->fac2[i];
04240 #endif
04241 # ifndef NDEBUG
04242 bolflux1 += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*rfield.anu[i];
04243 if( gptr->DustZ <= -1 )
04244 bolflux1 += rfield.SummedCon[i]*gptr->cs_pdt[i]*rfield.anu[i];
04245 # endif
04246 }
04247 gptr->hots1 = hots1;
04248 gptr->bolflux1 = bolflux1;
04249 gptr->pe1 = pe1;
04250 }
04251 else
04252 {
04253 hots1 = gptr->hots1;
04254 bolflux1 = gptr->bolflux1;
04255 pe1 = gptr->pe1;
04256 }
04257
04258
04259
04260
04261
04262 if( ipLya < MIN2(gptr->ipThresInf,rfield.nflux) )
04263 {
04264 hla1 = rfield.otslin[ipLya]*gv.bin[nd]->dstab1[ipLya]*0.75;
04265 }
04266 else if( ipLya < rfield.nflux )
04267 {
04268
04269 hla1 = rfield.otslin[ipLya]*gptr->fac1[ipLya];
04270 }
04271 else
04272 {
04273 hla1 = 0.;
04274 }
04275
04276 ASSERT( hcon1 >= 0. && hots1 >= 0. && hla1 >= 0. && bolflux1 >= 0. && pe1 >= 0. );
04277
04278 *hcon += gptr->FracPop*hcon1;
04279 *hots += gptr->FracPop*hots1;
04280 *hla += gptr->FracPop*hla1;
04281 gv.bin[nd]->BolFlux += gptr->FracPop*bolflux1;
04282 if( gv.lgDHetOn )
04283 gv.bin[nd]->GasHeatPhotoEl += gptr->FracPop*pe1;
04284
04285 # ifndef NDEBUG
04286 if( trace.lgTrace && trace.lgDustBug )
04287 {
04288 fprintf( ioQQQ, " Zg %ld bolflux: %.4e\n", gptr->DustZ,
04289 gptr->FracPop*bolflux1*EN1RYD*gv.bin[nd]->cnv_H_pCM3 );
04290 }
04291 # endif
04292
04293
04294
04295
04296
04297
04298
04299 rate = gptr->FracPop*gptr->ThermRate*gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
04300
04301 EhatThermionic = 2.*BOLTZMANN*gv.bin[nd]->tedust + MAX2(gptr->PotSurf*EN1RYD,0.);
04302 gv.bin[nd]->GrainCoolTherm += rate * (EhatThermionic + gptr->ThresSurf*EN1RYD);
04303 gv.bin[nd]->thermionic += rate * (EhatThermionic - gptr->PotSurf*EN1RYD);
04304 }
04305
04306
04307 norm = EN1RYD*gv.bin[nd]->cnv_H_pCM3;
04308
04309
04310 *hcon *= norm;
04311
04312
04313 *hots *= norm;
04314
04315
04316 *hla *= norm;
04317
04318 gv.bin[nd]->BolFlux *= norm;
04319
04320
04321
04322
04323
04324
04325
04326 GrainCollHeating(nd,&dcheat,dccool);
04327
04328
04329 gv.bin[nd]->GasHeatPhotoEl *= norm;
04330
04331 if( gv.lgBakesPAH_heat )
04332 {
04333
04334
04335
04336
04337
04338 gv.bin[nd]->GasHeatPhotoEl = 1.e-24*hmi.UV_Cont_rel2_Habing_TH85_depth*
04339 dense.gas_phase[ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((hmi.UV_Cont_rel2_Habing_TH85_depth*
04340
04341 phycon.sqrte/dense.eden),0.73)) + 3.65e-2*pow(phycon.te/1.e4,0.7)/
04342 (1.+2.e-4*(hmi.UV_Cont_rel2_Habing_TH85_depth*phycon.sqrte/dense.eden)))/gv.bin.size();
04343
04344 }
04345
04346
04347
04348 gv.bin[nd]->GasHeatPhotoEl *= gv.GrainHeatScaleFactor;
04349
04350
04351
04352
04353
04354
04355
04356
04357
04358
04359
04360
04361
04362
04363
04364
04365 gv.bin[nd]->GrainHeat = *hcon + *hots + dcheat - gv.bin[nd]->GrainCoolTherm;
04366
04367
04368 gv.bin[nd]->GrainHeatColl = dcheat;
04369
04370
04371
04372
04373
04374 if( gv.bin[nd]->GrainHeat > 0. )
04375 {
04376 bool lgOutOfBounds;
04377
04378
04379 x = log(MAX2(DBL_MIN,gv.bin[nd]->GrainHeat*gv.bin[nd]->cnv_CM3_pH));
04380
04381 splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,x,&y,&lgOutOfBounds);
04382 gv.bin[nd]->tedust = (realnum)exp(y);
04383 }
04384 else
04385 {
04386 gv.bin[nd]->GrainHeat = -1.;
04387 gv.bin[nd]->tedust = -1.;
04388 }
04389
04390 if( thermal.ConstGrainTemp > 0. )
04391 {
04392 bool lgOutOfBounds;
04393
04394 gv.bin[nd]->tedust = thermal.ConstGrainTemp;
04395
04396 x = log(gv.bin[nd]->tedust);
04397 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgOutOfBounds);
04398 gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;
04399 }
04400
04401
04402 gv.bin[nd]->TeGrainMax = (realnum)MAX2(gv.bin[nd]->TeGrainMax,gv.bin[nd]->tedust);
04403
04404 if( trace.lgTrace && trace.lgDustBug )
04405 {
04406 fprintf( ioQQQ, " >GrainTemperature finds %s Tdst %.5e hcon %.4e ",
04407 gv.bin[nd]->chDstLab, gv.bin[nd]->tedust, *hcon);
04408 fprintf( ioQQQ, "hots %.4e dcheat %.4e GrainCoolTherm %.4e\n",
04409 *hots, dcheat, gv.bin[nd]->GrainCoolTherm );
04410 }
04411 return;
04412 }
04413
04414
04415
04416 STATIC void PE_init(size_t nd,
04417 long nz,
04418 long i,
04419 double *cs1,
04420 double *cs2,
04421 double *cs_tot,
04422 double *cool1,
04423 double *cool2,
04424 double *ehat1,
04425 double *ehat2)
04426 {
04427 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
04428 long ipLo1 = gptr->ipThresInfVal;
04429 long ipLo2 = gptr->ipThresInf;
04430
04431 DEBUG_ENTRY( "PE_init()" );
04432
04433
04434 ASSERT( nd < gv.bin.size() );
04435 ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
04436 ASSERT( i >= 0 && i < rfield.nflux );
04437
04440
04441 if( i >= ipLo1 )
04442 {
04443
04444 *cs1 = gv.bin[nd]->dstab1[i]*gptr->yhat[i];
04445
04446
04447 *ehat1 = gptr->ehat[i];
04448
04449
04450
04451
04452
04453
04454
04455
04456
04457
04458
04459
04460
04461
04462
04463
04464
04465
04466
04467
04468 if( gptr->DustZ <= -1 )
04469 *cool1 = gptr->ThresSurf + gptr->PotSurf + *ehat1;
04470 else
04471 *cool1 = gptr->ThresSurfVal + gptr->PotSurf + *ehat1;
04472
04473 ASSERT( *ehat1 > 0. && *cool1 > 0. );
04474 }
04475 else
04476 {
04477 *cs1 = 0.;
04478 *ehat1 = 0.;
04479 *cool1 = 0.;
04480 }
04481
04482
04483 if( gptr->DustZ <= -1 && i >= ipLo2 )
04484 {
04485
04486 *cs2 = gptr->cs_pdt[i];
04487
04488 *ehat2 = rfield.anu[i] - gptr->ThresSurf - gptr->PotSurf;
04489
04490 *cool2 = rfield.anu[i];
04491
04492 ASSERT( *ehat2 > 0. && *cool2 > 0. );
04493 }
04494 else
04495 {
04496 *cs2 = 0.;
04497 *ehat2 = 0.;
04498 *cool2 = 0.;
04499 }
04500
04501 *cs_tot = gv.bin[nd]->dstab1[i] + *cs2;
04502 return;
04503 }
04504
04505
04506
04507 STATIC void GrainCollHeating(size_t nd,
04508 realnum *dcheat,
04509 realnum *dccool)
04510 {
04511 long int ion,
04512 nelem,
04513 nz;
04514 H2_type ipH2;
04515 double Accommodation,
04516 CollisionRateElectr,
04517 CollisionRateMol,
04518 CollisionRateIon,
04519 CoolTot,
04520 CoolBounce,
04521 CoolEmitted,
04522 CoolElectrons,
04523 CoolMolecules,
04524 CoolPotential,
04525 CoolPotentialGas,
04526 eta,
04527 HeatTot,
04528 HeatBounce,
04529 HeatCollisions,
04530 HeatElectrons,
04531 HeatIons,
04532 HeatMolecules,
04533 HeatRecombination,
04534 HeatChem,
04535 HeatCor,
04536 Stick,
04537 ve,
04538 WeightMol,
04539 xi;
04540
04541
04542
04543 const double H2_FORMATION_GRAIN_HEATING[H2_TOP] = { 0.20, 0.4, 1.72 };
04544
04545 DEBUG_ENTRY( "GrainCollHeating()" );
04546
04547
04548
04549
04550
04551
04552
04553
04554
04555
04556
04557
04558 HeatTot = 0.;
04559 CoolTot = 0.;
04560
04561 HeatIons = 0.;
04562
04563 gv.bin[nd]->ChemEn = 0.;
04564
04565
04566 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
04567 {
04568 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
04569
04570
04571
04572 double Heat1 = 0.;
04573 double Cool1 = 0.;
04574 double ChemEn1 = 0.;
04575
04576
04577
04578
04579
04580 for( ion=0; ion <= LIMELM; ion++ )
04581 {
04582
04583
04584
04585
04586
04587
04588
04589 CollisionRateIon = 0.;
04590 CoolPotential = 0.;
04591 CoolPotentialGas = 0.;
04592 HeatRecombination = 0.;
04593 HeatChem = 0.;
04594
04595
04596
04597 GrainScreen(ion,nd,nz,&eta,&xi);
04598
04599 for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
04600 {
04601 if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. )
04602 {
04603 double CollisionRateOne;
04604
04605
04606
04607
04608 #if defined( IGNORE_GRAIN_ION_COLLISIONS )
04609 Stick = 0.;
04610 #elif defined( WD_TEST2 )
04611 Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
04612 0. : STICK_ION;
04613 #else
04614 Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
04615 gv.bin[nd]->AccomCoef[nelem] : STICK_ION;
04616 #endif
04617
04618
04619
04620 CollisionRateOne = Stick*dense.xIonDense[nelem][ion]*GetAveVelocity( dense.AtomicWeight[nelem] );
04621 CollisionRateIon += CollisionRateOne;
04622
04623
04624
04625
04626
04627
04628
04629
04630 if( ion >= gptr->RecomZ0[nelem][ion] )
04631 {
04632 CoolPotential += CollisionRateOne * (double)ion *
04633 gptr->PotSurf;
04634 CoolPotentialGas += CollisionRateOne *
04635 (double)gptr->RecomZ0[nelem][ion] *
04636 gptr->PotSurf;
04637 }
04638 else
04639 {
04640 CoolPotential += CollisionRateOne * (double)ion *
04641 gptr->PotSurfInc;
04642 CoolPotentialGas += CollisionRateOne *
04643 (double)gptr->RecomZ0[nelem][ion] *
04644 gptr->PotSurfInc;
04645 }
04646
04647
04648
04649
04650 HeatRecombination += CollisionRateOne *
04651 gptr->RecomEn[nelem][ion];
04652 HeatChem += CollisionRateOne * gptr->ChemEn[nelem][ion];
04653 }
04654 }
04655
04656
04657
04658
04659
04660
04661 HeatCollisions = CollisionRateIon * 2.*BOLTZMANN*phycon.te*xi;
04662
04663
04664
04665
04666 CoolPotential *= eta*EN1RYD;
04667 CoolPotentialGas *= eta*EN1RYD;
04668
04669 HeatRecombination *= eta*EN1RYD;
04670 HeatChem *= eta*EN1RYD;
04671
04672 CoolEmitted = CollisionRateIon * 2.*BOLTZMANN*gv.bin[nd]->tedust*eta;
04673
04674
04675 Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;
04676
04677
04678
04679
04680 Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;
04681
04682 ChemEn1 += HeatChem;
04683 }
04684
04685
04686 HeatIons += gptr->FracPop*Heat1;
04687
04688 if( trace.lgTrace && trace.lgDustBug )
04689 {
04690 fprintf( ioQQQ, " Zg %ld ions heat/cool: %.4e %.4e\n", gptr->DustZ,
04691 gptr->FracPop*Heat1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04692 gptr->FracPop*Cool1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
04693 }
04694
04695
04696
04697
04698 ion = -1;
04699 Stick = ( gptr->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
04700
04701
04702 ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);
04703
04704
04705 CollisionRateElectr = Stick*dense.eden*ve;
04706
04707
04708
04709 GrainScreen(ion,nd,nz,&eta,&xi);
04710
04711 if( gptr->DustZ > gv.bin[nd]->LowestZg )
04712 {
04713 HeatCollisions = CollisionRateElectr*2.*BOLTZMANN*phycon.te*xi;
04714
04715
04716 CoolPotential = CollisionRateElectr * (double)ion*gptr->PotSurfInc*eta*EN1RYD;
04717
04718 HeatRecombination = CollisionRateElectr * gptr->ThresSurfInc*eta*EN1RYD;
04719 HeatBounce = 0.;
04720 CoolBounce = 0.;
04721 }
04722 else
04723 {
04724 HeatCollisions = 0.;
04725 CoolPotential = 0.;
04726 HeatRecombination = 0.;
04727
04728
04729
04730
04731 HeatBounce = CollisionRateElectr * 2.*BOLTZMANN*phycon.te*xi;
04732
04733
04734 CoolBounce = CollisionRateElectr *
04735 (-gptr->ThresSurfInc-gptr->PotSurfInc)*EN1RYD*eta;
04736 CoolBounce = MAX2(CoolBounce,0.);
04737 }
04738
04739
04740
04741 HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
04742 Heat1 += HeatElectrons;
04743
04744 CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
04745 Cool1 += CoolElectrons;
04746
04747 if( trace.lgTrace && trace.lgDustBug )
04748 {
04749 fprintf( ioQQQ, " Zg %ld electrons heat/cool: %.4e %.4e\n", gptr->DustZ,
04750 gptr->FracPop*HeatElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04751 gptr->FracPop*CoolElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
04752 }
04753
04754
04755
04756
04757
04758
04759
04760
04761
04762
04763
04764
04765 gptr->HeatingRate2 = HeatElectrons*gv.bin[nd]->IntArea/4. -
04766 gv.bin[nd]->GrainCoolTherm*gv.bin[nd]->cnv_CM3_pH;
04767
04768
04769
04770
04771
04772 HeatTot += gptr->FracPop*Heat1;
04773
04774
04775 CoolTot += gptr->FracPop*Cool1;
04776
04777 gv.bin[nd]->ChemEn += gptr->FracPop*ChemEn1;
04778 }
04779
04780
04781
04782
04783
04784
04785
04786
04787 WeightMol = 2.*dense.AtomicWeight[ipHYDROGEN];
04788 Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
04789
04790 #ifndef IGNORE_GRAIN_ION_COLLISIONS
04791
04792 CollisionRateMol = Accommodation*hmi.H2_total*
04793 sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
04794
04795
04796 ipH2 = gv.which_H2distr[gv.bin[nd]->matType];
04797
04798
04799 gv.bin[nd]->ChemEnH2 = gv.bin[nd]->rate_h2_form_grains_used*dense.xIonDense[ipHYDROGEN][0]*
04800 H2_FORMATION_GRAIN_HEATING[ipH2]*EN1EV;
04801
04802 gv.bin[nd]->ChemEnH2 /= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
04803 #else
04804 CollisionRateMol = 0.;
04805 gv.bin[nd]->ChemEnH2 = 0.;
04806 #endif
04807
04808
04809 WeightMol = dense.AtomicWeight[ipCARBON] + dense.AtomicWeight[ipOXYGEN];
04810 Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
04811 #ifndef IGNORE_GRAIN_ION_COLLISIONS
04812 CollisionRateMol += Accommodation*findspecieslocal("CO")->den*
04813 sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
04814 #else
04815 CollisionRateMol = 0.;
04816 #endif
04817
04818
04819 HeatCollisions = CollisionRateMol * 2.*BOLTZMANN*phycon.te;
04820 CoolEmitted = CollisionRateMol * 2.*BOLTZMANN*gv.bin[nd]->tedust;
04821
04822 HeatMolecules = HeatCollisions - CoolEmitted + gv.bin[nd]->ChemEnH2;
04823 HeatTot += HeatMolecules;
04824
04825
04826 CoolMolecules = HeatCollisions - CoolEmitted;
04827 CoolTot += CoolMolecules;
04828
04829 gv.bin[nd]->RateUp = 0.;
04830 gv.bin[nd]->RateDn = 0.;
04831 HeatCor = 0.;
04832 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
04833 {
04834 double d[4];
04835 double rate_dn = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
04836 double rate_up = GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
04837
04838 gv.bin[nd]->RateUp += gv.bin[nd]->chrg[nz]->FracPop*rate_up;
04839 gv.bin[nd]->RateDn += gv.bin[nd]->chrg[nz]->FracPop*rate_dn;
04840
04842 HeatCor += (gv.bin[nd]->chrg[nz]->FracPop*rate_up*gv.bin[nd]->chrg[nz]->ThresSurf -
04843 gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->ThresSurfInc +
04844 gv.bin[nd]->chrg[nz]->FracPop*rate_up*gv.bin[nd]->chrg[nz]->PotSurf -
04845 gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->PotSurfInc)*EN1RYD;
04846 }
04847
04848
04849 HeatTot += HeatCor;
04850
04851 if( trace.lgTrace && trace.lgDustBug )
04852 {
04853 fprintf( ioQQQ, " molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
04854 HeatMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04855 CoolMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04856 HeatCor*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
04857 }
04858
04859 *dcheat = (realnum)(HeatTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3);
04860 *dccool = ( gv.lgDColOn ) ? (realnum)(CoolTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3) : 0.f;
04861
04862 gv.bin[nd]->ChemEn *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
04863 gv.bin[nd]->ChemEnH2 *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
04864
04865
04866
04867
04868
04869
04870
04871
04872
04873
04874
04875
04876
04877
04878
04879
04880 gv.bin[nd]->HeatingRate1 = (HeatMolecules+HeatIons+HeatCor)*gv.bin[nd]->IntArea/4.;
04881
04882
04883 return;
04884 }
04885
04886
04887
04888 void GrainDrift(void)
04889 {
04890 long int i,
04891 loop;
04892 double alam,
04893 corr,
04894 dmomen,
04895 fac,
04896 fdrag,
04897 g0,
04898 g2,
04899 phi2lm,
04900 psi,
04901 rdust,
04902 si,
04903 vdold,
04904 volmom;
04905
04906 DEBUG_ENTRY( "GrainDrift()" );
04907
04908 vector<realnum> help( rfield.nflux );
04909 for( i=0; i < rfield.nflux; i++ )
04910 {
04911 help[i] = (rfield.flux[0][i]+rfield.ConInterOut[i]+rfield.outlin[0][i]+rfield.outlin_noplot[i])*
04912 rfield.anu[i];
04913 }
04914
04915 for( size_t nd=0; nd < gv.bin.size(); nd++ )
04916 {
04917
04918 dmomen = 0.;
04919 for( i=0; i < rfield.nflux; i++ )
04920 {
04921
04922 dmomen += help[i]*(gv.bin[nd]->dstab1[i] + gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]);
04923 }
04924 ASSERT( dmomen >= 0. );
04925 dmomen *= EN1RYD*4./gv.bin[nd]->IntArea;
04926
04927
04928 fac = 2*BOLTZMANN*phycon.te;
04929
04930
04931
04932 psi = gv.bin[nd]->dstpot*TE1RYD/phycon.te;
04933 if( psi > 0. )
04934 {
04935 rdust = 1.e-6;
04936 alam = log(20.702/rdust/psi*phycon.sqrte/dense.SqrtEden);
04937 }
04938 else
04939 {
04940 alam = 0.;
04941 }
04942
04943 phi2lm = POW2(psi)*alam;
04944 corr = 2.;
04945
04946 for( loop = 0; loop < 50 && fabs(corr-1.) > 0.001; loop++ )
04947 {
04948 vdold = gv.bin[nd]->DustDftVel;
04949
04950
04951 si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
04952 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04953 g2 = si/(1.329 + POW3(si));
04954
04955
04956
04957 fdrag = fac*dense.xIonDense[ipHYDROGEN][1]*(g0 + phi2lm*g2);
04958
04959
04960 si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.816e-6;
04961 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04962 g2 = si/(1.329 + POW3(si));
04963 fdrag += fac*dense.eden*(g0 + phi2lm*g2);
04964
04965
04966 si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
04967 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04968 fdrag += fac*(dense.xIonDense[ipHYDROGEN][0] + 1.1*dense.xIonDense[ipHELIUM][0])*g0;
04969
04970
04971 si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.551e-4;
04972 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04973 g2 = si/(1.329 + POW3(si));
04974 fdrag += fac*dense.xIonDense[ipHELIUM][1]*(g0 + phi2lm*g2);
04975
04976
04977
04978
04979 volmom = dmomen/SPEEDLIGHT;
04980
04981 if( fdrag > 0. )
04982 {
04983 corr = sqrt(volmom/fdrag);
04984 gv.bin[nd]->DustDftVel = (realnum)(vdold*corr);
04985 }
04986 else
04987 {
04988 corr = 1.;
04989 gv.lgNegGrnDrg = true;
04990 gv.bin[nd]->DustDftVel = 0.;
04991 }
04992
04993 if( trace.lgTrace && trace.lgDustBug )
04994 {
04995 fprintf( ioQQQ, " %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n",
04996 loop, gv.bin[nd]->DustDftVel, volmom );
04997 }
04998 }
04999 }
05000 return;
05001 }
05002
05003
05004
05005
05006 STATIC double GrnVryDpth(
05007
05008
05009
05010
05011
05012 size_t nd)
05013 {
05014 DEBUG_ENTRY( "GrnVryDpth()" );
05015
05016 ASSERT( nd < gv.bin.size() );
05017
05018
05019
05020
05021
05022
05023
05024
05025
05026
05027 double GrnVryDpth_v = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
05028
05029
05030
05031
05032
05033 return min(max(1.e-10,GrnVryDpth_v),1.);
05034 }