cloudy  trunk
grains.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*grain main routine to converge grains thermal solution */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "atmdat.h"
7 #include "rfield.h"
8 #include "hmi.h"
9 #include "trace.h"
10 #include "conv.h"
11 #include "ionbal.h"
12 #include "thermal.h"
13 #include "phycon.h"
14 #include "doppvel.h"
15 #include "taulines.h"
16 #include "mole.h"
17 #include "heavy.h"
18 #include "thirdparty.h"
19 #include "dense.h"
20 #include "ipoint.h"
21 #include "elementnames.h"
22 #include "grainvar.h"
23 #include "grains.h"
24 #include "iso.h"
25 
26 /* the next three defines are for debugging purposes only, uncomment to activate */
27 /* #define WD_TEST2 1 */
28 /* #define IGNORE_GRAIN_ION_COLLISIONS 1 */
29 /* #define IGNORE_THERMIONIC 1 */
30 
34 inline double ASINH(double x)
35 {
36  if( abs(x) <= 8.e-3 )
37  return ((0.075*pow2(x) - 1./6.)*pow2(x) + 1.)*x;
38  else if( abs(x) <= 1./sqrt(DBL_EPSILON) )
39  {
40  if( x < 0. )
41  return -log(sqrt(1. + pow2(x)) - x);
42  else
43  return log(sqrt(1. + pow2(x)) + x);
44  }
45  else
46  {
47  if( x < 0. )
48  return -(log(-x)+LN_TWO);
49  else
50  return log(x)+LN_TWO;
51  }
52 }
53 
54 /* no parentheses around PTR needed since it needs to be an lvalue */
55 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
56 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
57 
58 static const long MAGIC_AUGER_DATA = 20060126L;
59 
60 static const bool INCL_TUNNEL = true;
61 static const bool NO_TUNNEL = false;
62 
63 static const bool ALL_STAGES = true;
64 /* static const bool NONZERO_STAGES = false; */
65 
66 /* counts how many times GrainDrive has been called, set to zero in GrainZero */
67 static long int nCalledGrainDrive;
68 
69 /*================================================================================*/
70 /* these are used for setting up grain emissivities in InitEmissivities() */
71 
72 /* NTOP is number of bins for temps between GRAIN_TMID and GRAIN_TMAX */
73 static const long NTOP = NDEMS/5;
74 
75 /*================================================================================*/
76 /* these are used when iterating the grain charge in GrainCharge() */
77 static const double TOLER = CONSERV_TOL/10.;
78 static const long BRACKET_MAX = 50L;
79 
80 /* >>chng 06 feb 07, increased CT_LOOP_MAX (10 -> 25), T_LOOP_MAX (30 -> 50), pah.in, PvH */
81 
82 /* maximum number of tries to converge charge/temperature in GrainChargeTemp() */
83 static const long CT_LOOP_MAX = 25L;
84 
85 /* maximum number of tries to converge grain temperature in GrainChargeTemp() */
86 static const long T_LOOP_MAX = 50L;
87 
88 /* these will become the new tolerance levels used throughout the code */
89 static double HEAT_TOLER = DBL_MAX;
90 static double HEAT_TOLER_BIN = DBL_MAX;
91 static double CHRG_TOLER = DBL_MAX;
92 /* static double CHRG_TOLER_BIN = DBL_MAX; */
93 
94 /*================================================================================*/
95 /* miscellaneous grain physics */
96 
97 /* a_0 thru a_2 constants for calculating IP_V and EA, in cm */
98 static const double AC0 = 3.e-9;
99 static const double AC1G = 4.e-8;
100 static const double AC2G = 7.e-8;
101 
102 /* constants needed to calculate energy distribution of secondary electrons */
103 static const double ETILDE = 2.*SQRT2/EVRYD; /* sqrt(8) eV */
104 
105 /* constant for thermionic emissions, 7.501e20 e/cm^2/s/K^2 */
106 static const double THERMCONST = PI4*ELECTRON_MASS*POW2(BOLTZMANN)/POW3(HPLANCK);
107 
108 /* sticking probabilities */
109 static const double STICK_ELEC = 0.5;
110 static const double STICK_ION = 1.0;
111 
113 inline double one_elec(long nd)
114 {
115  return ELEM_CHARGE/EVRYD/gv.bin[nd]->Capacity;
116 }
117 
119 inline double pot2chrg(double x,
120  long nd)
121 {
122  return x/one_elec(nd) - 1.;
123 }
124 
126 inline double chrg2pot(double x,
127  long nd)
128 {
129  return (x+1.)*one_elec(nd);
130 }
131 
133 inline double elec_esc_length(double e, // energy of electron in Ryd
134  long nd)
135 {
136  // calculate escape length in cm
137  if( e <= gv.bin[nd]->le_thres )
138  return 1.e-7;
139  else
140  return 3.e-6*gv.bin[nd]->eec*sqrt(pow3(e*EVRYD*1.e-3));
141 }
142 
143 /* read data for electron energy spectrum of Auger electrons */
144 STATIC void ReadAugerData();
145 /* initialize the Auger data for grain bin nd between index ipBegin <= i < ipEnd */
146 STATIC void InitBinAugerData(size_t,long,long);
147 /* read a single line of data from data file */
148 STATIC void GetNextLine(const char*, FILE*, char[]);
149 /* initialize grain emissivities */
150 STATIC void InitEmissivities(void);
151 /* PlanckIntegral compute total radiative cooling due to large grains */
152 STATIC double PlanckIntegral(double,size_t,long);
153 /* invalidate charge dependent data from previous iteration */
154 STATIC void NewChargeData(long);
155 /* GrnStdDpth returns the grain abundance as a function of depth into cloud */
156 STATIC double GrnStdDpth(long);
157 /* iterate grain charge and temperature */
158 STATIC void GrainChargeTemp(void);
159 /* GrainCharge compute grains charge */
160 STATIC void GrainCharge(size_t,/*@out@*/double*);
161 /* grain electron recombination rates for single charge state */
162 STATIC double GrainElecRecomb1(size_t,long,/*@out@*/double*,/*@out@*/double*);
163 /* grain electron emission rates for single charge state */
164 STATIC double GrainElecEmis1(size_t,long,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*);
165 /* correction factors for grain charge screening (including image potential
166  * to correct for polarization of the grain as charged particle approaches). */
167 STATIC void GrainScreen(long,size_t,long,double*,double*);
168 /* helper function for GrainScreen */
169 STATIC double ThetaNu(double);
170 /* update items that depend on grain potential */
171 STATIC void UpdatePot(size_t,long,long,/*@out@*/double[],/*@out@*/double[]);
172 /* calculate charge state populations */
173 STATIC void GetFracPop(size_t,long,/*@in@*/double[],/*@in@*/double[],/*@out@*/long*);
174 /* this routine updates all quantities that depend only on grain charge and radius */
175 STATIC void UpdatePot1(size_t,long,long,long);
176 /* this routine updates all quantities that depend on grain charge, radius and temperature */
177 STATIC void UpdatePot2(size_t,long);
178 /* Helper function to calculate primary and secondary yields and the average electron energy at infinity */
179 inline void Yfunc(long,long,double,double,double,double,double,/*@out@*/double*,/*@out@*/double*,
180  /*@out@*/double*,/*@out@*/double*);
181 /* This calculates the y0 function for band electrons (Sect. 4.1.3/4.1.4 of WDB06) */
182 STATIC double y0b(size_t,long,long);
183 /* This calculates the y0 function for band electrons (Eq. 16 of WD01) */
184 STATIC double y0b01(size_t,long,long);
185 /* This calculates the y0 function for primary/secondary and Auger electrons (Eq. 9 of WDB06) */
186 STATIC double y0psa(size_t,long,long,double);
187 /* This calculates the y1 function for primary/secondary and Auger electrons (Eq. 6 of WDB06) */
188 STATIC double y1psa(size_t,long,double);
189 /* This calculates the y2 function for primary and Auger electrons (Eq. 8 of WDB06) */
190 inline double y2pa(double,double,long,double*);
191 /* This calculates the y2 function for secondary electrons (Eq. 20-21 of WDB06) */
192 inline double y2s(double,double,long,double*);
193 /* find highest ionization stage with non-zero population */
194 STATIC long HighestIonStage(void);
195 /* determine charge Z0 ion recombines to upon impact on grain */
196 STATIC void UpdateRecomZ0(size_t,long,bool);
197 /* helper routine for UpdatePot */
198 STATIC void GetPotValues(size_t,long,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
199  /*@out@*/double*,/*@out@*/double*,/*@out@*/double*,bool);
200 /* given grain nd in charge state nz, and incoming ion (nelem,ion),
201  * detemine outgoing ion (nelem,Z0) and chemical energy ChEn released
202  * ChemEn is net contribution of ion recombination to grain heating */
203 STATIC void GrainIonColl(size_t,long,long,long,const double[],const double[],/*@out@*/long*,
204  /*@out@*/realnum*,/*@out@*/realnum*);
205 /* initialize ion recombination rates on grain species nd */
206 STATIC void GrainChrgTransferRates(long);
207 /* this routine updates all grain quantities that depend on radius, except gv.dstab and gv.dstsc */
208 STATIC void GrainUpdateRadius1(void);
209 /* this routine adds all the grain opacities in gv.dstab and gv.dstsc */
210 STATIC void GrainUpdateRadius2();
211 /* GrainTemperature computes grains temperature, and gas cooling */
212 STATIC void GrainTemperature(size_t,/*@out@*/realnum*,/*@out@*/double*,/*@out@*/double*,
213  /*@out@*/double*);
214 /* helper routine for initializing quantities related to the photo-electric effect */
215 STATIC void PE_init(size_t,long,long,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
216  /*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*);
217 /* GrainCollHeating computes grains collisional heating cooling */
218 STATIC void GrainCollHeating(size_t,/*@out@*/realnum*,/*@out@*/realnum*);
219 /* GrnVryDpth user supplied function for the grain abundance as a function of depth into cloud */
220 STATIC double GrnVryDpth(size_t);
221 
222 
223 void AEInfo::p_clear0()
224 {
225  nData.clear();
226  IonThres.clear();
227  AvNumber.clear();
228  Energy.clear();
229 }
230 
231 void AEInfo::p_clear1()
232 {
233  nSubShell = 0;
234 }
235 
236 void ShellData::p_clear0()
237 {
238  p.clear();
239  y01.clear();
240  AvNr.clear();
241  Ener.clear();
242  y01A.clear();
243 }
244 
245 void ShellData::p_clear1()
246 {
247  nelem = LONG_MIN;
248  ns = LONG_MIN;
249  ionPot = -DBL_MAX;
250  ipLo = LONG_MIN;
251  nData = 0;
252 }
253 
254 void ChargeBin::p_clear0()
255 {
256  yhat.clear();
257  yhat_primary.clear();
258  ehat.clear();
259  cs_pdt.clear();
260  fac1.clear();
261  fac2.clear();
262 }
263 
264 void ChargeBin::p_clear1()
265 {
266  DustZ = LONG_MIN;
267  nfill = 0;
268  FracPop = -DBL_MAX;
269  tedust = 1.f;
270 }
271 
272 void GrainBin::p_clear0()
273 {
274  dstab1.clear();
275  pure_sc1.clear();
276  asym.clear();
277  y0b06.clear();
278  inv_att_len.clear();
279 
280  for( unsigned int ns=0; ns < sd.size(); ns++ )
281  delete sd[ns];
282  sd.clear();
283 
284  for( int nz=0; nz < NCHS; nz++ )
285  {
286  delete chrg[nz];
287  chrg[nz] = NULL;
288  }
289 }
290 
291 void GrainBin::p_clear1()
292 {
293  nDustFunc = DF_STANDARD;
294  lgPAHsInIonizedRegion = false;
295  avDGRatio = 0.;
296  dstfactor = 1.f;
297  dstAbund = -FLT_MAX;
298  GrnDpth = 1.f;
299  cnv_H_pGR = -DBL_MAX;
300  cnv_H_pCM3 = -DBL_MAX;
301  cnv_CM3_pGR = -DBL_MAX;
302  cnv_CM3_pH = -DBL_MAX;
303  cnv_GR_pH = -DBL_MAX;
304  cnv_GR_pCM3 = -DBL_MAX;
305  /* used to check that the energy grid resolution scale factor in
306  * grains opacity files is the same as current cloudy scale */
307  RSFCheck = 0.;
308  memset( dstems, 0, NDEMS*sizeof(double) );
309  memset( dstslp, 0, NDEMS*sizeof(double) );
310  memset( dstslp2, 0, NDEMS*sizeof(double) );
311  lgTdustConverged = false;
312  /* >>chng 00 jun 19, tedust has to be greater than zero
313  * to prevent division by zero in GrainElecEmis and GrainCollHeating, PvH */
314  tedust = 1.f;
315  TeGrainMax = FLT_MAX;
316  avdust = 0.;
317  lgChrgConverged = false;
318  LowestZg = LONG_MIN;
319  nfill = 0;
320  sd.reserve(15);
321  AveDustZ = -DBL_MAX;
322  dstpot = -DBL_MAX;
323  dstpotsav = -DBL_MAX;
324  LowestPot = -DBL_MAX;
325  RateUp = -DBL_MAX;
326  RateDn = -DBL_MAX;
327  StickElecNeg = -DBL_MAX;
328  StickElecPos = -DBL_MAX;
329  avdpot = 0.;
330  le_thres = FLT_MAX;
331  BolFlux = -DBL_MAX;
332  GrainCoolTherm = -DBL_MAX;
333  GasHeatPhotoEl = -DBL_MAX;
334  GrainHeat = DBL_MAX/10.;
335  GrainHeatColl = -DBL_MAX;
336  GrainGasCool = DBL_MAX/10.;
337  ChemEn = -DBL_MAX;
338  ChemEnH2 = -DBL_MAX;
339  thermionic = -DBL_MAX;
340  lgQHeat = false;
341  lgUseQHeat = false;
342  lgEverQHeat = false;
343  lgQHTooWide = false;
344  QHeatFailures = 0;
345  qnflux = LONG_MAX;
346  qnflux2 = LONG_MAX;
347  qtmin = -DBL_MAX;
348  qtmin_zone1 = -DBL_MAX;
349  HeatingRate1 = -DBL_MAX;
350  memset( DustEnth, 0, NDEMS*sizeof(double) );
351  memset( EnthSlp, 0, NDEMS*sizeof(double) );
352  memset( EnthSlp2, 0, NDEMS*sizeof(double) );
353  rate_h2_form_grains_HM79 = 0.;
354  rate_h2_form_grains_CT02 = 0.;
355  /* >>chng 04 feb 05, zero this rate in case "no molecules" is set, will.in, PvH */
356  rate_h2_form_grains_used = 0.;
357  DustDftVel = 1.e3f;
358  avdft = 0.;
359  /* NB - this number should not be larger than NCHU */
360  nChrgOrg = gv.nChrgRequested;
361  nChrg = nChrgOrg;
362  for( int nz=0; nz < NCHS; nz++ )
363  chrg[nz] = NULL;
364 }
365 
366 void GrainVar::p_clear0()
367 {
368  for( size_t nd=0; nd < bin.size(); nd++ )
369  delete bin[nd];
370  bin.clear();
371 
372  for( int nelem=0; nelem < LIMELM; nelem++ )
373  {
374  delete AugerData[nelem];
375  AugerData[nelem] = NULL;
376  }
377 
378  ReadRecord.clear();
379  anumin.clear();
380  anumax.clear();
381  dstab.clear();
382  dstsc.clear();
383  GrainEmission.clear();
384  GraphiteEmission.clear();
385  SilicateEmission.clear();
386 }
387 
388 void GrainVar::p_clear1()
389 {
390  bin.reserve(50);
391 
392  for( int nelem=0; nelem < LIMELM; nelem++ )
393  AugerData[nelem] = NULL;
394 
395  lgAnyDustVary = false;
396  TotalEden = 0.;
397  dHeatdT = 0.;
398  lgQHeatAll = false;
399  /* lgGrainElectrons - should grain electron source/sink be included in overall electron sum?
400  * default is true, set false with no grain electrons command */
401  lgGrainElectrons = true;
402  lgQHeatOn = true;
403  lgDHetOn = true;
404  lgDColOn = true;
405  GrainMetal = 1.;
406  dstAbundThresholdNear = 1.e-6f;
407  dstAbundThresholdFar = 1.e-3f;
408  lgWD01 = false;
409  nChrgRequested = NCHRG_DEFAULT;
410  /* by default grains always reevaluated - command grains reevaluate off sets to false */
411  lgReevaluate = true;
412  /* flag saying neg grain drift vel found */
413  lgNegGrnDrg = false;
414 
415  /* counts how many times GrainDrive has been called */
416  nCalledGrainDrive = 0;
417 
418  /* this is sest true with "set PAH Bakes" command - must also turn off
419  * grain heating with "grains no heat" to only get their results */
420  lgBakesPAH_heat = false;
421 
422  /* this is option to turn off all grain physics while leaving
423  * the opacity in, set false with no grain physics command */
424  lgGrainPhysicsOn = true;
425 
426  /* scale factor set with SET GRAINS HEAT command to rescale grain photoelectric
427  * heating as per Allers et al. 2005 */
428  GrainHeatScaleFactor = 1.f;
429 
430  /* the following entries define the physical behavior of each type of grains
431  * (entropy function, expression for Zmin and ionization potential, etc) */
432  which_enth[MAT_CAR] = ENTH_CAR;
433  which_zmin[MAT_CAR] = ZMIN_CAR;
434  which_pot[MAT_CAR] = POT_CAR;
435  which_ial[MAT_CAR] = IAL_CAR;
436  which_pe[MAT_CAR] = PE_CAR;
437  which_strg[MAT_CAR] = STRG_CAR;
438  which_H2distr[MAT_CAR] = H2_CAR;
439 
440  which_enth[MAT_SIL] = ENTH_SIL;
441  which_zmin[MAT_SIL] = ZMIN_SIL;
442  which_pot[MAT_SIL] = POT_SIL;
443  which_ial[MAT_SIL] = IAL_SIL;
444  which_pe[MAT_SIL] = PE_SIL;
445  which_strg[MAT_SIL] = STRG_SIL;
446  which_H2distr[MAT_SIL] = H2_SIL;
447 
448  which_enth[MAT_PAH] = ENTH_PAH;
449  which_zmin[MAT_PAH] = ZMIN_CAR;
450  which_pot[MAT_PAH] = POT_CAR;
451  which_ial[MAT_PAH] = IAL_CAR;
452  which_pe[MAT_PAH] = PE_CAR;
453  which_strg[MAT_PAH] = STRG_CAR;
454  which_H2distr[MAT_PAH] = H2_CAR;
455 
456  which_enth[MAT_CAR2] = ENTH_CAR2;
457  which_zmin[MAT_CAR2] = ZMIN_CAR;
458  which_pot[MAT_CAR2] = POT_CAR;
459  which_ial[MAT_CAR2] = IAL_CAR;
460  which_pe[MAT_CAR2] = PE_CAR;
461  which_strg[MAT_CAR2] = STRG_CAR;
462  which_H2distr[MAT_CAR2] = H2_CAR;
463 
464  which_enth[MAT_SIL2] = ENTH_SIL2;
465  which_zmin[MAT_SIL2] = ZMIN_SIL;
466  which_pot[MAT_SIL2] = POT_SIL;
467  which_ial[MAT_SIL2] = IAL_SIL;
468  which_pe[MAT_SIL2] = PE_SIL;
469  which_strg[MAT_SIL2] = STRG_SIL;
470  which_H2distr[MAT_SIL2] = H2_SIL;
471 
472  which_enth[MAT_PAH2] = ENTH_PAH2;
473  which_zmin[MAT_PAH2] = ZMIN_CAR;
474  which_pot[MAT_PAH2] = POT_CAR;
475  which_ial[MAT_PAH2] = IAL_CAR;
476  which_pe[MAT_PAH2] = PE_CAR;
477  which_strg[MAT_PAH2] = STRG_CAR;
478  which_H2distr[MAT_PAH2] = H2_CAR;
479 
480  for( int nelem=0; nelem < LIMELM; nelem++ )
481  {
482  for( int ion=0; ion <= nelem+1; ion++ )
483  {
484  for( int ion_to=0; ion_to <= nelem+1; ion_to++ )
485  {
486  GrainChTrRate[nelem][ion][ion_to] = 0.f;
487  }
488  }
489  }
490 
491  /* this sets the default abundance dependence for PAHs,
492  * proportional to n(H0) / n(Htot)
493  * changed with SET PAH command */
494  chPAH_abundance = "H";
495 }
496 
497 
498 /* this routine is called by zero(), so it should contain initializations
499  * that need to be done every time before the input lines are parsed */
500 void GrainZero(void)
501 {
502  DEBUG_ENTRY( "GrainZero()" );
503 
504  /* >>>chng 01 may 08, return memory possibly allocated in previous calls to cloudy(), PvH
505  * this routine MUST be called before ParseCommands() so that grain commands find a clean slate */
506  gv.clear();
507  return;
508 }
509 
510 
511 /* this routine is called by IterStart(), so anything that needs to be reset before each
512  * iteration starts should be put here; typically variables that are integrated over radius */
513 void GrainStartIter(void)
514 {
515  DEBUG_ENTRY( "GrainStartIter()" );
516 
517  if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
518  {
519  gv.lgNegGrnDrg = false;
520  gv.TotalDustHeat = 0.;
521  gv.GrnElecDonateMax = 0.;
522  gv.GrnElecHoldMax = 0.;
523  gv.dphmax = 0.f;
524  gv.dclmax = 0.f;
525 
526  for( size_t nd=0; nd < gv.bin.size(); nd++ )
527  {
528  /* >>chng 97 jul 5, save and reset this
529  * save grain potential */
530  gv.bin[nd]->dstpotsav = gv.bin[nd]->dstpot;
531  gv.bin[nd]->qtmin = ( gv.bin[nd]->qtmin_zone1 > 0. ) ?
532  gv.bin[nd]->qtmin_zone1 : DBL_MAX;
533  gv.bin[nd]->avdust = 0.;
534  gv.bin[nd]->avdpot = 0.;
535  gv.bin[nd]->avdft = 0.;
536  gv.bin[nd]->avDGRatio = 0.;
537  gv.bin[nd]->TeGrainMax = -1.f;
538  gv.bin[nd]->lgEverQHeat = false;
539  gv.bin[nd]->QHeatFailures = 0L;
540  gv.bin[nd]->lgQHTooWide = false;
541  gv.bin[nd]->lgPAHsInIonizedRegion = false;
542  gv.bin[nd]->nChrgOrg = gv.bin[nd]->nChrg;
543  }
544  }
545  return;
546 }
547 
548 
549 /* this routine is called by IterRestart(), so anything that needs to be
550  * reset or saved after an iteration is finished should be put here */
551 void GrainRestartIter(void)
552 {
553  DEBUG_ENTRY( "GrainRestartIter()" );
554 
555  if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
556  {
557  for( size_t nd=0; nd < gv.bin.size(); nd++ )
558  {
559  /* >>chng 97 jul 5, reset grain potential
560  * reset grain to pential to initial value from previous iteration */
561  gv.bin[nd]->dstpot = gv.bin[nd]->dstpotsav;
562  gv.bin[nd]->nChrg = gv.bin[nd]->nChrgOrg;
563  }
564  }
565  return;
566 }
567 
568 
569 /* this routine is called by ParseSet() */
570 void SetNChrgStates(long nChrg)
571 {
572  DEBUG_ENTRY( "SetNChrgStates()" );
573 
574  ASSERT( nChrg >= 2 && nChrg <= NCHU );
575  gv.nChrgRequested = nChrg;
576  return;
577 }
578 
579 
580 /*GrainsInit, called one time by opacitycreateall at initialization of calculation,
581  * called after commands have been parsed,
582  * not after every iteration or every model */
583 void GrainsInit(void)
584 {
585  long int i,
586  nelem;
587  unsigned int ns;
588 
589  DEBUG_ENTRY( "GrainsInit()" );
590 
591  if( trace.lgTrace && trace.lgDustBug )
592  {
593  fprintf( ioQQQ, " GrainsInit called.\n" );
594  }
595 
596  gv.anumin.resize( rfield.nupper );
597  gv.anumax.resize( rfield.nupper );
598  gv.dstab.resize( rfield.nupper );
599  gv.dstsc.resize( rfield.nupper );
600  gv.GrainEmission.resize( rfield.nupper );
601  gv.GraphiteEmission.resize( rfield.nupper );
602  gv.SilicateEmission.resize( rfield.nupper );
603 
604  /* >>chng 02 jan 15, initialize to zero in case grains are not used, needed in IonIron(), PvH */
605  for( nelem=0; nelem < LIMELM; nelem++ )
606  {
607  gv.elmSumAbund[nelem] = 0.f;
608  }
609 
610  for( i=0; i < rfield.nupper; i++ )
611  {
612  gv.dstab[i] = 0.;
613  gv.dstsc[i] = 0.;
614  /* >>chng 01 sep 12, moved next three initializations from GrainZero(), PvH */
615  gv.GrainEmission[i] = 0.;
616  gv.SilicateEmission[i] = 0.;
617  gv.GraphiteEmission[i] = 0.;
618  }
619 
620  if( !gv.lgDustOn() )
621  {
622  /* grains are not on, set all heating/cooling agents to zero */
623  gv.GrainHeatInc = 0.;
624  gv.GrainHeatDif = 0.;
625  gv.GrainHeatLya = 0.;
626  gv.GrainHeatCollSum = 0.;
627  gv.GrainHeatSum = 0.;
628  gv.GasCoolColl = 0.;
629  thermal.heating[0][13] = 0.;
630  thermal.heating[0][14] = 0.;
631  thermal.heating[0][25] = 0.;
632 
633  if( trace.lgTrace && trace.lgDustBug )
634  {
635  fprintf( ioQQQ, " GrainsInit exits.\n" );
636  }
637  return;
638  }
639 
640 #ifdef WD_TEST2
641  gv.lgWD01 = true;
642 #endif
643 
644  HEAT_TOLER = conv.HeatCoolRelErrorAllowed / 3.;
645  HEAT_TOLER_BIN = HEAT_TOLER / sqrt((double)gv.bin.size());
646  CHRG_TOLER = conv.EdenErrorAllowed / 3.;
647  /* CHRG_TOLER_BIN = CHRG_TOLER / sqrt(gv.bin.size()); */
648 
649  gv.anumin[0] = 0.f;
650  for( i=1; i < rfield.nupper; i++ )
651  gv.anumax[i-1] = gv.anumin[i] = (realnum)sqrt(rfield.anu[i-1]*rfield.anu[i]);
652  gv.anumax[rfield.nupper-1] = FLT_MAX;
653 
654  ReadAugerData();
655 
656  for( size_t nd=0; nd < gv.bin.size(); nd++ )
657  {
658  double help,atoms,p_rad,ThresInf,ThresInfVal,Emin,d[5];
659  long low1,low2,low3,lowm;
660 
661  /* sanity checks */
662  ASSERT( gv.bin[nd] != NULL );
663  ASSERT( gv.bin[nd]->nChrg >= 2 && gv.bin[nd]->nChrg <= NCHU );
664 
665  if( gv.bin[nd]->DustWorkFcn < rfield.anu[0] || gv.bin[nd]->DustWorkFcn > rfield.anu[rfield.nupper] )
666  {
667  fprintf( ioQQQ, " Grain work function for %s has insane value: %.4e\n",
668  gv.bin[nd]->chDstLab,gv.bin[nd]->DustWorkFcn );
669  cdEXIT(EXIT_FAILURE);
670  }
671 
672  /* this is QHEAT ALL command */
673  if( gv.lgQHeatAll )
674  {
675  gv.bin[nd]->lgQHeat = true;
676  }
677 
678  /* this is NO GRAIN QHEAT command, always takes precedence */
679  if( !gv.lgQHeatOn )
680  {
681  gv.bin[nd]->lgQHeat = false;
682  }
683 
684  /* >>chng 04 jun 01, disable quantum heating when constant grain temperature is used, PvH */
685  if( thermal.ConstGrainTemp > 0. )
686  {
687  gv.bin[nd]->lgQHeat = false;
688  }
689 
690 #ifndef IGNORE_QUANTUM_HEATING
691  gv.bin[nd]->lgQHTooWide = false;
692  gv.bin[nd]->qtmin = DBL_MAX;
693 #endif
694 
695  if( gv.bin[nd]->nDustFunc>DF_STANDARD || gv.bin[nd]->matType == MAT_PAH ||
696  gv.bin[nd]->matType == MAT_PAH2 )
697  gv.lgAnyDustVary = true;
698 
699  /* grain abundance may depend on radius,
700  * invalidate for now; GrainUpdateRadius1() will set correct value */
701  gv.bin[nd]->dstAbund = -FLT_MAX;
702 
703  gv.bin[nd]->GrnDpth = 1.f;
704 
705  gv.bin[nd]->qtmin_zone1 = -1.;
706 
707  /* this is threshold in Ryd above which to use X-ray prescription for electron escape length */
708  gv.bin[nd]->le_thres = gv.lgWD01 ? FLT_MAX :
709  (realnum)(pow(pow((double)gv.bin[nd]->dustp[0],0.85)/30.,2./3.)*1.e3/EVRYD);
710 
711  for( long nz=0; nz < NCHS; nz++ )
712  {
713  ASSERT( gv.bin[nd]->chrg[nz] == NULL );
714  gv.bin[nd]->chrg[nz] = new ChargeBin;
715  }
716 
717  /* >>chng 00 jun 19, this value is absolute lower limit for the grain
718  * potential, electrons cannot be bound for lower values..., PvH */
719  zmin_type zcase = gv.which_zmin[gv.bin[nd]->matType];
720  switch( zcase )
721  {
722  case ZMIN_CAR:
723  // this is Eq. 23a + 24 of WD01
724  help = gv.bin[nd]->AvRadius*1.e7;
725  help = ceil(-(1.2*POW2(help)+3.9*help+0.2)/1.44);
726  break;
727  case ZMIN_SIL:
728  // this is Eq. 23b + 24 of WD01
729  help = gv.bin[nd]->AvRadius*1.e7;
730  help = ceil(-(0.7*POW2(help)+2.5*help+0.8)/1.44);
731  break;
732  default:
733  fprintf( ioQQQ, " GrainsInit detected unknown Zmin type: %d\n" , zcase );
734  cdEXIT(EXIT_FAILURE);
735  }
736 
737  /* this is to assure that gv.bin[nd]->LowestZg > LONG_MIN */
738  ASSERT( help > (double)(LONG_MIN+1) );
739  low1 = nint(help);
740 
741  /* >>chng 01 apr 20, iterate to get LowestPot such that the exponent in the thermionic
742  * rate never becomes positive; the value can be derived by equating ThresInf >= 0;
743  * the new expression for Emin (see GetPotValues) cannot be inverted analytically,
744  * hence it is necessary to iterate for LowestPot. this also automatically assures that
745  * the expressions for ThresInf and LowestPot are consistent with each other, PvH */
746  low2 = low1;
747  GetPotValues(nd,low2,&ThresInf,&d[0],&d[1],&d[2],&d[3],&d[4],INCL_TUNNEL);
748  if( ThresInf < 0. )
749  {
750  low3 = 0;
751  /* do a bisection search for the lowest charge such that
752  * ThresInf >= 0, the end result will eventually be in low3 */
753  while( low3-low2 > 1 )
754  {
755  lowm = (low2+low3)/2;
756  GetPotValues(nd,lowm,&ThresInf,&d[0],&d[1],&d[2],&d[3],&d[4],INCL_TUNNEL);
757  if( ThresInf < 0. )
758  low2 = lowm;
759  else
760  low3 = lowm;
761  }
762  low2 = low3;
763  }
764 
765  /* the first term implements the minimum charge due to autoionization
766  * the second term assures that the exponent in the thermionic rate never
767  * becomes positive; the expression was derived by equating ThresInf >= 0 */
768  gv.bin[nd]->LowestZg = MAX2(low1,low2);
769  gv.bin[nd]->LowestPot = chrg2pot(gv.bin[nd]->LowestZg,nd);
770 
771  ns = 0;
772 
773  ASSERT( gv.bin[nd]->sd.size() == 0 );
774  gv.bin[nd]->sd.push_back( new ShellData );
775 
776  /* this is data for valence band */
777  gv.bin[nd]->sd[ns]->nelem = -1;
778  gv.bin[nd]->sd[ns]->ns = -1;
779  gv.bin[nd]->sd[ns]->ionPot = gv.bin[nd]->DustWorkFcn;
780 
781  /* now add data for inner shell photoionization */
782  for( nelem=ipLITHIUM; nelem < LIMELM && !gv.lgWD01; nelem++ )
783  {
784  if( gv.bin[nd]->elmAbund[nelem] > 0. )
785  {
786  if( gv.AugerData[nelem] == NULL )
787  {
788  fprintf( ioQQQ, " Grain Auger data are missing for element %s\n",
789  elementnames.chElementName[nelem] );
790  fprintf( ioQQQ, " Please include the NO GRAIN X-RAY TREATMENT command "
791  "to disable the Auger treatment in grains.\n" );
792  cdEXIT(EXIT_FAILURE);
793  }
794 
795  for( unsigned int j=0; j < gv.AugerData[nelem]->nSubShell; j++ )
796  {
797  ++ns;
798 
799  gv.bin[nd]->sd.push_back( new ShellData );
800 
801  gv.bin[nd]->sd[ns]->nelem = nelem;
802  gv.bin[nd]->sd[ns]->ns = j;
803  gv.bin[nd]->sd[ns]->ionPot = gv.AugerData[nelem]->IonThres[j];
804  }
805  }
806  }
807 
808  GetPotValues(nd,gv.bin[nd]->LowestZg,&d[0],&ThresInfVal,&d[1],&d[2],&d[3],&Emin,INCL_TUNNEL);
809 
810  for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
811  {
812  long ipLo;
813  double Ethres = ( ns == 0 ) ? ThresInfVal : gv.bin[nd]->sd[ns]->ionPot;
814  ShellData *sptr = gv.bin[nd]->sd[ns];
815 
816  sptr->ipLo = hunt_bisect( rfield.anu, rfield.nupper, Ethres ) + 1;
817 
818  ipLo = sptr->ipLo;
819  // allow 10 elements room for adjustment of rfield.nflux later on
820  // if the adjustment is larger, flex_arr will copy the store, so no problem
821  long len = rfield.nflux + 10 - ipLo;
822 
823  sptr->p.reserve( len );
824  sptr->p.alloc( ipLo, rfield.nflux );
825 
826  sptr->y01.reserve( len );
827  sptr->y01.alloc( ipLo, rfield.nflux );
828 
829  /* there are no Auger electrons from the band structure */
830  if( ns > 0 )
831  {
832  sptr->nData = gv.AugerData[sptr->nelem]->nData[sptr->ns];
833  sptr->AvNr.resize( sptr->nData );
834  sptr->Ener.resize( sptr->nData );
835  sptr->y01A.resize( sptr->nData );
836 
837  for( long n=0; n < sptr->nData; n++ )
838  {
839  sptr->AvNr[n] = gv.AugerData[sptr->nelem]->AvNumber[sptr->ns][n];
840  sptr->Ener[n] = gv.AugerData[sptr->nelem]->Energy[sptr->ns][n];
841 
842  sptr->y01A[n].reserve( len );
843  sptr->y01A[n].alloc( ipLo, rfield.nflux );
844  }
845  }
846  }
847 
848  gv.bin[nd]->y0b06.resize( rfield.nupper );
849 
850  InitBinAugerData( nd, 0, rfield.nflux );
851 
852  gv.bin[nd]->nfill = rfield.nflux;
853 
854  /* >>chng 00 jul 13, new sticking probability for electrons */
855  /* the second term is chance that electron passes through grain,
856  * 1-p_rad is chance that electron is ejected before grain settles
857  * see discussion in
858  * >>refer grain physics Weingartner & Draine, 2001, ApJS, 134, 263 */
860  gv.bin[nd]->StickElecPos = STICK_ELEC*(1. - exp(-gv.bin[nd]->AvRadius/elec_esc_length(0.,nd)));
861  atoms = gv.bin[nd]->AvVol*gv.bin[nd]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[nd]->atomWeight;
862  p_rad = 1./(1.+exp(20.-atoms));
863  gv.bin[nd]->StickElecNeg = gv.bin[nd]->StickElecPos*p_rad;
864 
865  /* >>chng 02 feb 15, these quantities depend on radius and are normally set
866  * in GrainUpdateRadius1(), however, it is necessary to initialize them here
867  * as well so that they are valid the first time hmole is called. */
868  gv.bin[nd]->GrnDpth = (realnum)GrnStdDpth(nd);
869  gv.bin[nd]->dstAbund = (realnum)(gv.bin[nd]->dstfactor*gv.GrainMetal*gv.bin[nd]->GrnDpth);
870  ASSERT( gv.bin[nd]->dstAbund > 0.f );
871  /* grain unit conversion, <unit>/H (default depl) -> <unit>/cm^3 (actual depl) */
872  gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
873  gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
874  /* grain unit conversion, <unit>/cm^3 (actual depl) -> <unit>/grain */
875  gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
876  gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;
877  }
878 
879  /* >>chng 02 dec 19, these quantities depend on radius and are normally set
880  * in GrainUpdateRadius1(), however, it is necessary to initialize them here
881  * as well so that they are valid for the initial printout in Cloudy, PvH */
882  /* calculate the summed grain abundances, these are valid at the inner radius;
883  * these numbers depend on radius and are updated in GrainUpdateRadius1() */
884  for( nelem=0; nelem < LIMELM; nelem++ )
885  {
886  gv.elmSumAbund[nelem] = 0.f;
887  }
888 
889  for( size_t nd=0; nd < gv.bin.size(); nd++ )
890  {
891  for( nelem=0; nelem < LIMELM; nelem++ )
892  {
893  gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(realnum)gv.bin[nd]->cnv_H_pCM3;
894  }
895  }
896 
897  gv.nzone = -1;
898  gv.GrnRecomTe = -1.;
899 
900  /* >>chng 01 nov 21, total grain opacities depend on charge and therefore on radius,
901  * invalidate for now; GrainUpdateRadius2() will set correct values */
902  for( i=0; i < rfield.nupper; i++ )
903  {
904  /* these are total absorption and scattering cross sections,
905  * the latter should contain the asymmetry factor (1-g) */
906  gv.dstab[i] = -DBL_MAX;
907  gv.dstsc[i] = -DBL_MAX;
908  }
909 
910  InitEmissivities();
911  InitEnthalpy();
912 
913  if( trace.lgDustBug && trace.lgTrace )
914  {
915  fprintf( ioQQQ, " There are %ld grain types turned on.\n", (unsigned long)gv.bin.size() );
916 
917  fprintf( ioQQQ, " grain depletion factors, dstfactor*GrainMetal=" );
918  for( size_t nd=0; nd < gv.bin.size(); nd++ )
919  {
920  fprintf( ioQQQ, "%10.2e", gv.bin[nd]->dstfactor*gv.GrainMetal );
921  }
922  fprintf( ioQQQ, "\n" );
923 
924  fprintf( ioQQQ, " nChrg =" );
925  for( size_t nd=0; nd < gv.bin.size(); nd++ )
926  {
927  fprintf( ioQQQ, " %ld", gv.bin[nd]->nChrg );
928  }
929  fprintf( ioQQQ, "\n" );
930 
931  fprintf( ioQQQ, " lowest charge (e) =" );
932  for( size_t nd=0; nd < gv.bin.size(); nd++ )
933  {
934  fprintf( ioQQQ, "%10.2e", pot2chrg(gv.bin[nd]->LowestPot,nd) );
935  }
936  fprintf( ioQQQ, "\n" );
937 
938  fprintf( ioQQQ, " nDustFunc flag for user requested custom depth dependence:" );
939  for( size_t nd=0; nd < gv.bin.size(); nd++ )
940  {
941  fprintf( ioQQQ, "%2i", gv.bin[nd]->nDustFunc );
942  }
943  fprintf( ioQQQ, "\n" );
944 
945  fprintf( ioQQQ, " Quantum heating flag:" );
946  for( size_t nd=0; nd < gv.bin.size(); nd++ )
947  {
948  fprintf( ioQQQ, "%2c", TorF(gv.bin[nd]->lgQHeat) );
949  }
950  fprintf( ioQQQ, "\n" );
951 
952  /* >>chng 01 nov 21, removed total abs and sct cross sections, they are invalid */
953  fprintf( ioQQQ, " NU(Ryd), Abs cross sec per proton\n" );
954 
955  fprintf( ioQQQ, " Ryd " );
956  for( size_t nd=0; nd < gv.bin.size(); nd++ )
957  {
958  fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
959  }
960  fprintf( ioQQQ, "\n" );
961 
962  for( i=0; i < rfield.nupper; i += 40 )
963  {
964  fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
965  for( size_t nd=0; nd < gv.bin.size(); nd++ )
966  {
967  fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->dstab1[i] );
968  }
969  fprintf( ioQQQ, "\n" );
970  }
971 
972  fprintf( ioQQQ, " NU(Ryd), Sct cross sec per proton\n" );
973 
974  fprintf( ioQQQ, " Ryd " );
975  for( size_t nd=0; nd < gv.bin.size(); nd++ )
976  {
977  fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
978  }
979  fprintf( ioQQQ, "\n" );
980 
981  for( i=0; i < rfield.nupper; i += 40 )
982  {
983  fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
984  for( size_t nd=0; nd < gv.bin.size(); nd++ )
985  {
986  fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->pure_sc1[i] );
987  }
988  fprintf( ioQQQ, "\n" );
989  }
990 
991  fprintf( ioQQQ, " NU(Ryd), Q abs\n" );
992 
993  fprintf( ioQQQ, " Ryd " );
994  for( size_t nd=0; nd < gv.bin.size(); nd++ )
995  {
996  fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
997  }
998  fprintf( ioQQQ, "\n" );
999 
1000  for( i=0; i < rfield.nupper; i += 40 )
1001  {
1002  fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
1003  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1004  {
1005  fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->dstab1[i]*4./gv.bin[nd]->IntArea );
1006  }
1007  fprintf( ioQQQ, "\n" );
1008  }
1009 
1010  fprintf( ioQQQ, " NU(Ryd), Q sct\n" );
1011 
1012  fprintf( ioQQQ, " Ryd " );
1013  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1014  {
1015  fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
1016  }
1017  fprintf( ioQQQ, "\n" );
1018 
1019  for( i=0; i < rfield.nupper; i += 40 )
1020  {
1021  fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
1022  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1023  {
1024  fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->pure_sc1[i]*4./gv.bin[nd]->IntArea );
1025  }
1026  fprintf( ioQQQ, "\n" );
1027  }
1028 
1029  fprintf( ioQQQ, " NU(Ryd), asymmetry factor\n" );
1030 
1031  fprintf( ioQQQ, " Ryd " );
1032  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1033  {
1034  fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
1035  }
1036  fprintf( ioQQQ, "\n" );
1037 
1038  for( i=0; i < rfield.nupper; i += 40 )
1039  {
1040  fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
1041  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1042  {
1043  fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->asym[i] );
1044  }
1045  fprintf( ioQQQ, "\n" );
1046  }
1047 
1048  fprintf( ioQQQ, " GrainsInit exits.\n" );
1049  }
1050  return;
1051 }
1052 
1053 /* read data for electron energy spectrum of Auger electrons */
1054 STATIC void ReadAugerData()
1055 {
1056  char chString[FILENAME_PATH_LENGTH_2];
1057  long version;
1058  FILE *fdes;
1059 
1060  DEBUG_ENTRY( "ReadAugerData()" );
1061 
1062  static const char chFile[] = "auger_spectrum.dat";
1063  fdes = open_data( chFile, "r" );
1064 
1065  GetNextLine( chFile, fdes, chString );
1066  sscanf( chString, "%ld", &version );
1067  if( version != MAGIC_AUGER_DATA )
1068  {
1069  fprintf( ioQQQ, " File %s has wrong version number\n", chFile );
1070  fprintf( ioQQQ, " please update your installation...\n" );
1071  cdEXIT(EXIT_FAILURE);
1072  }
1073 
1074  while( true )
1075  {
1076  int res;
1077  long nelem;
1078  unsigned int ns;
1079  AEInfo *ptr;
1080 
1081  GetNextLine( chFile, fdes, chString );
1082  res = sscanf( chString, "%ld", &nelem );
1083  ASSERT( res == 1 );
1084 
1085  if( nelem < 0 )
1086  break;
1087 
1088  ASSERT( nelem < LIMELM );
1089 
1090  ptr = new AEInfo;
1091 
1092  GetNextLine( chFile, fdes, chString );
1093  res = sscanf( chString, "%u", &ptr->nSubShell );
1094  ASSERT( res == 1 && ptr->nSubShell > 0 );
1095 
1096  ptr->nData.resize( ptr->nSubShell );
1097  ptr->IonThres.resize( ptr->nSubShell );
1098  ptr->Energy.resize( ptr->nSubShell );
1099  ptr->AvNumber.resize( ptr->nSubShell );
1100 
1101  for( ns=0; ns < ptr->nSubShell; ns++ )
1102  {
1103  unsigned int ss;
1104 
1105  GetNextLine( chFile, fdes, chString );
1106  res = sscanf( chString, "%u", &ss );
1107  ASSERT( res == 1 && ns == ss );
1108 
1109  GetNextLine( chFile, fdes, chString );
1110  res = sscanf( chString, "%le", &ptr->IonThres[ns] );
1111  ASSERT( res == 1 );
1112  ptr->IonThres[ns] /= EVRYD;
1113 
1114  GetNextLine( chFile, fdes, chString );
1115  res = sscanf( chString, "%u", &ptr->nData[ns] );
1116  ASSERT( res == 1 && ptr->nData[ns] > 0 );
1117 
1118  ptr->Energy[ns].resize( ptr->nData[ns] );
1119  ptr->AvNumber[ns].resize( ptr->nData[ns] );
1120 
1121  for( unsigned int n=0; n < ptr->nData[ns]; n++ )
1122  {
1123  GetNextLine( chFile, fdes, chString );
1124  res = sscanf(chString,"%le %le",&ptr->Energy[ns][n],&ptr->AvNumber[ns][n]);
1125  ASSERT( res == 2 );
1126  ptr->Energy[ns][n] /= EVRYD;
1127  ASSERT( ptr->Energy[ns][n] < ptr->IonThres[ns] );
1128  }
1129  }
1130 
1131  ASSERT( gv.AugerData[nelem] == NULL );
1132  gv.AugerData[nelem] = ptr;
1133  }
1134 
1135  fclose( fdes );
1136 }
1137 
1139 STATIC void InitBinAugerData(size_t nd,
1140  long ipBegin,
1141  long ipEnd)
1142 {
1143  DEBUG_ENTRY( "InitBinAugerData()" );
1144 
1145  long i, ipLo, nelem;
1146  unsigned int ns;
1147 
1148  flex_arr<realnum> temp( ipBegin, ipEnd );
1149  temp.zero();
1150 
1151  /* this converts gv.bin[nd]->elmAbund[nelem] to particle density inside the grain */
1152  double norm = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->AvVol;
1153 
1154  /* this loop calculates the probability that photoionization occurs in a given shell */
1155  for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
1156  {
1157  ipLo = max( gv.bin[nd]->sd[ns]->ipLo, ipBegin );
1158 
1159  gv.bin[nd]->sd[ns]->p.realloc( ipEnd );
1160 
1163  for( i=ipLo; i < ipEnd; i++ )
1164  {
1165  long nel,nsh;
1166  double phot_ev,cs;
1167 
1168  phot_ev = rfield.anu[i]*EVRYD;
1169 
1170  if( ns == 0 )
1171  {
1172  /* this is the valence band, defined as the sum of any
1173  * subshell not treated explicitly as an inner shell below */
1174  gv.bin[nd]->sd[ns]->p[i] = 0.;
1175 
1176  for( nelem=ipHYDROGEN; nelem < LIMELM && !gv.lgWD01; nelem++ )
1177  {
1178  if( gv.bin[nd]->elmAbund[nelem] == 0. )
1179  continue;
1180 
1181  long nshmax = Heavy.nsShells[nelem][0];
1182 
1183  for( nsh = gv.AugerData[nelem]->nSubShell; nsh < nshmax; nsh++ )
1184  {
1185  nel = nelem+1;
1186  cs = t_ADfA::Inst().phfit(nelem+1,nel,nsh+1,phot_ev);
1187  gv.bin[nd]->sd[ns]->p[i] +=
1188  (realnum)(norm*gv.bin[nd]->elmAbund[nelem]*cs*1e-18);
1189  }
1190  }
1191 
1192  temp[i] += gv.bin[nd]->sd[ns]->p[i];
1193  }
1194  else
1195  {
1196  /* this is photoionization from inner shells */
1197  nelem = gv.bin[nd]->sd[ns]->nelem;
1198  nel = nelem+1;
1199  nsh = gv.bin[nd]->sd[ns]->ns;
1200  cs = t_ADfA::Inst().phfit(nelem+1,nel,nsh+1,phot_ev);
1201  gv.bin[nd]->sd[ns]->p[i] =
1202  (realnum)(norm*gv.bin[nd]->elmAbund[nelem]*cs*1e-18);
1203  temp[i] += gv.bin[nd]->sd[ns]->p[i];
1204  }
1205  }
1206  }
1207 
1208  for( i=ipBegin; i < ipEnd && !gv.lgWD01; i++ )
1209  {
1210  /* this is Eq. 10 of WDB06 */
1211  if( rfield.anu[i] > 20./EVRYD )
1212  gv.bin[nd]->inv_att_len[i] = temp[i];
1213  }
1214 
1215  for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
1216  {
1217  ipLo = max( gv.bin[nd]->sd[ns]->ipLo, ipBegin );
1218  /* renormalize so that sum of probabilities is 1 */
1219  for( i=ipLo; i < ipEnd; i++ )
1220  {
1221  if( temp[i] > 0. )
1222  gv.bin[nd]->sd[ns]->p[i] /= temp[i];
1223  else
1224  gv.bin[nd]->sd[ns]->p[i] = ( ns == 0 ) ? 1.f : 0.f;
1225  }
1226  }
1227 
1228  temp.clear();
1229 
1230  for( ns=0; ns < gv.bin[nd]->sd.size(); ns++ )
1231  {
1232  long n;
1233  ShellData *sptr = gv.bin[nd]->sd[ns];
1234 
1235  ipLo = max( sptr->ipLo, ipBegin );
1236 
1237  /* initialize the yield for primary electrons */
1238  sptr->y01.realloc( ipEnd );
1239 
1240  for( i=ipLo; i < ipEnd; i++ )
1241  {
1242  double elec_en,yzero,yone;
1243 
1244  elec_en = MAX2(rfield.anu[i] - sptr->ionPot,0.);
1245  yzero = y0psa( nd, ns, i, elec_en );
1246 
1247  /* this is size-dependent geometrical yield enhancement
1248  * defined in Weingartner & Draine, 2001; modified in WDB06 */
1249  yone = y1psa( nd, i, elec_en );
1250 
1251  if( ns == 0 )
1252  {
1253  gv.bin[nd]->y0b06[i] = (realnum)yzero;
1254  sptr->y01[i] = (realnum)yone;
1255  }
1256  else
1257  {
1258  sptr->y01[i] = (realnum)(yzero*yone);
1259  }
1260  }
1261 
1262  /* there are no Auger electrons from the band structure */
1263  if( ns > 0 )
1264  {
1265  /* initialize the yield for Auger electrons */
1266  for( n=0; n < sptr->nData; n++ )
1267  {
1268  sptr->y01A[n].realloc( ipEnd );
1269 
1270  for( i=ipLo; i < ipEnd; i++ )
1271  {
1272  double yzero = sptr->AvNr[n] * y0psa( nd, ns, i, sptr->Ener[n] );
1273 
1274  /* this is size-dependent geometrical yield enhancement
1275  * defined in Weingartner & Draine, 2001; modified in WDB06 */
1276  double yone = y1psa( nd, i, sptr->Ener[n] );
1277 
1278  sptr->y01A[n][i] = (realnum)(yzero*yone);
1279  }
1280  }
1281  }
1282  }
1283 }
1284 
1285 /* read a single line of data from data file */
1286 STATIC void GetNextLine(const char *chFile,
1287  FILE *io,
1288  char chLine[]) /* chLine[FILENAME_PATH_LENGTH_2] */
1289 {
1290  char *str;
1291 
1292  DEBUG_ENTRY( "GetNextLine()" );
1293 
1294  do
1295  {
1296  if( read_whole_line( chLine, FILENAME_PATH_LENGTH_2, io ) == NULL )
1297  {
1298  fprintf( ioQQQ, " Could not read from %s\n", chFile );
1299  if( feof(io) )
1300  fprintf( ioQQQ, " EOF reached\n");
1301  cdEXIT(EXIT_FAILURE);
1302  }
1303  }
1304  while( chLine[0] == '#' );
1305 
1306  /* erase comment part of the line */
1307  str = strstr_s(chLine,"#");
1308  if( str != NULL )
1309  *str = '\0';
1310  return;
1311 }
1312 
1313 STATIC void InitEmissivities(void)
1314 {
1315  double fac,
1316  fac2,
1317  mul,
1318  tdust;
1319  long int i;
1320 
1321  DEBUG_ENTRY( "InitEmissivities()" );
1322 
1323  if( trace.lgTrace && trace.lgDustBug )
1324  {
1325  fprintf( ioQQQ, " InitEmissivities starts\n" );
1326  fprintf( ioQQQ, " ND Tdust Emis BB Check 4pi*a^2*<Q>\n" );
1327  }
1328 
1329  ASSERT( NTOP >= 2 && NDEMS > 2*NTOP );
1330  fac = exp(log(GRAIN_TMID/GRAIN_TMIN)/(double)(NDEMS-NTOP));
1331  tdust = GRAIN_TMIN;
1332  for( i=0; i < NDEMS-NTOP; i++ )
1333  {
1334  gv.dsttmp[i] = log(tdust);
1335  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1336  {
1337  gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
1338  }
1339  tdust *= fac;
1340  }
1341 
1342  /* temperatures above GRAIN_TMID are unrealistic -> make grid gradually coarser */
1343  fac2 = exp(log(GRAIN_TMAX/GRAIN_TMID/powi(fac,NTOP-1))/(double)((NTOP-1)*NTOP/2));
1344  for( i=NDEMS-NTOP; i < NDEMS; i++ )
1345  {
1346  gv.dsttmp[i] = log(tdust);
1347  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1348  {
1349  gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
1350  }
1351  fac *= fac2;
1352  tdust *= fac;
1353  }
1354 
1355  /* sanity checks */
1356  mul = 1.;
1357  ASSERT( fabs(gv.dsttmp[0] - log(GRAIN_TMIN)) < 10.*mul*DBL_EPSILON );
1358  mul = sqrt((double)(NDEMS-NTOP));
1359  ASSERT( fabs(gv.dsttmp[NDEMS-NTOP] - log(GRAIN_TMID)) < 10.*mul*DBL_EPSILON );
1360  mul = (double)NTOP + sqrt((double)NDEMS);
1361  ASSERT( fabs(gv.dsttmp[NDEMS-1] - log(GRAIN_TMAX)) < 10.*mul*DBL_EPSILON );
1362 
1363  /* now find slopes form spline fit */
1364  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1365  {
1366  /* set up coefficients for spline */
1367  spline(gv.bin[nd]->dstems,gv.dsttmp,NDEMS,2e31,2e31,gv.bin[nd]->dstslp);
1368  spline(gv.dsttmp,gv.bin[nd]->dstems,NDEMS,2e31,2e31,gv.bin[nd]->dstslp2);
1369  }
1370 
1371 # if 0
1372  /* test the dstems interpolation */
1373  nd = nint(fudge(0));
1374  ASSERT( nd >= 0 && nd < gv.bin.size() );
1375  for( i=0; i < 2000; i++ )
1376  {
1377  double x,y,z;
1378  z = pow(10.,-40. + (double)i/50.);
1379  splint(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,log(z),&y);
1380  if( exp(y) > GRAIN_TMIN && exp(y) < GRAIN_TMAX )
1381  {
1382  x = PlanckIntegral(exp(y),nd,3);
1383  printf(" input %.6e temp %.6e output %.6e rel. diff. %.6e\n",z,exp(y),x,(x-z)/z);
1384  }
1385  }
1386  cdEXIT(EXIT_SUCCESS);
1387 # endif
1388  return;
1389 }
1390 
1391 
1392 /* PlanckIntegral compute total radiative cooling due to grains */
1393 STATIC double PlanckIntegral(double tdust,
1394  size_t nd,
1395  long int ip)
1396 {
1397  long int i;
1398 
1399  double arg,
1400  ExpM1,
1401  integral1 = 0., /* integral(Planck) */
1402  integral2 = 0., /* integral(Planck*abs_cs) */
1403  Planck1,
1404  Planck2,
1405  TDustRyg,
1406  x;
1407 
1408  DEBUG_ENTRY( "PlanckIntegral()" );
1409 
1410  /******************************************************************
1411  *
1412  * >>>chng 99 mar 12, this sub rewritten following Peter van Hoof
1413  * comments. Original coding was in single precision, and for
1414  * very low temperature the exponential was set to zero. As
1415  * a result Q was far too large for grain temperatures below 10K
1416  *
1417  ******************************************************************/
1418 
1419  /* Boltzmann factors for Planck integration */
1420  TDustRyg = TE1RYD/tdust;
1421 
1422  x = 0.999*log(DBL_MAX);
1423 
1424  for( i=0; i < rfield.nupper; i++ )
1425  {
1426  /* this is hnu/kT for grain at this temp and photon energy */
1427  arg = TDustRyg*rfield.anu[i];
1428 
1429  /* want the number exp(hnu/kT) - 1, two expansions */
1430  if( arg < 1.e-5 )
1431  {
1432  /* for small arg expand exp(hnu/kT) - 1 to second order */
1433  ExpM1 = arg*(1. + arg/2.);
1434  }
1435  else
1436  {
1437  /* for large arg, evaluate the full expansion */
1438  ExpM1 = exp(MIN2(x,arg)) - 1.;
1439  }
1440 
1441  Planck1 = PI4*2.*HPLANCK/POW2(SPEEDLIGHT)*POW2(FR1RYD)*POW2(FR1RYD)*
1442  rfield.anu3[i]/ExpM1*rfield.widflx[i];
1443  Planck2 = Planck1*gv.bin[nd]->dstab1[i];
1444 
1445  /* add integral over RJ tail, maybe useful for extreme low temps */
1446  if( i == 0 )
1447  {
1448  integral1 = Planck1/rfield.widflx[0]*rfield.anu[0]/3.;
1449  integral2 = Planck2/rfield.widflx[0]*rfield.anu[0]/5.;
1450  }
1451  /* if we are in the Wien tail - exit */
1452  if( Planck1/integral1 < DBL_EPSILON && Planck2/integral2 < DBL_EPSILON )
1453  break;
1454 
1455  integral1 += Planck1;
1456  integral2 += Planck2;
1457  }
1458 
1459  /* this is an option to print out every few steps, when 'trace grains' is set */
1460  if( trace.lgDustBug && trace.lgTrace && ip%10 == 0 )
1461  {
1462  fprintf( ioQQQ, " %4ld %11.4e %11.4e %11.4e %11.4e\n", (unsigned long)nd, tdust,
1463  integral2, integral1/4./5.67051e-5/powi(tdust,4), integral2*4./integral1 );
1464  }
1465 
1466  ASSERT( integral2 > 0. );
1467  return integral2;
1468 }
1469 
1470 
1471 /* invalidate charge dependent data from previous iteration */
1472 STATIC void NewChargeData(long nd)
1473 {
1474  long nz;
1475 
1476  DEBUG_ENTRY( "NewChargeData()" );
1477 
1478  for( nz=0; nz < NCHS; nz++ )
1479  {
1480  gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
1481  gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
1482  gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
1483  gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
1484  gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
1485 
1487  gv.bin[nd]->chrg[nz]->ThermRate = -DBL_MAX;
1488  gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
1489  gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
1490 
1491  gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
1492  gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
1493  gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
1494  }
1495 
1496  if( !fp_equal(phycon.te,gv.GrnRecomTe) )
1497  {
1498  for( nz=0; nz < NCHS; nz++ )
1499  {
1500  memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
1501  memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );
1502  }
1503  }
1504 
1505  if( nzone != gv.nzone )
1506  {
1507  for( nz=0; nz < NCHS; nz++ )
1508  {
1509  gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
1510  }
1511  }
1512  return;
1513 }
1514 
1515 
1516 /* GrnStdDpth sets the standard behavior of the grain abundance as a function
1517  * of depth into cloud - user-define code should go in GrnVryDpth */
1518 STATIC double GrnStdDpth(long int nd)
1519 {
1520  double GrnStdDpth_v;
1521 
1522  DEBUG_ENTRY( "GrnStdDpth()" );
1523 
1524  /* NB NB - this routine defines the standard depth dependence of the grain abundance,
1525  * to implement user-defined behavior of the abundance (invoked with the "function"
1526  * keyword on the command line), modify the routine GrnVryDpth at the end of this file,
1527  * DO NOT MODIFY THIS ROUTINE */
1528 
1529  if( gv.bin[nd]->nDustFunc == DF_STANDARD )
1530  {
1531  if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
1532  {
1533  /* default behavior for PAH's */
1534  if( gv.chPAH_abundance == "H" )
1535  {
1536  /* the scale factor is the hydrogen atomic fraction, small when gas is ionized
1537  * or molecular and unity when atomic. This function is observed for PAHs
1538  * across the Orion Bar, the PAHs are strong near the ionization front and
1539  * weak in the ionized and molecular gas */
1540  /* >>chng 04 sep 28, propto atomic fraction */
1541  GrnStdDpth_v = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
1542  }
1543  else if( gv.chPAH_abundance == "H,H2" )
1544  {
1545  /* the scale factor is the total of the hydrogen atomic and molecular fractions,
1546  * small when gas is ionized and unity when atomic or molecular. This function is
1547  * observed for PAHs across the Orion Bar, the PAHs are strong near the ionization
1548  * front and weak in the ionized and molecular gas */
1549  /* >>chng 04 sep 28, propto atomic fraction */
1550  GrnStdDpth_v = (dense.xIonDense[ipHYDROGEN][0]+2*hmi.H2_total)/dense.gas_phase[ipHYDROGEN];
1551  }
1552  else if( gv.chPAH_abundance == "CON" )
1553  {
1554  /* constant abundance - unphysical, used only for testing */
1555  GrnStdDpth_v = 1.;
1556  }
1557  else
1558  {
1559  fprintf(ioQQQ,"Invalid argument to SET PAH: %s\n",gv.chPAH_abundance.c_str());
1560  TotalInsanity();
1561  }
1562  }
1563  else
1564  {
1565  /* default behavior for all other types of grains */
1566  GrnStdDpth_v = 1.;
1567  }
1568  }
1569  else if( gv.bin[nd]->nDustFunc == DF_USER_FUNCTION )
1570  {
1571  GrnStdDpth_v = GrnVryDpth(nd);
1572  }
1573  else if( gv.bin[nd]->nDustFunc == DF_SUBLIMATION )
1574  {
1575  // abundance depends on temperature relative to sublimation
1576  // "grain function sublimation" command
1577  GrnStdDpth_v = sexp( pow3( gv.bin[nd]->tedust / gv.bin[nd]->Tsublimat ) );
1578  }
1579  else
1580  {
1581  TotalInsanity();
1582  }
1583 
1584  GrnStdDpth_v = max(1.e-10,GrnStdDpth_v);
1585 
1586  return GrnStdDpth_v;
1587 }
1588 
1589 
1590 /* this is the main routine that drives the grain physics */
1591 void GrainDrive(void)
1592 {
1593  DEBUG_ENTRY( "GrainDrive()" );
1594 
1595  /* gv.lgGrainPhysicsOn set false with no grain physics command */
1596  if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
1597  {
1598  static double tesave = -1.;
1599  static long int nzonesave = -1;
1600 
1601  /* option to only reevaluate grain physics if something has changed.
1602  * gv.lgReevaluate is set false with keyword no reevaluate on grains command
1603  * option to force constant reevaluation of grain physics -
1604  * by default is true
1605  * usually reevaluate grains at all times, but NO REEVALUATE will
1606  * save some time but may affect stability */
1607  if( gv.lgReevaluate || conv.lgSearch || nzonesave != nzone ||
1608  /* need to reevaluate the grains when temp changes since */
1609  ! fp_equal( phycon.te, tesave ) ||
1610  /* >>chng 03 dec 30, check that electrons locked in grains are not important,
1611  * if they are, then reevaluate */
1612  fabs(gv.TotalEden)/dense.eden > conv.EdenErrorAllowed/5. ||
1613  /* >>chng 04 aug 06, always reevaluate when thermal effects of grains are important,
1614  * first is collisional energy exchange with gas, second is grain photoionization */
1615  (fabs( gv.GasCoolColl ) + fabs( thermal.heating[0][13] ))/SDIV(thermal.ctot)>0.1 )
1616  {
1617  nzonesave = nzone;
1618  tesave = phycon.te;
1619 
1620  if( trace.nTrConvg >= 5 )
1621  {
1622  fprintf( ioQQQ, " grain5 calling GrainChargeTemp\n");
1623  }
1624  /* find dust charge and temperature - this must be called at least once per zone
1625  * since grain abundances, set here, may change with depth */
1626  GrainChargeTemp();
1627 
1628  /* >>chng 04 jan 31, moved call to GrainDrift to ConvPresTempEdenIoniz(), PvH */
1629  }
1630  }
1631  else if( gv.lgDustOn() && !gv.lgGrainPhysicsOn )
1632  {
1633  /* very minimalistic treatment of grains; only extinction of continuum is considered
1634  * however, the absorbed energy is not reradiated, so this creates thermal imbalance! */
1635  if( nCalledGrainDrive == 0 )
1636  {
1637  long nelem, ion, ion_to;
1638 
1639  /* when not doing grain physics still want some exported quantities
1640  * to be reasonable, grain temperature used for H2 formation */
1641  gv.GasHeatPhotoEl = 0.;
1642  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1643  {
1644  long nz;
1645 
1646  /* this disables warnings about PAHs in the ionized region */
1647  gv.bin[nd]->lgPAHsInIonizedRegion = false;
1648 
1649  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
1650  {
1651  gv.bin[nd]->chrg[nz]->DustZ = nz;
1652  gv.bin[nd]->chrg[nz]->FracPop = ( nz == 0 ) ? 1. : 0.;
1653  gv.bin[nd]->chrg[nz]->nfill = 0;
1654  gv.bin[nd]->chrg[nz]->tedust = 100.f;
1655  }
1656 
1657  gv.bin[nd]->AveDustZ = 0.;
1658  gv.bin[nd]->dstpot = chrg2pot(0.,nd);
1659 
1660  gv.bin[nd]->tedust = 100.f;
1661  gv.bin[nd]->TeGrainMax = 100.;
1662 
1663  /* set all heating/cooling agents to zero */
1664  gv.bin[nd]->BolFlux = 0.;
1665  gv.bin[nd]->GrainCoolTherm = 0.;
1666  gv.bin[nd]->GasHeatPhotoEl = 0.;
1667  gv.bin[nd]->GrainHeat = 0.;
1668  gv.bin[nd]->GrainHeatColl = 0.;
1669  gv.bin[nd]->ChemEn = 0.;
1670  gv.bin[nd]->ChemEnH2 = 0.;
1671  gv.bin[nd]->thermionic = 0.;
1672 
1673  gv.bin[nd]->lgUseQHeat = false;
1674  gv.bin[nd]->lgEverQHeat = false;
1675  gv.bin[nd]->QHeatFailures = 0;
1676 
1677  gv.bin[nd]->DustDftVel = 0.;
1678 
1679  gv.bin[nd]->avdust = gv.bin[nd]->tedust;
1680  gv.bin[nd]->avdft = 0.f;
1681  gv.bin[nd]->avdpot = (realnum)(gv.bin[nd]->dstpot*EVRYD);
1682  gv.bin[nd]->avDGRatio = -1.f;
1683 
1684  /* >>chng 06 jul 21, add this here as well as in GrainTemperature so that can
1685  * get fake heating when grain physics is turned off */
1686  if( 0 && gv.lgBakesPAH_heat )
1687  {
1688  /* this is a dirty hack to get BT94 PE heating rate
1689  * for PAH's included, for Lorentz Center 2004 PDR meeting, PvH */
1690  /*>>>refer PAH heating Bakes, E.L.O., & Tielens, A.G.G.M. 1994, ApJ, 427, 822 */
1691  /* >>chng 05 aug 12, change from +=, which added additional heating to what exists already,
1692  * to simply = to set the heat, this equation gives total heating */
1693  gv.bin[nd]->GasHeatPhotoEl = 1.e-24*hmi.UV_Cont_rel2_Habing_TH85_depth*
1694  dense.gas_phase[ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((hmi.UV_Cont_rel2_Habing_TH85_depth*
1695  sqrt(phycon.te)/dense.eden),0.73)) + 3.65e-2*pow(phycon.te/1.e4,0.7)/
1696  (1.+2.e-4*(hmi.UV_Cont_rel2_Habing_TH85_depth*sqrt(phycon.te)/dense.eden)))/gv.bin.size() *
1697  gv.GrainHeatScaleFactor;
1698  gv.GasHeatPhotoEl += gv.bin[nd]->GasHeatPhotoEl;
1699  }
1700  }
1701 
1702  gv.TotalEden = 0.;
1703  gv.GrnElecDonateMax = 0.f;
1704  gv.GrnElecHoldMax = 0.f;
1705 
1706  for( nelem=0; nelem < LIMELM; nelem++ )
1707  {
1708  for( ion=0; ion <= nelem+1; ion++ )
1709  {
1710  for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1711  {
1712  gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
1713  }
1714  }
1715  }
1716 
1717  /* set all heating/cooling agents to zero */
1718  gv.GrainHeatInc = 0.;
1719  gv.GrainHeatDif = 0.;
1720  gv.GrainHeatLya = 0.;
1721  gv.GrainHeatCollSum = 0.;
1722  gv.GrainHeatSum = 0.;
1723  gv.GrainHeatChem = 0.;
1724  gv.GasCoolColl = 0.;
1725  gv.TotalDustHeat = 0.f;
1726  gv.dphmax = 0.f;
1727  gv.dclmax = 0.f;
1728 
1729  thermal.heating[0][13] = 0.;
1730  thermal.heating[0][14] = 0.;
1731  thermal.heating[0][25] = 0.;
1732  }
1733 
1734  if( nCalledGrainDrive == 0 || gv.lgAnyDustVary )
1735  {
1736  GrainUpdateRadius1();
1737  GrainUpdateRadius2();
1738  }
1739  }
1740 
1741  ++nCalledGrainDrive;
1742  return;
1743 }
1744 
1745 /* iterate grain charge and temperature */
1746 STATIC void GrainChargeTemp(void)
1747 {
1748  long int i,
1749  ion,
1750  ion_to,
1751  nelem,
1752  nz;
1753  realnum dccool = FLT_MAX;
1754  double delta,
1755  GasHeatNet,
1756  hcon = DBL_MAX,
1757  hla = DBL_MAX,
1758  hots = DBL_MAX,
1759  oldtemp,
1760  oldTotalEden,
1761  ratio,
1762  ThermRatio;
1763 
1764  static long int oldZone = -1;
1765  static double oldTe = -DBL_MAX,
1766  oldHeat = -DBL_MAX;
1767 
1768  DEBUG_ENTRY( "GrainChargeTemp()" );
1769 
1770  if( trace.lgTrace && trace.lgDustBug )
1771  {
1772  fprintf( ioQQQ, "\n GrainChargeTemp called lgSearch%2c\n\n", TorF(conv.lgSearch) );
1773  }
1774 
1775  oldTotalEden = gv.TotalEden;
1776 
1777  /* these will sum heating agents over grain populations */
1778  gv.GrainHeatInc = 0.;
1779  gv.GrainHeatDif = 0.;
1780  gv.GrainHeatLya = 0.;
1781  gv.GrainHeatCollSum = 0.;
1782  gv.GrainHeatSum = 0.;
1783  gv.GrainHeatChem = 0.;
1784 
1785  gv.GasCoolColl = 0.;
1786  gv.GasHeatPhotoEl = 0.;
1787  gv.GasHeatTherm = 0.;
1788 
1789  gv.TotalEden = 0.;
1790 
1791  for( nelem=0; nelem < LIMELM; nelem++ )
1792  {
1793  for( ion=0; ion <= nelem+1; ion++ )
1794  {
1795  for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1796  {
1797  gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
1798  }
1799  }
1800  }
1801 
1802  gv.HighestIon = HighestIonStage();
1803 
1804  /* this sets dstAbund and conversion factors, but not gv.dstab and gv.dstsc! */
1805  GrainUpdateRadius1();
1806 
1807  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1808  {
1809  double one;
1810  double ChTdBracketLo = 0., ChTdBracketHi = -DBL_MAX;
1811  long relax = ( conv.lgSearch ) ? 3 : 1;
1812 
1813  /* >>chng 02 nov 11, added test for the presence of PAHs in the ionized region, PvH */
1814  if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
1815  {
1816  if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.50 )
1817  {
1818  gv.bin[nd]->lgPAHsInIonizedRegion = true;
1819  }
1820  }
1821 
1822  /* >>chng 01 sep 13, dynamically allocate backup store, remove ncell dependence, PvH */
1823  /* allocate data inside loop to avoid accidental spillover to next iteration */
1824  /* >>chng 04 jan 18, no longer delete and reallocate data inside loop to speed up the code, PvH */
1825  NewChargeData(nd);
1826 
1827  if( trace.lgTrace && trace.lgDustBug )
1828  {
1829  fprintf( ioQQQ, " >>GrainChargeTemp starting grain %s\n",
1830  gv.bin[nd]->chDstLab );
1831  }
1832 
1833  delta = 2.*TOLER;
1834  /* >>chng 01 nov 29, relax max no. of iterations during initial search */
1835  for( i=0; i < relax*CT_LOOP_MAX && delta > TOLER; ++i )
1836  {
1837  string which;
1838  long j;
1839  double TdBracketLo = 0., TdBracketHi = -DBL_MAX;
1840  double ThresEst = 0.;
1841  oldtemp = gv.bin[nd]->tedust;
1842 
1843  /* solve for charge using previous estimate for grain temp
1844  * grain temp only influences thermionic emissions
1845  * Thermratio is fraction thermionic emissions contribute
1846  * to the total electron loss rate of the grain */
1847  GrainCharge(nd,&ThermRatio);
1848 
1849  ASSERT( gv.bin[nd]->GrainHeat > 0. );
1850  ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
1851 
1852  /* >>chng 04 may 31, in conditions where collisions become an important
1853  * heating/cooling source (e.g. gas that is predominantly heated by cosmic
1854  * rays), the heating rate depends strongly on the assumed dust temperature.
1855  * hence it is necessary to iterate for the dust temperature. PvH */
1856  gv.bin[nd]->lgTdustConverged = false;
1857  for( j=0; j < relax*T_LOOP_MAX; ++j )
1858  {
1859  double oldTemp2 = gv.bin[nd]->tedust;
1860  double oldHeat2 = gv.bin[nd]->GrainHeat;
1861  double oldCool = gv.bin[nd]->GrainGasCool;
1862 
1863  /* now solve grain temp using new value for grain potential */
1864  GrainTemperature(nd,&dccool,&hcon,&hots,&hla);
1865 
1866  gv.bin[nd]->GrainGasCool = dccool;
1867 
1868  if( trace.lgTrace && trace.lgDustBug )
1869  {
1870  fprintf( ioQQQ, " >>loop %ld BracketLo %.6e BracketHi %.6e",
1871  j, TdBracketLo, TdBracketHi );
1872  }
1873 
1874  /* this test assures that convergence can only happen if GrainHeat > 0
1875  * and therefore the value of tedust is guaranteed to be valid as well */
1876  /* >>chng 04 aug 05, test that gas cooling is converged as well,
1877  * in deep PDRs gas cooling depends critically on grain temperature, PvH */
1878  if( fabs(gv.bin[nd]->GrainHeat-oldHeat2) < HEAT_TOLER*gv.bin[nd]->GrainHeat &&
1879  fabs(gv.bin[nd]->GrainGasCool-oldCool) < HEAT_TOLER_BIN*thermal.ctot )
1880  {
1881  gv.bin[nd]->lgTdustConverged = true;
1882  if( trace.lgTrace && trace.lgDustBug )
1883  fprintf( ioQQQ, " converged\n" );
1884  break;
1885  }
1886 
1887  /* update the bracket for the solution */
1888  if( gv.bin[nd]->tedust < oldTemp2 )
1889  TdBracketHi = oldTemp2;
1890  else
1891  TdBracketLo = oldTemp2;
1892 
1893  /* GrainTemperature yields a new estimate for tedust, and initially
1894  * that estimate will be used. In most zones this will converge quickly.
1895  * However, sometimes the solution will oscillate and converge very
1896  * slowly. So, as soon as j >= 2 and the bracket is set up, we will
1897  * force convergence by using a bisection search within the bracket */
1900  /* this test assures that TdBracketHi is initialized */
1901  if( TdBracketHi > TdBracketLo )
1902  {
1903  /* if j >= 2, the solution is converging too slowly
1904  * so force convergence by doing a bisection search */
1905  if( ( j >= 2 && TdBracketLo > 0. ) ||
1906  gv.bin[nd]->tedust <= TdBracketLo ||
1907  gv.bin[nd]->tedust >= TdBracketHi )
1908  {
1909  gv.bin[nd]->tedust = (realnum)(0.5*(TdBracketLo + TdBracketHi));
1910  if( trace.lgTrace && trace.lgDustBug )
1911  fprintf( ioQQQ, " bisection\n" );
1912  }
1913  else
1914  {
1915  if( trace.lgTrace && trace.lgDustBug )
1916  fprintf( ioQQQ, " iteration\n" );
1917  }
1918  }
1919  else
1920  {
1921  if( trace.lgTrace && trace.lgDustBug )
1922  fprintf( ioQQQ, " iteration\n" );
1923  }
1924 
1925  ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
1926  }
1927 
1928  if( gv.bin[nd]->lgTdustConverged )
1929  {
1930  /* update the bracket for the solution */
1931  if( gv.bin[nd]->tedust < oldtemp )
1932  ChTdBracketHi = oldtemp;
1933  else
1934  ChTdBracketLo = oldtemp;
1935  }
1936  else
1937  {
1938  bool lgBoundErr;
1939  double y, x = log(gv.bin[nd]->tedust);
1940  /* make sure GrainHeat is consistent with value of tedust */
1941  splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgBoundErr);
1942  gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;
1943 
1944  fprintf( ioQQQ," PROBLEM temperature of grain species %s (Tg=%.3eK) not converged\n",
1945  gv.bin[nd]->chDstLab , gv.bin[nd]->tedust );
1946  ConvFail("grai","");
1947  if( lgAbort )
1948  {
1949  return;
1950  }
1951  }
1952 
1953  ASSERT( gv.bin[nd]->GrainHeat > 0. );
1954  ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
1955 
1956  /* delta estimates relative change in electron emission rate
1957  * due to the update in the grain temperature, if it is small
1958  * we won't bother to iterate (which is usually the case)
1959  * the formula assumes that thermionic emission is the only
1960  * process that depends on grain temperature */
1962  ratio = gv.bin[nd]->tedust/oldtemp;
1963  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
1964  {
1965  ThresEst += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->ThresInf;
1966  }
1967  delta = ThresEst*TE1RYD/gv.bin[nd]->tedust*(ratio - 1.);
1969  delta = ( delta < 0.9*log(DBL_MAX) ) ?
1970  ThermRatio*fabs(POW2(ratio)*exp(delta)-1.) : DBL_MAX;
1971 
1972  /* >>chng 06 feb 07, bracket grain temperature to force convergence when oscillating, PvH */
1973  if( delta > TOLER )
1974  {
1975  if( trace.lgTrace && trace.lgDustBug )
1976  which = "iteration";
1977 
1978  /* The loop above yields a new estimate for tedust, and initially that
1979  * estimate will be used. In most zones this will converge very quickly.
1980  * However, sometimes the solution will oscillate and converge very
1981  * slowly. So, as soon as i >= 2 and the bracket is set up, we will
1982  * force convergence by using a bisection search within the bracket */
1985  /* this test assures that ChTdBracketHi is initialized */
1986  if( ChTdBracketHi > ChTdBracketLo )
1987  {
1988  /* if i >= 2, the solution is converging too slowly
1989  * so force convergence by doing a bisection search */
1990  if( ( i >= 2 && ChTdBracketLo > 0. ) ||
1991  gv.bin[nd]->tedust <= ChTdBracketLo ||
1992  gv.bin[nd]->tedust >= ChTdBracketHi )
1993  {
1994  gv.bin[nd]->tedust = (realnum)(0.5*(ChTdBracketLo + ChTdBracketHi));
1995  if( trace.lgTrace && trace.lgDustBug )
1996  which = "bisection";
1997  }
1998  }
1999  }
2000 
2001  if( trace.lgTrace && trace.lgDustBug )
2002  {
2003  fprintf( ioQQQ, " >>GrainChargeTemp finds delta=%.4e, ", delta );
2004  fprintf( ioQQQ, " old/new temp=%.5e %.5e, ", oldtemp, gv.bin[nd]->tedust );
2005  if( delta > TOLER )
2006  fprintf( ioQQQ, "doing another %s\n", which.c_str() );
2007  else
2008  fprintf( ioQQQ, "converged\n" );
2009  }
2010  }
2011  if( delta > TOLER )
2012  {
2013  fprintf( ioQQQ, " PROBLEM charge/temperature not converged for %s zone %.2f\n",
2014  gv.bin[nd]->chDstLab , fnzone );
2015  ConvFail("grai","");
2016  }
2017 
2018  /* add in ion recombination rates on this grain species */
2019  /* ionbal.lgGrainIonRecom is 1 by default, set to 0 with
2020  * no grain neutralization command */
2021  if( ionbal.lgGrainIonRecom )
2022  GrainChrgTransferRates(nd);
2023 
2024  /* >>chng 04 jan 31, moved call to UpdateRadius2 outside loop, PvH */
2025 
2026  /* following used to keep track of heating agents in printout
2027  * no physics done with GrainHeatInc
2028  * dust heating by incident continuum, and elec friction before ejection */
2029  gv.GrainHeatInc += hcon;
2030  /* remember total heating by diffuse fields, for printout (includes Lya) */
2031  gv.GrainHeatDif += hots;
2032  /* GrainHeatLya - total heating by LA in this zone, erg cm-3 s-1, only here
2033  * for eventual printout, hots is total ots line heating */
2034  gv.GrainHeatLya += hla;
2035 
2036  /* this will be total collisional heating, for printing in lines */
2037  gv.GrainHeatCollSum += gv.bin[nd]->GrainHeatColl;
2038 
2039  /* GrainHeatSum is total heating of all grain types in this zone,
2040  * will be carried by total cooling, only used in lines to print tot heat
2041  * printed as entry "GraT 0 " */
2042  gv.GrainHeatSum += gv.bin[nd]->GrainHeat;
2043 
2044  /* net amount of chemical energy donated by recombining ions and molecule formation */
2045  gv.GrainHeatChem += gv.bin[nd]->ChemEn + gv.bin[nd]->ChemEnH2;
2046 
2047  /* dccool is gas cooling due to collisions with grains - negative if net heating
2048  * zero if NO GRAIN GAS COLLISIONAL EXCHANGE command included */
2049  gv.GasCoolColl += dccool;
2050  gv.GasHeatPhotoEl += gv.bin[nd]->GasHeatPhotoEl;
2051  gv.GasHeatTherm += gv.bin[nd]->thermionic;
2052 
2053  /* this is grain charge in e/cm^3, positive number means grain supplied free electrons */
2054  /* >>chng 01 mar 24, changed DustZ+1 to DustZ, PvH */
2055  one = 0.;
2056  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2057  {
2058  one += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
2059  gv.bin[nd]->cnv_GR_pCM3;
2060  }
2061  /* electron density contributed by grains, cm-3 */
2062  gv.TotalEden += one;
2063  {
2064  /*@-redef@*/
2065  enum {DEBUG_LOC=false};
2066  /*@+redef@*/
2067  if( DEBUG_LOC )
2068  {
2069  fprintf(ioQQQ," DEBUG grn chr nz\t%.2f\teden\t%.3e\tnd\t%li",
2070  fnzone,
2071  dense.eden,
2072  (unsigned long)nd);
2073  fprintf(ioQQQ,"\tne\t%.2e\tAveDustZ\t%.2e\t%.2e\t%.2e\t%.2e",
2074  one,
2075  gv.bin[nd]->AveDustZ,
2076  gv.bin[nd]->chrg[0]->FracPop,(double)gv.bin[nd]->chrg[0]->DustZ,
2077  gv.bin[nd]->cnv_GR_pCM3);
2078  fprintf(ioQQQ,"\n");
2079  }
2080  }
2081 
2082  if( trace.lgTrace && trace.lgDustBug )
2083  {
2084  fprintf(ioQQQ," %s Pot %.5e Thermal %.5e GasCoolColl %.5e" ,
2085  gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot, gv.bin[nd]->GrainHeat, dccool );
2086  fprintf(ioQQQ," GasPEHeat %.5e GasThermHeat %.5e ChemHeat %.5e\n\n" ,
2087  gv.bin[nd]->GasHeatPhotoEl, gv.bin[nd]->thermionic, gv.bin[nd]->ChemEn );
2088  }
2089  }
2090 
2091  /* >>chng 04 aug 06, added test of convergence of the net gas heating/cooling, PvH */
2092  GasHeatNet = gv.GasHeatPhotoEl + gv.GasHeatTherm - gv.GasCoolColl;
2093 
2094  if( !fp_equal(phycon.te,gv.GrnRecomTe) )
2095  {
2096  oldZone = gv.nzone;
2097  oldTe = gv.GrnRecomTe;
2098  oldHeat = gv.GasHeatNet;
2099  }
2100 
2101  /* >>chng 04 aug 07, added estimate for heating derivative, PvH */
2102  if( nzone == oldZone && !fp_equal(phycon.te,oldTe) )
2103  {
2104  gv.dHeatdT = (GasHeatNet-oldHeat)/(phycon.te-oldTe);
2105  }
2106 
2107  /* >>chng 04 sep 15, add test for convergence of gv.TotalEden, PvH */
2108  if( nzone != gv.nzone || !fp_equal(phycon.te,gv.GrnRecomTe) ||
2109  fabs(gv.GasHeatNet-GasHeatNet) > HEAT_TOLER*thermal.ctot ||
2110  fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
2111  {
2112  /* >>chng 04 aug 07, add test whether eden on grain converged */
2113  /* flag that change in eden was too large */
2114  /*conv.lgConvEden = false;*/
2115  if( fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
2116  {
2117  conv.setConvIonizFail( "grn eden chg" , oldTotalEden, gv.TotalEden);
2118  }
2119  else if( fabs(gv.GasHeatNet-GasHeatNet) > HEAT_TOLER*thermal.ctot )
2120  {
2121  conv.setConvIonizFail( "grn het chg" , gv.GasHeatNet, GasHeatNet);
2122  }
2123  else if( !fp_equal(phycon.te,gv.GrnRecomTe) )
2124  {
2125  conv.setConvIonizFail( "grn ter chg" , gv.GrnRecomTe, phycon.te);
2126  }
2127  else if( nzone != gv.nzone )
2128  {
2129  conv.setConvIonizFail( "grn zon chg" , gv.nzone, nzone );
2130  }
2131  else
2132  TotalInsanity();
2133  }
2134 
2135  /* printf( "DEBUG GasHeatNet %.6e -> %.6e TotalEden %e -> %e conv.lgConvIoniz %c\n",
2136  gv.GasHeatNet, GasHeatNet, gv.TotalEden, oldTotalEden, TorF(conv.lgConvIoniz) ); */
2137  /* printf( "DEBUG %.2f %e %e\n", fnzone, phycon.te, dense.eden ); */
2138 
2139  /* update total grain opacities in gv.dstab and gv.dstsc,
2140  * they depend on grain charge and may depend on depth
2141  * also add in the photo-dissociation cs in gv.dstab */
2142  GrainUpdateRadius2();
2143 
2144  gv.nzone = nzone;
2145  gv.GrnRecomTe = phycon.te;
2146  gv.GasHeatNet = GasHeatNet;
2147 
2148 #ifdef WD_TEST2
2149  printf("wd test: proton fraction %.5e Total DustZ %.6f heating/cooling rate %.5e %.5e\n",
2150  dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
2151  gv.bin[0]->AveDustZ,gv.GasHeatPhotoEl/dense.gas_phase[ipHYDROGEN]/fudge(0),
2152  gv.GasCoolColl/dense.gas_phase[ipHYDROGEN]/fudge(0));
2153 #endif
2154 
2155  if( trace.lgTrace )
2156  {
2157  /*@-redef@*/
2158  enum {DEBUG_LOC=false};
2159  /*@+redef@*/
2160  if( DEBUG_LOC )
2161  {
2162  fprintf( ioQQQ, " %.2f Grain surface charge transfer rates\n", fnzone );
2163 
2164  for( nelem=0; nelem < LIMELM; nelem++ )
2165  {
2166  if( dense.lgElmtOn[nelem] )
2167  {
2168  printf( " %s:", elementnames.chElementSym[nelem] );
2169  for( ion=dense.IonLow[nelem]; ion <= dense.IonHigh[nelem]; ++ion )
2170  {
2171  for( ion_to=0; ion_to <= nelem+1; ion_to++ )
2172  {
2173  if( gv.GrainChTrRate[nelem][ion][ion_to] > 0.f )
2174  {
2175  printf( " %ld->%ld %.2e", ion, ion_to,
2176  gv.GrainChTrRate[nelem][ion][ion_to] );
2177  }
2178  }
2179  }
2180  printf( "\n" );
2181  }
2182  }
2183  }
2184 
2185  fprintf( ioQQQ, " %.2f Grain contribution to electron density %.2e\n",
2186  fnzone , gv.TotalEden );
2187 
2188  fprintf( ioQQQ, " Grain electons: " );
2189  for( size_t nd=0; nd < gv.bin.size(); nd++ )
2190  {
2191  double sum = 0.;
2192  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2193  {
2194  sum += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
2195  gv.bin[nd]->cnv_GR_pCM3;
2196  }
2197  fprintf( ioQQQ, " %.2e", sum );
2198  }
2199  fprintf( ioQQQ, "\n" );
2200 
2201  fprintf( ioQQQ, " Grain potentials:" );
2202  for( size_t nd=0; nd < gv.bin.size(); nd++ )
2203  {
2204  fprintf( ioQQQ, " %.2e", gv.bin[nd]->dstpot );
2205  }
2206  fprintf( ioQQQ, "\n" );
2207 
2208  fprintf( ioQQQ, " Grain temperatures:" );
2209  for( size_t nd=0; nd < gv.bin.size(); nd++ )
2210  {
2211  fprintf( ioQQQ, " %.2e", gv.bin[nd]->tedust );
2212  }
2213  fprintf( ioQQQ, "\n" );
2214 
2215  fprintf( ioQQQ, " GrainCollCool: %.6e\n", gv.GasCoolColl );
2216  }
2217 
2218  /*if( nzone > 900)
2219  fprintf(ioQQQ,"DEBUG cool\t%.2f\t%e\t%e\t%e\n",
2220  fnzone,
2221  phycon.te ,
2222  dense.eden,
2223  gv.GasCoolColl );*/
2224  return;
2225 }
2226 
2227 
2228 STATIC void GrainCharge(size_t nd,
2229  /*@out@*/double *ThermRatio) /* ratio of thermionic to total rate */
2230 {
2231  bool lgBigError,
2232  lgInitial;
2233  long backup,
2234  i,
2235  loopMax,
2236  newZlo,
2237  nz,
2238  power,
2239  stride,
2240  stride0,
2241  Zlo;
2242  double crate,
2243  csum1,
2244  csum1a,
2245  csum1b,
2246  csum2,
2247  csum3,
2248  netloss0 = -DBL_MAX,
2249  netloss1 = -DBL_MAX,
2250  rate_up[NCHU],
2251  rate_dn[NCHU],
2252  step;
2253 
2254  DEBUG_ENTRY( "GrainCharge()" );
2255 
2256  /* find dust charge */
2257  if( trace.lgTrace && trace.lgDustBug )
2258  {
2259  fprintf( ioQQQ, " Charge loop, search %c,", TorF(conv.lgSearch) );
2260  }
2261 
2262  ASSERT( nd < gv.bin.size() );
2263 
2264  for( nz=0; nz < NCHS; nz++ )
2265  {
2266  gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
2267  }
2268 
2269  /* this algorithm determines the value of Zlo and the charge state populations
2270  * in the n-charge state model as described in:
2271  *
2272  * >>refer grain physics van Hoof P.A.M., Weingartner J.C., et al., 2004, MNRAS, 350, 1330
2273  *
2274  * the algorithm first uses the n charge states to bracket the solution by
2275  * separating the charge states with a stride that is an integral power of
2276  * n-1, e.g. like this for an n == 4 calculation:
2277  *
2278  * +gain +gain -gain -gain
2279  * | | | |
2280  * -8 1 10 19
2281  *
2282  * for each of the charge states the total electron emission and recombination
2283  * rates are calculated. the solution is considered properly bracketed if the
2284  * sign of the net gain rate (emission-recombination) is different for the first
2285  * and the last of the n charge states. then the algorithm searches for two
2286  * adjacent charge states where the net gain rate changes sign and divides
2287  * that space into n-1 equal parts, like here
2288  *
2289  * net gain + + + -
2290  * | | | |
2291  * Zg 1 4 7 10
2292  *
2293  * note that the first and the last charge state can be retained from the
2294  * previous iteration and do not need to be recalculated (UpdatePot as well
2295  * as GrainElecEmis1 and GrainElecRecomb1 have mechanisms for re-using
2296  * previously calculated data, so GrainCharge doesn't need to worry about
2297  * this). The dividing algorithm is repeated until the stride equals 1, like
2298  * here
2299  *
2300  * net gain +---
2301  * ||||
2302  * Zg 7 10
2303  *
2304  * finally, the bracket may need to be shifted in order for the criterion
2305  * J1 x J2 <= 0 to be fulfilled (see the paper quoted above for a detailed
2306  * explanation). in the example here one shift might be necessary:
2307  *
2308  * net gain ++--
2309  * ||||
2310  * Zg 6 9
2311  *
2312  * for n == 2, the algorithm would have to be slightly different. In order
2313  * to avoid complicating the code, we force the code to use n == 3 in the
2314  * first two steps of the process, reverting back to n == 2 upon the last
2315  * step. this should not produce any noticeable additional CPU overhead */
2316 
2317  lgInitial = ( gv.bin[nd]->chrg[0]->DustZ == LONG_MIN );
2318 
2319  backup = gv.bin[nd]->nChrg;
2320  gv.bin[nd]->nChrg = MAX2(gv.bin[nd]->nChrg,3);
2321 
2322  stride0 = gv.bin[nd]->nChrg-1;
2323 
2324  /* set up initial bracket for grain charge, will be checked below */
2325  if( lgInitial )
2326  {
2327  double xxx;
2328  step = MAX2((double)(-gv.bin[nd]->LowestZg),1.);
2329  power = (int)(log(step)/log((double)stride0));
2330  power = MAX2(power,0);
2331  xxx = powi((double)stride0,power);
2332  stride = nint(xxx);
2333  Zlo = gv.bin[nd]->LowestZg;
2334  }
2335  else
2336  {
2337  /* the previous solution is the best choice here */
2338  stride = 1;
2339  Zlo = gv.bin[nd]->chrg[0]->DustZ;
2340  }
2341  UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2342 
2343  /* check if the grain charge is correctly bracketed */
2344  for( i=0; i < BRACKET_MAX; i++ )
2345  {
2346  netloss0 = rate_up[0] - rate_dn[0];
2347  netloss1 = rate_up[gv.bin[nd]->nChrg-1] - rate_dn[gv.bin[nd]->nChrg-1];
2348 
2349  if( netloss0*netloss1 <= 0. )
2350  break;
2351 
2352  if( netloss1 > 0. )
2353  Zlo += (gv.bin[nd]->nChrg-1)*stride;
2354 
2355  if( i > 0 )
2356  stride *= stride0;
2357 
2358  if( netloss1 < 0. )
2359  Zlo -= (gv.bin[nd]->nChrg-1)*stride;
2360 
2361  Zlo = MAX2(Zlo,gv.bin[nd]->LowestZg);
2362  UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2363  }
2364 
2365  if( netloss0*netloss1 > 0. ) {
2366  fprintf( ioQQQ, " insanity: could not bracket grain charge for %s\n", gv.bin[nd]->chDstLab );
2367  ShowMe();
2368  cdEXIT(EXIT_FAILURE);
2369  }
2370 
2371  /* home in on the charge */
2372  while( stride > 1 )
2373  {
2374  stride /= stride0;
2375 
2376  netloss0 = rate_up[0] - rate_dn[0];
2377  for( nz=0; nz < gv.bin[nd]->nChrg-1; nz++ )
2378  {
2379  netloss1 = rate_up[nz+1] - rate_dn[nz+1];
2380 
2381  if( netloss0*netloss1 <= 0. )
2382  {
2383  Zlo = gv.bin[nd]->chrg[nz]->DustZ;
2384  break;
2385  }
2386 
2387  netloss0 = netloss1;
2388  }
2389  UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2390  }
2391 
2392  ASSERT( netloss0*netloss1 <= 0. );
2393 
2394  gv.bin[nd]->nChrg = backup;
2395 
2396  /* >>chng 04 feb 15, relax upper limit on initial search when nChrg is much too large, PvH */
2397  loopMax = ( lgInitial ) ? 4*gv.bin[nd]->nChrg : 2*gv.bin[nd]->nChrg;
2398 
2399  lgBigError = true;
2400  for( i=0; i < loopMax; i++ )
2401  {
2402  GetFracPop( nd, Zlo, rate_up, rate_dn, &newZlo );
2403 
2404  if( newZlo == Zlo )
2405  {
2406  lgBigError = false;
2407  break;
2408  }
2409 
2410  Zlo = newZlo;
2411  UpdatePot( nd, Zlo, 1, rate_up, rate_dn );
2412  }
2413 
2414  if( lgBigError ) {
2415  fprintf( ioQQQ, " insanity: could not converge grain charge for %s\n", gv.bin[nd]->chDstLab );
2416  ShowMe();
2417  cdEXIT(EXIT_FAILURE);
2418  }
2419 
2422  gv.bin[nd]->lgChrgConverged = true;
2423 
2424  gv.bin[nd]->AveDustZ = 0.;
2425  crate = csum3 = 0.;
2426  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2427  {
2428  double d[4];
2429 
2430  gv.bin[nd]->AveDustZ += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->DustZ;
2431 
2432  crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
2433  csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[3];
2434  }
2435  gv.bin[nd]->dstpot = chrg2pot(gv.bin[nd]->AveDustZ,nd);
2436  *ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;
2437 
2438  ASSERT( *ThermRatio >= 0. );
2439 
2440  if( trace.lgTrace && trace.lgDustBug )
2441  {
2442  double d[4];
2443 
2444  fprintf( ioQQQ, "\n" );
2445 
2446  crate = csum1a = csum1b = csum2 = csum3 = 0.;
2447  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2448  {
2449  crate += gv.bin[nd]->chrg[nz]->FracPop*
2450  GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
2451  csum1a += gv.bin[nd]->chrg[nz]->FracPop*d[0];
2452  csum1b += gv.bin[nd]->chrg[nz]->FracPop*d[1];
2453  csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[2];
2454  csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[3];
2455  }
2456 
2457  fprintf( ioQQQ, " ElecEm rate1a=%.4e, rate1b=%.4e, ", csum1a, csum1b );
2458  fprintf( ioQQQ, "rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum2, csum3, crate );
2459  if( crate > 0. )
2460  {
2461  fprintf( ioQQQ, " rate1a/sum=%.4e, rate1b/sum=%.4e, ", csum1a/crate, csum1b/crate );
2462  fprintf( ioQQQ, "rate2/sum=%.4e, rate3/sum=%.4e\n", csum2/crate, csum3/crate );
2463  }
2464 
2465  crate = csum1 = csum2 = 0.;
2466  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2467  {
2468  crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecRecomb1(nd,nz,&d[0],&d[1]);
2469  csum1 += gv.bin[nd]->chrg[nz]->FracPop*d[0];
2470  csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[1];
2471  }
2472 
2473  fprintf( ioQQQ, " ElecRc rate1=%.4e, rate2=%.4e, sum=%.4e\n", csum1, csum2, crate );
2474  if( crate > 0. )
2475  {
2476  fprintf( ioQQQ, " rate1/sum=%.4e, rate2/sum=%.4e\n", csum1/crate, csum2/crate );
2477  }
2478 
2479  fprintf( ioQQQ, " Charging rates:" );
2480  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2481  {
2482  fprintf( ioQQQ, " Zg %ld up %.4e dn %.4e",
2483  gv.bin[nd]->chrg[nz]->DustZ, rate_up[nz], rate_dn[nz] );
2484  }
2485  fprintf( ioQQQ, "\n" );
2486 
2487  fprintf( ioQQQ, " FracPop: fnzone %.2f nd %ld AveZg %.5e",
2488  fnzone, (unsigned long)nd, gv.bin[nd]->AveDustZ );
2489  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2490  {
2491  fprintf( ioQQQ, " Zg %ld %.5f", Zlo+nz, gv.bin[nd]->chrg[nz]->FracPop );
2492  }
2493  fprintf( ioQQQ, "\n" );
2494 
2495  fprintf( ioQQQ, " >Grain potential:%12.12s %.6fV",
2496  gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot*EVRYD );
2497  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2498  {
2499  fprintf( ioQQQ, " Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V",
2500  gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInf*EVRYD,
2501  gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInfVal*EVRYD );
2502  }
2503  fprintf( ioQQQ, "\n" );
2504  }
2505  return;
2506 }
2507 
2508 
2509 /* grain electron recombination rates for single charge state */
2510 STATIC double GrainElecRecomb1(size_t nd,
2511  long nz,
2512  /*@out@*/ double *sum1,
2513  /*@out@*/ double *sum2)
2514 {
2515  long ion,
2516  nelem;
2517  double eta,
2518  rate,
2519  Stick,
2520  ve,
2521  xi;
2522 
2523  DEBUG_ENTRY( "GrainElecRecomb1()" );
2524 
2525  ASSERT( nd < gv.bin.size() );
2526  ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
2527 
2528  /* >>chng 01 may 31, try to find cached results first */
2529  /* >>chng 04 feb 11, moved cached data to ChargeBin, PvH */
2530  if( gv.bin[nd]->chrg[nz]->RSum1 >= 0. )
2531  {
2532  *sum1 = gv.bin[nd]->chrg[nz]->RSum1;
2533  *sum2 = gv.bin[nd]->chrg[nz]->RSum2;
2534  rate = *sum1 + *sum2;
2535  return rate;
2536  }
2537 
2538  /* -1 makes psi correct for impact by electrons */
2539  ion = -1;
2540  /* VE is mean (not RMS) electron velocity */
2541  /*ve = TePowers.sqrte*6.2124e5;*/
2542  ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);
2543 
2544  Stick = ( gv.bin[nd]->chrg[nz]->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
2545  /* >>chng 00 jul 19, replace classical results with results including image potential
2546  * to correct for polarization of the grain as charged particle approaches. */
2547  GrainScreen(ion,nd,nz,&eta,&xi);
2548  /* this is grain surface recomb rate for electrons */
2549  *sum1 = ( gv.bin[nd]->chrg[nz]->DustZ > gv.bin[nd]->LowestZg ) ? Stick*dense.eden*ve*eta : 0.;
2550 
2551  /* >>chng 00 jul 13, add in gain rate from atoms and ions, PvH */
2552  *sum2 = 0.;
2553 
2554 #ifndef IGNORE_GRAIN_ION_COLLISIONS
2555  for( ion=0; ion <= LIMELM; ion++ )
2556  {
2557  double CollisionRateAll = 0.;
2558 
2559  for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
2560  {
2561  if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
2562  gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] > ion )
2563  {
2564  /* this is rate with which charged ion strikes grain */
2565  CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*GetAveVelocity( dense.AtomicWeight[nelem] )*
2566  (double)(gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]-ion);
2567  }
2568  }
2569 
2570  if( CollisionRateAll > 0. )
2571  {
2572  /* >>chng 00 jul 19, replace classical results with results
2573  * including image potential to correct for polarization
2574  * of the grain as charged particle approaches. */
2575  GrainScreen(ion,nd,nz,&eta,&xi);
2576  *sum2 += CollisionRateAll*eta;
2577  }
2578  }
2579 #endif
2580 
2581  rate = *sum1 + *sum2;
2582 
2583  /* >>chng 01 may 31, store results so that they may be used agian */
2584  gv.bin[nd]->chrg[nz]->RSum1 = *sum1;
2585  gv.bin[nd]->chrg[nz]->RSum2 = *sum2;
2586 
2587  ASSERT( *sum1 >= 0. && *sum2 >= 0. );
2588  return rate;
2589 }
2590 
2591 
2592 /* grain electron emission rates for single charge state */
2593 STATIC double GrainElecEmis1(size_t nd,
2594  long nz,
2595  /*@out@*/ double *sum1a,
2596  /*@out@*/ double *sum1b,
2597  /*@out@*/ double *sum2,
2598  /*@out@*/ double *sum3)
2599 {
2600  long int i,
2601  ion,
2602  ipLo,
2603  nelem;
2604  double eta,
2605  rate,
2606  xi;
2607 
2608  DEBUG_ENTRY( "GrainElecEmis1()" );
2609 
2610  ASSERT( nd < gv.bin.size() );
2611  ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
2612 
2613  /* >>chng 01 may 31, try to find cached results first */
2614  /* >>chng 04 feb 11, moved cached data to ChargeBin, PvH */
2615  if( gv.bin[nd]->chrg[nz]->ESum1a >= 0. )
2616  {
2617  *sum1a = gv.bin[nd]->chrg[nz]->ESum1a;
2618  *sum1b = gv.bin[nd]->chrg[nz]->ESum1b;
2619  *sum2 = gv.bin[nd]->chrg[nz]->ESum2;
2620  /* don't cache thermionic rates as they depend on grain temp */
2621  *sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;
2622  rate = *sum1a + *sum1b + *sum2 + *sum3;
2623  return rate;
2624  }
2625 
2626  /* this is the loss rate due to photo-electric effect (includes Auger and secondary electrons) */
2627  /* >>chng 01 dec 18, added code for modeling secondary electron emissions, PvH */
2628  /* this code does a crude correction for the Auger effect in grains,
2629  * it is roughly valid for neutral and negative grains, but overestimates
2630  * the effect for positively charged grains. Many of the Auger electrons have
2631  * rather low energies and will not make it out of the potential well for
2632  * high grain potentials typical of AGN conditions, see Table 4.1 & 4.2 of
2633  * >>refer grain physics Dwek E. & Smith R.K., 1996, ApJ, 459, 686 */
2634  /* >>chng 06 jan 31, this code has been completely rewritten following
2635  * >>refer grain physics Weingartner J.C., Draine B.T, Barr D.K., 2006, ApJ, 645, 1188 */
2636 
2637  *sum1a = 0.;
2638  ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
2639  for( i=ipLo; i < rfield.nflux; i++ )
2640  {
2641 # ifdef WD_TEST2
2642  *sum1a += rfield.flux[0][i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i];
2643 # else
2644  *sum1a += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i];
2645 # endif
2646  }
2647  /* normalize to rates per cm^2 of projected grain area */
2648  *sum1a /= gv.bin[nd]->IntArea/4.;
2649 
2650  *sum1b = 0.;
2651  if( gv.bin[nd]->chrg[nz]->DustZ <= -1 )
2652  {
2653  ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
2654  for( i=ipLo; i < rfield.nflux; i++ )
2655  {
2656  /* >>chng 00 jul 17, use description of Weingartner & Draine, 2001 */
2657 # ifdef WD_TEST2
2658  *sum1b += rfield.flux[0][i]*gv.bin[nd]->chrg[nz]->cs_pdt[i];
2659 # else
2660  *sum1b += rfield.SummedCon[i]*gv.bin[nd]->chrg[nz]->cs_pdt[i];
2661 # endif
2662  }
2663  *sum1b /= gv.bin[nd]->IntArea/4.;
2664  }
2665 
2666  /* >>chng 00 jun 19, add in loss rate due to recombinations with ions, PvH */
2667  *sum2 = 0.;
2668 # ifndef IGNORE_GRAIN_ION_COLLISIONS
2669  for( ion=0; ion <= LIMELM; ion++ )
2670  {
2671  double CollisionRateAll = 0.;
2672 
2673  for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
2674  {
2675  if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
2676  ion > gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] )
2677  {
2678  /* this is rate with which charged ion strikes grain */
2679  CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*GetAveVelocity( dense.AtomicWeight[nelem] )*
2680  (double)(ion-gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]);
2681  }
2682  }
2683 
2684  if( CollisionRateAll > 0. )
2685  {
2686  /* >>chng 00 jul 19, replace classical results with results
2687  * including image potential to correct for polarization
2688  * of the grain as charged particle approaches. */
2689  GrainScreen(ion,nd,nz,&eta,&xi);
2690  *sum2 += CollisionRateAll*eta;
2691  }
2692  }
2693 # endif
2694 
2695  /* >>chng 01 may 30, moved calculation of ThermRate to UpdatePot */
2696  /* >>chng 01 jan 19, multiply by 4 since thermionic emissions scale with total grain
2697  * surface area while the above two processes scale with projected grain surface area, PvH */
2698  *sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;
2699 
2702  rate = *sum1a + *sum1b + *sum2 + *sum3;
2703 
2704  /* >>chng 01 may 31, store results so that they may be used agian */
2705  gv.bin[nd]->chrg[nz]->ESum1a = *sum1a;
2706  gv.bin[nd]->chrg[nz]->ESum1b = *sum1b;
2707  gv.bin[nd]->chrg[nz]->ESum2 = *sum2;
2708 
2709  ASSERT( *sum1a >= 0. && *sum1b >= 0. && *sum2 >= 0. && *sum3 >= 0. );
2710  return rate;
2711 }
2712 
2713 
2714 /* correction factors for grain charge screening (including image potential
2715  * to correct for polarization of the grain as charged particle approaches). */
2716 STATIC void GrainScreen(long ion,
2717  size_t nd,
2718  long nz,
2719  /*@out@*/ double *eta,
2720  /*@out@*/ double *xi)
2721 {
2722 
2723  /* >>chng 04 jan 31, start caching eta, xi results, PvH */
2724 
2725  /* add 1 to allow for electron charge ion = -1 */
2726  long ind = ion+1;
2727 
2728  DEBUG_ENTRY( "GrainScreen()" );
2729 
2730  ASSERT( ind >= 0 && ind < LIMELM+2 );
2731 
2732  if( gv.bin[nd]->chrg[nz]->eta[ind] > 0. )
2733  {
2734  *eta = gv.bin[nd]->chrg[nz]->eta[ind];
2735  *xi = gv.bin[nd]->chrg[nz]->xi[ind];
2736  return;
2737  }
2738 
2739  /* >>refer grain physics Draine & Sutin, 1987, ApJ, 320, 803
2740  * eta = J-tilde (eq. 3.3 thru 3.5), xi = Lambda-tilde/2. (eq. 3.8 thru 3.10) */
2741  if( ion == 0 )
2742  {
2743  *eta = 1.;
2744  *xi = 1.;
2745  }
2746  else
2747  {
2748  /* >>chng 01 jan 03, assume that grain charge is distributed in two states just below and
2749  * above the average charge, instead of the delta function we assume elsewhere. by averaging
2750  * over the distribution we smooth out the discontinuities of the the Draine & Sutin expressions
2751  * around nu == 0. they were the cause of temperature instabilities in globule.in. PvH */
2752  /* >>chng 01 may 07, revert back to single charge state since only integral charge states are
2753  * fed into this routine now, making the two-charge state approximation obsolete, PvH */
2754  double nu = (double)gv.bin[nd]->chrg[nz]->DustZ/(double)ion;
2755  double tau = gv.bin[nd]->Capacity*BOLTZMANN*phycon.te*1.e-7/POW2((double)ion*ELEM_CHARGE);
2756  if( nu < 0. )
2757  {
2758  *eta = (1. - nu/tau)*(1. + sqrt(2./(tau - 2.*nu)));
2759  *xi = (1. - nu/(2.*tau))*(1. + 1./sqrt(tau - nu));
2760  }
2761  else if( nu == 0. )
2762  {
2763  *eta = 1. + sqrt(PI/(2.*tau));
2764  *xi = 1. + 0.75*sqrt(PI/(2.*tau));
2765  }
2766  else
2767  {
2768  double theta_nu = ThetaNu(nu);
2769  /* >>>chng 00 jul 27, avoid passing functions to macro so set to local temp var */
2770  double xxx = 1. + 1./sqrt(4.*tau+3.*nu);
2771  *eta = POW2(xxx)*exp(-theta_nu/tau);
2772 # ifdef WD_TEST2
2773  *xi = (1. + nu/(2.*tau))*(1. + 1./sqrt(3./(2.*tau)+3.*nu))*exp(-theta_nu/tau);
2774 # else
2775  /* >>chng 01 jan 24, use new expression for xi which only contains the excess
2776  * energy above the potential barrier of the incoming particle (accurate to
2777  * 2% or better), and then add in potential barrier separately, PvH */
2778  xxx = 0.25*pow(nu/tau,0.75)/(pow(nu/tau,0.75) + pow((25.+3.*nu)/5.,0.75)) +
2779  (1. + 0.75*sqrt(PI/(2.*tau)))/(1. + sqrt(PI/(2.*tau)));
2780  *xi = (MIN2(xxx,1.) + theta_nu/(2.*tau))*(*eta);
2781 # endif
2782  }
2783 
2784  ASSERT( *eta >= 0. && *xi >= 0. );
2785  }
2786 
2787  gv.bin[nd]->chrg[nz]->eta[ind] = *eta;
2788  gv.bin[nd]->chrg[nz]->xi[ind] = *xi;
2789 
2790  return;
2791 }
2792 
2793 
2794 STATIC double ThetaNu(double nu)
2795 {
2796  double theta_nu;
2797 
2798  DEBUG_ENTRY( "ThetaNu()" );
2799 
2800  if( nu > 0. )
2801  {
2802  double xxx;
2803  const double REL_TOLER = 10.*DBL_EPSILON;
2804 
2805  /* >>chng 01 jan 24, get first estimate for xi_nu and iteratively refine, PvH */
2806  double xi_nu = 1. + 1./sqrt(3.*nu);
2807  double xi_nu2 = POW2(xi_nu);
2808  do
2809  {
2810  double old = xi_nu;
2811  /* >>chng 04 feb 01, use Newton-Raphson to speed up convergence, PvH */
2812  /* xi_nu = sqrt(1.+sqrt((2.*POW2(xi_nu)-1.)/(nu*xi_nu))); */
2813  double fnu = 2.*xi_nu2 - 1. - nu*xi_nu*POW2(xi_nu2 - 1.);
2814  double dfdxi = 4.*xi_nu - nu*((5.*xi_nu2 - 6.)*xi_nu2 + 1.);
2815  xi_nu -= fnu/dfdxi;
2816  xi_nu2 = POW2(xi_nu);
2817  xxx = fabs(old-xi_nu);
2818  } while( xxx > REL_TOLER*xi_nu );
2819 
2820  theta_nu = nu/xi_nu - 1./(2.*xi_nu2*(xi_nu2-1.));
2821  }
2822  else
2823  {
2824  theta_nu = 0.;
2825  }
2826  return theta_nu;
2827 }
2828 
2829 
2830 /* update items that depend on grain potential */
2831 STATIC void UpdatePot(size_t nd,
2832  long Zlo,
2833  long stride,
2834  /*@out@*/ double rate_up[], /* rate_up[NCHU] */
2835  /*@out@*/ double rate_dn[]) /* rate_dn[NCHU] */
2836 {
2837  long nz,
2838  Zg;
2839  double BoltzFac,
2840  HighEnergy;
2841 
2842  DEBUG_ENTRY( "UpdatePot()" );
2843 
2844  ASSERT( nd < gv.bin.size() );
2845  ASSERT( Zlo >= gv.bin[nd]->LowestZg );
2846  ASSERT( stride >= 1 );
2847 
2848  /* >>chng 00 jul 17, use description of Weingartner & Draine, 2001 */
2849  /* >>chng 01 mar 21, assume that grain charge is distributed in two states just below and
2850  * above the average charge. */
2851  /* >>chng 01 may 07, this routine now completely supports the hybrid grain
2852  * charge model, and the average charge state is not used anywhere anymore, PvH */
2853  /* >>chng 01 may 30, reorganize code such that all relevant data can be copied
2854  * when a valid set of data is available from a previous call, this saves CPU time, PvH */
2855  /* >>chng 04 jan 17, reorganized code to use pointers to the charge bins, PvH */
2856 
2857  if( trace.lgTrace && trace.lgDustBug )
2858  {
2859  fprintf( ioQQQ, " %ld/%ld", Zlo, stride );
2860  }
2861 
2862  if( gv.bin[nd]->nfill < rfield.nflux )
2863  {
2864  InitBinAugerData( nd, gv.bin[nd]->nfill, rfield.nflux );
2865  gv.bin[nd]->nfill = rfield.nflux;
2866  }
2867 
2868  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2869  {
2870  long ind, zz;
2871  double d[4];
2872  ChargeBin *ptr;
2873 
2874  Zg = Zlo + nz*stride;
2875 
2876  /* check if charge state is already present */
2877  ind = NCHS-1;
2878  for( zz=0; zz < NCHS-1; zz++ )
2879  {
2880  if( gv.bin[nd]->chrg[zz]->DustZ == Zg )
2881  {
2882  ind = zz;
2883  break;
2884  }
2885  }
2886 
2887  /* in the code below the old charge bins are shifted to the back in order to assure
2888  * that the most recently used ones are upfront; they are more likely to be reused */
2889  ptr = gv.bin[nd]->chrg[ind];
2890 
2891  /* make room to swap in charge state */
2892  for( zz=ind-1; zz >= nz; zz-- )
2893  gv.bin[nd]->chrg[zz+1] = gv.bin[nd]->chrg[zz];
2894 
2895  gv.bin[nd]->chrg[nz] = ptr;
2896 
2897  if( gv.bin[nd]->chrg[nz]->DustZ != Zg )
2898  UpdatePot1(nd,nz,Zg,0);
2899  else if( gv.bin[nd]->chrg[nz]->nfill < rfield.nflux )
2900  UpdatePot1(nd,nz,Zg,gv.bin[nd]->chrg[nz]->nfill);
2901 
2902  UpdatePot2(nd,nz);
2903 
2904  rate_up[nz] = GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
2905  rate_dn[nz] = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
2906 
2907  /* sanity checks */
2908  ASSERT( gv.bin[nd]->chrg[nz]->DustZ == Zg );
2909  ASSERT( gv.bin[nd]->chrg[nz]->nfill >= rfield.nflux );
2910  ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
2911  }
2912 
2913  /* determine highest energy to be considered by quantum heating routines.
2914  * since the Boltzmann distribution is resolved, the upper limit has to be
2915  * high enough that a negligible amount of energy is in the omitted tail */
2916  /* >>chng 03 jan 26, moved this code from GrainChargeTemp to UpdatePot
2917  * since the new code depends on grain potential, HTT91.in, PvH */
2918  BoltzFac = (-log(CONSERV_TOL) + 8.)*BOLTZMANN/EN1RYD;
2919  HighEnergy = 0.;
2920  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
2921  {
2922  /* >>chng 04 jan 21, changed phycon.te -> MAX2(phycon.te,gv.bin[nd]->tedust), PvH */
2923  HighEnergy = MAX2(HighEnergy,
2924  MAX2(gv.bin[nd]->chrg[nz]->ThresInfInc,0.) + BoltzFac*MAX2(phycon.te,gv.bin[nd]->tedust));
2925  }
2926  HighEnergy = MIN2(HighEnergy,rfield.anu[rfield.nupper-1]);
2927  gv.bin[nd]->qnflux2 = ipoint(HighEnergy);
2928  gv.bin[nd]->qnflux = MAX2(rfield.nflux,gv.bin[nd]->qnflux2);
2929 
2930  ASSERT( gv.bin[nd]->qnflux <= rfield.nupper-1 );
2931  return;
2932 }
2933 
2934 
2935 /* calculate charge state populations */
2936 STATIC void GetFracPop(size_t nd,
2937  long Zlo,
2938  /*@in@*/ double rate_up[], /* rate_up[NCHU] */
2939  /*@in@*/ double rate_dn[], /* rate_dn[NCHU] */
2940  /*@out@*/ long *newZlo)
2941 {
2942  bool lgRedo;
2943  long i,
2944  nz;
2945  double netloss[2],
2946  pop[2][NCHU-1];
2947 
2948 
2949  DEBUG_ENTRY( "GetFracPop()" );
2950 
2951  ASSERT( nd < gv.bin.size() );
2952  ASSERT( Zlo >= gv.bin[nd]->LowestZg );
2953 
2954  /* solve level populations for levels 0..nChrg-2 (i == 0) and
2955  * levels 1..nChrg-1 (i == 1), and determine net electron loss rate
2956  * for each of those subsystems. Next we demand that
2957  * netloss[1] <= 0 <= netloss[0]
2958  * and determine FracPop by linearly adding the subsystems such that
2959  * 0 == frac*netloss[0] + (1-frac)*netloss[1]
2960  * this assures that all charge state populations are positive */
2961  do
2962  {
2963  for( i=0; i < 2; i++ )
2964  {
2965  long j, k;
2966  double sum;
2967 
2968  sum = pop[i][0] = 1.;
2969  for( j=1; j < gv.bin[nd]->nChrg-1; j++ )
2970  {
2971  nz = i + j;
2972  if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
2973  {
2974  pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
2975  sum += pop[i][j];
2976  }
2977  else
2978  {
2979  for( k=0; k < j; k++ )
2980  {
2981  pop[i][k] = 0.;
2982  }
2983  pop[i][j] = 1.;
2984  sum = pop[i][j];
2985  }
2986  /* guard against overflow */
2987  if( pop[i][j] > sqrt(DBL_MAX) )
2988  {
2989  for( k=0; k <= j; k++ )
2990  {
2991  pop[i][k] /= DBL_MAX/10.;
2992  }
2993  sum /= DBL_MAX/10.;
2994  }
2995  }
2996  netloss[i] = 0.;
2997  for( j=0; j < gv.bin[nd]->nChrg-1; j++ )
2998  {
2999  nz = i + j;
3000  pop[i][j] /= sum;
3001  netloss[i] += pop[i][j]*(rate_up[nz] - rate_dn[nz]);
3002  }
3003  }
3004 
3005  /* ascertain that the choice of Zlo was correct, this is to ensure positive
3006  * level populations and continuous emission and recombination rates */
3007  if( netloss[0]*netloss[1] > 0. )
3008  *newZlo = ( netloss[1] > 0. ) ? Zlo + 1 : Zlo - 1;
3009  else
3010  *newZlo = Zlo;
3011 
3012  /* >>chng 04 feb 15, add protection for roundoff error when nChrg is much too large;
3013  * netloss[0/1] can be almost zero, but may have arbitrary sign due to roundoff error;
3014  * this can lead to a spurious lowering of Zlo below LowestZg, orion_pdr10.in,
3015  * since newZlo cannot be lowered, lower nChrg instead and recalculate, PvH */
3016  /* >>chng 04 feb 15, also lower nChrg if population of highest charge state is marginal */
3017  if( gv.bin[nd]->nChrg > 2 &&
3018  ( *newZlo < gv.bin[nd]->LowestZg ||
3019  ( *newZlo == Zlo && pop[1][gv.bin[nd]->nChrg-2] < DBL_EPSILON ) ) )
3020  {
3021  gv.bin[nd]->nChrg--;
3022  lgRedo = true;
3023  }
3024  else
3025  {
3026  lgRedo = false;
3027  }
3028 
3029 # if 0
3030  printf( " fnzone %.2f nd %ld Zlo %ld newZlo %ld netloss %.4e %.4e nChrg %ld lgRedo %d\n",
3031  fnzone, nd, Zlo, *newZlo, netloss[0], netloss[1], gv.bin[nd]->nChrg, lgRedo );
3032 # endif
3033  }
3034  while( lgRedo );
3035 
3036  if( *newZlo < gv.bin[nd]->LowestZg )
3037  {
3038  fprintf( ioQQQ, " could not converge charge state populations for %s\n", gv.bin[nd]->chDstLab );
3039  ShowMe();
3040  cdEXIT(EXIT_FAILURE);
3041  }
3042 
3043  if( *newZlo == Zlo )
3044  {
3045  double frac0, frac1;
3046 # ifndef NDEBUG
3047  double test1, test2, test3, x1, x2;
3048 # endif
3049 
3050  /* frac1 == 1-frac0, but calculating it like this avoids cancellation errors */
3051  frac0 = netloss[1]/(netloss[1]-netloss[0]);
3052  frac1 = -netloss[0]/(netloss[1]-netloss[0]);
3053 
3054  gv.bin[nd]->chrg[0]->FracPop = frac0*pop[0][0];
3055  gv.bin[nd]->chrg[gv.bin[nd]->nChrg-1]->FracPop = frac1*pop[1][gv.bin[nd]->nChrg-2];
3056  for( nz=1; nz < gv.bin[nd]->nChrg-1; nz++ )
3057  {
3058  gv.bin[nd]->chrg[nz]->FracPop = frac0*pop[0][nz] + frac1*pop[1][nz-1];
3059  }
3060 
3061 # ifndef NDEBUG
3062  test1 = test2 = test3 = 0.;
3063  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
3064  {
3065  ASSERT( gv.bin[nd]->chrg[nz]->FracPop >= 0. );
3066  test1 += gv.bin[nd]->chrg[nz]->FracPop;
3067  test2 += gv.bin[nd]->chrg[nz]->FracPop*rate_up[nz];
3068  test3 += gv.bin[nd]->chrg[nz]->FracPop*(rate_up[nz]-rate_dn[nz]);
3069  }
3070  x1 = fabs(test1-1.);
3071  x2 = fabs( safe_div(test3, test2, 0.) );
3072  ASSERT( MAX2(x1,x2) < 10.*sqrt((double)gv.bin[nd]->nChrg)*DBL_EPSILON );
3073 # endif
3074  }
3075  return;
3076 }
3077 
3078 
3079 /* this routine updates all quantities that depend only on grain charge and radius
3080  *
3081  * NB NB - All global data in grain.c and grainvar.h that are charge dependent should
3082  * be calculated here or in UpdatePot2 (except gv.bin[nd]->chrg[nz]->FracPop
3083  * which is special).
3084  *
3085  * NB NB - the code assumes that the data calculated here remain valid throughout
3086  * the model, i.e. do NOT depend on grain temperature, etc.
3087  */
3088 STATIC void UpdatePot1(size_t nd,
3089  long nz,
3090  long Zg,
3091  long ipStart)
3092 {
3093  long i,
3094  ipLo,
3095  nfill;
3096  double c1,
3097  cnv_GR_pH,
3098  d[2],
3099  DustWorkFcn,
3100  Elo,
3101  Emin,
3102  ThresInf,
3103  ThresInfVal;
3104 
3105  double *anu = rfield.anu;
3106 
3107  /* constants in the expression for the photodetachment cross section
3108  * taken from
3109  * >>refer grain physics Weingartner & Draine, ApJS, 2001, 134, 263 */
3110  const double INV_DELTA_E = EVRYD/3.;
3111  const double CS_PDT = 1.2e-17;
3112 
3113  DEBUG_ENTRY( "UpdatePot1()" );
3114 
3115  /* >>chng 04 jan 23, replaced rfield.nflux with rfield.nupper so
3116  * that data remains valid throughout the model, PvH */
3117  /* >>chng 04 jan 26, start using ipStart to solve the validity problem
3118  * ipStart == 0 means that a full initialization needs to be done
3119  * ipStart > 0 means that rfield.nflux was updated and yhat(_primary), ehat,
3120  * cs_pdt, fac1, and fac2 need to be reallocated, PvH */
3121  if( ipStart == 0 )
3122  {
3123  gv.bin[nd]->chrg[nz]->DustZ = Zg;
3124 
3125  /* invalidate eta and xi storage */
3126  memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
3127  memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );
3128 
3129  GetPotValues(nd,Zg,&gv.bin[nd]->chrg[nz]->ThresInf,&gv.bin[nd]->chrg[nz]->ThresInfVal,
3130  &gv.bin[nd]->chrg[nz]->ThresSurf,&gv.bin[nd]->chrg[nz]->ThresSurfVal,
3131  &gv.bin[nd]->chrg[nz]->PotSurf,&gv.bin[nd]->chrg[nz]->Emin,INCL_TUNNEL);
3132 
3133  /* >>chng 01 may 09, do not use tunneling corrections for incoming electrons */
3134  /* >>chng 01 nov 25, add gv.bin[nd]->chrg[nz]->ThresInfInc, PvH */
3135  GetPotValues(nd,Zg-1,&gv.bin[nd]->chrg[nz]->ThresInfInc,&d[0],&gv.bin[nd]->chrg[nz]->ThresSurfInc,
3136  &d[1],&gv.bin[nd]->chrg[nz]->PotSurfInc,&gv.bin[nd]->chrg[nz]->EminInc,NO_TUNNEL);
3137 
3138  gv.bin[nd]->chrg[nz]->ipThresInfVal =
3139  hunt_bisect( rfield.anu, rfield.nupper, gv.bin[nd]->chrg[nz]->ThresInfVal ) + 1;
3140  }
3141 
3142  ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
3143 
3144  /* remember how far the yhat(_primary), ehat, cs_pdt, fac1, and fac2 arrays were filled in */
3145  gv.bin[nd]->chrg[nz]->nfill = rfield.nflux;
3146  nfill = gv.bin[nd]->chrg[nz]->nfill;
3147 
3148  /* >>chng 04 feb 07, only allocate arrays from ipLo to nfill to save memory, PvH */
3149  long len = nfill + 10 - ipLo;
3150  if( ipStart == 0 )
3151  {
3152  gv.bin[nd]->chrg[nz]->yhat.reserve( len );
3153  gv.bin[nd]->chrg[nz]->yhat_primary.reserve( len );
3154  gv.bin[nd]->chrg[nz]->ehat.reserve( len );
3155  gv.bin[nd]->chrg[nz]->yhat.alloc( ipLo, nfill );
3156  gv.bin[nd]->chrg[nz]->yhat_primary.alloc( ipLo, nfill );
3157  gv.bin[nd]->chrg[nz]->ehat.alloc( ipLo, nfill );
3158  }
3159  else
3160  {
3161  gv.bin[nd]->chrg[nz]->yhat.realloc( nfill );
3162  gv.bin[nd]->chrg[nz]->yhat_primary.realloc( nfill );
3163  gv.bin[nd]->chrg[nz]->ehat.realloc( nfill );
3164  }
3165 
3166  double GrainPot = chrg2pot(Zg,nd);
3167 
3168  if( nfill > ipLo )
3169  {
3170  DustWorkFcn = gv.bin[nd]->DustWorkFcn;
3171  Elo = -gv.bin[nd]->chrg[nz]->PotSurf;
3172  ThresInfVal = gv.bin[nd]->chrg[nz]->ThresInfVal;
3173  Emin = gv.bin[nd]->chrg[nz]->Emin;
3174 
3178  ASSERT( gv.bin[nd]->sd[0]->ipLo <= ipLo );
3179 
3180  for( i=max(ipLo,ipStart); i < nfill; i++ )
3181  {
3182  long n;
3183  unsigned int ns=0;
3184  double Yp1,Ys1,Ehp,Ehs,yp,ya,ys,eyp,eya,eys;
3185  double yzero,yone,Ehi,Ehi_band,Wcorr,Eel;
3186  ShellData *sptr = gv.bin[nd]->sd[ns];
3187 
3188  yp = ya = ys = 0.;
3189  eyp = eya = eys = 0.;
3190 
3191  /* calculate yield for band structure */
3192  Ehi = Ehi_band = anu[i] - ThresInfVal - Emin;
3193  Wcorr = ThresInfVal + Emin - GrainPot;
3194  Eel = anu[i] - DustWorkFcn;
3195  yzero = y0b( nd, nz, i );
3196  yone = sptr->y01[i];
3197  Yfunc(nd,nz,yzero*yone,sptr->p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3198  yp += Yp1;
3199  ys += Ys1;
3200  eyp += Yp1*Ehp;
3201  eys += Ys1*Ehs;
3202 
3203  /* add in yields for inner shell ionization */
3204  unsigned int nsmax = gv.bin[nd]->sd.size();
3205  for( ns=1; ns < nsmax; ns++ )
3206  {
3207  sptr = gv.bin[nd]->sd[ns];
3208 
3209  if( i < sptr->ipLo )
3210  continue;
3211 
3212  Ehi = Ehi_band + Wcorr - sptr->ionPot;
3213  Eel = anu[i] - sptr->ionPot;
3214  Yfunc(nd,nz,sptr->y01[i],sptr->p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3215  yp += Yp1;
3216  ys += Ys1;
3217  eyp += Yp1*Ehp;
3218  eys += Ys1*Ehs;
3219 
3220  /* add in Auger yields */
3221  long nmax = sptr->nData;
3222  for( n=0; n < nmax; n++ )
3223  {
3224  double max = sptr->AvNr[n]*sptr->p[i];
3225  Ehi = sptr->Ener[n] - GrainPot;
3226  Eel = sptr->Ener[n];
3227  Yfunc(nd,nz,sptr->y01A[n][i],max,Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3228  ya += Yp1;
3229  ys += Ys1;
3230  eya += Yp1*Ehp;
3231  eys += Ys1*Ehs;
3232  }
3233  }
3234  /* add up primary, Auger, and secondary yields */
3235  gv.bin[nd]->chrg[nz]->yhat[i] = (realnum)(yp + ya + ys);
3236  gv.bin[nd]->chrg[nz]->yhat_primary[i] = min((realnum)yp,1.f);
3237  gv.bin[nd]->chrg[nz]->ehat[i] = ( gv.bin[nd]->chrg[nz]->yhat[i] > 0. ) ?
3238  (realnum)((eyp + eya + eys)/gv.bin[nd]->chrg[nz]->yhat[i]) : 0.f;
3239 
3240  ASSERT( yp <= 1.00001 );
3241  ASSERT( gv.bin[nd]->chrg[nz]->ehat[i] >= 0.f );
3242  }
3243  }
3244 
3245  if( ipStart == 0 )
3246  {
3247  /* >>chng 01 jan 08, ThresInf[nz] and ThresInfVal[nz] may become zero in
3248  * initial stages of grain potential search, PvH */
3249  /* >>chng 01 oct 10, use bisection search to find ipThresInf, ipThresInfVal. On C scale now */
3250  gv.bin[nd]->chrg[nz]->ipThresInf =
3251  hunt_bisect( rfield.anu, rfield.nupper, gv.bin[nd]->chrg[nz]->ThresInf ) + 1;
3252  }
3253 
3254  ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
3255 
3256  len = nfill + 10 - ipLo;
3257 
3258  if( Zg <= -1 )
3259  {
3260  /* >>chng 04 feb 07, only allocate arrays from ipLo to nfill to save memory, PvH */
3261  if( ipStart == 0 )
3262  {
3263  gv.bin[nd]->chrg[nz]->cs_pdt.reserve( len );
3264  gv.bin[nd]->chrg[nz]->cs_pdt.alloc( ipLo, nfill );
3265  }
3266  else
3267  {
3268  gv.bin[nd]->chrg[nz]->cs_pdt.realloc( nfill );
3269  }
3270 
3271  if( nfill > ipLo )
3272  {
3273  c1 = -CS_PDT*(double)Zg;
3274  ThresInf = gv.bin[nd]->chrg[nz]->ThresInf;
3275  cnv_GR_pH = gv.bin[nd]->cnv_GR_pH;
3276 
3277  for( i=max(ipLo,ipStart); i < nfill; i++ )
3278  {
3279  double x = (anu[i] - ThresInf)*INV_DELTA_E;
3280  double cs = c1*x/POW2(1.+(1./3.)*POW2(x));
3281  /* this cross section must be at default depletion for consistency
3282  * with dstab1, hence no factor dstAbund should be applied here */
3283  gv.bin[nd]->chrg[nz]->cs_pdt[i] = MAX2(cs,0.)*cnv_GR_pH;
3284  }
3285  }
3286  }
3287 
3288  /* >>chng 04 feb 07, added fac1, fac2 to optimize loop for photoelectric heating, PvH */
3289  if( ipStart == 0 )
3290  {
3291  gv.bin[nd]->chrg[nz]->fac1.reserve( len );
3292  gv.bin[nd]->chrg[nz]->fac2.reserve( len );
3293  gv.bin[nd]->chrg[nz]->fac1.alloc( ipLo, nfill );
3294  gv.bin[nd]->chrg[nz]->fac2.alloc( ipLo, nfill );
3295  }
3296  else
3297  {
3298  gv.bin[nd]->chrg[nz]->fac1.realloc( nfill );
3299  gv.bin[nd]->chrg[nz]->fac2.realloc( nfill );
3300  }
3301 
3302  if( nfill > ipLo )
3303  {
3304  for( i=max(ipLo,ipStart); i < nfill; i++ )
3305  {
3306  double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;
3307 
3308  /* >>chng 04 jan 25, created separate routine PE_init, PvH */
3309  PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);
3310 
3311  gv.bin[nd]->chrg[nz]->fac1[i] = cs_tot*anu[i]-cs1*cool1-cs2*cool2;
3312  gv.bin[nd]->chrg[nz]->fac2[i] = cs1*ehat1 + cs2*ehat2;
3313 
3314  ASSERT( gv.bin[nd]->chrg[nz]->fac1[i] >= 0. && gv.bin[nd]->chrg[nz]->fac2[i] >= 0. );
3315  }
3316  }
3317 
3318  if( ipStart == 0 )
3319  {
3320  /* >>chng 00 jul 05, determine ionization stage Z0 the ion recombines to */
3321  /* >>chng 04 jan 20, use ALL_STAGES here so that result remains valid throughout the model */
3322  UpdateRecomZ0(nd,nz,ALL_STAGES);
3323  }
3324 
3325  /* invalidate the remaining fields */
3326  gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
3327 
3328  gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
3329  gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
3330  gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
3331  gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
3332  gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
3333 
3334  gv.bin[nd]->chrg[nz]->tedust = 1.f;
3335 
3336  gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
3337  gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
3338  gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
3339  gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
3340 
3341  gv.bin[nd]->chrg[nz]->BolFlux = -DBL_MAX;
3342  gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
3343  gv.bin[nd]->chrg[nz]->GrainHeatColl = -DBL_MAX;
3344  gv.bin[nd]->chrg[nz]->GasHeatPhotoEl = -DBL_MAX;
3345  gv.bin[nd]->chrg[nz]->GasHeatTherm = -DBL_MAX;
3346  gv.bin[nd]->chrg[nz]->GrainCoolTherm = -DBL_MAX;
3347  gv.bin[nd]->chrg[nz]->ChemEnIon = -DBL_MAX;
3348  gv.bin[nd]->chrg[nz]->ChemEnH2 = -DBL_MAX;
3349 
3350  gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
3351 
3352  /* sanity check */
3353  ASSERT( gv.bin[nd]->chrg[nz]->ipThresInf <= gv.bin[nd]->chrg[nz]->ipThresInfVal );
3354  return;
3355 }
3356 
3357 
3358 /* this routine updates all quantities that depend on grain charge, radius and temperature
3359  *
3360  * NB NB - All global data in grain.c and grainvar.h that are charge dependent should
3361  * be calculated here or in UpdatePot1 (except gv.bin[nd]->chrg[nz]->FracPop
3362  * which is special).
3363  *
3364  * NB NB - the code assumes that the data calculated here may vary throughout the model,
3365  * e.g. because of a dependence on grain temperature
3366  */
3367 STATIC void UpdatePot2(size_t nd,
3368  long nz)
3369 {
3370  double ThermExp;
3371 
3372  DEBUG_ENTRY( "UpdatePot2()" );
3373 
3374  /* >>chng 00 jun 19, add in loss rate due to thermionic emission of electrons, PvH */
3375  ThermExp = gv.bin[nd]->chrg[nz]->ThresInf*TE1RYD/gv.bin[nd]->tedust;
3376  /* ThermExp is guaranteed to be >= 0. */
3377  gv.bin[nd]->chrg[nz]->ThermRate = THERMCONST*gv.bin[nd]->ThermEff*POW2(gv.bin[nd]->tedust)*exp(-ThermExp);
3378 # if defined( WD_TEST2 ) || defined( IGNORE_THERMIONIC )
3379  gv.bin[nd]->chrg[nz]->ThermRate = 0.;
3380 # endif
3381  return;
3382 }
3383 
3384 
3385 /* Helper function to calculate primary and secondary yields and the average electron energy at infinity */
3386 inline void Yfunc(long nd,
3387  long nz,
3388  double y01,
3389  double maxval,
3390  double Elo,
3391  double Ehi,
3392  double Eel,
3393  /*@out@*/ double *Yp,
3394  /*@out@*/ double *Ys,
3395  /*@out@*/ double *Ehp,
3396  /*@out@*/ double *Ehs)
3397 {
3398  DEBUG_ENTRY( "Yfunc()" );
3399 
3400  long Zg = gv.bin[nd]->chrg[nz]->DustZ;
3401  double y2pr, y2sec;
3402 
3403  ASSERT( Ehi >= Elo );
3404 
3405  y2pr = y2pa( Elo, Ehi, Zg, Ehp );
3406 
3407  if( y2pr > 0. )
3408  {
3409  pe_type pcase = gv.which_pe[gv.bin[nd]->matType];
3410  double eps, f3;
3411 
3412  *Yp = y2pr*min(y01,maxval);
3413 
3414  y2sec = y2s( Elo, Ehi, Zg, Ehs );
3415  if( pcase == PE_CAR )
3416  eps = 117./EVRYD;
3417  else if( pcase == PE_SIL )
3418  eps = 155./EVRYD;
3419  else
3420  {
3421  fprintf( ioQQQ, " Yfunc: unknown type for PE effect: %d\n" , pcase );
3422  cdEXIT(EXIT_FAILURE);
3423  }
3424  /* this is Eq. 18 of WDB06 */
3425  /* Eel may be negative near threshold -> set yield to zero */
3426  f3 = max(Eel,0.)/(eps*elec_esc_length(Eel,nd)*gv.bin[nd]->eyc);
3427  *Ys = y2sec*f3*min(y01,maxval);
3428  }
3429  else
3430  {
3431  *Yp = 0.;
3432  *Ys = 0.;
3433  *Ehp = 0.;
3434  *Ehs = 0.;
3435  }
3436  return;
3437 }
3438 
3439 
3440 /* This calculates the y0 function for band electrons (Sect. 4.1.3/4.1.4 of WDB06) */
3441 STATIC double y0b(size_t nd,
3442  long nz,
3443  long i) /* incident photon energy is anu[i] */
3444 {
3445  double yzero;
3446 
3447  DEBUG_ENTRY( "y0b()" );
3448 
3449  if( gv.lgWD01 )
3450  yzero = y0b01( nd, nz, i );
3451  else
3452  {
3453  double Eph = rfield.anu[i];
3454 
3455  if( Eph <= 20./EVRYD )
3456  yzero = y0b01( nd, nz, i );
3457  else if( Eph < 50./EVRYD )
3458  {
3459  double y0a = y0b01( nd, nz, i );
3460  double y0b = gv.bin[nd]->y0b06[i];
3461  /* constant 1.09135666... is 1./log(50./20.) */
3462  double frac = log(Eph*(EVRYD/20.))*1.0913566679372915;
3463 
3464  yzero = y0a * exp(log(y0b/y0a)*frac);
3465  }
3466  else
3467  yzero = gv.bin[nd]->y0b06[i];
3468  }
3469 
3470  ASSERT( yzero > 0. );
3471  return yzero;
3472 }
3473 
3474 
3475 /* This calculates the y0 function for band electrons (Eq. 16 of WD01) */
3476 STATIC double y0b01(size_t nd,
3477  long nz,
3478  long i) /* incident photon energy is anu[i] */
3479 {
3480  pe_type pcase = gv.which_pe[gv.bin[nd]->matType];
3481  double xv, yzero;
3482 
3483  DEBUG_ENTRY( "y0b01()" );
3484 
3485  xv = MAX2((rfield.anu[i] - gv.bin[nd]->chrg[nz]->ThresSurfVal)/gv.bin[nd]->DustWorkFcn,0.);
3486 
3487  switch( pcase )
3488  {
3489  case PE_CAR:
3490  /* >>refer grain physics Bakes & Tielens, 1994, ApJ, 427, 822 */
3491  xv = POW2(xv)*POW3(xv);
3492  yzero = xv/((1./9.e-3) + (3.7e-2/9.e-3)*xv);
3493  break;
3494  case PE_SIL:
3495  /* >>refer grain physics Weingartner & Draine, 2001 */
3496  yzero = xv/(2.+10.*xv);
3497  break;
3498  default:
3499  fprintf( ioQQQ, " y0b01: unknown type for PE effect: %d\n" , pcase );
3500  cdEXIT(EXIT_FAILURE);
3501  }
3502 
3503  ASSERT( yzero > 0. );
3504  return yzero;
3505 }
3506 
3507 
3508 /* This calculates the y0 function for primary/secondary and Auger electrons (Eq. 9 of WDB06) */
3509 STATIC double y0psa(size_t nd,
3510  long ns, /* shell number */
3511  long i, /* incident photon energy is anu[i] */
3512  double Eel) /* emitted electron energy */
3513 {
3514  double yzero, leola;
3515 
3516  DEBUG_ENTRY( "y0psa()" );
3517 
3518  ASSERT( i >= gv.bin[nd]->sd[ns]->ipLo );
3519 
3520  /* this is l_e/l_a */
3521  leola = elec_esc_length(Eel,nd)*gv.bin[nd]->inv_att_len[i];
3522 
3523  ASSERT( leola > 0. );
3524 
3525  /* this is Eq. 9 of WDB06 */
3526  if( leola < 1.e4 )
3527  yzero = gv.bin[nd]->sd[ns]->p[i]*leola*(1. - leola*log(1.+1./leola));
3528  else
3529  {
3530  double x = 1./leola;
3531  yzero = gv.bin[nd]->sd[ns]->p[i]*(((-1./5.*x+1./4.)*x-1./3.)*x+1./2.);
3532  }
3533 
3534  ASSERT( yzero > 0. );
3535  return yzero;
3536 }
3537 
3538 
3539 /* This calculates the y1 function for primary/secondary and Auger electrons (Eq. 6 of WDB06) */
3540 STATIC double y1psa(size_t nd,
3541  long i, /* incident photon energy is anu[i] */
3542  double Eel) /* emitted electron energy */
3543 {
3544  double alpha, beta, af, bf, yone;
3545 
3546  DEBUG_ENTRY( "y1psa()" );
3547 
3548  beta = gv.bin[nd]->AvRadius*gv.bin[nd]->inv_att_len[i];
3549  if( beta > 1.e-4 )
3550  bf = pow2(beta) - 2.*beta + 2. - 2.*exp(-beta);
3551  else
3552  bf = ((1./60.*beta - 1./12.)*beta + 1./3.)*pow3(beta);
3553 
3554  alpha = beta + gv.bin[nd]->AvRadius/elec_esc_length(Eel,nd);
3555  if( alpha > 1.e-4 )
3556  af = pow2(alpha) - 2.*alpha + 2. - 2.*exp(-alpha);
3557  else
3558  af = ((1./60.*alpha - 1./12.)*alpha + 1./3.)*pow3(alpha);
3559 
3560  yone = pow2(beta/alpha)*af/bf;
3561 
3562  ASSERT( yone > 0. );
3563  return yone;
3564 }
3565 
3566 
3567 /* This calculates the y2 function for primary and Auger electrons (Eq. 8 of WDB06) */
3568 inline double y2pa(double Elo,
3569  double Ehi,
3570  long Zg,
3571  double *Ehp)
3572 {
3573  DEBUG_ENTRY( "y2pa()" );
3574 
3575  double ytwo;
3576 
3577  if( Zg > -1 )
3578  {
3579  if( Ehi > 0. )
3580  {
3581  double x = Elo/Ehi;
3582  *Ehp = 0.5*Ehi*(1.-2.*x)/(1.-3.*x);
3583  // use Taylor expansion for small arguments to avoid blowing assert
3584  ytwo = ( abs(x) > 1e-4 ) ? (1.-3.*x)/pow3(1.-x) : 1. - (3. + 8.*x)*x*x;
3585  ASSERT( *Ehp > 0. && *Ehp <= Ehi && ytwo > 0. && ytwo <= 1. );
3586  }
3587  else
3588  {
3589  *Ehp = 0.;
3590  ytwo = 0.;
3591  }
3592  }
3593  else
3594  {
3595  if( Ehi > Elo )
3596  {
3597  *Ehp = 0.5*(Elo+Ehi);
3598  ytwo = 1.;
3599  ASSERT( *Ehp >= Elo && *Ehp <= Ehi );
3600  }
3601  else
3602  {
3603  *Ehp = 0.;
3604  ytwo = 0.;
3605  }
3606  }
3607  return ytwo;
3608 }
3609 
3610 
3611 /* This calculates the y2 function for secondary electrons (Eqs. 20-21 of WDB06) */
3612 inline double y2s(double Elo,
3613  double Ehi,
3614  long Zg,
3615  double *Ehs)
3616 {
3617  DEBUG_ENTRY( "y2s()" );
3618 
3619  double ytwo;
3620 
3621  if( Zg > -1 )
3622  {
3623  if( !gv.lgWD01 && Ehi > 0. )
3624  {
3625  double yl = Elo/ETILDE;
3626  double yh = Ehi/ETILDE;
3627  double x = yh - yl;
3628  double E0, N0;
3629  if( x < 0.01 )
3630  {
3631  // use series expansions to avoid cancellation error
3632  double x2 = x*x, x3 = x2*x, x4 = x3*x, x5 = x4*x;
3633  double yh2 = yh*yh, yh3 = yh2*yh, yh4 = yh3*yh, yh5 = yh4*yh;
3634  double help1 = 2.*x-yh;
3635  double help2 = (6.*x3-15.*yh*x2+12.*yh2*x-3.*yh3)/4.;
3636  double help3 = (22.*x5-95.*yh*x4+164.*yh2*x3-141.*yh3*x2+60.*yh4*x-10.*yh5)/16.;
3637  N0 = yh*(help1 - help2 + help3)/x2;
3638 
3639  help1 = (3.*x-yh)/3.;
3640  help2 = (15.*x3-25.*yh*x2+15.*yh2*x-3.*yh3)/20.;
3641  help3 = (1155.*x5-3325.*yh*x4+4305.*yh2*x3-2961.*yh3*x2+1050.*yh4*x-150.*yh5)/1680.;
3642  E0 = ETILDE*yh2*(help1 - help2 + help3)/x2;
3643  }
3644  else
3645  {
3646  double sR0 = (1. + yl*yl);
3647  double sqR0 = sqrt(sR0);
3648  double sqRh = sqrt(1. + x*x);
3649  double alpha = sqRh/(sqRh - 1.);
3650  if( yh/sqR0 < 0.01 )
3651  {
3652  // use series expansions to avoid cancellation error
3653  double z = yh*(yh - 2.*yl)/sR0;
3654  N0 = ((((7./256.*z-5./128.)*z+1./16.)*z-1./8.)*z+1./2.)*z/(sqRh-1.);
3655 
3656  double yl2 = yl*yl, yl3 = yl2*yl, yl4 = yl3*yl;
3657  double help1 = yl/2.;
3658  double help2 = (2.*yl2-1.)/3.;
3659  double help3 = (6.*yl3-9.*yl)/8.;
3660  double help4 = (8.*yl4-24.*yl2+3.)/10.;
3661  double h = yh/sR0;
3662  E0 = -alpha*Ehi*(((help4*h + help3)*h + help2)*h + help1)*h/sqR0;
3663  }
3664  else
3665  {
3666  N0 = alpha*(1./sqR0 - 1./sqRh);
3667  E0 = alpha*ETILDE*(ASINH(x*sqR0 + yl*sqRh) - yh/sqRh);
3668  }
3669  }
3670  ASSERT( N0 > 0. && N0 <= 1. );
3671 
3672  *Ehs = E0/N0;
3673 
3674  ASSERT( *Ehs > 0. && *Ehs <= Ehi );
3675 
3676  ytwo = N0;
3677  }
3678  else
3679  {
3680  *Ehs = 0.;
3681  ytwo = 0.;
3682  }
3683  }
3684  else
3685  {
3686  if( !gv.lgWD01 && Ehi > Elo )
3687  {
3688  double yl = Elo/ETILDE;
3689  double yh = Ehi/ETILDE;
3690  double x = yh - yl;
3691  double x2 = x*x;
3692  if( x > 0.025 )
3693  {
3694  double sqRh = sqrt(1. + x2);
3695  double alpha = sqRh/(sqRh - 1.);
3696  *Ehs = alpha*ETILDE*(ASINH(x) - yh/sqRh + yl);
3697  }
3698  else
3699  {
3700  // use series expansion to avoid cancellation error
3701  *Ehs = Ehi - (Ehi-Elo)*((-37./840.*x2 + 1./10.)*x2 + 1./3.);
3702  }
3703 
3704  ASSERT( *Ehs >= Elo && *Ehs <= Ehi );
3705 
3706  ytwo = 1.;
3707  }
3708  else
3709  {
3710  *Ehs = 0.;
3711  ytwo = 0.;
3712  }
3713  }
3714  return ytwo;
3715 }
3716 
3717 
3718 /* find highest ionization stage with non-zero population */
3719 STATIC long HighestIonStage(void)
3720 {
3721  long high,
3722  ion,
3723  nelem;
3724 
3725  DEBUG_ENTRY( "HighestIonStage()" );
3726 
3727  high = 0;
3728  for( nelem=LIMELM-1; nelem >= 0; nelem-- )
3729  {
3730  if( dense.lgElmtOn[nelem] )
3731  {
3732  for( ion=nelem+1; ion >= 0; ion-- )
3733  {
3734  if( ion == high || dense.xIonDense[nelem][ion] > 0. )
3735  break;
3736  }
3737  high = MAX2(high,ion);
3738  }
3739  if( nelem <= high )
3740  break;
3741  }
3742  return high;
3743 }
3744 
3745 
3746 STATIC void UpdateRecomZ0(size_t nd,
3747  long nz,
3748  bool lgAllIonStages)
3749 {
3750  long hi_ion,
3751  i,
3752  ion,
3753  nelem,
3754  Zg;
3755  double d[5],
3756  phi_s_up[LIMELM+1],
3757  phi_s_dn[2];
3758 
3759  DEBUG_ENTRY( "UpdateRecomZ0()" );
3760 
3761  Zg = gv.bin[nd]->chrg[nz]->DustZ;
3762 
3763  hi_ion = ( lgAllIonStages ) ? LIMELM : gv.HighestIon;
3764 
3765  phi_s_up[0] = gv.bin[nd]->chrg[nz]->ThresSurf;
3766  for( i=1; i <= LIMELM; i++ )
3767  {
3768  if( i <= hi_ion )
3769  GetPotValues(nd,Zg+i,&d[0],&d[1],&phi_s_up[i],&d[2],&d[3],&d[4],INCL_TUNNEL);
3770  else
3771  phi_s_up[i] = -DBL_MAX;
3772  }
3773  phi_s_dn[0] = gv.bin[nd]->chrg[nz]->ThresSurfInc;
3774  GetPotValues(nd,Zg-2,&d[0],&d[1],&phi_s_dn[1],&d[2],&d[3],&d[4],NO_TUNNEL);
3775 
3776  /* >>chng 01 may 09, use GrainIonColl which properly tracks step-by-step charge changes */
3777  for( nelem=0; nelem < LIMELM; nelem++ )
3778  {
3779  if( dense.lgElmtOn[nelem] )
3780  {
3781  for( ion=0; ion <= nelem+1; ion++ )
3782  {
3783  if( lgAllIonStages || dense.xIonDense[nelem][ion] > 0. )
3784  {
3785  GrainIonColl(nd,nz,nelem,ion,phi_s_up,phi_s_dn,
3786  &gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion],
3787  &gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion],
3788  &gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion]);
3789  }
3790  else
3791  {
3792  gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] = ion;
3793  gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion] = 0.f;
3794  gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion] = 0.f;
3795  }
3796  }
3797  }
3798  }
3799  return;
3800 }
3801 
3802 STATIC void GetPotValues(size_t nd,
3803  long Zg,
3804  /*@out@*/ double *ThresInf,
3805  /*@out@*/ double *ThresInfVal,
3806  /*@out@*/ double *ThresSurf,
3807  /*@out@*/ double *ThresSurfVal,
3808  /*@out@*/ double *PotSurf,
3809  /*@out@*/ double *Emin,
3810  bool lgUseTunnelCorr)
3811 {
3812  double dstpot,
3813  dZg = (double)Zg,
3814  IP_v;
3815 
3816  DEBUG_ENTRY( "GetPotValues()" );
3817 
3818  /* >>chng 01 may 07, this routine now completely supports the hybrid grain charge model,
3819  * the change for this routine is that now it is only fed integral charge states; calculation
3820  * of IP has also been changed in accordance with Weingartner & Draine, 2001, PvH */
3821 
3822  /* this is average grain potential in Rydberg */
3823  dstpot = chrg2pot(dZg,nd);
3824 
3825  /* >>chng 01 mar 20, include O(a^-2) correction terms in ionization potential */
3826  /* these are correction terms for the ionization potential that are
3827  * important for small grains. See Weingartner & Draine, 2001, Eq. 2 */
3828  IP_v = gv.bin[nd]->DustWorkFcn + dstpot - 0.5*one_elec(nd) + (dZg+2.)*AC0/gv.bin[nd]->AvRadius*one_elec(nd);
3829 
3830  /* >>chng 01 mar 01, use new expresssion for ThresInfVal, ThresSurfVal following the discussion
3831  * with Joe Weingartner. Also include the Schottky effect (see
3832  * >>refer grain physics Spitzer, 1948, ApJ, 107, 6,
3833  * >>refer grain physics Draine & Sutin, 1987, ApJ, 320, 803), PvH */
3834  if( Zg <= -1 )
3835  {
3836  pot_type pcase = gv.which_pot[gv.bin[nd]->matType];
3837  double IP;
3838 
3839  IP = gv.bin[nd]->DustWorkFcn - gv.bin[nd]->BandGap + dstpot - 0.5*one_elec(nd);
3840  switch( pcase )
3841  {
3842  case POT_CAR:
3843  IP -= AC1G/(gv.bin[nd]->AvRadius+AC2G)*one_elec(nd);
3844  break;
3845  case POT_SIL:
3846  /* do nothing */
3847  break;
3848  default:
3849  fprintf( ioQQQ, " GetPotValues detected unknown type for ionization pot: %d\n" , pcase );
3850  cdEXIT(EXIT_FAILURE);
3851  }
3852 
3853  /* prevent valence electron from becoming less bound than attached electron; this
3854  * can happen for very negative, large graphitic grains and is not physical, PvH */
3855  IP_v = MAX2(IP,IP_v);
3856 
3857  if( Zg < -1 )
3858  {
3859  /* >>chng 01 apr 20, use improved expression for tunneling effect, PvH */
3860  double help = fabs(dZg+1);
3861  /* this is the barrier height solely due to the Schottky effect */
3862  *Emin = -ThetaNu(help)*one_elec(nd);
3863  if( lgUseTunnelCorr )
3864  {
3865  /* this is the barrier height corrected for tunneling effects */
3866  *Emin *= 1. - 2.124e-4/(pow(gv.bin[nd]->AvRadius,(realnum)0.45)*pow(help,0.26));
3867  }
3868  }
3869  else
3870  {
3871  *Emin = 0.;
3872  }
3873 
3874  *ThresInf = IP - *Emin;
3875  *ThresInfVal = IP_v - *Emin;
3876  *ThresSurf = *ThresInf;
3877  *ThresSurfVal = *ThresInfVal;
3878  *PotSurf = *Emin;
3879  }
3880  else
3881  {
3882  *ThresInf = IP_v;
3883  *ThresInfVal = IP_v;
3884  *ThresSurf = *ThresInf - dstpot;
3885  *ThresSurfVal = *ThresInfVal - dstpot;
3886  *PotSurf = dstpot;
3887  *Emin = 0.;
3888  }
3889  return;
3890 }
3891 
3892 
3893 /* given grain nd in charge state nz, and incoming ion (nelem,ion),
3894  * detemine outgoing ion (nelem,Z0) and chemical energy ChEn released
3895  * ChemEn is net contribution of ion recombination to grain heating */
3896 STATIC void GrainIonColl(size_t nd,
3897  long int nz,
3898  long int nelem,
3899  long int ion,
3900  const double phi_s_up[], /* phi_s_up[LIMELM+1] */
3901  const double phi_s_dn[], /* phi_s_dn[2] */
3902  /*@out@*/long *Z0,
3903  /*@out@*/realnum *ChEn,
3904  /*@out@*/realnum *ChemEn)
3905 {
3906  long Zg;
3907  double d[5];
3908  double phi_s;
3909 
3910  long save = ion;
3911 
3912  DEBUG_ENTRY( "GrainIonColl()" );
3913  if( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (realnum)phi_s_up[0] )
3914  {
3915  /* ion will get electron(s) */
3916  *ChEn = 0.f;
3917  *ChemEn = 0.f;
3918  Zg = gv.bin[nd]->chrg[nz]->DustZ;
3919  phi_s = phi_s_up[0];
3920  do
3921  {
3922  *ChEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] - (realnum)phi_s;
3923  *ChemEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1];
3924  /* this is a correction for the imperfections in the n-charge state model:
3925  * since the transfer gets modeled as n single-electron transfers, instead of one
3926  * n-electron transfer, a correction for the difference in binding energy is needed */
3927  *ChemEn -= (realnum)(phi_s - phi_s_up[0]);
3928  --ion;
3929  ++Zg;
3930  phi_s = phi_s_up[save-ion];
3931  } while( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (realnum)phi_s );
3932 
3933  *Z0 = ion;
3934  }
3935  else if( ion <= nelem && gv.bin[nd]->chrg[nz]->DustZ > gv.bin[nd]->LowestZg &&
3936  rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (realnum)phi_s_dn[0] )
3937  {
3938  /* grain will get electron(s) */
3939  *ChEn = 0.f;
3940  *ChemEn = 0.f;
3941  Zg = gv.bin[nd]->chrg[nz]->DustZ;
3942  phi_s = phi_s_dn[0];
3943  do
3944  {
3945  *ChEn += (realnum)phi_s - rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
3946  *ChemEn -= rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
3947  /* this is a correction for the imperfections in the n-charge state model:
3948  * since the transfer gets modeled as n single-electron transfers, instead of one
3949  * n-electron transfer, a correction for the difference in binding energy is needed */
3950  *ChemEn += (realnum)(phi_s - phi_s_dn[0]);
3951  ++ion;
3952  --Zg;
3953 
3954  if( ion-save < 2 )
3955  phi_s = phi_s_dn[ion-save];
3956  else
3957  GetPotValues(nd,Zg-1,&d[0],&d[1],&phi_s,&d[2],&d[3],&d[4],NO_TUNNEL);
3958 
3959  } while( ion <= nelem && Zg > gv.bin[nd]->LowestZg &&
3960  rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (realnum)phi_s );
3961  *Z0 = ion;
3962  }
3963  else
3964  {
3965  /* nothing happens */
3966  *ChEn = 0.f;
3967  *ChemEn = 0.f;
3968  *Z0 = ion;
3969  }
3970 /* printf(" GrainIonColl: nelem %ld ion %ld -> %ld, ChEn %.6f\n",nelem,save,*Z0,*ChEn); */
3971  return;
3972 }
3973 
3974 
3975 /* initialize grain-ion charge transfer rates on grain species nd */
3976 STATIC void GrainChrgTransferRates(long nd)
3977 {
3978  long nz;
3979  double fac0 = STICK_ION*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
3980 
3981  DEBUG_ENTRY( "GrainChrgTransferRates()" );
3982 
3983 # ifndef IGNORE_GRAIN_ION_COLLISIONS
3984 
3985  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
3986  {
3987  long ion;
3988  ChargeBin *gptr = gv.bin[nd]->chrg[nz];
3989  double fac1 = gptr->FracPop*fac0;
3990 
3991  if( fac1 == 0. )
3992  continue;
3993 
3994  for( ion=0; ion <= LIMELM; ion++ )
3995  {
3996  long nelem;
3997  double eta, fac2, xi;
3998 
3999  /* >>chng 00 jul 19, replace classical results with results including image potential
4000  * to correct for polarization of the grain as charged particle approaches. */
4001  GrainScreen(ion,nd,nz,&eta,&xi);
4002 
4003  fac2 = eta*fac1;
4004 
4005  if( fac2 == 0. )
4006  continue;
4007 
4008  for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
4009  {
4010  if( dense.lgElmtOn[nelem] && ion != gptr->RecomZ0[nelem][ion] )
4011  {
4012  gv.GrainChTrRate[nelem][ion][gptr->RecomZ0[nelem][ion]] +=
4013  (realnum)(fac2*GetAveVelocity( dense.AtomicWeight[nelem] )*atmdat.lgCTOn);
4014  }
4015  }
4016  }
4017  }
4018 # endif
4019  return;
4020 }
4021 
4022 
4023 /* this routine updates all grain quantities that depend on radius,
4024  * except gv.dstab and gv.dstsc which are updated in GrainUpdateRadius2() */
4025 STATIC void GrainUpdateRadius1(void)
4026 {
4027  long nelem;
4028 
4029  DEBUG_ENTRY( "GrainUpdateRadius1()" );
4030 
4031  for( nelem=0; nelem < LIMELM; nelem++ )
4032  {
4033  gv.elmSumAbund[nelem] = 0.f;
4034  }
4035 
4036  /* grain abundance may be a function of depth */
4037  for( size_t nd=0; nd < gv.bin.size(); nd++ )
4038  {
4039  gv.bin[nd]->GrnDpth = (realnum)GrnStdDpth(nd);
4040  gv.bin[nd]->dstAbund = (realnum)(gv.bin[nd]->dstfactor*gv.GrainMetal*gv.bin[nd]->GrnDpth);
4041  ASSERT( gv.bin[nd]->dstAbund > 0.f );
4042 
4043  /* grain unit conversion, <unit>/H (default depl) -> <unit>/cm^3 (actual depl) */
4044  gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
4045  gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
4046  /* grain unit conversion, <unit>/cm^3 (actual depl) -> <unit>/grain */
4047  gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
4048  gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;
4049 
4050  /* >>chng 01 dec 05, calculate the number density of each element locked in grains,
4051  * summed over all grain bins. this number uses the actual depletion of the grains
4052  * and is already multiplied with hden, units cm^-3. */
4053  for( nelem=0; nelem < LIMELM; nelem++ )
4054  {
4055  gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(realnum)gv.bin[nd]->cnv_H_pCM3;
4056  }
4057  }
4058  return;
4059 }
4060 
4061 
4062 /* this routine adds all the grain opacities in gv.dstab and gv.dstsc, this could not be
4063  * done in GrainUpdateRadius1 since charge and FracPop must be converged first */
4064 STATIC void GrainUpdateRadius2()
4065 {
4066  DEBUG_ENTRY( "GrainUpdateRadius2()" );
4067 
4068  for( long i=0; i < rfield.nupper; i++ )
4069  {
4070  gv.dstab[i] = 0.;
4071  gv.dstsc[i] = 0.;
4072  }
4073 
4074  /* >>chng 06 oct 05 rjrw, reorder loops */
4075  /* >>chng 11 dec 12 reorder loops so they can be vectorized, PvH */
4076  for( size_t nd=0; nd < gv.bin.size(); nd++ )
4077  {
4078  realnum dstAbund = gv.bin[nd]->dstAbund;
4079 
4080  /* >>chng 01 mar 26, from nupper to nflux */
4081  for( long i=0; i < rfield.nflux; i++ )
4082  {
4083  /* these are total absorption and scattering cross sections,
4084  * the latter should contain the asymmetry factor (1-g) */
4085  /* this is effective area per proton, scaled by depletion
4086  * dareff(nd) = darea(nd) * dstAbund(nd) */
4087  /* grain abundance may be a function of depth */
4088  /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
4089  gv.dstab[i] += gv.bin[nd]->dstab1[i]*dstAbund;
4090  gv.dstsc[i] += gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]*dstAbund;
4091  }
4092 
4093  for( long nz=0; nz < gv.bin[nd]->nChrg; nz++ )
4094  {
4095  ChargeBin *gptr = gv.bin[nd]->chrg[nz];
4096  if( gptr->DustZ <= -1 )
4097  {
4098  double FracPop = gptr->FracPop;
4099 
4100  for( long i=gptr->ipThresInf; i < rfield.nflux; i++ )
4101  gv.dstab[i] += FracPop*gptr->cs_pdt[i]*dstAbund;
4102  }
4103  }
4104  }
4105 
4106  for( long i=0; i < rfield.nflux; i++ )
4107  {
4108  /* this must be positive, zero in case of uncontrolled underflow */
4109  ASSERT( gv.dstab[i] > 0. && gv.dstsc[i] > 0. );
4110  }
4111  return;
4112 }
4113 
4114 
4115 /* GrainTemperature computes grains temperature, and gas cooling */
4116 STATIC void GrainTemperature(size_t nd,
4117  /*@out@*/ realnum *dccool,
4118  /*@out@*/ double *hcon,
4119  /*@out@*/ double *hots,
4120  /*@out@*/ double *hla)
4121 {
4122  long int i,
4123  ipLya,
4124  nz;
4125  double EhatThermionic,
4126  norm,
4127  rate,
4128  x,
4129  y;
4130  realnum dcheat;
4131 
4132  DEBUG_ENTRY( "GrainTemperature()" );
4133 
4134  /* sanity checks */
4135  ASSERT( nd < gv.bin.size() );
4136 
4137  if( trace.lgTrace && trace.lgDustBug )
4138  {
4139  fprintf( ioQQQ, " GrainTemperature starts for grain %s\n", gv.bin[nd]->chDstLab );
4140  }
4141 
4142  /* >>chng 01 may 07, this routine now completely supports the hybrid grain
4143  * charge model, and the average charge state is not used anywhere anymore, PvH */
4144 
4145  /* direct heating by incident continuum (all energies) */
4146  *hcon = 0.;
4147  /* heating by diffuse ots fields */
4148  *hots = 0.;
4149  /* heating by Ly alpha alone, for output only, is already included in hots */
4150  *hla = 0.;
4151 
4152  ipLya = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont() - 1;
4153 
4154  /* integrate over ionizing continuum; energy goes to dust and gas
4155  * GasHeatPhotoEl is what heats the gas */
4156  gv.bin[nd]->GasHeatPhotoEl = 0.;
4157 
4158  gv.bin[nd]->GrainCoolTherm = 0.;
4159  gv.bin[nd]->thermionic = 0.;
4160 
4161  dcheat = 0.f;
4162  *dccool = 0.f;
4163 
4164  gv.bin[nd]->BolFlux = 0.;
4165 
4166  /* >>chng 04 jan 25, moved initialization of phiTilde to qheat_init(), PvH */
4167 
4168  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
4169  {
4170  ChargeBin *gptr = gv.bin[nd]->chrg[nz];
4171 
4172  double hcon1 = 0.;
4173  double hots1 = 0.;
4174  double hla1 = 0.;
4175  double bolflux1 = 0.;
4176  double pe1 = 0.;
4177 
4178  /* >>chng 04 may 31, introduced lgReEvaluate2 to save time when iterating Tdust, PvH */
4179  bool lgReEvaluate1 = gptr->hcon1 < 0.;
4180  bool lgReEvaluate2 = gptr->hots1 < 0.;
4181  bool lgReEvaluate = lgReEvaluate1 || lgReEvaluate2;
4182 
4183  /* integrate over incident continuum for non-ionizing energies */
4184  if( lgReEvaluate )
4185  {
4186  long loopmax = MIN2(gptr->ipThresInf,rfield.nflux);
4187  for( i=0; i < loopmax; i++ )
4188  {
4189  double fac = gv.bin[nd]->dstab1[i]*rfield.anu[i];
4190 
4191  if( lgReEvaluate1 )
4192  hcon1 += rfield.flux[0][i]*fac;
4193 
4194  if( lgReEvaluate2 )
4195  {
4196  hots1 += rfield.SummedDif[i]*fac;
4197 # ifndef NDEBUG
4198  bolflux1 += rfield.SummedCon[i]*fac;
4199 # endif
4200  }
4201  }
4202  }
4203 
4204  /* >>chng 01 mar 02, use new expresssions for grain cooling and absorption
4205  * cross sections following the discussion with Joe Weingartner, PvH */
4206  /* >>chng 04 feb 07, use fac1, fac2 to optimize this loop, PvH */
4207  /* >>chng 06 nov 21 rjrw, factor logic out of loops */
4208 
4209  /* this is heating by incident radiation field */
4210  if( lgReEvaluate1 )
4211  {
4212  for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
4213  {
4214  hcon1 += rfield.flux[0][i]*gptr->fac1[i];
4215  }
4216  /* >>chng 04 feb 07, remember hcon1 for possible later use, PvH */
4217  gptr->hcon1 = hcon1;
4218  }
4219  else
4220  {
4221  hcon1 = gptr->hcon1;
4222  }
4223 
4224  if( lgReEvaluate2 )
4225  {
4226  for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
4227  {
4228  /* this is heating by all diffuse fields:
4229  * SummedDif has all continua and lines */
4230  hots1 += rfield.SummedDif[i]*gptr->fac1[i];
4231  /* GasHeatPhotoEl is rate grain photoionization heats the gas */
4232 #ifdef WD_TEST2
4233  pe1 += rfield.flux[0][i]*gptr->fac2[i];
4234 #else
4235  pe1 += rfield.SummedCon[i]*gptr->fac2[i];
4236 #endif
4237 # ifndef NDEBUG
4238  bolflux1 += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*rfield.anu[i];
4239  if( gptr->DustZ <= -1 )
4240  bolflux1 += rfield.SummedCon[i]*gptr->cs_pdt[i]*rfield.anu[i];
4241 # endif
4242  }
4243  gptr->hots1 = hots1;
4244  gptr->bolflux1 = bolflux1;
4245  gptr->pe1 = pe1;
4246  }
4247  else
4248  {
4249  hots1 = gptr->hots1;
4250  bolflux1 = gptr->bolflux1;
4251  pe1 = gptr->pe1;
4252  }
4253 
4254  /* heating by Ly A on dust in this zone,
4255  * only used for printout; Ly-a is already in OTS fields */
4256  /* >>chng 00 apr 18, moved calculation of hla, by PvH */
4257  /* >>chng 04 feb 01, moved calculation of hla1 outside loop for optimization, PvH */
4258  if( ipLya < MIN2(gptr->ipThresInf,rfield.nflux) )
4259  {
4260  hla1 = rfield.otslin[ipLya]*gv.bin[nd]->dstab1[ipLya]*0.75;
4261  }
4262  else if( ipLya < rfield.nflux )
4263  {
4264  /* >>chng 00 apr 18, include photo-electric effect, by PvH */
4265  hla1 = rfield.otslin[ipLya]*gptr->fac1[ipLya];
4266  }
4267  else
4268  {
4269  hla1 = 0.;
4270  }
4271 
4272  ASSERT( hcon1 >= 0. && hots1 >= 0. && hla1 >= 0. && bolflux1 >= 0. && pe1 >= 0. );
4273 
4274  *hcon += gptr->FracPop*hcon1;
4275  *hots += gptr->FracPop*hots1;
4276  *hla += gptr->FracPop*hla1;
4277  gv.bin[nd]->BolFlux += gptr->FracPop*bolflux1;
4278  if( gv.lgDHetOn )
4279  gv.bin[nd]->GasHeatPhotoEl += gptr->FracPop*pe1;
4280 
4281 # ifndef NDEBUG
4282  if( trace.lgTrace && trace.lgDustBug )
4283  {
4284  fprintf( ioQQQ, " Zg %ld bolflux: %.4e\n", gptr->DustZ,
4285  gptr->FracPop*bolflux1*EN1RYD*gv.bin[nd]->cnv_H_pCM3 );
4286  }
4287 # endif
4288 
4289  /* add in thermionic emissions (thermal evaporation of electrons), it gives a cooling
4290  * term for the grain. thermionic emissions will not be treated separately in quantum
4291  * heating since they are only important when grains are heated to near-sublimation
4292  * temperatures; under those conditions quantum heating effects will never be important.
4293  * in order to maintain energy balance they will be added to the ion contribution though */
4294  /* ThermRate is normalized per cm^2 of grain surface area, scales with total grain area */
4295  rate = gptr->FracPop*gptr->ThermRate*gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
4296  /* >>chng 01 mar 02, PotSurf[nz] term was incorrectly taken into account, PvH */
4297  EhatThermionic = 2.*BOLTZMANN*gv.bin[nd]->tedust + MAX2(gptr->PotSurf*EN1RYD,0.);
4298  gv.bin[nd]->GrainCoolTherm += rate * (EhatThermionic + gptr->ThresSurf*EN1RYD);
4299  gv.bin[nd]->thermionic += rate * (EhatThermionic - gptr->PotSurf*EN1RYD);
4300  }
4301 
4302  /* norm is used to convert all heating rates to erg/cm^3/s */
4303  norm = EN1RYD*gv.bin[nd]->cnv_H_pCM3;
4304 
4305  /* hcon is radiative heating by incident radiation field */
4306  *hcon *= norm;
4307 
4308  /* hots is total heating of the grain by diffuse fields */
4309  *hots *= norm;
4310 
4311  /* heating by Ly alpha alone, for output only, is already included in hots */
4312  *hla *= norm;
4313 
4314  gv.bin[nd]->BolFlux *= norm;
4315 
4316  /* heating by thermal collisions with gas does work
4317  * DCHEAT is grain collisional heating by gas
4318  * DCCOOL is gas cooling due to collisions with grains
4319  * they are different since grain surface recombinations
4320  * heat the grains, but do not cool the gas ! */
4321  /* >>chng 03 nov 06, moved call after renorm of BolFlux, so that GrainCollHeating can look at it, PvH */
4322  GrainCollHeating(nd,&dcheat,dccool);
4323 
4324  /* GasHeatPhotoEl is what heats the gas */
4325  gv.bin[nd]->GasHeatPhotoEl *= norm;
4326 
4327  if( gv.lgBakesPAH_heat )
4328  {
4329  /* this is a dirty hack to get BT94 PE heating rate
4330  * for PAH's included, for Lorentz Center 2004 PDR meeting, PvH */
4331  /*>>>refer PAH heating Bakes, E.L.O., & Tielens, A.G.G.M. 1994, ApJ, 427, 822 */
4332  /* >>chng 05 aug 12, change from +=, which added additional heating to what exists already,
4333  * to simply = to set the heat, this equation gives total heating */
4334  gv.bin[nd]->GasHeatPhotoEl = 1.e-24*hmi.UV_Cont_rel2_Habing_TH85_depth*
4335  dense.gas_phase[ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((hmi.UV_Cont_rel2_Habing_TH85_depth*
4336  /*>>chng 06 jul 21, use phycon.sqrte in next two lines */
4337  phycon.sqrte/dense.eden),0.73)) + 3.65e-2*pow(phycon.te/1.e4,0.7)/
4338  (1.+2.e-4*(hmi.UV_Cont_rel2_Habing_TH85_depth*phycon.sqrte/dense.eden)))/gv.bin.size();
4339 
4340  }
4341 
4342  /* >>chng 06 jun 01, add optional scale factor, set with command
4343  * set grains heat, to rescale PE heating as per Allers et al. 2005 */
4344  gv.bin[nd]->GasHeatPhotoEl *= gv.GrainHeatScaleFactor;
4345 
4346  /* >>chng 01 nov 29, removed next statement, PvH */
4347  /* dust often hotter than gas during initial TE search */
4348  /* if( nzone <= 2 ) */
4349  /* dcheat = MAX2(0.f,dcheat); */
4350 
4351  /* find power absorbed by dust and resulting temperature
4352  *
4353  * hcon is heating from incident continuum (all energies)
4354  * hots is heating from ots continua and lines
4355  * dcheat is net grain collisional and chemical heating by
4356  * particle collisions and recombinations
4357  * GrainCoolTherm is grain cooling by thermionic emissions
4358  *
4359  * GrainHeat is net heating of this grain type,
4360  * to be balanced by radiative cooling */
4361  gv.bin[nd]->GrainHeat = *hcon + *hots + dcheat - gv.bin[nd]->GrainCoolTherm;
4362 
4363  /* remember collisional heating for this grain species */
4364  gv.bin[nd]->GrainHeatColl = dcheat;
4365 
4366  /* >>chng 04 may 31, replace ASSERT of GrainHeat > 0. with if-statement and let
4367  * GrainChargeTemp sort out the consquences of GrainHeat becoming negative, PvH */
4368  /* in case where the thermionic rates become very large,
4369  * or collisional cooling dominates, this may become negative */
4370  if( gv.bin[nd]->GrainHeat > 0. )
4371  {
4372  bool lgOutOfBounds;
4373  /* now find temperature, GrainHeat is sum of total heating of grain
4374  * >>chng 97 jul 17, divide by abundance here */
4375  x = log(MAX2(DBL_MIN,gv.bin[nd]->GrainHeat*gv.bin[nd]->cnv_CM3_pH));
4376  /* >>chng 96 apr 27, as per Peter van Hoof comment */
4377  splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,x,&y,&lgOutOfBounds);
4378  gv.bin[nd]->tedust = (realnum)exp(y);
4379  }
4380  else
4381  {
4382  gv.bin[nd]->GrainHeat = -1.;
4383  gv.bin[nd]->tedust = -1.;
4384  }
4385 
4386  if( thermal.ConstGrainTemp > 0. )
4387  {
4388  bool lgOutOfBounds;
4389  /* use temperature set with constant grain temperature command */
4390  gv.bin[nd]->tedust = thermal.ConstGrainTemp;
4391  /* >>chng 04 jun 01, make sure GrainHeat is consistent with value of tedust, PvH */
4392  x = log(gv.bin[nd]->tedust);
4393  splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgOutOfBounds);
4394  gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;
4395  }
4396 
4397  /* save for later possible printout */
4398  gv.bin[nd]->TeGrainMax = (realnum)MAX2(gv.bin[nd]->TeGrainMax,gv.bin[nd]->tedust);
4399 
4400  if( trace.lgTrace && trace.lgDustBug )
4401  {
4402  fprintf( ioQQQ, " >GrainTemperature finds %s Tdst %.5e hcon %.4e ",
4403  gv.bin[nd]->chDstLab, gv.bin[nd]->tedust, *hcon);
4404  fprintf( ioQQQ, "hots %.4e dcheat %.4e GrainCoolTherm %.4e\n",
4405  *hots, dcheat, gv.bin[nd]->GrainCoolTherm );
4406  }
4407  return;
4408 }
4409 
4410 
4411 /* helper routine for initializing quantities related to the photo-electric effect */
4412 STATIC void PE_init(size_t nd,
4413  long nz,
4414  long i,
4415  /*@out@*/ double *cs1,
4416  /*@out@*/ double *cs2,
4417  /*@out@*/ double *cs_tot,
4418  /*@out@*/ double *cool1,
4419  /*@out@*/ double *cool2,
4420  /*@out@*/ double *ehat1,
4421  /*@out@*/ double *ehat2)
4422 {
4423  ChargeBin *gptr = gv.bin[nd]->chrg[nz];
4424  long ipLo1 = gptr->ipThresInfVal;
4425  long ipLo2 = gptr->ipThresInf;
4426 
4427  DEBUG_ENTRY( "PE_init()" );
4428 
4429  /* sanity checks */
4430  ASSERT( nd < gv.bin.size() );
4431  ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
4432  ASSERT( i >= 0 && i < rfield.nflux );
4433 
4436  /* contribution from valence band */
4437  if( i >= ipLo1 )
4438  {
4439  /* effective cross section for photo-ejection */
4440  *cs1 = gv.bin[nd]->dstab1[i]*gptr->yhat[i];
4441  /* >>chng 00 jul 17, use description of Weingartner & Draine, 2001 */
4442  /* ehat1 is the average energy of the escaping electron at infinity */
4443  *ehat1 = gptr->ehat[i];
4444 
4445  /* >>chng 01 nov 27, changed de-excitation energy to conserve energy,
4446  * this branch treats valence band ionizations, but for negative grains an
4447  * electron from the conduction band will de-excite into the hole in the
4448  * valence band, reducing the amount of potential energy lost. It is assumed
4449  * that no photons are ejected in the process. PvH */
4450  /* >>chng 06 mar 19, reorganized this routine in the wake of the introduction
4451  * of the WDB06 X-ray physics. The basic functionality is still the same, but
4452  * the meaning is not. A single X-ray photon can eject multiple electrons from
4453  * either the conduction band, valence band or an inner shell. In the WDB06
4454  * approximation all these electrons are assumed to be ejected from a grain
4455  * with the same charge. After the primary ejection, Auger cascades will fill
4456  * up any inner shell holes, so energetically it is as if all these electrons
4457  * come from the outermost band (either conduction or valence band, depending
4458  * on the grain charge). Recombination will also be into the outermost band
4459  * so that way energy conservation is assured. It is assumed that these Auger
4460  * cascades are so fast that they can be treated as a single event as far as
4461  * quantum heating is concerned. */
4462 
4463  /* cool1 is the amount by which photo-ejection cools the grain */
4464  if( gptr->DustZ <= -1 )
4465  *cool1 = gptr->ThresSurf + gptr->PotSurf + *ehat1;
4466  else
4467  *cool1 = gptr->ThresSurfVal + gptr->PotSurf + *ehat1;
4468 
4469  ASSERT( *ehat1 > 0. && *cool1 > 0. );
4470  }
4471  else
4472  {
4473  *cs1 = 0.;
4474  *ehat1 = 0.;
4475  *cool1 = 0.;
4476  }
4477 
4478  /* contribution from conduction band */
4479  if( gptr->DustZ <= -1 && i >= ipLo2 )
4480  {
4481  /* effective cross section for photo-detechment */
4482  *cs2 = gptr->cs_pdt[i];
4483  /* ehat2 is the average energy of the escaping electron at infinity */
4484  *ehat2 = rfield.anu[i] - gptr->ThresSurf - gptr->PotSurf;
4485  /* cool2 is the amount by which photo-detechment cools the grain */
4486  *cool2 = rfield.anu[i];
4487 
4488  ASSERT( *ehat2 > 0. && *cool2 > 0. );
4489  }
4490  else
4491  {
4492  *cs2 = 0.;
4493  *ehat2 = 0.;
4494  *cool2 = 0.;
4495  }
4496 
4497  *cs_tot = gv.bin[nd]->dstab1[i] + *cs2;
4498  return;
4499 }
4500 
4501 
4502 /* GrainCollHeating compute grains collisional heating cooling */
4503 STATIC void GrainCollHeating(size_t nd,
4504  /*@out@*/ realnum *dcheat,
4505  /*@out@*/ realnum *dccool)
4506 {
4507  long int ion,
4508  nelem,
4509  nz;
4510  H2_type ipH2;
4511  double Accommodation,
4512  CollisionRateElectr, /* rate electrons strike grains */
4513  CollisionRateMol, /* rate molecules strike grains */
4514  CollisionRateIon, /* rate ions strike grains */
4515  CoolTot,
4516  CoolBounce,
4517  CoolEmitted,
4518  CoolElectrons,
4519  CoolMolecules,
4520  CoolPotential,
4521  CoolPotentialGas,
4522  eta,
4523  HeatTot,
4524  HeatBounce,
4525  HeatCollisions,
4526  HeatElectrons,
4527  HeatIons,
4528  HeatMolecules,
4529  HeatRecombination, /* sum of abundances of ions times velocity times ionization potential times eta */
4530  HeatChem,
4531  HeatCor,
4532  Stick,
4533  ve,
4534  WeightMol,
4535  xi;
4536 
4537  /* energy deposited into grain by formation of a single H2 molecule, in eV,
4538  * >>refer grain physics Takahashi J., Uehara H., 2001, ApJ, 561, 843 */
4539  const double H2_FORMATION_GRAIN_HEATING[H2_TOP] = { 0.20, 0.4, 1.72 };
4540 
4541  DEBUG_ENTRY( "GrainCollHeating()" );
4542 
4543 
4544  /* >>chng 01 may 07, this routine now completely supports the hybrid grain
4545  * charge model, and the average charge state is not used anywhere anymore, PvH */
4546 
4547  /* this subroutine evaluates the gas heating-cooling rate
4548  * (erg cm^-3 s^-1) due to grain gas collisions.
4549  * the net effect can be positive or negative,
4550  * depending on whether the grains or gas are hotter
4551  * the physics is described in
4552  * >>refer grain physics Baldwin, Ferland, Martin et al., 1991, ApJ 374, 580 */
4553 
4554  HeatTot = 0.;
4555  CoolTot = 0.;
4556 
4557  HeatIons = 0.;
4558 
4559  gv.bin[nd]->ChemEn = 0.;
4560 
4561  /* loop over the charge states */
4562  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
4563  {
4564  ChargeBin *gptr = gv.bin[nd]->chrg[nz];
4565 
4566  /* HEAT1 will be rate collisions heat the grain
4567  * COOL1 will be rate collisions cool the gas kinetics */
4568  double Heat1 = 0.;
4569  double Cool1 = 0.;
4570  double ChemEn1 = 0.;
4571 
4572  /* ============================================================================= */
4573  /* heating/cooling due to neutrals and positive ions */
4574 
4575  /* loop over all stages of ionization */
4576  for( ion=0; ion <= LIMELM; ion++ )
4577  {
4578  /* this is heating of grains due to recombination energy of species,
4579  * and assumes that every ion is fully neutralized upon striking the grain surface.
4580  * all radiation produced in the recombination process is absorbed within the grain
4581  *
4582  * ion=0 are neutrals, ion=1 are single ions, etc
4583  * each population is weighted by the AVERAGE velocity
4584  * */
4585  CollisionRateIon = 0.;
4586  CoolPotential = 0.;
4587  CoolPotentialGas = 0.;
4588  HeatRecombination = 0.;
4589  HeatChem = 0.;
4590 
4591  /* >>chng 00 jul 19, replace classical results with results including image potential
4592  * to correct for polarization of the grain as charged particle approaches. */
4593  GrainScreen(ion,nd,nz,&eta,&xi);
4594 
4595  for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
4596  {
4597  if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. )
4598  {
4599  double CollisionRateOne;
4600 
4601  /* >>chng 00 apr 05, use correct accomodation coefficient, by PvH
4602  * the coefficient is defined at the end of appendix A.10 of BFM
4603  * assume ion sticking prob is unity */
4604 #if defined( IGNORE_GRAIN_ION_COLLISIONS )
4605  Stick = 0.;
4606 #elif defined( WD_TEST2 )
4607  Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
4608  0. : STICK_ION;
4609 #else
4610  Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
4611  gv.bin[nd]->AccomCoef[nelem] : STICK_ION;
4612 #endif
4613  /* this is rate with which charged ion strikes grain */
4614  /* >>chng 00 may 02, this had left 2./SQRTPI off */
4615  /* >>chng 00 may 05, use average speed instead of 2./SQRTPI*Doppler, PvH */
4616  CollisionRateOne = Stick*dense.xIonDense[nelem][ion]*GetAveVelocity( dense.AtomicWeight[nelem] );
4617  CollisionRateIon += CollisionRateOne;
4618  /* >>chng 01 nov 26, use PotSurfInc when appropriate:
4619  * the values for the surface potential used here make it
4620  * consistent with the rest of the code and preserve energy.
4621  * NOTE: For incoming particles one should use PotSurfInc with
4622  * Schottky effect for positive ion, for outgoing particles
4623  * one should use PotSurf for Zg+ion-Z_0-1 (-1 because PotSurf
4624  * assumes electron going out), these corrections are small
4625  * and will be neglected for now, PvH */
4626  if( ion >= gptr->RecomZ0[nelem][ion] )
4627  {
4628  CoolPotential += CollisionRateOne * (double)ion *
4629  gptr->PotSurf;
4630  CoolPotentialGas += CollisionRateOne *
4631  (double)gptr->RecomZ0[nelem][ion] *
4632  gptr->PotSurf;
4633  }
4634  else
4635  {
4636  CoolPotential += CollisionRateOne * (double)ion *
4637  gptr->PotSurfInc;
4638  CoolPotentialGas += CollisionRateOne *
4639  (double)gptr->RecomZ0[nelem][ion] *
4640  gptr->PotSurfInc;
4641  }
4642  /* this is sum of all energy liberated as ion recombines to Z0 in grain */
4643  /* >>chng 00 jul 05, subtract energy needed to get
4644  * electron out of grain potential well, PvH */
4645  /* >>chng 01 may 09, chemical energy now calculated in GrainIonColl, PvH */
4646  HeatRecombination += CollisionRateOne *
4647  gptr->RecomEn[nelem][ion];
4648  HeatChem += CollisionRateOne * gptr->ChemEn[nelem][ion];
4649  }
4650  }
4651 
4652  /* >>chng 00 may 01, Boltzmann factor had multiplied all of factor instead
4653  * of only first and last term. pvh */
4654 
4655  /* equation 29 from Balwin et al 91 */
4656  /* this is direct collision rate, 2kT * xi, first term in eq 29 */
4657  HeatCollisions = CollisionRateIon * 2.*BOLTZMANN*phycon.te*xi;
4658  /* this is change in energy due to charge acceleration within grain's potential
4659  * this is exactly balanced by deceleration of incoming electrons and accelaration
4660  * of outgoing photo-electrons and thermionic emissions; all these terms should
4661  * add up to zero (total charge of grain should remain constant) */
4662  CoolPotential *= eta*EN1RYD;
4663  CoolPotentialGas *= eta*EN1RYD;
4664  /* this is recombination energy released within grain */
4665  HeatRecombination *= eta*EN1RYD;
4666  HeatChem *= eta*EN1RYD;
4667  /* energy carried away by neutrals after recombination, so a cooling term */
4668  CoolEmitted = CollisionRateIon * 2.*BOLTZMANN*gv.bin[nd]->tedust*eta;
4669 
4670  /* total GraC 0 in the emission line output */
4671  Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;
4672 
4673  /* rate kinetic energy lost from gas - gas cooling - eq 32 in BFM */
4674  /* this GrGC 0 in the main output */
4675  /* >>chng 00 may 05, reversed sign of gas cooling contribution */
4676  Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;
4677 
4678  ChemEn1 += HeatChem;
4679  }
4680 
4681  /* remember grain heating by ion collisions for quantum heating treatment */
4682  HeatIons += gptr->FracPop*Heat1;
4683 
4684  if( trace.lgTrace && trace.lgDustBug )
4685  {
4686  fprintf( ioQQQ, " Zg %ld ions heat/cool: %.4e %.4e\n", gptr->DustZ,
4687  gptr->FracPop*Heat1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
4688  gptr->FracPop*Cool1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
4689  }
4690 
4691  /* ============================================================================= */
4692  /* heating/cooling due to electrons */
4693 
4694  ion = -1;
4695  Stick = ( gptr->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
4696  /* VE is mean (not RMS) electron velocity */
4697  /*ve = TePowers.sqrte*6.2124e5;*/
4698  ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);
4699 
4700  /* electron arrival rate - eqn 29 again */
4701  CollisionRateElectr = Stick*dense.eden*ve;
4702 
4703  /* >>chng 00 jul 19, replace classical results with results including image potential
4704  * to correct for polarization of the grain as charged particle approaches. */
4705  GrainScreen(ion,nd,nz,&eta,&xi);
4706 
4707  if( gptr->DustZ > gv.bin[nd]->LowestZg )
4708  {
4709  HeatCollisions = CollisionRateElectr*2.*BOLTZMANN*phycon.te*xi;
4710  /* this is change in energy due to charge acceleration within grain's potential
4711  * this term (perhaps) adds up to zero when summed over all charged particles */
4712  CoolPotential = CollisionRateElectr * (double)ion*gptr->PotSurfInc*eta*EN1RYD;
4713  /* >>chng 00 jul 05, this is term for energy released due to recombination, PvH */
4714  HeatRecombination = CollisionRateElectr * gptr->ThresSurfInc*eta*EN1RYD;
4715  HeatBounce = 0.;
4716  CoolBounce = 0.;
4717  }
4718  else
4719  {
4720  HeatCollisions = 0.;
4721  CoolPotential = 0.;
4722  HeatRecombination = 0.;
4723  /* >>chng 00 jul 05, add in terms for electrons that bounce off grain, PvH */
4724  /* >>chng 01 mar 09, remove these terms, their contribution is negligible, and replace
4725  * them with similar terms that describe electrons that are captured by grains at Z_min,
4726  * these electrons are not in a bound state and the grain will quickly autoionize, PvH */
4727  HeatBounce = CollisionRateElectr * 2.*BOLTZMANN*phycon.te*xi;
4728  /* >>chng 01 mar 14, replace (2kT_g - phi_g) term with -EA; for autoionizing states EA is
4729  * usually higher than phi_g, so more energy is released back into the electron gas, PvH */
4730  CoolBounce = CollisionRateElectr *
4731  (-gptr->ThresSurfInc-gptr->PotSurfInc)*EN1RYD*eta;
4732  CoolBounce = MAX2(CoolBounce,0.);
4733  }
4734 
4735  /* >>chng 00 may 02, CoolPotential had not been included */
4736  /* >>chng 00 jul 05, HeatRecombination had not been included */
4737  HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
4738  Heat1 += HeatElectrons;
4739 
4740  CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
4741  Cool1 += CoolElectrons;
4742 
4743  if( trace.lgTrace && trace.lgDustBug )
4744  {
4745  fprintf( ioQQQ, " Zg %ld electrons heat/cool: %.4e %.4e\n", gptr->DustZ,
4746  gptr->FracPop*HeatElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
4747  gptr->FracPop*CoolElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
4748  }
4749 
4750  /* add quantum heating due to recombination of electrons, subtract thermionic cooling */
4751 
4752  /* calculate net heating rate in erg/H/s at standard depl
4753  * include contributions for recombining electrons, autoionizing electrons
4754  * and subtract thermionic emissions here since it is inverse process
4755  *
4756  * NB - in extreme conditions this rate may become negative (if there
4757  * is an intense radiation field leading to very hot grains, but no ionizing
4758  * photons, hence very few free electrons). we assume that the photon rates
4759  * are high enough under those circumstances to avoid phiTilde becoming negative,
4760  * but we will check that in qheat1 anyway. */
4761  gptr->HeatingRate2 = HeatElectrons*gv.bin[nd]->IntArea/4. -
4762  gv.bin[nd]->GrainCoolTherm*gv.bin[nd]->cnv_CM3_pH;
4763 
4764  /* >>chng 04 jan 25, moved inclusion into phitilde to qheat_init(), PvH */
4765 
4766  /* heating/cooling above is in erg/s/cm^2 -> multiply with projected grain area per cm^3 */
4767  /* GraC 0 is integral of dcheat, the total collisional heating of the grain */
4768  HeatTot += gptr->FracPop*Heat1;
4769 
4770  /* GrGC 0 total cooling of gas integrated */
4771  CoolTot += gptr->FracPop*Cool1;
4772 
4773  gv.bin[nd]->ChemEn += gptr->FracPop*ChemEn1;
4774  }
4775 
4776  /* ============================================================================= */
4777  /* heating/cooling due to molecules */
4778 
4779  /* these rates do not depend on charge, hence they are outside of nz loop */
4780 
4781  /* sticking prob for H2 onto grain,
4782  * estimated from accomodation coefficient defined at end of A.10 in BFM */
4783  WeightMol = 2.*dense.AtomicWeight[ipHYDROGEN];
4784  Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
4785  /* molecular hydrogen onto grains */
4786 #ifndef IGNORE_GRAIN_ION_COLLISIONS
4787  /*CollisionRateMol = Accommodation*findspecies("H2")->den* */
4788  CollisionRateMol = Accommodation*hmi.H2_total*
4789  sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
4790  /* >>chng 03 feb 12, added grain heating by H2 formation on the surface, PvH
4791  * >>refer grain H2 heat Takahashi & Uehara, ApJ, 561, 843 */
4792  ipH2 = gv.which_H2distr[gv.bin[nd]->matType];
4793  /* this is rate in erg/cm^3/s */
4794  /* >>chng 04 may 26, changed dense.gas_phase[ipHYDROGEN] -> dense.xIonDense[ipHYDROGEN][0], PvH */
4795  gv.bin[nd]->ChemEnH2 = gv.bin[nd]->rate_h2_form_grains_used*dense.xIonDense[ipHYDROGEN][0]*
4796  H2_FORMATION_GRAIN_HEATING[ipH2]*EN1EV;
4797  /* convert to rate per cm^2 of projected grain surface area used here */
4798  gv.bin[nd]->ChemEnH2 /= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
4799 #else
4800  CollisionRateMol = 0.;
4801  gv.bin[nd]->ChemEnH2 = 0.;
4802 #endif
4803 
4804  /* now add in CO */
4805  WeightMol = dense.AtomicWeight[ipCARBON] + dense.AtomicWeight[ipOXYGEN];
4806  Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
4807 #ifndef IGNORE_GRAIN_ION_COLLISIONS
4808  CollisionRateMol += Accommodation*findspecieslocal("CO")->den*
4809  sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
4810 #else
4811  CollisionRateMol = 0.;
4812 #endif
4813 
4814  /* xi and eta are unity for neutrals and so ignored */
4815  HeatCollisions = CollisionRateMol * 2.*BOLTZMANN*phycon.te;
4816  CoolEmitted = CollisionRateMol * 2.*BOLTZMANN*gv.bin[nd]->tedust;
4817 
4818  HeatMolecules = HeatCollisions - CoolEmitted + gv.bin[nd]->ChemEnH2;
4819  HeatTot += HeatMolecules;
4820 
4821  /* >>chng 00 may 05, reversed sign of gas cooling contribution */
4822  CoolMolecules = HeatCollisions - CoolEmitted;
4823  CoolTot += CoolMolecules;
4824 
4825  gv.bin[nd]->RateUp = 0.;
4826  gv.bin[nd]->RateDn = 0.;
4827  HeatCor = 0.;
4828  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
4829  {
4830  double d[4];
4831  double rate_dn = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
4832  double rate_up = GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
4833 
4834  gv.bin[nd]->RateUp += gv.bin[nd]->chrg[nz]->FracPop*rate_up;
4835  gv.bin[nd]->RateDn += gv.bin[nd]->chrg[nz]->FracPop*rate_dn;
4836 
4838  HeatCor += (gv.bin[nd]->chrg[nz]->FracPop*rate_up*gv.bin[nd]->chrg[nz]->ThresSurf -
4839  gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->ThresSurfInc +
4840  gv.bin[nd]->chrg[nz]->FracPop*rate_up*gv.bin[nd]->chrg[nz]->PotSurf -
4841  gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->PotSurfInc)*EN1RYD;
4842  }
4843  /* >>chng 01 nov 24, correct for imperfections in the n-charge state model,
4844  * these corrections should add up to zero, but are actually small but non-zero, PvH */
4845  HeatTot += HeatCor;
4846 
4847  if( trace.lgTrace && trace.lgDustBug )
4848  {
4849  fprintf( ioQQQ, " molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
4850  HeatMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
4851  CoolMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
4852  HeatCor*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
4853  }
4854 
4855  *dcheat = (realnum)(HeatTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3);
4856  *dccool = ( gv.lgDColOn ) ? (realnum)(CoolTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3) : 0.f;
4857 
4858  gv.bin[nd]->ChemEn *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
4859  gv.bin[nd]->ChemEnH2 *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
4860 
4861  /* add quantum heating due to molecule/ion collisions */
4862 
4863  /* calculate heating rate in erg/H/s at standard depl
4864  * include contributions from molecules/neutral atoms and recombining ions
4865  *
4866  * in fully ionized conditions electron heating rates will be much higher
4867  * than ion and molecule rates since electrons are so much faster and grains
4868  * tend to be positive. in non-ionized conditions the main contribution will
4869  * come from neutral atoms and molecules, so it is appropriate to treat both
4870  * the same. in fully ionized conditions we don't care since unimportant.
4871  *
4872  * NB - if grains are hotter than ambient gas, the heating rate may become negative.
4873  * if photon rates are not high enough to prevent phiTilde from becoming negative,
4874  * we will raise a flag while calculating the quantum heating in qheat1 */
4875  /* >>chng 01 nov 26, add in HeatCor as well, otherwise energy imbalance will result, PvH */
4876  gv.bin[nd]->HeatingRate1 = (HeatMolecules+HeatIons+HeatCor)*gv.bin[nd]->IntArea/4.;
4877 
4878  /* >>chng 04 jan 25, moved inclusion into phiTilde to qheat_init(), PvH */
4879  return;
4880 }
4881 
4882 
4883 /* GrainDrift computes grains drift velocity */
4884 void GrainDrift(void)
4885 {
4886  long int i,
4887  loop;
4888  double alam,
4889  corr,
4890  dmomen,
4891  fac,
4892  fdrag,
4893  g0,
4894  g2,
4895  phi2lm,
4896  psi,
4897  rdust,
4898  si,
4899  vdold,
4900  volmom;
4901 
4902  DEBUG_ENTRY( "GrainDrift()" );
4903 
4904  vector<realnum> help( rfield.nflux );
4905  for( i=0; i < rfield.nflux; i++ )
4906  {
4907  help[i] = (rfield.flux[0][i]+rfield.ConInterOut[i]+rfield.outlin[0][i]+rfield.outlin_noplot[i])*
4908  rfield.anu[i];
4909  }
4910 
4911  for( size_t nd=0; nd < gv.bin.size(); nd++ )
4912  {
4913  /* find momentum absorbed by grain */
4914  dmomen = 0.;
4915  for( i=0; i < rfield.nflux; i++ )
4916  {
4917  /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
4918  dmomen += help[i]*(gv.bin[nd]->dstab1[i] + gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]);
4919  }
4920  ASSERT( dmomen >= 0. );
4921  dmomen *= EN1RYD*4./gv.bin[nd]->IntArea;
4922 
4923  /* now find force on grain, and drift velocity */
4924  fac = 2*BOLTZMANN*phycon.te;
4925 
4926  /* now PSI defined by
4927  * >>refer grain physics Draine and Salpeter 79 Ap.J. 231, 77 (1979) */
4928  psi = gv.bin[nd]->dstpot*TE1RYD/phycon.te;
4929  if( psi > 0. )
4930  {
4931  rdust = 1.e-6;
4932  alam = log(20.702/rdust/psi*phycon.sqrte/dense.SqrtEden);
4933  }
4934  else
4935  {
4936  alam = 0.;
4937  }
4938 
4939  phi2lm = POW2(psi)*alam;
4940  corr = 2.;
4941  /* >>chng 04 jan 31, increased loop limit 10 -> 50, precision -> 0.001, PvH */
4942  for( loop = 0; loop < 50 && fabs(corr-1.) > 0.001; loop++ )
4943  {
4944  vdold = gv.bin[nd]->DustDftVel;
4945 
4946  /* interactions with protons */
4947  si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
4948  g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4949  g2 = si/(1.329 + POW3(si));
4950 
4951  /* drag force due to protons, both linear and square in velocity
4952  * equation 4 from D+S Ap.J. 231, p77. */
4953  fdrag = fac*dense.xIonDense[ipHYDROGEN][1]*(g0 + phi2lm*g2);
4954 
4955  /* drag force due to interactions with electrons */
4956  si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.816e-6;
4957  g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4958  g2 = si/(1.329 + POW3(si));
4959  fdrag += fac*dense.eden*(g0 + phi2lm*g2);
4960 
4961  /* drag force due to collisions with hydrogen and helium atoms */
4962  si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
4963  g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4964  fdrag += fac*(dense.xIonDense[ipHYDROGEN][0] + 1.1*dense.xIonDense[ipHELIUM][0])*g0;
4965 
4966  /* drag force due to interactions with helium ions */
4967  si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.551e-4;
4968  g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4969  g2 = si/(1.329 + POW3(si));
4970  fdrag += fac*dense.xIonDense[ipHELIUM][1]*(g0 + phi2lm*g2);
4971 
4972  /* this term does not work
4973  * 2 HEIII*(G0+4.*PSI**2*(ALAM-0.693)*G2) )
4974  * this is total momentum absorbed by dust per unit vol */
4975  volmom = dmomen/SPEEDLIGHT;
4976 
4977  if( fdrag > 0. )
4978  {
4979  corr = sqrt(volmom/fdrag);
4980  gv.bin[nd]->DustDftVel = (realnum)(vdold*corr);
4981  }
4982  else
4983  {
4984  corr = 1.;
4985  gv.lgNegGrnDrg = true;
4986  gv.bin[nd]->DustDftVel = 0.;
4987  }
4988 
4989  if( trace.lgTrace && trace.lgDustBug )
4990  {
4991  fprintf( ioQQQ, " %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n",
4992  loop, gv.bin[nd]->DustDftVel, volmom );
4993  }
4994  }
4995  }
4996  return;
4997 }
4998 
4999 /* GrnVryDpth sets the grain abundance as a function of depth into cloud
5000  * this is intended as a playpen where the user can alter things at will
5001  * standard, documented, code should go in GrnStdDpth */
5002 STATIC double GrnVryDpth(
5003 
5004 /* nd is the number of the grain bin. The values are listed in the Cloudy output,
5005  * under "Average Grain Properties", and can easily be obtained by doing a trial
5006  * run without varying the grain abundance and setting stop zone to 1 */
5007 
5008  size_t nd)
5009 {
5010  DEBUG_ENTRY( "GrnVryDpth()" );
5011 
5012  ASSERT( nd < gv.bin.size() );
5013 
5014  /* --- DEFINE THE USER-SUPPLIED GRAIN ABUNDANCE FUNCTION HERE --- */
5015 
5016  /* This is the code that gets activated by the keyword "function" on the command line */
5017 
5018  /* NB some quantities may still be undefined on the first call to this routine. */
5019 
5020  /* the scale factor is the hydrogen atomic fraction, small when gas is ionized or molecular and
5021  * unity when atomic. This function is observed for PAHs across the Orion Bar, the PAHs are
5022  * strong near the ionization front and weak in the ionized and molecular gas */
5023  double GrnVryDpth_v = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
5024 
5025  /* This routine must return a scale factor >= 1.e-10 for the grain abundance at this position.
5026  * See Section A.3 in Hazy for more details on how grain abundances are calculated. The function
5027  * A_rel(r) mentioned there is this function times the multiplication factor set with the METALS
5028  * command (usually 1) and the multiplication factor set with the GRAINS command */
5029  return max(1.e-10,GrnVryDpth_v);
5030 }
t_elementnames::chElementName
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
GrainVar::GrnRecomTe
double GrnRecomTe
Definition: grainvar.h:532
BRACKET_MAX
static const long BRACKET_MAX
Definition: grains.cpp:78
zmin_type
zmin_type
Definition: grainvar.h:55
PE_CAR
@ PE_CAR
Definition: grainvar.h:74
GrainVar::GrainEmission
vector< realnum > GrainEmission
Definition: grainvar.h:578
thermal.h
GrainBin::nChrg
long nChrg
Definition: grainvar.h:438
FR1RYD
const UNUSED double FR1RYD
Definition: physconst.h:195
GrainBin::rate_h2_form_grains_CT02
double rate_h2_form_grains_CT02
Definition: grainvar.h:429
GrainVar::GrainHeatDif
double GrainHeatDif
Definition: grainvar.h:550
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
GrainCharge
STATIC void GrainCharge(size_t, double *)
Definition: grains.cpp:2228
TorF
char TorF(bool l)
Definition: cddefines.h:710
GrainVar::which_ial
ial_type which_ial[MAT_TOP]
Definition: grainvar.h:514
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
AC0
static const double AC0
Definition: grains.cpp:98
PE_SIL
@ PE_SIL
Definition: grainvar.h:75
ELEM_CHARGE
const UNUSED double ELEM_CHARGE
Definition: physconst.h:112
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
GrainBin::dstpot
double dstpot
Definition: grainvar.h:388
DF_SUBLIMATION
@ DF_SUBLIMATION
Definition: grainvar.h:41
ZMIN_CAR
@ ZMIN_CAR
Definition: grainvar.h:56
splint
void splint(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)
Definition: thirdparty.h:140
t_dense::eden
double eden
Definition: dense.h:190
THERMCONST
static const double THERMCONST
Definition: grains.cpp:106
flex_arr< realnum >
ChargeBin::ThresSurf
double ThresSurf
Definition: grainvar.h:232
one_elec
double one_elec(long nd)
Definition: grains.cpp:113
dense
t_dense dense
Definition: dense.cpp:24
AEInfo::Energy
vector< vector< double > > Energy
Definition: grainvar.h:182
STRG_SIL
@ STRG_SIL
Definition: grainvar.h:81
elementnames.h
GrainBin::rate_h2_form_grains_HM79
double rate_h2_form_grains_HM79
Definition: grainvar.h:428
GrainVar::which_H2distr
H2_type which_H2distr[MAT_TOP]
Definition: grainvar.h:517
GrainBin::GasHeatPhotoEl
double GasHeatPhotoEl
Definition: grainvar.h:403
GrainUpdateRadius2
STATIC void GrainUpdateRadius2()
Definition: grains.cpp:4064
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:175
AEInfo::nSubShell
unsigned int nSubShell
Definition: grainvar.h:177
GrainBin::dstpotsav
double dstpotsav
Definition: grainvar.h:389
GrainVar::lgBakesPAH_heat
bool lgBakesPAH_heat
Definition: grainvar.h:481
rfield
t_rfield rfield
Definition: rfield.cpp:8
t_rfield::flux
realnum ** flux
Definition: rfield.h:86
ShellData::ionPot
double ionPot
Definition: grainvar.h:143
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
ChargeBin::bolflux1
double bolflux1
Definition: grainvar.h:256
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
GrainVar::anumin
vector< realnum > anumin
Definition: grainvar.h:519
POT_SIL
@ POT_SIL
Definition: grainvar.h:63
GrainVar::GasHeatTherm
double GasHeatTherm
Definition: grainvar.h:546
GrainVar::GrainHeatLya
double GrainHeatLya
Definition: grainvar.h:549
GrainBin::avdpot
realnum avdpot
Definition: grainvar.h:395
GrainVar::lgQHeatAll
bool lgQHeatAll
Definition: grainvar.h:569
ENTH_CAR2
@ ENTH_CAR2
Definition: grainvar.h:47
GrainChargeTemp
STATIC void GrainChargeTemp(void)
Definition: grains.cpp:1746
ChargeBin::ehat
flex_arr< realnum > ehat
Definition: grainvar.h:238
AC2G
static const double AC2G
Definition: grains.cpp:100
realnum
float realnum
Definition: cddefines.h:103
conv.h
rfield.h
GrainVar::dHeatdT
double dHeatdT
Definition: grainvar.h:555
t_conv::EdenErrorAllowed
double EdenErrorAllowed
Definition: conv.h:267
STATIC
#define STATIC
Definition: cddefines.h:97
ipLITHIUM
const int ipLITHIUM
Definition: cddefines.h:307
GrainVar::TotalEden
double TotalEden
Definition: grainvar.h:528
ChargeBin::ThresSurfVal
double ThresSurfVal
Definition: grainvar.h:234
GetPotValues
STATIC void GetPotValues(size_t, long, double *, double *, double *, double *, double *, double *, bool)
Definition: grains.cpp:3802
GrainVar::lgReevaluate
bool lgReevaluate
Definition: grainvar.h:477
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
spline
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
Definition: thirdparty.h:117
mole.h
GrainBin::p_clear0
void p_clear0()
Definition: grains.cpp:272
GrainVar::lgDHetOn
bool lgDHetOn
Definition: grainvar.h:485
GrainBin::nChrgOrg
long nChrgOrg
Definition: grainvar.h:437
GrnStdDpth
STATIC double GrnStdDpth(long)
IAL_SIL
@ IAL_SIL
Definition: grainvar.h:69
GrainVar::GrainHeatSum
double GrainHeatSum
Definition: grainvar.h:548
y0psa
STATIC double y0psa(size_t, long, long, double)
Definition: grains.cpp:3509
ENTH_SIL
@ ENTH_SIL
Definition: grainvar.h:48
t_conv::setConvIonizFail
void setConvIonizFail(const char *reason, double oldval, double newval)
Definition: conv.h:107
ShellData::p_clear0
void p_clear0()
Definition: grains.cpp:236
y0b
STATIC double y0b(size_t, long, long)
Definition: grains.cpp:3441
MAT_CAR
@ MAT_CAR
Definition: grainvar.h:101
GrainBin::GrnDpth
realnum GrnDpth
Definition: grainvar.h:347
GrainVar::GasHeatPhotoEl
double GasHeatPhotoEl
Definition: grainvar.h:545
GRAIN_TMIN
const double GRAIN_TMIN
Definition: grainvar.h:17
GrainsInit
void GrainsInit(void)
Definition: grains.cpp:583
ChargeBin::PotSurf
double PotSurf
Definition: grainvar.h:227
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
GetFracPop
STATIC void GetFracPop(size_t, long, double[], double[], long *)
Definition: grains.cpp:2936
ChargeBin::hots1
double hots1
Definition: grainvar.h:255
ChargeBin
Definition: grainvar.h:199
thirdparty.h
t_dense::SqrtEden
double SqrtEden
Definition: dense.h:212
GrainVar::GraphiteEmission
vector< realnum > GraphiteEmission
Definition: grainvar.h:579
phycon
t_phycon phycon
Definition: phycon.cpp:6
GrainStartIter
void GrainStartIter(void)
Definition: grains.cpp:513
trace.h
GrainBin::BolFlux
double BolFlux
Definition: grainvar.h:401
flex_arr::clear
void clear()
Definition: container_classes.h:1928
NDEMS
const int NDEMS
Definition: grainvar.h:14
ipoint.h
T_LOOP_MAX
static const long T_LOOP_MAX
Definition: grains.cpp:86
t_rfield::outlin
realnum ** outlin
Definition: rfield.h:199
ReadAugerData
STATIC void ReadAugerData()
Definition: grains.cpp:1054
STICK_ELEC
static const double STICK_ELEC
Definition: grains.cpp:109
ChargeBin::cs_pdt
flex_arr< double > cs_pdt
Definition: grainvar.h:239
NewChargeData
STATIC void NewChargeData(long)
Definition: grains.cpp:1472
y0b01
STATIC double y0b01(size_t, long, long)
Definition: grains.cpp:3476
GrainVar::lgGrainPhysicsOn
bool lgGrainPhysicsOn
Definition: grainvar.h:479
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
SQRT2
const UNUSED double SQRT2
Definition: physconst.h:41
GrainBin::EnthSlp
double EnthSlp[NDEMS]
Definition: grainvar.h:424
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
ShellData::p_clear1
void p_clear1()
Definition: grains.cpp:245
ChargeBin::HeatingRate2
double HeatingRate2
Definition: grainvar.h:275
GrainBin::cnv_CM3_pH
double cnv_CM3_pH
Definition: grainvar.h:352
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:450
SetNChrgStates
void SetNChrgStates(long nChrg)
Definition: grains.cpp:570
ENTH_CAR
@ ENTH_CAR
Definition: grainvar.h:46
GrainUpdateRadius1
STATIC void GrainUpdateRadius1(void)
Definition: grains.cpp:4025
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
GrainBin::pure_sc1
vector< double > pure_sc1
Definition: grainvar.h:362
strstr_s
const char * strstr_s(const char *haystack, const char *needle)
Definition: cddefines.h:1429
ChargeBin::FracPop
double FracPop
Definition: grainvar.h:224
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
GrainVar::GasCoolColl
double GasCoolColl
Definition: grainvar.h:544
t_hmi::H2_total
double H2_total
Definition: hmi.h:16
GrainVar::which_pe
pe_type which_pe[MAT_TOP]
Definition: grainvar.h:515
ShellData::ipLo
long ipLo
Definition: grainvar.h:144
UpdatePot2
STATIC void UpdatePot2(size_t, long)
Definition: grains.cpp:3367
GrainVar::which_strg
strg_type which_strg[MAT_TOP]
Definition: grainvar.h:516
GrainVar::GrnElecDonateMax
realnum GrnElecDonateMax
Definition: grainvar.h:530
iso.h
PI4
const UNUSED double PI4
Definition: physconst.h:35
GrainCollHeating
STATIC void GrainCollHeating(size_t, realnum *, realnum *)
Definition: grains.cpp:4503
GrainBin::lgQHTooWide
bool lgQHTooWide
Definition: grainvar.h:415
flex_arr::realloc
void realloc(size_type end)
Definition: container_classes.h:1999
GrainVar::p_clear1
void p_clear1()
Definition: grains.cpp:388
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
t_thermal::ctot
double ctot
Definition: thermal.h:112
ShellData::nelem
long nelem
Definition: grainvar.h:141
GrainBin::nfill
long nfill
Definition: grainvar.h:383
ChargeBin::ThresSurfInc
double ThresSurfInc
Definition: grainvar.h:233
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
atmdat.h
GrainVar::clear
void clear()
Definition: grainvar.h:464
GrainBin::ChemEnH2
double ChemEnH2
Definition: grainvar.h:408
ShellData::p
flex_arr< realnum > p
Definition: grainvar.h:145
GetNextLine
STATIC void GetNextLine(const char *, FILE *, char[])
Definition: grains.cpp:1286
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:16
x1
static double x1[83]
Definition: atmdat_3body.cpp:28
POW2
#define POW2
Definition: cddefines.h:929
ShellData
Definition: grainvar.h:121
GrainBin::GrainHeat
double GrainHeat
Definition: grainvar.h:404
MIN2
#define MIN2
Definition: cddefines.h:761
MAT_CAR2
@ MAT_CAR2
Definition: grainvar.h:104
GrainBin::lgEverQHeat
bool lgEverQHeat
Definition: grainvar.h:414
ATOMIC_MASS_UNIT
const UNUSED double ATOMIC_MASS_UNIT
Definition: physconst.h:88
INCL_TUNNEL
static const bool INCL_TUNNEL
Definition: grains.cpp:60
GrainBin::RateDn
double RateDn
Definition: grainvar.h:392
nzone
long int nzone
Definition: cddefines.cpp:14
GrainBin::qtmin
double qtmin
Definition: grainvar.h:420
GrainBin::avdft
realnum avdft
Definition: grainvar.h:435
GrainVar::which_pot
pot_type which_pot[MAT_TOP]
Definition: grainvar.h:513
NCHU
const int NCHU
Definition: grainvar.h:25
GrainVar::lgGrainElectrons
bool lgGrainElectrons
Definition: grainvar.h:494
GrainVar::dstAbundThresholdNear
realnum dstAbundThresholdNear
Definition: grainvar.h:567
MAT_SIL2
@ MAT_SIL2
Definition: grainvar.h:105
PI
const UNUSED double PI
Definition: physconst.h:29
STRG_CAR
@ STRG_CAR
Definition: grainvar.h:80
GrainDrift
void GrainDrift(void)
Definition: grains.cpp:4884
hunt_bisect
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:270
CT_LOOP_MAX
static const long CT_LOOP_MAX
Definition: grains.cpp:83
GrainVar::nChrgRequested
long nChrgRequested
Definition: grainvar.h:534
AC1G
static const double AC1G
Definition: grains.cpp:99
ChargeBin::yhat_primary
flex_arr< realnum > yhat_primary
Definition: grainvar.h:237
GrainBin::lgPAHsInIonizedRegion
bool lgPAHsInIonizedRegion
Definition: grainvar.h:314
UpdateRecomZ0
STATIC void UpdateRecomZ0(size_t, long, bool)
Definition: grains.cpp:3746
ALL_STAGES
static const bool ALL_STAGES
Definition: grains.cpp:63
t_ADfA::phfit
double phfit(long int nz, long int ne, long int is, double e)
Definition: atmdat_adfa.cpp:269
AEInfo::p_clear1
void p_clear1()
Definition: grains.cpp:231
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
GrainVar::elmSumAbund
realnum elmSumAbund[LIMELM]
Definition: grainvar.h:507
t_hmi::UV_Cont_rel2_Habing_TH85_depth
realnum UV_Cont_rel2_Habing_TH85_depth
Definition: hmi.h:64
GrainBin::p_clear1
void p_clear1()
Definition: grains.cpp:291
y2s
double y2s(double, double, long, double *)
Definition: grains.cpp:3612
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
ChargeBin::yhat
flex_arr< realnum > yhat
Definition: grainvar.h:236
ChargeBin::ipThresInf
long ipThresInf
Definition: grainvar.h:221
GrainBin::DustDftVel
realnum DustDftVel
Definition: grainvar.h:434
pot_type
pot_type
Definition: grainvar.h:61
STICK_ION
static const double STICK_ION
Definition: grains.cpp:110
dense.h
AEInfo::AvNumber
vector< vector< double > > AvNumber
Definition: grainvar.h:180
H2_type
H2_type
Definition: grainvar.h:85
PE_init
STATIC void PE_init(size_t, long, long, double *, double *, double *, double *, double *, double *, double *)
Definition: grains.cpp:4412
GrainVar::which_enth
enth_type which_enth[MAT_TOP]
Definition: grainvar.h:511
ChargeBin::PotSurfInc
double PotSurfInc
Definition: grainvar.h:228
NCHRG_DEFAULT
const int NCHRG_DEFAULT
Definition: grainvar.h:27
GrainBin::cnv_H_pCM3
double cnv_H_pCM3
Definition: grainvar.h:350
GrainVar::chPAH_abundance
string chPAH_abundance
Definition: grainvar.h:498
GrainBin::dstAbund
realnum dstAbund
Definition: grainvar.h:346
trace
t_trace trace
Definition: trace.cpp:5
y2pa
double y2pa(double, double, long, double *)
Definition: grains.cpp:3568
cddefines.h
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:961
GrainBin::qnflux2
long qnflux2
Definition: grainvar.h:418
GrainVar::GrainChTrRate
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
Definition: grainvar.h:541
EN1EV
const UNUSED double EN1EV
Definition: physconst.h:192
MAGIC_AUGER_DATA
static const long MAGIC_AUGER_DATA
Definition: grains.cpp:58
pot2chrg
double pot2chrg(double x, long nd)
Definition: grains.cpp:119
atoms
t_atoms atoms
Definition: atoms.cpp:5
GrainBin::chrg
ChargeBin * chrg[NCHS]
Definition: grainvar.h:439
ipH2
@ ipH2
Definition: collision.h:17
GrainVar::lgQHeatOn
bool lgQHeatOn
Definition: grainvar.h:486
thermal
t_thermal thermal
Definition: thermal.cpp:5
GrainVar::which_zmin
zmin_type which_zmin[MAT_TOP]
Definition: grainvar.h:512
POW3
#define POW3
Definition: cddefines.h:936
ChargeBin::p_clear1
void p_clear1()
Definition: grains.cpp:264
flex_arr::zero
void zero()
Definition: container_classes.h:1933
t_trace::nTrConvg
int nTrConvg
Definition: trace.h:27
GrainVar::anumax
vector< realnum > anumax
Definition: grainvar.h:520
t_thermal::ConstGrainTemp
realnum ConstGrainTemp
Definition: thermal.h:47
GrainChrgTransferRates
STATIC void GrainChrgTransferRates(long)
Definition: grains.cpp:3976
GrainBin::dstslp
double dstslp[NDEMS]
Definition: grainvar.h:367
UpdatePot
STATIC void UpdatePot(size_t, long, long, double[], double[])
Definition: grains.cpp:2831
GrainBin::GrainGasCool
double GrainGasCool
Definition: grainvar.h:406
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
GrainBin::QHeatFailures
long QHeatFailures
Definition: grainvar.h:416
pe_type
pe_type
Definition: grainvar.h:73
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
ThetaNu
STATIC double ThetaNu(double)
Definition: grains.cpp:2794
GrainBin::AveDustZ
double AveDustZ
Definition: grainvar.h:386
GrainVar::AugerData
AEInfo * AugerData[LIMELM]
Definition: grainvar.h:536
ChargeBin::DustZ
long DustZ
Definition: grainvar.h:220
ChargeBin::fac2
flex_arr< double > fac2
Definition: grainvar.h:259
nCalledGrainDrive
static long int nCalledGrainDrive
Definition: grains.cpp:67
GrainVar::GrainHeatScaleFactor
realnum GrainHeatScaleFactor
Definition: grainvar.h:557
GrainBin::ChemEn
double ChemEn
Definition: grainvar.h:407
nint
long nint(double x)
Definition: cddefines.h:719
ShellData::ns
long ns
Definition: grainvar.h:142
GrainVar::lgDColOn
bool lgDColOn
Definition: grainvar.h:490
t_thermal::heating
double heating[LIMELM][LIMELM]
Definition: thermal.h:158
GrainBin::EnthSlp2
double EnthSlp2[NDEMS]
Definition: grainvar.h:425
MAT_PAH
@ MAT_PAH
Definition: grainvar.h:103
t_rfield::SummedDif
realnum * SummedDif
Definition: rfield.h:172
GrainBin::rate_h2_form_grains_used
double rate_h2_form_grains_used
Definition: grainvar.h:430
grains.h
GrainVar::GasHeatNet
double GasHeatNet
Definition: grainvar.h:547
GrainVar::HighestIon
long HighestIon
Definition: grainvar.h:533
splint_safe
void splint_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
Definition: thirdparty.h:166
ChargeBin::p_clear0
void p_clear0()
Definition: grains.cpp:254
AEInfo::nData
vector< unsigned int > nData
Definition: grainvar.h:178
ENTH_PAH
@ ENTH_PAH
Definition: grainvar.h:50
GrainBin::lgUseQHeat
bool lgUseQHeat
Definition: grainvar.h:413
ZMIN_SIL
@ ZMIN_SIL
Definition: grainvar.h:57
t_rfield::nflux
long int nflux
Definition: rfield.h:43
SPEEDLIGHT
const UNUSED double SPEEDLIGHT
Definition: physconst.h:100
heavy.h
GrainBin::y0b06
vector< realnum > y0b06
Definition: grainvar.h:385
t_trace::lgDustBug
bool lgDustBug
Definition: trace.h:79
HEAT_TOLER
static double HEAT_TOLER
Definition: grains.cpp:89
hmi.h
H2_TOP
@ H2_TOP
Definition: grainvar.h:93
ELECTRON_MASS
const UNUSED double ELECTRON_MASS
Definition: physconst.h:91
Yfunc
void Yfunc(long, long, double, double, double, double, double, double *, double *, double *, double *)
Definition: grains.cpp:3386
LN_TWO
const UNUSED double LN_TWO
Definition: physconst.h:50
H2_CAR
@ H2_CAR
Definition: grainvar.h:88
t_rfield::SummedCon
double * SummedCon
Definition: rfield.h:171
MAX2
#define MAX2
Definition: cddefines.h:782
GrainBin::sd
vector< ShellData * > sd
Definition: grainvar.h:384
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
GrainBin::cnv_CM3_pGR
double cnv_CM3_pGR
Definition: grainvar.h:351
GrnVryDpth
STATIC double GrnVryDpth(size_t)
Definition: grains.cpp:5002
LIMELM
const int LIMELM
Definition: cddefines.h:258
pow2
T pow2(T a)
Definition: cddefines.h:931
y1psa
STATIC double y1psa(size_t, long, double)
Definition: grains.cpp:3540
ShellData::AvNr
vector< double > AvNr
Definition: grainvar.h:148
GrainElecRecomb1
STATIC double GrainElecRecomb1(size_t, long, double *, double *)
Definition: grains.cpp:2510
HPLANCK
const UNUSED double HPLANCK
Definition: physconst.h:103
InitBinAugerData
STATIC void InitBinAugerData(size_t, long, long)
Definition: grains.cpp:1139
GrainVar::ReadRecord
vector< string > ReadRecord
Definition: grainvar.h:496
GrainRestartIter
void GrainRestartIter(void)
Definition: grains.cpp:551
t_dense::IonLow
long int IonLow[LIMELM+1]
Definition: dense.h:119
GRAIN_TMAX
const double GRAIN_TMAX
Definition: grainvar.h:19
MAT_PAH2
@ MAT_PAH2
Definition: grainvar.h:106
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
GrainVar::TotalDustHeat
realnum TotalDustHeat
Definition: grainvar.h:559
GrainBin::cnv_GR_pH
double cnv_GR_pH
Definition: grainvar.h:353
GrainVar::dstsc
vector< double > dstsc
Definition: grainvar.h:525
GrainVar::GrainHeatCollSum
double GrainHeatCollSum
Definition: grainvar.h:552
GrainBin::tedust
realnum tedust
Definition: grainvar.h:371
GrainBin::inv_att_len
vector< realnum > inv_att_len
Definition: grainvar.h:397
x2
static double x2[63]
Definition: atmdat_3body.cpp:20
t_rfield::nupper
long int nupper
Definition: rfield.h:46
GrainBin::le_thres
realnum le_thres
Definition: grainvar.h:396
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
GrainBin::dstems
double dstems[NDEMS]
Definition: grainvar.h:366
GetAveVelocity
realnum GetAveVelocity(realnum massAMU)
Definition: temp_change.cpp:530
AEInfo::IonThres
vector< double > IonThres
Definition: grainvar.h:179
Energy
Definition: energy.h:7
t_rfield::otslin
realnum * otslin
Definition: rfield.h:193
HighestIonStage
STATIC long HighestIonStage(void)
Definition: grains.cpp:3719
GrainBin::cnv_H_pGR
double cnv_H_pGR
Definition: grainvar.h:349
GrainBin::StickElecNeg
double StickElecNeg
Definition: grainvar.h:393
ChargeBin::pe1
double pe1
Definition: grainvar.h:257
GrainBin::dstab1
vector< double > dstab1
Definition: grainvar.h:361
UpdatePot1
STATIC void UpdatePot1(size_t, long, long, long)
Definition: grains.cpp:3088
fudge
double fudge(long int ipnt)
Definition: service.cpp:481
NO_TUNNEL
static const bool NO_TUNNEL
Definition: grains.cpp:61
ChargeBin::ThermRate
double ThermRate
Definition: grainvar.h:235
doppvel.h
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
GrainVar::dstab
vector< double > dstab
Definition: grainvar.h:524
IAL_CAR
@ IAL_CAR
Definition: grainvar.h:68
GrainBin::qnflux
long qnflux
Definition: grainvar.h:417
ipH2p
const int ipH2p
Definition: iso.h:29
grainvar.h
GrainBin::dstslp2
double dstslp2[NDEMS]
Definition: grainvar.h:368
ASINH
double ASINH(double x)
Definition: grains.cpp:34
powi
double powi(double, long int)
Definition: service.cpp:604
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
GrainBin::LowestZg
long LowestZg
Definition: grainvar.h:382
TOLER
static const double TOLER
Definition: grains.cpp:77
t_rfield::anu
double * anu
Definition: rfield.h:58
PlanckIntegral
STATIC double PlanckIntegral(double, size_t, long)
DF_USER_FUNCTION
@ DF_USER_FUNCTION
Definition: grainvar.h:40
ionbal.h
GrainVar::lgNegGrnDrg
bool lgNegGrnDrg
Definition: grainvar.h:482
GrainBin::lgQHeat
bool lgQHeat
Definition: grainvar.h:412
GrainBin::thermionic
double thermionic
Definition: grainvar.h:409
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
GrainIonColl
STATIC void GrainIonColl(size_t, long, long, long, const double[], const double[], long *, realnum *, realnum *)
min
long min(int a, long b)
Definition: cddefines.h:723
GrainDrive
void GrainDrive(void)
Definition: grains.cpp:1591
ShellData::y01A
vector< flex_arr< realnum > > y01A
Definition: grainvar.h:150
hmi
t_hmi hmi
Definition: hmi.cpp:5
GrainVar::dclmax
realnum dclmax
Definition: grainvar.h:561
GrainVar::lgWD01
bool lgWD01
Definition: grainvar.h:475
physconst.h
ChargeBin::hcon1
double hcon1
Definition: grainvar.h:254
GrainBin::lgChrgConverged
bool lgChrgConverged
Definition: grainvar.h:381
GrainBin::GrainCoolTherm
double GrainCoolTherm
Definition: grainvar.h:402
conv
t_conv conv
Definition: conv.cpp:5
gv
GrainVar gv
Definition: grainvar.cpp:5
ChargeBin::ChemEn
realnum ChemEn[LIMELM][LIMELM+1]
Definition: grainvar.h:262
GrainBin::cnv_GR_pCM3
double cnv_GR_pCM3
Definition: grainvar.h:354
ChargeBin::ipThresInfVal
long ipThresInfVal
Definition: grainvar.h:222
GrainBin::nDustFunc
df_type nDustFunc
Definition: grainvar.h:313
NCHS
const int NCHS
Definition: grainvar.h:22
GrainBin::GrainHeatColl
double GrainHeatColl
Definition: grainvar.h:405
t_rfield::widflx
realnum * widflx
Definition: rfield.h:65
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:164
NTOP
static const long NTOP
Definition: grains.cpp:73
t_rfield::outlin_noplot
realnum * outlin_noplot
Definition: rfield.h:200
InitEnthalpy
void InitEnthalpy(void)
Definition: grains_qheat.cpp:2482
GrainBin::RSFCheck
double RSFCheck
Definition: grainvar.h:357
ChargeBin::RecomEn
realnum RecomEn[LIMELM][LIMELM+1]
Definition: grainvar.h:261
GrainBin::avdust
realnum avdust
Definition: grainvar.h:373
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
ENTH_PAH2
@ ENTH_PAH2
Definition: grainvar.h:51
taulines.h
ShellData::y01
flex_arr< realnum > y01
Definition: grainvar.h:146
GRAIN_TMID
const double GRAIN_TMID
Definition: grainvar.h:18
t_rfield::anu3
realnum * anu3
Definition: rfield.h:80
ChargeBin::RecomZ0
long RecomZ0[LIMELM][LIMELM+1]
Definition: grainvar.h:241
t_conv::HeatCoolRelErrorAllowed
realnum HeatCoolRelErrorAllowed
Definition: conv.h:278
phycon.h
GrainVar::dsttmp
double dsttmp[NDEMS]
Definition: grainvar.h:564
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
ShowMe
void ShowMe(void)
Definition: service.cpp:181
ETILDE
static const double ETILDE
Definition: grains.cpp:103
AEInfo
Definition: grainvar.h:157
ChargeBin::nfill
long nfill
Definition: grainvar.h:223
ipH1s
const int ipH1s
Definition: iso.h:27
t_atmdat::lgCTOn
bool lgCTOn
Definition: atmdat.h:177
t_phycon::sqrte
double sqrte
Definition: phycon.h:48
GrainBin::avDGRatio
realnum avDGRatio
Definition: grainvar.h:334
TE1RYD
const UNUSED double TE1RYD
Definition: physconst.h:183
GrainBin::qtmin_zone1
double qtmin_zone1
Definition: grainvar.h:421
GrainBin::asym
vector< double > asym
Definition: grainvar.h:363
elec_esc_length
double elec_esc_length(double e, long nd)
Definition: grains.cpp:133
ChargeBin::tedust
realnum tedust
Definition: grainvar.h:253
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
GrainTemperature
STATIC void GrainTemperature(size_t, realnum *, double *, double *, double *)
Definition: grains.cpp:4116
ENTH_SIL2
@ ENTH_SIL2
Definition: grainvar.h:49
GrainVar::GrainHeatChem
double GrainHeatChem
Definition: grainvar.h:553
t_ionbal::lgGrainIonRecom
int lgGrainIonRecom
Definition: ionbal.h:228
GrainVar::GrainMetal
realnum GrainMetal
Definition: grainvar.h:506
GrainBin::DustEnth
double DustEnth[NDEMS]
Definition: grainvar.h:423
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
fnzone
double fnzone
Definition: cddefines.cpp:15
GrainVar::SilicateEmission
vector< realnum > SilicateEmission
Definition: grainvar.h:580
GrainBin::LowestPot
double LowestPot
Definition: grainvar.h:390
GrainBin::TeGrainMax
realnum TeGrainMax
Definition: grainvar.h:372
GrainElecEmis1
STATIC double GrainElecEmis1(size_t, long, double *, double *, double *, double *)
Definition: grains.cpp:2593
molezone::den
double den
Definition: mole.h:358
ShellData::Ener
vector< double > Ener
Definition: grainvar.h:149
POT_CAR
@ POT_CAR
Definition: grainvar.h:62
t_conv::lgSearch
bool lgSearch
Definition: conv.h:175
BOLTZMANN
const UNUSED double BOLTZMANN
Definition: physconst.h:97
t_phycon::te
double te
Definition: phycon.h:11
GrainZero
void GrainZero(void)
Definition: grains.cpp:500
GrainBin::RateUp
double RateUp
Definition: grainvar.h:391
DF_STANDARD
@ DF_STANDARD
Definition: grainvar.h:39
t_Heavy::ipHeavy
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
GrainBin::dstfactor
realnum dstfactor
Definition: grainvar.h:345
ChargeBin::fac1
flex_arr< double > fac1
Definition: grainvar.h:258
GrainVar::lgAnyDustVary
bool lgAnyDustVary
Definition: grainvar.h:480
GrainVar::dstAbundThresholdFar
realnum dstAbundThresholdFar
Definition: grainvar.h:568
GrainVar::lgDustOn
bool lgDustOn() const
Definition: grainvar.h:471
pow3
T pow3(T a)
Definition: cddefines.h:938
GrainVar::GrnElecHoldMax
realnum GrnElecHoldMax
Definition: grainvar.h:531
EN1RYD
const UNUSED double EN1RYD
Definition: physconst.h:179
GrainVar::nzone
long nzone
Definition: grainvar.h:523
MAT_SIL
@ MAT_SIL
Definition: grainvar.h:102
CONSERV_TOL
const double CONSERV_TOL
Definition: grainvar.h:35
CHRG_TOLER
static double CHRG_TOLER
Definition: grains.cpp:91
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
GrainBin::HeatingRate1
double HeatingRate1
Definition: grainvar.h:422
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
save
t_save save
Definition: save.cpp:5
GrainVar::p_clear0
void p_clear0()
Definition: grains.cpp:366
GrainScreen
STATIC void GrainScreen(long, size_t, long, double *, double *)
Definition: grains.cpp:2716
H2_SIL
@ H2_SIL
Definition: grainvar.h:87
GrainBin::StickElecPos
double StickElecPos
Definition: grainvar.h:394
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
ConvFail
void ConvFail(const char chMode[], const char chDetail[])
Definition: conv_fail.cpp:18
GrainVar::dphmax
realnum dphmax
Definition: grainvar.h:560
ShellData::nData
long nData
Definition: grainvar.h:147
chrg2pot
double chrg2pot(double x, long nd)
Definition: grains.cpp:126
GrainBin::lgTdustConverged
bool lgTdustConverged
Definition: grainvar.h:370
HEAT_TOLER_BIN
static double HEAT_TOLER_BIN
Definition: grains.cpp:90
GrainVar::GrainHeatInc
double GrainHeatInc
Definition: grainvar.h:551
InitEmissivities
STATIC void InitEmissivities(void)
Definition: grains.cpp:1313
AEInfo::p_clear0
void p_clear0()
Definition: grains.cpp:223
Generated by   doxygen 1.8.17