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