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