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