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

Generated for cloudy by doxygen 1.8.1.1