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