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