cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
grains_qheat.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*GrainMakeDiffuse main routine for generating the grain diffuse emission, called by RT_diffuse */
4 #include "cddefines.h"
5 #include "rfield.h"
6 #include "phycon.h"
7 #include "dense.h"
8 #include "hmi.h"
9 #include "thermal.h"
10 #include "trace.h"
11 #include "iterations.h"
12 #include "vectorize.h"
13 #include "grainvar.h"
14 #include "grains.h"
15 
16 inline double no_atoms(size_t nd)
17 {
18  return gv.bin[nd]->AvVol*gv.bin[nd]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[nd]->atomWeight;
19 }
20 
21 /* NB NB -- the sequence below has been carefully chosen and should NEVER be
22  * altered unless you really know what you are doing !! */
23 /* >>chng 03 jan 16, introduced QH_THIGH_FAIL and started using splint_safe and spldrv_safe
24  * throughout the code; this solves the parispn_a031_sil.in bug, PvH */
25 /* >>chng 03 jan 16, rescheduled QH_STEP_FAIL as non-fatal; it can be triggered at very low temps, PvH */
26 typedef enum {
27  /* the following are OK */
28  /* 0 1 2 3 */
29  QH_OK, QH_ANALYTIC, QH_ANALYTIC_RELAX, QH_DELTA,
30  /* the following are mild errors we already recovered from */
31  /* 4 5 6 7 */
32  QH_NEGRATE_FAIL, QH_LOOP_FAIL, QH_ARRAY_FAIL, QH_THIGH_FAIL,
33  /* any of these errors will prompt qheat to do another iteration */
34  /* 8 9 10 11 */
35  QH_RETRY, QH_CONV_FAIL, QH_BOUND_FAIL, QH_DELTA_FAIL,
36  /* any error larger than QH_NO_REBIN will cause GetProbDistr_LowLimit to return
37  * before even attempting to rebin the results; we may be able to recover though... */
38  /* 12 13 14 15 */
39  QH_NO_REBIN, QH_LOW_FAIL, QH_HIGH_FAIL, QH_STEP_FAIL,
40  /* any case larger than QH_FATAL is truly pathological; there is no hope of recovery */
41  /* 16 17 18 19 */
42  QH_FATAL, QH_WIDE_FAIL, QH_NBIN_FAIL, QH_REBIN_FAIL
43 } QH_Code;
44 
45 /*================================================================================*/
46 /* definitions relating to quantum heating routines */
47 
48 /* this is the minimum number of bins for quantum heating calculation to be valid */
49 static const long NQMIN = 10L;
50 
51 /* this is the lowest value for dPdlnT that should be included in the modeling */
52 static const double PROB_CUTOFF_LO = 1.e-15;
53 static const double PROB_CUTOFF_HI = 1.e-20;
54 static const double SAFETY = 1.e+8;
55 
56 /* during the first NSTARTUP steps, use special algorithm to calculate stepsize */
57 static const long NSTARTUP = 5L;
58 
59 /* if the average number of multiple events is above this number
60  * don't try full quantum heating treatment. */
61 static const double MAX_EVENTS = 150.;
62 
63 /* maximum number of tries for quantum heating routine */
64 /* >>chng 02 jan 30, changed LOOP_MAX from 10 -> 20, PvH */
65 static const long LOOP_MAX = 20L;
66 
67 /* if all else fails, divide temp estimate by this factor */
68 static const double DEF_FAC = 3.;
69 
70 /* total probability for all bins should not deviate more than this from 1. */
71 static const double PROB_TOL = 0.02;
72 
73 /* after NQTEST steps, make an estimate if prob distr will converge in NQGRID steps */
74 /* >>chng 02 jan 30, change NQTEST from 1000 -> 500, PvH */
75 static const long NQTEST = 500L;
76 
77 /* if the ratio fwhm/Umax is lower than this number
78  * don't try full quantum heating treatment. */
79 static const double FWHM_RATIO = 0.1;
80 /* if the ratio fwhm/Umax is lower than this number
81  * don't even try analytic approximation, use delta function instead */
82 static const double FWHM_RATIO2 = 0.007;
83 
84 /* maximum number of steps for integrating quantum heating probability distribution */
85 static const long MAX_LOOP = 2*NQGRID;
86 
87 /* this is the tolerance used while integrating the temperature distribution of the grains */
88 static const double QHEAT_TOL = 5.e-3;
89 
90 /* maximum number of tries before we declare that probability distribution simply won't fit */
91 static const long WIDE_FAIL_MAX = 3;
92 
93 /* multipliers for PROB_TOL used in GetProbDistr_HighLimit */
94 static const double STRICT = 1.;
95 static const double RELAX = 15.;
96 
97 /* when rebinning quantum heating results, make ln(qtemp[i]/qtemp[i-1]) == QT_RATIO */
98 /* this constant determines the accuracy with which the Wien tail of the grain emission is
99  * calculated; if x = h*nu/k*T_gr, d = QT_RATIO-1., and p(T) = p0 * T_gr^C, then the
100  * relative accuracy of the flux in the Wien tail is
101  * rel. acc. = d^2/24*fabs(C^2 + (2x-1)*C + x^2 -2*x)
102  * A typical value of x would be x = 12 (this corresponds to Bnu ~ 1e-2*Bnu(peak)), and C
103  * varies around -4 near the highest temperatures, so QT_RATIO = 1.074 should converge the
104  * spectrum to 1% at x = 12 and better closer to the peak of the spectrum. Since C can
105  * vary, the actual error should be between 0 and 2% at x=12. */
106 static const double QT_RATIO = 1.074;
107 
108 
109 /*================================================================================*/
110 /* global variables */
111 
112 /* these data define the enthalpy function for silicates
113  * derived from:
114  * >>refer grain physics Guhathakurta & Draine, 1989, ApJ, 345, 230
115  * coefficients converted to rydberg energy units, per unit atom
116  * assuming a density of 3.3 g/cm^3 and pure MgSiFeO4.
117  * this is not right, but the result is correct because number
118  * of atoms will be calculated using the same assumption. */
119 
120 /* this is the specific density of silicate in g/cm^3 */
121 static const double DEN_SIL = 3.30;
122 
123 /* these are the mean molecular weights per atom for MgSiFeO4 and SiO2 in amu */
124 static const double MW_SIL = 24.6051;
125 /*static const double MW_SIO2 = 20.0283;*/
126 
127 static const double tlim[5]={0.,50.,150.,500.,DBL_MAX};
128 static const double ppower[4]={2.00,1.30,0.68,0.00};
129 static const double cval[4]={
130  1.40e3/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
131  2.20e4/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
132  4.80e5/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
133  3.41e7/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD};
134 
135 
136 /* initialize phiTilde */
137 STATIC void qheat_init(size_t,/*@out@*/vector<double>&,/*@out@*/double*);
138 /* worker routine, this implements the algorithm of Guhathakurtha & Draine */
139 STATIC void GetProbDistr_LowLimit(size_t,double,double,double,/*@in@*/const vector<double>&,
140  /*@in@*/const vector<double>&,/*@out@*/vector<double>&,
141  /*@out@*/vector<double>&,/*@out@*/vector<double>&,
142  /*@out@*/long*,/*@out@*/double*,long*,/*@out@*/QH_Code*);
143 /* try two consecutive integration steps using stepsize "step/2." (yielding p[k]),
144  * and also one double integration step using stepsize "step" (yielding p2k). */
145 STATIC double TryDoubleStep(vector<double>&,vector<double>&,vector<double>&,vector<double>&,
146  vector<double>&,const vector<double>&,const vector<double>&,double,
147  /*@out@*/double*,double,long,size_t,/*@out@*/bool*);
148 /* calculate logarithmic integral from (x1,y1) to (x2,y2) */
149 STATIC double log_integral(double,double,double,double,double,double,double,double);
150 /* scan the probability distribution for valid range */
151 STATIC void ScanProbDistr(const vector<double>&,const vector<double>&,long,double,long,double,
152  /*@out@*/long*,/*@out@*/long*,/*@out@*/long*,long*,QH_Code*);
153 /* rebin the quantum heating results to speed up RT_diffuse */
154 STATIC long RebinQHeatResults(size_t,long,long,vector<double>&,vector<double>&,vector<double>&,
155  vector<double>&,vector<double>&,vector<double>&,vector<double>&,QH_Code*);
156 /* calculate approximate probability distribution in high intensity limit */
157 STATIC void GetProbDistr_HighLimit(long,double,double,double,/*@out@*/vector<double>&,/*@out@*/vector<double>&,
158  /*@out@*/vector<double>&,/*@out@*/double*,
159  /*@out@*/long*,/*@out@*/double*,/*@out@*/QH_Code*);
160 /* derivative of the enthalpy function dU/dT */
161 STATIC double uderiv(double,size_t);
162 /* enthalpy function */
163 STATIC double ufunct(double,size_t,/*@out@*/bool*);
164 /* inverse enthalpy function */
165 STATIC double inv_ufunct(double,size_t,/*@out@*/bool*);
166 /* helper function for calculating enthalpy, uses Debye theory */
167 STATIC double DebyeDeriv(double,long);
168 
169 /* >>chng 01 oct 29, introduced gv.bin[nd]->cnv_H_pGR, cnv_GR_pH, etc. PvH */
170 
171 
172 #ifndef __AVX__
173 
174 /* this assures 1e-4 relative precision in the evaluation of exp(x)-1 below */
175 static const realnum LIM1 = 2.e-4f;
176 static const realnum LIM2 = sqrtf(6.e-4f);
177 static const realnum LIM3 = cbrtf(24.e-4f);
178 
179 // fast version of expm1() with reduced precision, it implicitly assumes x >= 0!
180 inline realnum fast_expm1(realnum x)
181 {
182  if( x < LIM1 )
183  return x;
184  else if( x < LIM2 )
185  return (x/2.f + 1.f)*x;
186  else if( x < LIM3 )
187  return ((x/6.f + 0.5f)*x +1.f)*x;
188  else
189  return expf(x) - 1.f;
190 }
191 
192 #endif
193 
194 STATIC long GrainMakeDiffuseSingle(double Tgrain, double fracpop, avx_ptr<realnum>& flux, long nflux)
195 {
196  avx_ptr<realnum> arg(nflux), val(nflux);
197  const realnum x = log(0.999f*FLT_MAX);
198 
199  double hokT = TE1RYD/Tgrain;
200  // NB NB -- the loops in this routine consume lots of CPU time, they should be well optimized!!
201  for( long i=0; i < nflux; i++ )
202  {
203  arg[i] = hokT*rfield.anu(i);
204  if( arg[i] > x )
205  {
206  nflux = i;
207  break;
208  }
209  }
210 #ifdef __AVX__
211  vexpm1( arg.ptr0(), val.ptr0(), 0, nflux );
212  for( long i=0; i < nflux; i++ )
213  flux[i] += fracpop/val[i];
214 #else
215  for( long i=0; i < nflux; i++ )
216  flux[i] += fracpop/fast_expm1(arg[i]);
217 #endif
218  return nflux;
219 }
220 
221 /* main routine for generating the grain diffuse emission, called by RT_diffuse */
222 void GrainMakeDiffuse()
223 {
224  DEBUG_ENTRY( "GrainMakeDiffuse()" );
225 
226  const double factor = 2.*PI4*pow2(FR1RYD/SPEEDLIGHT)*FR1RYD;
227 
228  /* save grain emission per unit volume */
229  for( long i=0; i < rfield.nflux; i++ )
230  {
231  /* used in RT_diffuse to derive total emission */
232  gv.GrainEmission[i] = 0.;
233  gv.SilicateEmission[i] = 0.;
234  gv.GraphiteEmission[i] = 0.;
235  }
236 
237  vector<double> qtemp(NQGRID);
238  vector<double> qprob(NQGRID);
239  avx_ptr<realnum> flux(rfield.nflux);
240 
241  for( size_t nd=0; nd < gv.bin.size(); nd++ )
242  {
243  /* this local copy is necessary to keep lint happy */
244  bool lgLocalQHeat = gv.bin[nd]->lgQHeat;
245  /* >>chng 04 nov 09, do not evaluate quantum heating if abundance is negligible, PvH
246  * this prevents PAH's deep inside molecular regions from failing if GrnVryDepth is used */
247  /* >>chng 04 dec 31, introduced separate thresholds near I-front and in molecular region, PvH */
248  realnum threshold = ( dense.xIonDense[ipHYDROGEN][0]+dense.xIonDense[ipHYDROGEN][1] > hmi.H2_total ) ?
249  gv.dstAbundThresholdNear : gv.dstAbundThresholdFar;
250  long qnbin=-200;
251 
252  if( lgLocalQHeat && gv.bin[nd]->dstAbund >= threshold )
253  {
254  qheat(qtemp,qprob,&qnbin,nd);
255 
256  if( gv.bin[nd]->lgUseQHeat )
257  {
258  ASSERT( qnbin > 0 );
259  }
260  }
261  else
262  {
263  /* >> chng 04 dec 31, repaired bug lgUseQHeat not set when abundance below threshold, PvH */
264  gv.bin[nd]->lgUseQHeat = false;
265  }
266 
267  long loopmax = 0;
268  memset( flux.data(), 0, size_t(rfield.nflux*sizeof(flux[0])) );
269 
270  if( lgLocalQHeat && gv.bin[nd]->lgUseQHeat )
271  {
272  for( long j=0; j < qnbin; j++ )
273  {
274  long maxi = GrainMakeDiffuseSingle(qtemp[j], qprob[j], flux, rfield.nflux);
275  loopmax = max(loopmax, maxi);
276  }
277  }
278  else
279  {
280  loopmax = GrainMakeDiffuseSingle(gv.bin[nd]->tedust, 1., flux, rfield.nflux);
281  }
282 
283  realnum* gt_ptr;
284  /* unit emission for each different grain type */
285  strg_type scase = gv.which_strg[gv.bin[nd]->matType];
286  switch( scase )
287  {
288  case STRG_SIL:
289  gt_ptr = get_ptr(gv.SilicateEmission);
290  break;
291  case STRG_CAR:
292  gt_ptr = get_ptr(gv.GraphiteEmission);
293  break;
294  default:
295  TotalInsanity();
296  }
297 
298  double fac = factor*gv.bin[nd]->cnv_H_pCM3;
299  // use two separate loops so that they can be vectorized
300  for( long i=0; i < loopmax; i++ )
301  flux[i] *= realnum(fac*gv.bin[nd]->dstab1[i]*rfield.anu2(i)*rfield.widflx(i));
302  for( long i=0; i < loopmax; i++ )
303  {
304  /* remember local emission -- these are zeroed out on each zone
305  * above, and now incremented so is unit emission from this zone */
306  gv.GrainEmission[i] += flux[i];
307  gt_ptr[i] += flux[i];
308  }
309  }
310 
311 # ifndef NDEBUG
312  /*********************************************************************************
313  *
314  * Following are three checks on energy and charge conservation by the grain code.
315  * Their primary function is as an internal consistency check, so that coding
316  * errors get caught as early as possible. This is why the program terminates as
317  * soon as any one of the checks fails.
318  *
319  * NB NB - despite appearances, these checks do NOT guarantee overall energy
320  * conservation in the Cloudy model to the asserted tolerance, see note 1B !
321  *
322  * Note 1: there are two sources for energy imbalance in the grain code (see A & B).
323  * A: Interpolation in dstems. The code calculates how much energy the grains
324  * emit in thermal radiation (gv.bin[nd]->GrainHeat), and converts that into
325  * an (average) grain temperature by reverse interpolation in dstems. If
326  * quantum heating is not used, that temperature is used directly to generate
327  * the local diffuse emission. Hence the finite resolution of the dstems grid
328  * can lead to small errors in flux. This is tested in Check 1. The maximum
329  * error of interpolating in dstems scales with NDEMS^-3. The same problem
330  * can also occur when quantum heating is used, however, the fact that many
331  * bins are used will probably lead to a cancellation effect of the errors.
332  * B: RT_OTS_Update gets called AFTER grain() has completed, so the grain heating
333  * was calculated using a different version of the OTS fields than the one
334  * that gets fed into the RT routines (where the actual attenuation of the
335  * radiation fields by the grain opacity is done). This can lead to an energy
336  * imbalance, depending on how accurate the convergence of the OTS fields is.
337  * This is outside the control of the grain code and is therefore NOT checked.
338  * Rather, the grain code remembers the contribution from the old OTS fields
339  * (through gv.bin[nd]->BolFlux) and uses that in Check 3. In most models the
340  * difference will be well below 0.1%, but in AGN type models where OTS continua
341  * are important, the energy imbalance can be of the order of 0.5% of the grain
342  * heating (status nov 2001). On 04 jan 25 the initialization of phiTilde has
343  * been moved to qheat, implying that phiTilde now uses the updated version of
344  * the OTS fields. The total amount of radiated energy however is still based
345  * on gv.bin[nd]->GrainHeat which uses the old version of the OTS fields.
346  * C: Energy conservation for collisional processes is guaranteed by adding in
347  * (very small) correction terms. These corrections are needed to cover up
348  * small imperfection in the theory, and cannot be avoided without making the
349  * already very complex theory even more complex.
350  * D: Photo-electric heating and collisional cooling can have an important effect
351  * on the total heating balance of the gas. Both processes depend strongly on
352  * the grain charge, so assuring proper charge balance is important as well.
353  * This is tested in Check 2.
354  *
355  * Note 2: for quantum heating it is important to resolve the Maxwell distribution
356  * of the electrons and ions far enough into the tail that the total amount of
357  * energy contained in the remaining tail is negligible. If this is not the case,
358  * the assert at the beginning of the qheat() routine will fail. This is because
359  * the code filling up the phiTilde array in GrainCollHeating() assumes a value for
360  * the average particle energy based on a Maxwell distribution going to infinity.
361  * If the maximum energy used is too low, the assumed average energy would be
362  * incorrect.
363  *
364  *********************************************************************************/
365 
366  bool lgNoTdustFailures = true;
367  for( size_t nd=0; nd < gv.bin.size(); nd++ )
368  {
369  if( !gv.bin[nd]->lgTdustConverged )
370  {
371  lgNoTdustFailures = false;
372  break;
373  }
374  }
375 
376  /* CHECK 1: does the grain thermal emission conserve energy ? */
377  double BolFlux = 0.;
378  for( long i=0; i < rfield.nflux; i++ )
379  {
380  BolFlux += gv.GrainEmission[i]*rfield.anu(i)*EN1RYD;
381  }
382  double Comparison1 = 0.;
383  for( size_t nd=0; nd < gv.bin.size(); nd++ )
384  {
385  if( gv.bin[nd]->tedust < gv.bin[nd]->Tsublimat )
386  Comparison1 += CONSERV_TOL*gv.bin[nd]->GrainHeat;
387  else
388  /* for high temperatures the interpolation in dstems
389  * is less accurate, so we have to be more lenient */
390  Comparison1 += 10.*CONSERV_TOL*gv.bin[nd]->GrainHeat;
391  }
392 
393  /* >>chng 04 mar 11, add constant grain temperature to pass assert */
394  /* >>chng 04 jun 01, deleted test for constant grain temperature, PvH */
395  ASSERT( fabs(BolFlux-gv.GrainHeatSum) < Comparison1 );
396 
397  /* CHECK 2: assert charging balance */
398  for( size_t nd=0; nd < gv.bin.size(); nd++ )
399  {
400  double ave = 0.5*(gv.bin[nd]->RateDn+gv.bin[nd]->RateUp);
401  /* >>chng 04 dec 16, add lgAbort - do not throw assert if we are in
402  * process of aborting */
403  ASSERT( lgAbort || fabs(gv.bin[nd]->RateDn-gv.bin[nd]->RateUp) < CONSERV_TOL*ave );
404  }
405 
406  if( lgNoTdustFailures && gv.lgDHetOn && gv.lgDColOn && thermal.ConstGrainTemp == 0. )
407  {
408  /* CHECK 3: calculate the total energy donated to grains, must be balanced by
409  * the energy emitted in thermal radiation plus various forms of gas heating */
410  Comparison1 = 0.;
411  for( size_t nd=0; nd < gv.bin.size(); nd++ )
412  {
413  Comparison1 += gv.bin[nd]->BolFlux;
414  }
415  /* add in collisional heating of grains by plasma (if positive) */
416  Comparison1 += MAX2(gv.GasCoolColl,0.);
417  /* add in net amount of chemical energy donated by recombining ions and molecule formation */
418  Comparison1 += gv.GrainHeatChem;
419 
420  /* thermal emis PE effect gas heating by coll thermionic emis */
421  double Comparison2 = gv.GrainHeatSum+thermal.heating(0,13)+thermal.heating(0,14)+thermal.heating(0,25);
422 
423  /* >>chng 06 jun 02, add test on gv.GrainHeatScaleFactor so that assert not thrown
424  * when set grain heat command is used */
425  ASSERT( lgAbort || gv.GrainHeatScaleFactor != 1.f || gv.lgBakesPAH_heat ||
426  fabs(Comparison1-Comparison2)/Comparison2 < CONSERV_TOL );
427  }
428 # endif
429  return;
430 }
431 
432 
433 /****************************************************************************
434  *
435  * qheat: driver routine for non-equilibrium heating of grains
436  *
437  * This routine calls GetProbDistr_LowLimit, GetProbDistr_HighLimit
438  * (which do the actual non-equilibrium calculations), and does the
439  * subsequent quality control.
440  *
441  * Written by Peter van Hoof (UK, CITA, QUB).
442  *
443  ****************************************************************************/
444 
445 /* this is the new version of the quantum heating code, used starting Cloudy 96 beta 3 */
446 
447 void qheat(/*@out@*/ vector<double>& qtemp, /* qtemp[NQGRID] */
448  /*@out@*/ vector<double>& qprob, /* qprob[NQGRID] */
449  /*@out@*/ long int *qnbin,
450  size_t nd)
451 {
452  bool lgBoundErr,
453  lgDelta,
454  lgNegRate,
455  lgOK,
456  lgTried;
457  long int i,
458  nWideFail;
459  QH_Code ErrorCode,
460  ErrorCode2,
461  ErrorStart;
462  double c0,
463  c1,
464  c2,
465  check,
466  DefFac,
467  deriv,
468  fwhm,
469  FwhmRatio,
470  integral,
471  minBracket,
472  maxBracket,
473  new_tmin,
474  NumEvents,
475  rel_tol,
476  Tmax,
477  tol,
478  Umax,
479  xx,
480  y;
481 
482  DEBUG_ENTRY( "qheat()" );
483 
484  /* sanity checks */
485  ASSERT( gv.bin[nd]->lgQHeat );
486  ASSERT( nd < gv.bin.size() );
487 
488  if( trace.lgTrace && trace.lgDustBug )
489  {
490  fprintf( ioQQQ, "\n >>>>qheat called for grain %s\n", gv.bin[nd]->chDstLab );
491  }
492 
493  /* >>chng 01 aug 22, allocate space */
494  /* phiTilde is continuum corrected for photo-electric effect, in events/H/s/cell, default depl */
495  vector<double> phiTilde(rfield.nflux_with_check);
496  vector<double> Phi(rfield.nflux_with_check);
497  vector<double> PhiDrv(rfield.nflux_with_check);
498  vector<double> dPdlnT(NQGRID);
499 
500  qheat_init( nd, phiTilde, &check );
501 
502  check += gv.bin[nd]->GrainHeatColl-gv.bin[nd]->GrainCoolTherm;
503 
504  xx = integral = 0.;
505  c0 = c1 = c2 = 0.;
506  lgNegRate = false;
507  tol = 1.;
508 
509  /* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
510  /* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
511  for( i=gv.bin[nd]->qnflux-1; i >= 0; i-- )
512  {
513  /* >>chng 97 jul 17, to summed continuum */
514  /* >>chng 00 mar 30, to phiTilde, to keep track of photo-electric effect and collisions, by PvH */
515  /* >>chng 01 oct 10, use trapezoidal rule for integrating Phi, reverse direction of integration. */
516  /* >>chng 01 oct 30, change normalization of Phi, PhiDrv from <unit>/cm^3 -> <unit>/grain, PvH */
517  /* phiTilde has units events/H/s, PhiDrv[i] has units events/grain/s/Ryd */
518  /* there are minus signs here because we are integrating from infinity downwards */
519  y = -phiTilde[i]*gv.bin[nd]->cnv_H_pGR;
520  PhiDrv[i] = y/rfield.widflx(i);
521  xx -= y;
522  /* Phi[i] is integral from exactly rfield.anumin(i) to infinity to second-order precision, PvH */
523  /* Phi[i] has units events/grain/s */
524  Phi[i] = xx;
525 
526 # ifndef NDEBUG
527  /* trapezoidal rule is not needed for integral, this is also second-order correct */
528  integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu(i)*EN1RYD;
529 # endif
530 
531  /* c<n> has units Ryd^(n+1)/grain/s */
532  c0 += Phi[i]*rfield.widflx(i);
533  c1 += Phi[i]*rfield.anu(i)*rfield.widflx(i);
534  c2 += Phi[i]*rfield.anu2(i)*rfield.widflx(i);
535 
536  lgNegRate = lgNegRate || ( phiTilde[i] < 0. );
537  }
538 
539  /* sanity check */
540  ASSERT( fabs(check-integral)/check <= CONSERV_TOL );
541 
542 # if 0
543  {
544  char fnam[50];
545  FILE *file;
546 
547  sprintf(fnam,"Lambda_%2.2ld.asc",nd);
548  file = open_data(fnam,"w");
549  for( i=0; i < NDEMS; ++i )
550  fprintf(file,"%e %e %e\n",
551  exp(gv.dsttmp[i]),
552  ufunct(exp(gv.dsttmp[i]),nd,&lgBoundErr),
553  exp(gv.bin[nd]->dstems[i])*gv.bin[nd]->cnv_H_pGR/EN1RYD);
554  fclose(file);
555 
556  sprintf(fnam,"Phi_%2.2ld.asc",nd);
557  file = open_data(fnam,"w");
558  for( i=0; i < gv.bin[nd]->qnflux; ++i )
559  fprintf(file,"%e %e\n", rfield.anu(i),Phi[i]);
560  fclose(file);
561  }
562 # endif
563 
564  /* Tmax is where p(U) will peak (at least in high intensity limit) */
565  Tmax = gv.bin[nd]->tedust;
566  /* grain enthalpy at peak of p(U), in Ryd */
567  Umax = ufunct(Tmax,nd,&lgBoundErr);
568  ASSERT( !lgBoundErr ); /* this should never happen */
569  /* y is dln(Lambda)/dlnT at Tmax */
570  spldrv_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(Tmax),&y,&lgBoundErr);
571  ASSERT( !lgBoundErr ); /* this should never happen */
572  /* deriv is dLambda/dU at Umax, in 1/grain/s */
573  deriv = y*c0/(uderiv(Tmax,nd)*Tmax);
574  /* estimated FWHM of probability distribution, in Ryd */
575  fwhm = sqrt(8.*LN_TWO*c1/deriv);
576 
577  NumEvents = pow2(fwhm)*c0/(4.*LN_TWO*c2);
578  FwhmRatio = fwhm/Umax;
579 
580  /* >>chng 01 nov 15, change ( NumEvents > MAX_EVENTS2 ) --> ( FwhmRatio < FWHM_RATIO2 ), PvH */
581  lgDelta = ( FwhmRatio < FWHM_RATIO2 );
582  /* high intensity case is always OK since we will use equilibrium treatment */
583  lgOK = lgDelta;
584 
585  ErrorStart = QH_OK;
586 
587  if( lgDelta )
588  {
589  /* in this case we ignore negative rates, equilibrium treatment is good enough */
590  lgNegRate = false;
591  ErrorStart = MAX2(ErrorStart,QH_DELTA);
592  }
593 
594  if( lgNegRate )
595  ErrorStart = MAX2(ErrorStart,QH_NEGRATE_FAIL);
596 
597  ErrorCode = ErrorStart;
598 
599  if( trace.lgTrace && trace.lgDustBug )
600  {
601  double Rate2 = 0.;
602  for( int nz=0; nz < gv.bin[nd]->nChrg; nz++ )
603  Rate2 += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->HeatingRate2;
604 
605  fprintf( ioQQQ, " grain heating: %.4e, integral %.4e, total rate %.4e lgNegRate %c\n",
606  gv.bin[nd]->GrainHeat,integral,Phi[0],TorF(lgNegRate));
607  fprintf( ioQQQ, " av grain temp %.4e av grain enthalpy (Ryd) %.4e\n",
608  gv.bin[nd]->tedust,Umax);
609  fprintf( ioQQQ, " fwhm^2/(4ln2*c2/c0): %.4e fwhm (Ryd) %.4e fwhm/Umax %.4e\n",
610  NumEvents,fwhm,FwhmRatio );
611  fprintf( ioQQQ, " HeatingRate1 %.4e HeatingRate2 %.4e lgQHTooWide %c\n",
612  gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3, Rate2*gv.bin[nd]->cnv_H_pCM3,
613  TorF(gv.bin[nd]->lgQHTooWide) );
614  }
615 
616  /* these two variables will bracket qtmin, they should only be needed during the initial search phase */
617  minBracket = GRAIN_TMIN;
618  maxBracket = gv.bin[nd]->tedust;
619 
620  /* >>chng 02 jan 30, introduced lgTried to avoid running GetProbDistr_HighLimit over and over..., PvH */
621  lgTried = false;
622  /* >>chng 02 aug 06, introduced QH_WIDE_FAIL and nWideFail, PvH */
623  nWideFail = 0;
624  /* >>chng 03 jan 27, introduced DefFac to increase factor for repeated LOW_FAIL's, PvH */
625  DefFac = DEF_FAC;
626  /* >>chng 04 nov 10, introduced rel_tol to increase precision in case of repeated CONV_FAIL's, PvH */
627  rel_tol = 1.;
628 
629  /* if average number of multiple photon events is too high, lgOK is set to true */
630  /* >>chng 02 aug 12, added gv.bin[nd]->lgQHTooWide to prevent unnecessary looping here.
631  * In general the number of integration steps needed to integrate the probability distribution
632  * will increase monotonically with depth into the cloud. Hence, once the distribution becomes
633  * too wide to fit into NQGRID steps (which only happens for extremely cold grains in deeply
634  * shielded conditions) there is no hope of ever converging GetProbDistr_LowLimit and the code
635  * will waste a lot of CPU time establishing this for every zone again. So once the distribution
636  * becomes too wide we immediately skip to the analytic approximation to save time, PvH */
637  for( i=0; i < LOOP_MAX && !lgOK && !gv.bin[nd]->lgQHTooWide; i++ )
638  {
639  if( gv.bin[nd]->qtmin >= gv.bin[nd]->tedust )
640  {
641  /* >>chng 02 jul 31, was gv.bin[nd]->qtmin = 0.7*gv.bin[nd]->tedust */
642  /* >>chng 03 nov 10, changed Umax/exp(+... to Umax*exp(-... to avoid overflow, PvH */
643  double Ulo = Umax*exp(-sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO))*fwhm/Umax);
644  double MinEnth = exp(gv.bin[nd]->DustEnth[0]);
645  Ulo = MAX2(Ulo,MinEnth);
646  gv.bin[nd]->qtmin = inv_ufunct(Ulo,nd,&lgBoundErr);
647  ASSERT( !lgBoundErr ); /* this should never happen */
648  /* >>chng 02 jul 30, added this test; prevents problems with ASSERT below, PvH */
649  if( gv.bin[nd]->qtmin <= minBracket || gv.bin[nd]->qtmin >= maxBracket )
650  gv.bin[nd]->qtmin = sqrt(minBracket*maxBracket);
651  }
652  gv.bin[nd]->qtmin = MAX2(gv.bin[nd]->qtmin,GRAIN_TMIN);
653 
654  ASSERT( minBracket <= gv.bin[nd]->qtmin && gv.bin[nd]->qtmin <= maxBracket );
655 
656  ErrorCode = ErrorStart;
657 
658  /* >>chng 01 nov 15, add ( FwhmRatio >= FWHM_RATIO ), PvH */
659  if( FwhmRatio >= FWHM_RATIO && NumEvents <= MAX_EVENTS )
660  {
661  GetProbDistr_LowLimit(nd,rel_tol,Umax,fwhm,Phi,PhiDrv,qtemp,qprob,dPdlnT,qnbin,
662  &new_tmin,&nWideFail,&ErrorCode);
663 
664  /* >>chng 02 jan 07, various changes to improve convergence for very cold grains, PvH */
665  if( ErrorCode == QH_DELTA_FAIL && fwhm < Umax && !lgTried )
666  {
667  double dummy;
668 
669  /* this situation can mean two things: either the photon rate is so high that
670  * the code needs too many steps to integrate the probability distribution,
671  * or alternatively, tmin is far too low and the code needs too many steps
672  * to overcome the rising side of the probability distribution.
673  * So we call GetProbDistr_HighLimit first to determine if the former is the
674  * case; if that fails then the latter must be true and we reset QH_DELTA_FAIL */
675  ErrorCode = MAX2(ErrorStart,QH_ANALYTIC);
676  /* use dummy to avoid losing estimate for new_tmin from GetProbDistr_LowLimit */
677  /* >>chng 02 aug 06, introduced STRICT and RELAX, PvH */
678  GetProbDistr_HighLimit(nd,STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
679  &ErrorCode);
680 
681  if( ErrorCode >= QH_RETRY )
682  {
683  ErrorCode = QH_DELTA_FAIL;
684  lgTried = true;
685  }
686  }
687 
688  /* >>chng 02 aug 07 major rewrite of the logic below */
689  if( ErrorCode < QH_NO_REBIN )
690  {
691  if( new_tmin < minBracket || new_tmin > maxBracket )
692  ++nWideFail;
693 
694  if( nWideFail < WIDE_FAIL_MAX )
695  {
696  if( new_tmin <= minBracket )
697  new_tmin = sqrt(gv.bin[nd]->qtmin*minBracket);
698  if( new_tmin >= maxBracket )
699  new_tmin = sqrt(gv.bin[nd]->qtmin*maxBracket);
700  }
701  else
702  {
703  ErrorCode = MAX2(ErrorCode,QH_WIDE_FAIL);
704  }
705 
706  if( ErrorCode == QH_CONV_FAIL )
707  {
708  rel_tol *= 0.9;
709  }
710  }
711  else if( ErrorCode == QH_LOW_FAIL )
712  {
713  double help1 = gv.bin[nd]->qtmin*sqrt(DefFac);
714  double help2 = sqrt(gv.bin[nd]->qtmin*maxBracket);
715  minBracket = gv.bin[nd]->qtmin;
716  new_tmin = MIN2(help1,help2);
717  /* increase factor in case we get repeated LOW_FAIL's */
718  DefFac += 1.5;
719  }
720  else if( ErrorCode == QH_HIGH_FAIL )
721  {
722  double help = sqrt(gv.bin[nd]->qtmin*minBracket);
723  maxBracket = gv.bin[nd]->qtmin;
724  new_tmin = MAX2(gv.bin[nd]->qtmin/DEF_FAC,help);
725  }
726  else
727  {
728  new_tmin = sqrt(minBracket*maxBracket);
729  }
730  }
731  else
732  {
733  GetProbDistr_HighLimit(nd,STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&new_tmin,
734  &ErrorCode);
735  }
736 
737  gv.bin[nd]->qtmin = new_tmin;
738 
739  lgOK = ( ErrorCode < QH_RETRY );
740 
741  if( ErrorCode >= QH_FATAL )
742  break;
743 
744  if( ErrorCode != QH_LOW_FAIL )
745  DefFac = DEF_FAC;
746 
747  if( trace.lgTrace && trace.lgDustBug )
748  {
749  fprintf( ioQQQ, " GetProbDistr returns code %d\n", ErrorCode );
750  if( !lgOK )
751  {
752  fprintf( ioQQQ, " >>>>qheat trying another iteration, qtmin bracket %.4e %.4e",
753  minBracket,maxBracket );
754  fprintf( ioQQQ, " nWideFail %ld\n", nWideFail );
755  }
756  }
757  }
758 
759  if( ErrorCode == QH_WIDE_FAIL )
760  gv.bin[nd]->lgQHTooWide = true;
761 
762  /* >>chng 03 jan 17, added test for !lgDelta, PvH */
763  /* if( gv.bin[nd]->lgQHTooWide ) */
764  if( gv.bin[nd]->lgQHTooWide && !lgDelta )
765  ErrorCode = MAX2(ErrorCode,QH_WIDE_FAIL);
766 
767 /* if( ErrorCode >= QH_RETRY ) */
768 /* printf( " *** PROBLEM loop not converged, errorcode %d\n",ErrorCode ); */
769 
770  /* The quantum heating code tends to run into trouble when it goes deep into the neutral zone,
771  * especially if the original spectrum was very hard, as is the case in high excitation PNe or AGN.
772  * You then get a bipartition in the spectrum where most of the photons have low energy, while
773  * there still are hard X-ray photons left. The fact that the average energy per photon is low
774  * forces the quantum code to take tiny little steps when integrating the probability distribution,
775  * while the fact that X-ray photons are still present implies that big temperature spikes still
776  * occur and hence the temperature distribution is very wide. Therefore the code needs a zillion
777  * steps to integrate the probability distribution and simply runs out of room. As a last resort
778  * try the analytic approximation with relaxed constraints used below. */
779  /* >>chng 02 oct 03, vary Umax and fwhm to force fit with fwhm/Umax remaining constant */
780  /* >>chng 03 jan 17, changed test so that last resort always gets tried when lgOK = lgDelta = false, PvH */
781  /* if( !lgOK && FwhmRatio >= FWHM_RATIO && NumEvents <= MAX_EVENTS ) */
782  if( !lgOK && !lgDelta )
783  {
784  double Umax2 = Umax*sqrt(tol);
785  double fwhm2 = fwhm*sqrt(tol);
786 
787  for( i=0; i < LOOP_MAX; ++i )
788  {
789  double dummy;
790 
791  ErrorCode2 = MAX2(ErrorStart,QH_ANALYTIC);
792  GetProbDistr_HighLimit(nd,RELAX,Umax2,fwhm2,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
793  &ErrorCode2);
794 
795  lgOK = ( ErrorCode2 < QH_RETRY );
796  if( lgOK )
797  {
798  gv.bin[nd]->qtmin = dummy;
799  ErrorCode = ErrorCode2;
800  break;
801  }
802  else
803  {
804  Umax2 *= sqrt(tol);
805  fwhm2 *= sqrt(tol);
806  }
807  }
808  }
809 
810  if( nzone == 1 )
811  gv.bin[nd]->qtmin_zone1 = gv.bin[nd]->qtmin;
812 
813  gv.bin[nd]->lgUseQHeat = ( lgOK && !lgDelta );
814  gv.bin[nd]->lgEverQHeat = ( gv.bin[nd]->lgEverQHeat || gv.bin[nd]->lgUseQHeat );
815 
816  if( lgOK )
817  {
818  if( trace.lgTrace && trace.lgDustBug )
819  fprintf( ioQQQ, " >>>>qheat converged with code: %d\n", ErrorCode );
820  }
821  else
822  {
823  *qnbin = 0;
824  ++gv.bin[nd]->QHeatFailures;
825  fprintf( ioQQQ, " PROBLEM qheat did not converge grain %s in zone %ld, error code = %d\n",
826  gv.bin[nd]->chDstLab,nzone,ErrorCode );
827  }
828 
829  if( gv.QHSaveFile != NULL && ( iterations.lgLastIt || !gv.lgQHPunLast ) )
830  {
831  fprintf( gv.QHSaveFile, "\nDust Temperature Distribution: grain %s zone %ld\n",
832  gv.bin[nd]->chDstLab,nzone );
833 
834  fprintf( gv.QHSaveFile, "Equilibrium temperature: %.2f\n", gv.bin[nd]->tedust );
835 
836  if( gv.bin[nd]->lgUseQHeat )
837  {
838  /* >>chng 01 oct 09, remove qprob from output, it depends on step size, PvH */
839  fprintf( gv.QHSaveFile, "Number of bins: %ld\n", *qnbin );
840  fprintf( gv.QHSaveFile, " Tgrain dP/dlnT\n" );
841  for( i=0; i < *qnbin; i++ )
842  {
843  fprintf( gv.QHSaveFile, "%.4e %.4e\n", qtemp[i],dPdlnT[i] );
844  }
845  }
846  else
847  {
848  fprintf( gv.QHSaveFile, "**** quantum heating was not used\n" );
849  }
850  }
851  return;
852 }
853 
854 
855 /* initialize phiTilde */
856 STATIC void qheat_init(size_t nd,
857  /*@out@*/ vector<double>& phiTilde, /* phiTilde[rfield.nflux_with_check] */
858  /*@out@*/ double *check)
859 {
860  long i,
861  nz;
862  double sum = 0.;
863 
864  /*@-redef@*/
865  enum {DEBUG_LOC=false};
866  /*@+redef@*/
867 
868  DEBUG_ENTRY( "qheat_init()" );
869 
870  /* sanity checks */
871  ASSERT( gv.bin[nd]->lgQHeat );
872  ASSERT( nd < gv.bin.size() );
873 
874  *check = 0.;
875 
876  /* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
877  /* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
878  for( i=0; i < gv.bin[nd]->qnflux; i++ )
879  {
880  phiTilde[i] = 0.;
881  }
882 
883  /* fill in auxiliary array for quantum heating routine
884  * it reshuffles the input spectrum times the cross section to take
885  * the photo-electric effect into account. this prevents the quantum
886  * heating routine from having to calculate this effect over and over
887  * again; it can do a straight integration instead, making the code
888  * a lot simpler and faster. this initializes the array for non-ionizing
889  * energies, the reshuffling for higher energies is done in the next loop
890  * phiTilde has units events/H/s/cell at default depletion */
891 
892  double NegHeatingRate = 0.;
893 
894  for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
895  {
896  double check1 = 0.;
897  ChargeBin *gptr = gv.bin[nd]->chrg[nz];
898 
899  // rfield.nPositive may have increased since the last call to GrainDrive()
900  // if so, arrays like gptr->fac1 would not be initialized up to nPositive
901  long limit = min( rfield.nPositive, gptr->nfill );
902 
903  /* integrate over incident continuum for non-ionizing energies */
904  for( i=0; i < min(gptr->ipThresInf,limit); i++ )
905  {
906  check1 += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*rfield.anu(i);
907  phiTilde[i] += gptr->FracPop*rfield.SummedCon[i]*gv.bin[nd]->dstab1[i];
908  }
909 
910  /* >>chng 01 mar 02, use new expressions for grain cooling and absorption
911  * cross sections following the discussion with Joe Weingartner, PvH */
912  for( i=gptr->ipThresInf; i < limit; i++ )
913  {
914  long ipLo2 = gptr->ipThresInfVal;
915  double cs1 = ( i >= ipLo2 ) ? gv.bin[nd]->dstab1[i]*gptr->yhat_primary[i] : 0.;
916 
917  check1 += rfield.SummedCon[i]*gptr->fac1[i];
918  /* this accounts for the photons that are fully absorbed by grain */
919  phiTilde[i] += gptr->FracPop*rfield.SummedCon[i]*MAX2(gv.bin[nd]->dstab1[i]-cs1,0.);
920 
921  /* >>chng 01 oct 10, use bisection search to find ip. On C scale now */
922 
923  /* this accounts for photons that eject an electron from the valence band */
924  if( cs1 > 0. )
925  {
926  /* we treat photo-ejection and all subsequent de-excitation cascades
927  * from the conduction/valence band as one simultaneous event */
928  /* the condition cs1 > 0. assures i >= ipLo2 */
929  /* ratio is number of ejected electrons per primary ionization */
930  double ratio = ( gv.lgWD01 ) ? 1. : gptr->yhat[i]/gptr->yhat_primary[i];
931  /* ehat is average energy of ejected electron at infinity */
932  double ehat = gptr->ehat[i];
933  double cool1, sign = 1.;
934  realnum xx;
935 
936  if( gptr->DustZ <= -1 )
937  cool1 = gptr->ThresSurf + gptr->PotSurf + ehat;
938  else
939  cool1 = gptr->ThresSurfVal + gptr->PotSurf + ehat;
940  /* secondary electrons are assumed to have the same Elo and Ehi as the
941  * primary electrons that excite them. This neglects the threshold for
942  * the ejection of the secondary electron and can cause xx to become
943  * negative if Ehi is less than that threshold. To conserve energy we
944  * will simply assume a negative rate here. Since secondary electrons
945  * are generally not important this should have little impact on the
946  * overall temperature distribution */
947  xx = rfield.anu(i) - (realnum)(ratio*cool1);
948  if( xx < 0.f )
949  {
950  xx = -xx;
951  sign = -1.;
952  }
953  long ipLo = rfield.ipointC( max(xx,rfield.emm()) );
954  /* for grains in hard X-ray environments, the coarseness of the grid can
955  * lead to inaccuracies in the integral over phiTilde that would trip the
956  * sanity check in qheat(), here we correct for the energy mismatch */
957  double corr = xx/rfield.anu(ipLo);
958  phiTilde[ipLo] += sign*corr*gptr->FracPop*rfield.SummedCon[i]*cs1;
959  }
960 
961  /* no need to account for photons that eject an electron from the conduction band */
962  /* >>chng 01 dec 11, cool2 always equal to rfield.anu(i) -> no grain heating */
963  }
964 
965  *check += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
966 
967  sum += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
968 
969  if( DEBUG_LOC )
970  {
971  double integral = 0.;
972  for( i=0; i < gv.bin[nd]->qnflux; i++ )
973  {
974  integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu(i)*EN1RYD;
975  }
976  dprintf( ioQQQ, " integral test 1: integral %.6e %.6e\n", integral, sum );
977  }
978 
979  /* add quantum heating due to recombination of electrons, subtract thermionic cooling */
980 
981  /* gptr->HeatingRate2 is net heating rate in erg/H/s at standard depl
982  * includes contributions for recombining electrons, autoionizing electrons
983  * subtracted by thermionic emissions here since it is inverse process
984  *
985  * NB - in extreme conditions this rate may be negative (if there
986  * is an intense radiation field leading to very hot grains, but no ionizing
987  * photons, hence very few free electrons). we assume that the photon rates
988  * are high enough under those circumstances to avoid phiTilde becoming negative,
989  * but we will check that in qheat1 anyway. */
990 
991  /* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, pah_crash.in, PvH */
992  if( gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
993  {
994  double Sum,ESum,DSum,E_av2,Corr;
995  double fac = BOLTZMANN/EN1RYD*phycon.te;
996  /* E0 is barrier that electron needs to overcome, zero for positive grains */
997  /* >>chng 03 jan 23, added second term to correctly account for autoionizing states
998  * where ThresInfInc is negative, tested in small_grain.in, PvH */
999  double E0 = -(MIN2(gptr->PotSurfInc,0.) + MIN2(gptr->ThresInfInc,0.));
1000  /* >>chng 01 mar 02, this should be energy gap between top electron and infinity, PvH */
1001  /* >>chng 01 nov 21, use correct barrier: ThresInf[nz] -> ThresInfInc[nz], PvH */
1002  /* >>chng 03 jan 23, replaced -E0 with MIN2(PotSurfInc[nz],0.), PvH */
1003  double Einf = gptr->ThresInfInc + MIN2(gptr->PotSurfInc,0.);
1004  /* this is average energy deposited by one event, in erg
1005  * this value is derived from distribution assumed here, and may
1006  * not be the same as HeatElectrons/(CollisionRateElectr*eta) !! */
1007  /* >>chng 01 nov 21, use correct barrier: ThresInf[nz] -> ThresInfInc[nz], PvH */
1008  /* >>chng 03 jan 23, replaced ThresInfInc[nz] with MAX2(ThresInfInc[nz],0.), PvH */
1009  double E_av = MAX2(gptr->ThresInfInc,0.)*EN1RYD + 2.*BOLTZMANN*phycon.te;
1010  /* this is rate in events/H/s at standard depletion */
1011  double rate = gptr->HeatingRate2/E_av;
1012 
1013  double ylo = -exp(-E0/fac);
1014  /* this is highest kinetic energy of electron that can be represented in phiTilde */
1015  /* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
1016  /* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
1017  double Ehi = rfield.anumax(gv.bin[nd]->qnflux-1)-Einf;
1018  double yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
1019  /* renormalize rate so that integral over phiTilde*anu gives correct total energy */
1020  rate /= yhi-ylo;
1021 
1022  /* correct for fractional population of this charge state */
1023  rate *= gptr->FracPop;
1024 
1025  /* >>chng 03 jan 24, add code to correct for discretization errors, hotdust.in, PvH */
1026  vector<double> RateArr(gv.bin[nd]->qnflux);
1027  Sum = ESum = DSum = 0.;
1028 
1029  /* >>chng 04 jan 21, replaced gv.bin[nd]->qnflux -> gv.bin[nd]->qnflux2, PvH */
1030  for( i=0; i < gv.bin[nd]->qnflux2; i++ )
1031  {
1032  Ehi = rfield.anumax(i) - Einf;
1033  if( Ehi >= E0 )
1034  {
1035  /* Ehi is kinetic energy of electron at infinity */
1036  yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
1037  /* >>chng 01 mar 24, use MAX2 to protect against roundoff error, PvH */
1038  RateArr[i] = rate*MAX2(yhi-ylo,0.);
1039  Sum += RateArr[i];
1040  ESum += rfield.anu(i)*RateArr[i];
1041 # ifndef NDEBUG
1042  DSum += rfield.widflx(i)*RateArr[i];
1043 # endif
1044  ylo = yhi;
1045  }
1046  else
1047  {
1048  RateArr[i] = 0.;
1049  }
1050  }
1051  E_av2 = ESum/Sum*EN1RYD;
1052  ASSERT( fabs(E_av-E_av2) <= DSum/Sum*EN1RYD );
1053  Corr = E_av/E_av2;
1054 
1055  for( i=0; i < gv.bin[nd]->qnflux2; i++ )
1056  {
1057  phiTilde[i] += RateArr[i]*Corr;
1058  }
1059 
1060  sum += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
1061 
1062  if( DEBUG_LOC )
1063  {
1064  double integral = 0.;
1065  for( i=0; i < gv.bin[nd]->qnflux; i++ )
1066  {
1067  integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu(i)*EN1RYD;
1068  }
1069  dprintf( ioQQQ, " integral test 2: integral %.6e %.6e\n", integral, sum );
1070  }
1071  }
1072  else
1073  {
1074  NegHeatingRate += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
1075  }
1076  }
1077 
1078  /* ============================================================================= */
1079 
1080  /* add quantum heating due to molecule/ion collisions */
1081 
1082  /* gv.bin[nd]->HeatingRate1 is heating rate in erg/H/s at standard depl
1083  * includes contributions from molecules/neutral atoms and recombining ions
1084  *
1085  * in fully ionized conditions electron heating rates will be much higher
1086  * than ion and molecule rates since electrons are so much faster and grains
1087  * tend to be positive. in non-ionized conditions the main contribution will
1088  * come from neutral atoms and molecules, so it is appropriate to treat both
1089  * the same. in fully ionized conditions we don't care since unimportant.
1090  *
1091  * NB - if grains are hotter than ambient gas, the heating rate may become negative.
1092  * if photon rates are not high enough to prevent phiTilde from becoming negative,
1093  * we will raise a flag while calculating the quantum heating in qheat1 */
1094 
1095  /* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, PvH */
1096  if( gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
1097  {
1098  /* limits for Taylor expansion of (1+x)*exp(-x) */
1099  /* these choices will assure only 6 digits precision */
1100  const double LIM2 = cbrt(3.e-6);
1101  const double LIM3 = powpq(8.e-6,1,4);
1102  /* if gas temperature is higher than grain temperature we will
1103  * consider Maxwell-Boltzmann distribution of incoming particles
1104  * and ignore distribution of outgoing particles, if grains
1105  * are hotter than ambient gas, we use reverse treatment */
1106  double fac = BOLTZMANN/EN1RYD*MAX2(phycon.te,gv.bin[nd]->tedust);
1107  /* this is average energy deposited/extracted by one event, in erg */
1108  double E_av = 2.*BOLTZMANN*MAX2(phycon.te,gv.bin[nd]->tedust);
1109  /* this is rate in events/H/s at standard depletion */
1110  double rate = gv.bin[nd]->HeatingRate1/E_av;
1111 
1112  double ylo = -1.;
1113  /* this is highest energy of incoming/outgoing particle that can be represented in phiTilde */
1114  /* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
1115  /* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
1116  double Ehi = rfield.anumax(gv.bin[nd]->qnflux-1);
1117  double yhi = -(Ehi/fac+1.)*exp(-Ehi/fac);
1118  /* renormalize rate so that integral over phiTilde*anu gives correct total energy */
1119  rate /= yhi-ylo;
1120 
1121  for( i=0; i < gv.bin[nd]->qnflux2; i++ )
1122  {
1123  /* Ehi is kinetic energy of incoming/outgoing particle
1124  * we assume that Ehi-E0 is deposited/extracted from grain */
1125  /* Ehi = rfield.anumax(i); */
1126  double x = rfield.anumax(i)/fac;
1127  /* (1+x)*exp(-x) = 1 - 1/2*x^2 + 1/3*x^3 - 1/8*x^4 + O(x^5)
1128  * = 1 - Sum_n=2^infty (-x)^n/(n*(n-2)!) */
1129  if( x > LIM3 )
1130  yhi = -(x+1.)*exp(-x);
1131  else if( x > LIM2 )
1132  yhi = -(((1./3.)*x - 0.5)*x*x + 1.);
1133  else
1134  yhi = -(1. - 0.5*x*x);
1135 
1136  /* >>chng 01 mar 24, use MAX2 to protect against roundoff error, PvH */
1137  phiTilde[i] += rate*MAX2(yhi-ylo,0.);
1138  ylo = yhi;
1139  }
1140 
1141  sum += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
1142 
1143  if( DEBUG_LOC )
1144  {
1145  double integral = 0.;
1146  for( i=0; i < gv.bin[nd]->qnflux; i++ )
1147  {
1148  integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu(i)*EN1RYD;
1149  }
1150  dprintf( ioQQQ, " integral test 3: integral %.6e %.6e\n", integral, sum );
1151  }
1152  }
1153  else
1154  {
1155  NegHeatingRate += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
1156  }
1157 
1158  /* here we account for the negative heating rates, we simply do that by scaling the entire
1159  * phiTilde array down by a constant factor such that the total amount of energy is conserved
1160  * This treatment assures that phiTilde never goes negative, which avoids problems further on */
1161  if( NegHeatingRate < 0. )
1162  {
1163  double scale_fac = (sum+NegHeatingRate)/sum;
1164  for( i=0; i < gv.bin[nd]->qnflux; i++ )
1165  phiTilde[i] *= scale_fac;
1166 
1167  if( DEBUG_LOC )
1168  {
1169  sum += NegHeatingRate;
1170 
1171  double integral = 0.;
1172  for( i=0; i < gv.bin[nd]->qnflux; i++ )
1173  {
1174  integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu(i)*EN1RYD;
1175  }
1176  dprintf( ioQQQ, " integral test 4: integral %.6e %.6e\n", integral, sum );
1177  }
1178  }
1179 
1180  return;
1181 }
1182 
1183 
1184 /*******************************************************************************************
1185  *
1186  * GetProbDistr_LowLimit: main routine for calculating non-equilibrium heating of grains
1187  *
1188  * This routine implements the algorithm outlined in:
1189  * >>refer grain physics Guhathakurtha & Draine, 1989, ApJ, 345, 230
1190  *
1191  * The original (fortran) version of the code was written by Kevin Volk.
1192  *
1193  * Heavily modified and adapted for new style grains by Peter van Hoof.
1194  *
1195  *******************************************************************************************/
1196 
1197 STATIC void GetProbDistr_LowLimit(size_t nd,
1198  double rel_tol,
1199  double Umax,
1200  double fwhm,
1201  /*@in@*/ const vector<double>& Phi, /* Phi[NQGRID] */
1202  /*@in@*/ const vector<double>& PhiDrv, /* PhiDrv[NQGRID] */
1203  /*@out@*/ vector<double>& qtemp, /* qtemp[NQGRID] */
1204  /*@out@*/ vector<double>& qprob, /* qprob[NQGRID] */
1205  /*@out@*/ vector<double>& dPdlnT, /* dPdlnT[NQGRID] */
1206  /*@out@*/ long int *qnbin,
1207  /*@out@*/ double *new_tmin,
1208  long *nWideFail,
1209  /*@out@*/ QH_Code *ErrorCode)
1210 {
1211  bool lgAllNegSlope,
1212  lgBoundErr;
1213  long int j,
1214  k,
1215  l,
1216  nbad,
1217  nbin,
1218  nend=0,
1219  nmax,
1220  nok,
1221  nstart=0,
1222  nstart2=0;
1223  double dCool=0.,
1224  dlnLdlnT,
1225  dlnpdlnU,
1226  fac = 0.,
1227  maxVal,
1228  NextStep,
1229  qtmin1,
1230  RadCooling,
1231  sum,
1232  y;
1233  vector<double> delu(NQGRID);
1234  vector<double> Lambda(NQGRID);
1235  vector<double> p(NQGRID);
1236  vector<double> u1(NQGRID);
1237 
1238 
1239  DEBUG_ENTRY( "GetProbDistr_LowLimit()" );
1240 
1241  /* sanity checks */
1242  ASSERT( nd < gv.bin.size() );
1243 
1244  if( trace.lgTrace && trace.lgDustBug )
1245  {
1246  fprintf( ioQQQ, " GetProbDistr_LowLimit called for grain %s\n", gv.bin[nd]->chDstLab );
1247  fprintf( ioQQQ, " got qtmin1 %.4e\n", gv.bin[nd]->qtmin);
1248  }
1249 
1250  qtmin1 = gv.bin[nd]->qtmin;
1251  qtemp[0] = qtmin1;
1252  /* u1 holds enthalpy in Ryd/grain */
1253  u1[0] = ufunct(qtemp[0],nd,&lgBoundErr);
1254  ASSERT( !lgBoundErr ); /* this should never happen */
1255  /* >>chng 00 mar 22, taken out factor 4, factored in hden and dstAbund
1256  * interpolate in dstems array instead of integrating explicitly, by PvH */
1257  splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[0]),&y,&lgBoundErr);
1258  ASSERT( !lgBoundErr ); /* this should never happen */
1259  /* Lambda holds the radiated energy for grains in this bin, in Ryd/s/grain */
1260  Lambda[0] = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
1261  /* set up first step of integration */
1262  /* >>chng 01 nov 14, set to 2.*Lambda[0]/Phi[0] instead of u1[0],
1263  * this choice assures that p[1] doesn't make a large jump from p[0], PvH */
1264  delu[0] = 2.*Lambda[0]/Phi[0];
1265  p[0] = PROB_CUTOFF_LO;
1266  dPdlnT[0] = p[0]*qtemp[0]*uderiv(qtemp[0],nd);
1267  RadCooling = 0.5*p[0]*Lambda[0]*delu[0];
1268  NextStep = 0.01*Lambda[0]/Phi[0];
1269  /* >>chng 03 nov 10, added extra safeguard against stepsize too small, PvH */
1270  if( NextStep < 5.*DBL_EPSILON*u1[0] )
1271  {
1272  *ErrorCode = MAX2(*ErrorCode,QH_STEP_FAIL);
1273  return;
1274  }
1275 
1276  nbad = 0;
1277  k = 0;
1278 
1279  *qnbin = 0;
1280  *new_tmin = qtmin1;
1281  lgAllNegSlope = true;
1282  maxVal = dPdlnT[0];
1283  nmax = 0;
1284  double p_max = p[0];
1285 
1286  /* this test neglects a negative contribution which is impossible to calculate, so it may
1287  * fail to detect cases where the probability distribution starts dropping immediately,
1288  * we will use a second test using the variable lgAllNegSlope below to catch those cases, PvH */
1289  spldrv_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[0]),&dlnLdlnT,&lgBoundErr);
1290  ASSERT( !lgBoundErr ); /* this should never happen */
1291  dlnpdlnU = u1[0]*Phi[0]/Lambda[0] - dlnLdlnT*u1[0]/(qtemp[0]*uderiv(qtemp[0],nd));
1292  if( dlnpdlnU < 0. )
1293  {
1294  /* >>chng 03 nov 06, confirm this by integrating first step..., pah_crash.in, PvH */
1295  (void)TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,p_max,k,nd,&lgBoundErr);
1296  dPdlnT[2] = p[2]*qtemp[2]*uderiv(qtemp[2],nd);
1297 
1298  if( dPdlnT[2] < dPdlnT[0] )
1299  {
1300  /* if dPdlnT starts falling immediately,
1301  * qtmin1 was too high and convergence is impossible */
1302  *ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);
1303  return;
1304  }
1305  }
1306 
1307  /* NB NB -- every break in this loop should set *ErrorCode (except for regular stop criterion) !! */
1308  for( l=0; l < MAX_LOOP; ++l )
1309  {
1310  double rerr = TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,p_max,k,nd,&lgBoundErr);
1311 
1312  /* this happens if the grain temperature in qtemp becomes higher than GRAIN_TMAX
1313  * nothing that TryDoubleStep returns can be trusted, so this check should be first */
1314  if( lgBoundErr )
1315  {
1316  nbad += 2;
1317  *ErrorCode = MAX2(*ErrorCode,QH_THIGH_FAIL);
1318  break;
1319  }
1320 
1321  /* estimate new stepsize */
1322  if( rerr > rel_tol*QHEAT_TOL )
1323  {
1324  nbad += 2;
1325 
1326  /* step is rejected, decrease stepsize and try again */
1327  NextStep *= sqrt(0.9*rel_tol*QHEAT_TOL/rerr);
1328 
1329  /* stepsize too small, this can happen at extreme low temperatures */
1330  if( NextStep < 5.*DBL_EPSILON*u1[k] )
1331  {
1332  *ErrorCode = MAX2(*ErrorCode,QH_STEP_FAIL);
1333  break;
1334  }
1335 
1336  continue;
1337  }
1338  else
1339  {
1340  /* step was OK, adjust stepsize */
1341  k += 2;
1342 
1343  p_max = max(p_max,p[k-1]);
1344  p_max = max(p_max,p[k]);
1345 
1346  /* >>chng 03 nov 10, safeguard against division by zero, PvH */
1347  NextStep *= MIN2(cbrt(0.9*rel_tol*QHEAT_TOL/MAX2(rerr,1.e-50)),4.);
1348  NextStep = MIN2(NextStep,Lambda[k]/Phi[0]);
1349  }
1350 
1351  dPdlnT[k-1] = p[k-1]*qtemp[k-1]*uderiv(qtemp[k-1],nd);
1352  dPdlnT[k] = p[k]*qtemp[k]*uderiv(qtemp[k],nd);
1353 
1354  lgAllNegSlope = lgAllNegSlope && ( dPdlnT[k] < dPdlnT[k-2] );
1355 
1356  if( dPdlnT[k-1] > maxVal )
1357  {
1358  maxVal = dPdlnT[k-1];
1359  nmax = k-1;
1360  }
1361  if( dPdlnT[k] > maxVal )
1362  {
1363  maxVal = dPdlnT[k];
1364  nmax = k;
1365  }
1366 
1367  RadCooling += dCool;
1368 
1369 // if( nzone >= 24 && nd == 0 ) {
1370 // printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k-1,qtemp[k-1],u1[k-1],p[k-1],dPdlnT[k-1]);
1371 // printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k,qtemp[k],u1[k],p[k],dPdlnT[k]);
1372 // }
1373 
1374  /* if qtmin is far too low, p[k] can become extremely large, exceeding
1375  * even double precision range. the following check should prevent overflows */
1376  /* >>chng 01 nov 07, sqrt(DBL_MAX) -> sqrt(DBL_MAX/100.) so that sqrt(p[k]*p[k+1]) is safe */
1377  if( p[k] > sqrt(DBL_MAX/100.) )
1378  {
1379  *ErrorCode = MAX2(*ErrorCode,QH_LOW_FAIL);
1380  break;
1381  }
1382 
1383  /* this may catch a bug in the Compaq C compiler V6.3-025
1384  * if this gets triggered, try compiling with -ieee */
1385  ASSERT( p[k] > 0. && dPdlnT[k] > 0. && RadCooling > 0. );
1386 
1387  /* do a check for negative slope and if there will be enough room to store results */
1388  if( k > 0 && k%NQTEST == 0 )
1389  {
1390  double wid, avStep, factor;
1391  /* >>chng 02 jul 31, do additional test for HIGH_FAIL,
1392  * first test before main loop doesn't catch everything, PvH */
1393  if( lgAllNegSlope )
1394  {
1395  *ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);
1396  break;
1397  }
1398 
1399  /* this is a lower limit for the total width of the probability distr */
1400  /* >>chng 02 jan 30, corrected calculation of wid and avStep, PvH */
1401  wid = (sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO)) +
1402  sqrt(-log(PROB_CUTOFF_HI)/(4.*LN_TWO)))*fwhm/Umax;
1403  avStep = log(u1[k]/u1[0])/(double)k;
1404  /* make an estimate for the number of steps needed */
1405  /* >>chng 02 jan 30, included factor 1.5 because stepsize increases near peak, PvH */
1406  /* >>chng 02 aug 06, changed 1.5 to sliding scale because of laqheat2.in test, PvH */
1407  factor = 1.1 + 3.9*(1.0 - sqrt((double)k/(double)NQGRID));
1408  if( wid/avStep > factor*(double)NQGRID )
1409  {
1410  *ErrorCode = MAX2(*ErrorCode,QH_ARRAY_FAIL);
1411  break;
1412  }
1413  }
1414 
1415  /* if we run out of room to store results, do regular break
1416  * the code below will sort out if integration is valid or not */
1417  if( k >= NQGRID-2 )
1418  {
1419  *ErrorCode = MAX2(*ErrorCode,QH_ARRAY_FAIL);
1420  break;
1421  }
1422 
1423  /* force thermal equilibrium of the grains */
1424  fac = RadCooling*gv.bin[nd]->cnv_GR_pCM3*EN1RYD/gv.bin[nd]->GrainHeat;
1425 
1426  /* this is regular stop criterion */
1427  if( dPdlnT[k] < dPdlnT[k-1] && dPdlnT[k]/fac < PROB_CUTOFF_HI )
1428  {
1429  break;
1430  }
1431  }
1432 
1433  if( l == MAX_LOOP )
1434  *ErrorCode = MAX2(*ErrorCode,QH_LOOP_FAIL);
1435 
1436  nok = k;
1437  nbin = k+1;
1438 
1439  /* there are insufficient bins to attempt rebinning */
1440  if( *ErrorCode < QH_NO_REBIN && nbin < NQMIN )
1441  *ErrorCode = MAX2(*ErrorCode,QH_NBIN_FAIL);
1442 
1443  /* >>chng 02 aug 07, do some basic checks on the distribution first */
1444  if( *ErrorCode < QH_NO_REBIN )
1445  ScanProbDistr(u1,dPdlnT,nbin,maxVal,nmax,qtmin1,&nstart,&nstart2,&nend,nWideFail,ErrorCode);
1446 
1447  if( *ErrorCode >= QH_NO_REBIN )
1448  {
1449  return;
1450  }
1451 
1452  for( j=0; j < nbin; j++ )
1453  {
1454  p[j] /= fac;
1455  dPdlnT[j] /= fac;
1456  }
1457  RadCooling /= fac;
1458 
1459  /* >>chng 02 aug 08, moved RebinQHeatResults from here further down, this improves new_tmin estimate */
1460  *new_tmin = 0.;
1461  for( j=nstart; j < nbin; j++ )
1462  {
1463  if( dPdlnT[j] < PROB_CUTOFF_LO )
1464  {
1465  *new_tmin = qtemp[j];
1466  }
1467  else
1468  {
1469  if( j == nstart )
1470  {
1471  /* if dPdlnT[nstart] is too high, but qtmin1 is already close to GRAIN_TMIN,
1472  * then don't bother signaling a QH_BOUND_FAIL. grains below GRAIN_TMIN have the
1473  * peak of their thermal emission beyond 3 meter, so they really are irrelevant
1474  * since free-free emission from electrons will drown this grain emission... */
1475  if( dPdlnT[j] > SAFETY*PROB_CUTOFF_LO && qtmin1 > 1.2*GRAIN_TMIN )
1476  *ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
1477 
1478  /* >>chng 02 aug 11, use nstart2 for more reliable extrapolation */
1479  if( dPdlnT[nstart2] < 0.999*dPdlnT[nstart2+NSTARTUP] )
1480  {
1481  /* >>chng 02 aug 09, new formula for extrapolating new_tmin, PvH */
1482  /* this assumes that at low temperatures the behaviour
1483  * is as follows: dPdlnT(T) = C1 * exp( -C2/T**3 ) */
1484  double T1 = qtemp[nstart2];
1485  double T2 = qtemp[nstart2+NSTARTUP];
1486  double delta_y = log(dPdlnT[nstart2+NSTARTUP]/dPdlnT[nstart2]);
1487  double c2 = delta_y/(1./pow3(T1)-1./pow3(T2));
1488  double help = c2/pow3(T1) + log(dPdlnT[nstart2]/PROB_CUTOFF_LO);
1489  *new_tmin = cbrt(c2/help);
1490  }
1491 
1492  /* >>chng 04 nov 09, in case of lower bound failure, assure qtmin is lowered, PvH */
1493  if( dPdlnT[j] > SAFETY*PROB_CUTOFF_LO && *new_tmin >= qtmin1 )
1494  {
1495  double delta_x = log(qtemp[nstart2+NSTARTUP]/qtemp[nstart2]);
1496  double delta_y = log(dPdlnT[nstart2+NSTARTUP]/dPdlnT[nstart2]);
1497  delta_x *= log(PROB_CUTOFF_LO/dPdlnT[nstart2])/delta_y;
1498  *new_tmin = qtemp[nstart2]*exp(delta_x);
1499  if( *new_tmin < qtmin1 )
1500  /* in general this estimate will be too low -> use geometric mean */
1501  *new_tmin = sqrt( *new_tmin*qtmin1 );
1502  else
1503  /* last resort... */
1504  *new_tmin = qtmin1/DEF_FAC;
1505  }
1506  }
1507  break;
1508  }
1509  }
1510  *new_tmin = MAX3(*new_tmin,qtmin1/DEF_FAC,GRAIN_TMIN);
1511 
1512  ASSERT( *new_tmin < gv.bin[nd]->tedust );
1513 
1514  /* >>chng 02 jan 30, prevent excessive looping when prob distribution simply won't fit, PvH */
1515  if( dPdlnT[nbin-1] > SAFETY*PROB_CUTOFF_HI )
1516  {
1517  if( *ErrorCode == QH_ARRAY_FAIL || *ErrorCode == QH_LOOP_FAIL )
1518  {
1519  ++(*nWideFail);
1520 
1521  if( *nWideFail < WIDE_FAIL_MAX )
1522  {
1523  /* this indicates that low end was OK, but we ran out of room
1524  * to store the high end -> try GetProbDistr_HighLimit instead */
1525  *ErrorCode = MAX2(*ErrorCode,QH_DELTA_FAIL);
1526  }
1527  else
1528  {
1529  *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
1530  }
1531  }
1532  else
1533  {
1534  *ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
1535  }
1536  }
1537 
1538  /* >>chng 01 may 11, rebin the quantum heating results
1539  *
1540  * for grains in intense radiation fields, the code needs high resolution for stability
1541  * and therefore produces lots of small bins, even though the grains never make large
1542  * excurions from the equilibrium temperature; adding in the resulting spectra in RT_diffuse
1543  * takes up an excessive amount of CPU time where most CPU is spent on grains for which
1544  * the quantum treatment matters least, and moreover on temperature bins with very low
1545  * probability; rebinning the results on a coarser grid should help reduce the overhead */
1546  /* >>chng 02 aug 07, use nstart and nend while rebinning */
1547 
1548  nbin = RebinQHeatResults(nd,nstart,nend,p,qtemp,qprob,dPdlnT,u1,delu,Lambda,ErrorCode);
1549 
1550  /* >>chng 01 jul 13, add fail-safe for failure in RebinQHeatResults */
1551  if( nbin == 0 )
1552  {
1553  return;
1554  }
1555 
1556  *qnbin = nbin;
1557 
1558  sum = 0.;
1559  for( j=0; j < nbin; j++ )
1560  {
1561  sum += qprob[j];
1562  }
1563 
1564  /* the fact that the probability normalization fails may indicate that the distribution is
1565  * too close to a delta function to resolve, but another possibility is that the radiation
1566  * field is extremely diluted allowing a sizable fraction of the grains to cool below
1567  * GRAIN_TMIN. In the latter case we don't raise QH_CONV_FAIL since these very cool grains
1568  * only contribute at very long radio wavelengths (more than 1 meter) */
1569  if( fabs(sum-1.) > PROB_TOL && qtmin1 > 1.2*GRAIN_TMIN )
1570  *ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);
1571 
1572  if( trace.lgTrace && trace.lgDustBug )
1573  {
1574  fprintf( ioQQQ,
1575  " zone %ld %s nbin %ld nok %ld nbad %ld total prob %.4e rel_tol %.3e new_tmin %.4e\n",
1576  nzone,gv.bin[nd]->chDstLab,nbin,nok,nbad,sum,rel_tol,*new_tmin );
1577  }
1578  return;
1579 }
1580 
1581 
1582 /* try two consecutive integration steps using stepsize "step/2." (yielding p[k]),
1583  * and also one double integration step using stepsize "step" (yielding p2k).
1584  * the difference fabs(p2k-p[k])/(3.*p[k]) can be used to estimate the relative
1585  * accuracy of p[k] and will be used to adapt the stepsize to an optimal value */
1586 STATIC double TryDoubleStep(vector<double>& u1,
1587  vector<double>& delu,
1588  vector<double>& p,
1589  vector<double>& qtemp,
1590  vector<double>& Lambda,
1591  const vector<double>& Phi,
1592  const vector<double>& PhiDrv,
1593  double step,
1594  /*@out@*/ double *cooling,
1595  double p_max,
1596  long index,
1597  size_t nd,
1598  /*@out@*/ bool *lgBoundFail)
1599 {
1600  long i,
1601  j,
1602  jlo,
1603  k=-1000;
1604  double bval_jk,
1605  cooling2,
1606  p2k,
1607  RelErrCool,
1608  RelErrPk,
1609  sum,
1610  sum2 = -DBL_MAX,
1611  trap1,
1612  trap12 = -DBL_MAX,
1613  trap2,
1614  uhi,
1615  ulo,
1616  umin,
1617  y;
1618 
1619  DEBUG_ENTRY( "TryDoubleStep()" );
1620 
1621  /* sanity checks */
1622  ASSERT( index >= 0 && index < NQGRID-2 && nd < gv.bin.size() && step > 0. );
1623 
1624  ulo = rfield.anumin(0);
1625  /* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
1626  /* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
1627  uhi = rfield.anumax(gv.bin[nd]->qnflux-1);
1628 
1629  /* >>chng 01 nov 21, skip initial bins if they have very low probability */
1630  jlo = 0;
1631  while( p[jlo] < PROB_CUTOFF_LO*p_max )
1632  jlo++;
1633 
1634  for( i=1; i <= 2; i++ )
1635  {
1636  bool lgErr;
1637  long ipLo = 0;
1638  // using gv.bin[nd]->qnflux is OK since gv.bin[nd]->qnflux < gv.bin[nd]->nflux_with_check
1639  long ipHi = gv.bin[nd]->qnflux;
1640  k = index + i;
1641  delu[k] = step/2.;
1642  u1[k] = u1[k-1] + delu[k];
1643  qtemp[k] = inv_ufunct(u1[k],nd,lgBoundFail);
1644  splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[k]),&y,&lgErr);
1645  *lgBoundFail = *lgBoundFail || lgErr;
1646  Lambda[k] = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
1647 
1648  sum = sum2 = 0.;
1649  trap1 = trap2 = trap12 = 0.;
1650 
1651  /* this loop uses up a large fraction of the total CPU time, it should be well optimized */
1652  for( j=jlo; j < k; j++ )
1653  {
1654  umin = u1[k] - u1[j];
1655 
1656  if( umin >= uhi )
1657  {
1658  /* for increasing j, umin will decrease monotonically. If ( umin > uhi ) holds for
1659  * the current value of j, it must have held for the previous value as well. Hence
1660  * both trap1 and trap2 are zero at this point and we would only be adding zero
1661  * to sum. Therefore we skip this step to save CPU time */
1662  continue;
1663  }
1664  else if( umin > ulo )
1665  {
1666  /* do a bisection search such that rfield.anumin(ipLo) <= umin < rfield.anumin(ipHi)
1667  * explicit bisection search is faster, which is important here to save CPU time.
1668  * on the first iteration ipLo equals 0 and the first while loop will be skipped;
1669  * after that umin is monotonically decreasing, and ipHi is retained from the
1670  * previous iteration since it is a valid upper limit; ipLo will equal ipHi-1 */
1671  long ipStep = 1;
1672  /* >>chng 03 feb 03 rjrw: hunt for lower bracket */
1673  while( rfield.anumin(ipLo) > umin )
1674  {
1675  ipHi = ipLo;
1676  ipLo -= ipStep;
1677  if( ipLo <= 0 )
1678  {
1679  ipLo = 0;
1680  break;
1681  }
1682  ipStep *= 2;
1683  }
1684  /* now do regular bisection search */
1685  while( ipHi-ipLo > 1 )
1686  {
1687  long ipMd = (ipLo+ipHi)/2;
1688  if( rfield.anumin(ipMd) > umin )
1689  ipHi = ipMd;
1690  else
1691  ipLo = ipMd;
1692  }
1693  /* Phi[i] is integral of PhiDrv from exactly rfield.anumin(i) to infinity */
1694  bval_jk = Phi[ipLo] + (umin - rfield.anumin(ipLo))*PhiDrv[ipLo];
1695  }
1696  else
1697  {
1698  bval_jk = Phi[0];
1699  }
1700 
1701  /* these two quantities are needed to take double step from index -> index+2 */
1702  trap12 = trap1;
1703  sum2 = sum;
1704 
1705  /* bval_jk*gv.bin[nd]->cnv_CM3_pGR is the total excitation rate from j to k and
1706  * higher due to photon absorptions and particle collisions, it already implements
1707  * Eq. 2.17 of Guhathakurtha & Draine, in events/grain/s */
1708  /* >>chng 00 mar 27, factored in hden (in Phi), by PvH */
1709  /* >>chng 00 mar 29, add in contribution due to particle collisions, by PvH */
1710  /* >>chng 01 mar 30, moved multiplication of bval_jk with gv.bin[nd]->cnv_CM3_pGR
1711  * out of loop, PvH */
1712  trap2 = p[j]*bval_jk;
1713  /* Trapezoidal rule, multiplication with factor 0.5 is done outside loop */
1714  sum += (trap1 + trap2)*delu[j];
1715  trap1 = trap2;
1716  }
1717 
1718  /* >>chng 00 mar 27, multiplied with delu, by PvH */
1719  /* >>chng 00 apr 05, taken out Lambda[0], improves convergence at low end dramatically!, by PvH */
1720  /* qprob[k] = sum*gv.bin[nd]->cnv_CM3_pGR*delu[k]/(Lambda[k] - Lambda[0]); */
1721  /* this expression includes factor 0.5 from trapezoidal rule above */
1722  /* p[k] = 0.5*(sum + (trap1 + p[k]*Phi[0])*delu[k])/Lambda[k] */
1723  p[k] = (sum + trap1*delu[k])/(2.*Lambda[k] - Phi[0]*delu[k]);
1724 
1725  // total failure -> force next step to be smaller
1726  if( p[k] <= 0. )
1727  return 3.*QHEAT_TOL;
1728  }
1729 
1730  /* this is estimate for p[k] using one double step of size "step" */
1731  p2k = (sum2 + trap12*step)/(2.*Lambda[k] - Phi[0]*step);
1732 
1733  // total failure -> force next step to be smaller
1734  if( p2k <= 0. )
1735  return 3.*QHEAT_TOL;
1736 
1737  RelErrPk = fabs(p2k-p[k])/p[k];
1738 
1739  /* this is radiative cooling due to the two probability bins we just added */
1740  /* simple trapezoidal rule will not do here, RelErrCool would never converge */
1741  double z[8];
1742  vlog(z, u1[k-2], u1[k-1], u1[k], p[k-2]*Lambda[k-2], p[k-1]*Lambda[k-1], p[k]*Lambda[k], p2k*Lambda[k], 1.);
1743  *cooling = log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k-1],p[k-1]*Lambda[k-1],z[0],z[3],z[1],z[4]);
1744  *cooling += log_integral(u1[k-1],p[k-1]*Lambda[k-1],u1[k],p[k]*Lambda[k],z[1],z[4],z[2],z[5]);
1745 
1746  /* same as cooling, but now for double step of size "step" */
1747  cooling2 = log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k],p2k*Lambda[k],z[0],z[3],z[2],z[6]);
1748 
1749  /* p[0] is not reliable, so ignore convergence test on cooling on first step */
1750  RelErrCool = ( index > 0 ) ? fabs(cooling2-(*cooling))/(*cooling) : 0.;
1751 
1752 // dprintf( ioQQQ, " TryDoubleStep k %ld p[k-1] %.4e p[k] %.4e p2k %.4e\n",k,p[k-1],p[k],p2k );
1753  /* error scales as O(step^3), so this is relative accuracy of p[k] or cooling */
1754  return MAX2(RelErrPk,RelErrCool)/3.;
1755 }
1756 
1757 
1758 /* calculate logarithmic integral from (xx1,yy1) to (xx2,yy2) */
1759 STATIC double log_integral(double xx1,
1760  double yy1,
1761  double xx2,
1762  double yy2,
1763  double log_xx1,
1764  double log_yy1,
1765  double log_xx2,
1766  double log_yy2)
1767 {
1768  DEBUG_ENTRY( "log_integral()" );
1769 
1770  double xx = log_xx2 - log_xx1;
1771  double eps = xx + log_yy2 - log_yy1;
1772  if( fabs(eps) < 1.e-4 )
1773  {
1774  return xx1*yy1*xx*(((eps/24. + 1./6.)*eps + 0.5)*eps + 1.);
1775  }
1776  else
1777  {
1778  return (xx2*yy2 - xx1*yy1)*xx/eps;
1779  }
1780 }
1781 
1782 
1783 /* scan the probability distribution for valid range */
1784 STATIC void ScanProbDistr(const vector<double>& u1, /* u1[nbin] */
1785  const vector<double>& dPdlnT, /* dPdlnT[nbin] */
1786  long nbin,
1787  double maxVal,
1788  long nmax,
1789  double qtmin1,
1790  /*@out@*/long *nstart,
1791  /*@out@*/long *nstart2,
1792  /*@out@*/long *nend,
1793  long *nWideFail,
1794  QH_Code *ErrorCode)
1795 {
1796  bool lgSetLo,
1797  lgSetHi;
1798  long i;
1799  double deriv_max,
1800  minVal;
1801  const double MIN_FAC_LO = 1.e4;
1802  const double MIN_FAC_HI = 1.e4;
1803 
1804  DEBUG_ENTRY( "ScanProbDistr()" );
1805 
1806  /* sanity checks */
1807  ASSERT( nbin > 0 && nmax >= 0 && nmax < nbin && maxVal > 0. );
1808 
1809  /* sometimes the probability distribution will start falling before settling on
1810  * a rising slope. In such a case nstart will point to the first rising point,
1811  * while nstart2 will point to the point with the steepest derivative on the
1812  * rising slope. The code will use the distribution from nstart to nend as a
1813  * valid probability distribution, but the will use the points near nstart2
1814  * to extrapolate a new value for qtmin if needed */
1815  minVal = maxVal;
1816  *nstart = nmax;
1817  for( i=nmax; i >= 0; --i )
1818  {
1819  if( dPdlnT[i] < minVal )
1820  {
1821  *nstart = i;
1822  minVal = dPdlnT[i];
1823  }
1824  }
1825  deriv_max = 0.;
1826  *nstart2 = nmax;
1827  for( i=nmax; i > *nstart; --i )
1828  {
1829  double deriv = log(dPdlnT[i]/dPdlnT[i-1])/log(u1[i]/u1[i-1]);
1830  if( deriv > deriv_max )
1831  {
1832  *nstart2 = i-1;
1833  deriv_max = deriv;
1834  }
1835  }
1836  *nend = nbin-1;
1837 
1838  /* now do quality control; these checks are more stringent than the ones in GetProbDistr_LowLimit */
1839  lgSetLo = ( nmax >= *nend || maxVal/dPdlnT[*nend] < MIN_FAC_HI );
1840  /* >>chng 03 jan 22, prevent problems if both dPdlnT and its derivative are continuously rising,
1841  * in which case both lgSetLo and lgSetHi are set and QH_WIDE_FAIL is triggered;
1842  * this can happen if qtmin is far too low; triggered by pahtest.in, PvH */
1843  if( lgSetLo )
1844  /* use relaxed test if lgSetLo is already set */
1845  lgSetHi = ( nmax <= *nstart || maxVal/dPdlnT[*nstart] < MIN_FAC_LO );
1846  else
1847  lgSetHi = ( nmax <= *nstart2 || maxVal/dPdlnT[*nstart2] < MIN_FAC_LO );
1848 
1849  if( lgSetLo && lgSetHi )
1850  {
1851  ++(*nWideFail);
1852 
1853  if( *nWideFail >= WIDE_FAIL_MAX )
1854  *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
1855  }
1856 
1857  if( lgSetLo )
1858  *ErrorCode = MAX2(*ErrorCode,QH_LOW_FAIL);
1859 
1860  /* if dPdlnT[nstart] is too high, but qtmin1 is already close to GRAIN_TMIN,
1861  * then don't bother signaling a QH_HIGH_FAIL. grains below GRAIN_TMIN have the
1862  * peak of their thermal emission beyond 3 meter, so they really are irrelevant
1863  * since free-free emission from electrons will drown this grain emission... */
1864  if( lgSetHi && qtmin1 > 1.2*GRAIN_TMIN )
1865  *ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);
1866 
1867  /* there are insufficient bins to attempt rebinning */
1868  if( *ErrorCode < QH_NO_REBIN && (*nend - *nstart) < NQMIN )
1869  *ErrorCode = MAX2(*ErrorCode,QH_NBIN_FAIL);
1870 
1871  if( trace.lgTrace && trace.lgDustBug )
1872  {
1873  fprintf( ioQQQ, " ScanProbDistr nstart %ld nstart2 %ld nend %ld nmax %ld maxVal %.3e",
1874  *nstart,*nstart2,*nend,nmax,maxVal );
1875  fprintf( ioQQQ, " dPdlnT[nstart] %.3e dPdlnT[nstart2] %.3e dPdlnT[nend] %.3e code %d\n",
1876  dPdlnT[*nstart],dPdlnT[*nstart2],dPdlnT[*nend],*ErrorCode );
1877  }
1878 
1879  if( *ErrorCode >= QH_NO_REBIN )
1880  {
1881  *nstart = -1;
1882  *nstart2 = -1;
1883  *nend = -2;
1884  }
1885  return;
1886 }
1887 
1888 
1889 /* rebin the quantum heating results to speed up RT_diffuse */
1890 STATIC long RebinQHeatResults(size_t nd,
1891  long nstart,
1892  long nend,
1893  vector<double>& p,
1894  vector<double>& qtemp,
1895  vector<double>& qprob,
1896  vector<double>& dPdlnT,
1897  vector<double>& u1,
1898  vector<double>& delu,
1899  vector<double>& Lambda,
1900  QH_Code *ErrorCode)
1901 {
1902  long i,
1903  newnbin;
1904  double fac,
1905  help,
1906  mul_fac,
1907  PP1,
1908  PP2,
1909  RadCooling,
1910  T1,
1911  T2,
1912  Ucheck,
1913  uu1,
1914  uu2;
1915 
1916  DEBUG_ENTRY( "RebinQHeatResults()" );
1917 
1918  /* sanity checks */
1919  ASSERT( nd < gv.bin.size() );
1920  /* >>chng 02 aug 07, changed oldnbin -> nstart..nend */
1921  ASSERT( nstart >= 0 && nstart < nend && nend < NQGRID );
1922 
1923  /* leading entries may have become very small or zero -> skip */
1924  for( i=nstart; i <= nend && dPdlnT[i] < PROB_CUTOFF_LO; i++ ) {}
1925 
1926  /* >>chng 04 oct 17, add fail-safe to keep lint happy, but this should never happen... */
1927  if( i >= NQGRID )
1928  {
1929  *ErrorCode = MAX2(*ErrorCode,QH_REBIN_FAIL);
1930  return 0;
1931  }
1932 
1933  /* >>chng 02 aug 15, use malloc to prevent stack overflows */
1934  vector<double> new_delu(NQGRID);
1935  vector<double> new_dPdlnT(NQGRID);
1936  vector<double> new_Lambda(NQGRID);
1937  vector<double> new_p(NQGRID);
1938  vector<double> new_qprob(NQGRID);
1939  vector<double> new_qtemp(NQGRID);
1940  vector<double> new_u1(NQGRID);
1941  double z[8];
1942 
1943  newnbin = 0;
1944 
1945  T1 = qtemp[i];
1946  PP1 = p[i];
1947  uu1 = u1[i];
1948 
1949  /* >>chng 04 feb 01, change 2.*NQMIN -> 1.5*NQMIN, PvH */
1950  help = pow(qtemp[nend]/qtemp[i],1./(1.5*NQMIN));
1951  mul_fac = MIN2(QT_RATIO,help);
1952 
1953  Ucheck = u1[i];
1954  RadCooling = 0.;
1955 
1956  while( i < nend )
1957  {
1958  bool lgBoundErr;
1959  bool lgDone= false;
1960  double s0 = 0.;
1961  double s1 = 0.;
1962  double wid = 0.;
1963  double xx,y;
1964 
1965  T2 = T1*mul_fac;
1966 
1967  do
1968  {
1969  double p1,p2,L1,L2,frac,slope;
1970  if( qtemp[i] <= T1 && T1 <= qtemp[i+1] )
1971  {
1972  /* >>chng 01 nov 15, copy uu2 into uu1 instead, PvH */
1973  /* uu1 = ufunct(T1,nd); */
1974  double xrlog = log(qtemp[i+1]/qtemp[i]);
1975  if( xrlog > 0. )
1976  {
1977  frac = log(T1/qtemp[i]);
1978  slope = log(p[i+1]/p[i])/xrlog;
1979  p1 = p[i]*exp(frac*slope);
1980  slope = log(Lambda[i+1]/Lambda[i])/xrlog;
1981  L1 = Lambda[i]*exp(frac*slope);
1982  }
1983  else
1984  {
1985  /* pathological case where slope is extremely steep (daniela.in) */
1986  p1 = sqrt(p[i]*p[i+1]);
1987  L1 = sqrt(Lambda[i]*Lambda[i+1]);
1988  }
1989  }
1990  else
1991  {
1992  /* >>chng 01 nov 15, copy uu2 into uu1 instead, PvH */
1993  /* uu1 = u1[i]; */
1994  p1 = p[i];
1995  L1 = Lambda[i];
1996  }
1997  if( qtemp[i] <= T2 && T2 <= qtemp[i+1] )
1998  {
1999  /* >>chng 02 apr 30, make sure this doesn't point beyond valid range, PvH */
2000  help = ufunct(T2,nd,&lgBoundErr);
2001  uu2 = MIN2(help,u1[i+1]);
2002  ASSERT( !lgBoundErr ); /* this should never be triggered */
2003  double xrlog = log(qtemp[i+1]/qtemp[i]);
2004  if( xrlog > 0. )
2005  {
2006  frac = log(T2/qtemp[i]);
2007  slope = log(p[i+1]/p[i])/xrlog;
2008  p2 = p[i]*exp(frac*slope);
2009  slope = log(Lambda[i+1]/Lambda[i])/xrlog;
2010  L2 = Lambda[i]*exp(frac*slope);
2011  }
2012  else
2013  {
2014  /* pathological case where slope is extremely steep */
2015  p2 = sqrt(p[i]*p[i+1]);
2016  L2 = sqrt(Lambda[i]*Lambda[i+1]);
2017  }
2018  lgDone = true;
2019  }
2020  else
2021  {
2022  uu2 = u1[i+1];
2023  p2 = p[i+1];
2024  L2 = Lambda[i+1];
2025  /* >>chng 01 nov 15, this caps the range in p(U) integrated in one bin
2026  * it helps avoid spurious QH_BOUND_FAIL's when flank is very steep, PvH */
2027  if( MAX2(p2,PP1)/MIN2(p2,PP1) > sqrt(SAFETY) )
2028  {
2029  lgDone = true;
2030  T2 = qtemp[i+1];
2031  }
2032  ++i;
2033  }
2034  PP2 = p2;
2035  wid += uu2 - uu1;
2036  /* sanity check */
2037  ASSERT( wid >= 0. );
2038  vlog(z,uu1,uu2,p1,p2,p1*L1,p2*L2,1.,1.);
2039  s0 += log_integral(uu1,p1,uu2,p2,z[0],z[2],z[1],z[3]);
2040  s1 += log_integral(uu1,p1*L1,uu2,p2*L2,z[0],z[4],z[1],z[5]);
2041  uu1 = uu2;
2042 
2043  } while( i < nend && ! lgDone );
2044 
2045  /* >>chng 01 nov 14, if T2 == qtemp[oldnbin-1], the code will try another iteration
2046  * break here to avoid zero divide, the assert on Ucheck tests if we are really finished */
2047  /* >>chng 01 dec 04, change ( s0 == 0. ) to ( s0 <= 0. ), PvH */
2048  if( s0 <= 0. )
2049  {
2050  ASSERT( wid == 0. );
2051  break;
2052  }
2053 
2054  new_qprob[newnbin] = s0;
2055  new_Lambda[newnbin] = s1/s0;
2056  xx = log(new_Lambda[newnbin]*EN1RYD*gv.bin[nd]->cnv_GR_pH);
2057  splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,xx,&y,&lgBoundErr);
2058  ASSERT( !lgBoundErr ); /* this should never be triggered */
2059  new_qtemp[newnbin] = exp(y);
2060  new_u1[newnbin] = ufunct(new_qtemp[newnbin],nd,&lgBoundErr);
2061  ASSERT( !lgBoundErr ); /* this should never be triggered */
2062  new_delu[newnbin] = wid;
2063  new_p[newnbin] = new_qprob[newnbin]/new_delu[newnbin];
2064  new_dPdlnT[newnbin] = new_p[newnbin]*new_qtemp[newnbin]*uderiv(new_qtemp[newnbin],nd);
2065 
2066  Ucheck += wid;
2067  RadCooling += new_qprob[newnbin]*new_Lambda[newnbin];
2068 
2069  T1 = T2;
2070  PP1 = PP2;
2071  ++newnbin;
2072  }
2073 
2074  /* >>chng 01 jul 13, add fail-safe */
2075  if( newnbin < NQMIN )
2076  {
2077  *ErrorCode = MAX2(*ErrorCode,QH_REBIN_FAIL);
2078  return 0;
2079  }
2080 
2081  fac = RadCooling*EN1RYD*gv.bin[nd]->cnv_GR_pCM3/gv.bin[nd]->GrainHeat;
2082 
2083  if( trace.lgTrace && trace.lgDustBug )
2084  {
2085  fprintf( ioQQQ, " RebinQHeatResults found tol1 %.4e tol2 %.4e\n",
2086  fabs(u1[nend]/Ucheck-1.), fabs(fac-1.) );
2087  }
2088 
2089  /* do quality control */
2090  /* >>chng 02 apr 30, tighten up check, PvH */
2091  ASSERT( fabs(u1[nend]/Ucheck-1.) < 10.*sqrt((double)(nend-nstart+newnbin))*DBL_EPSILON );
2092 
2093  if( fabs(fac-1.) > CONSERV_TOL )
2094  *ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);
2095 
2096  for( i=0; i < newnbin; i++ )
2097  {
2098  /* renormalize the distribution to assure energy conservation */
2099  p[i] = new_p[i]/fac;
2100  qtemp[i] = new_qtemp[i];
2101  qprob[i] = new_qprob[i]/fac;
2102  dPdlnT[i] = new_dPdlnT[i]/fac;
2103  u1[i] = new_u1[i];
2104  delu[i] = new_delu[i];
2105  Lambda[i] = new_Lambda[i];
2106 
2107  /* sanity checks */
2108  ASSERT( qtemp[i] > 0. && qprob[i] > 0. );
2109 
2110 /* printf(" rk %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",i,qtemp[i],u1[i],p[i],dPdlnT[i]); */
2111  }
2112  return newnbin;
2113 }
2114 
2115 
2116 /* calculate approximate probability distribution in high intensity limit */
2117 STATIC void GetProbDistr_HighLimit(long nd,
2118  double TolFac,
2119  double Umax,
2120  double fwhm,
2121  /*@out@*/vector<double>& qtemp,
2122  /*@out@*/vector<double>& qprob,
2123  /*@out@*/vector<double>& dPdlnT,
2124  /*@out@*/double *tol,
2125  /*@out@*/long *qnbin,
2126  /*@out@*/double *new_tmin,
2127  /*@out@*/QH_Code *ErrorCode)
2128 {
2129  bool lgBoundErr,
2130  lgErr;
2131  long i,
2132  nbin;
2133  double c1,
2134  c2,
2135  delu[NQGRID],
2136  fac,
2137  fac1,
2138  fac2,
2139  help1,
2140  help2,
2141  L1,
2142  L2,
2143  Lambda[NQGRID],
2144  mul_fac,
2145  p[NQGRID],
2146  p1,
2147  p2,
2148  RadCooling,
2149  sum,
2150  T1,
2151  T2,
2152  Tlo,
2153  Thi,
2154  Ulo,
2155  Uhi,
2156  uu1,
2157  uu2,
2158  xx,
2159  y;
2160 
2161  DEBUG_ENTRY( "GetProbDistr_HighLimit()" );
2162 
2163  if( trace.lgTrace && trace.lgDustBug )
2164  {
2165  fprintf( ioQQQ, " GetProbDistr_HighLimit called for grain %s\n", gv.bin[nd]->chDstLab );
2166  }
2167 
2168  c1 = sqrt(4.*LN_TWO/PI)/fwhm*exp(-pow2(fwhm/Umax)/(16.*LN_TWO));
2169  c2 = 4.*LN_TWO*pow2(Umax/fwhm);
2170 
2171  fac1 = fwhm/Umax*sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO));
2172  /* >>chng 03 nov 10, safeguard against underflow, PvH */
2173  help1 = Umax*exp(-fac1);
2174  help2 = exp(gv.bin[nd]->DustEnth[0]);
2175  Ulo = MAX2(help1,help2);
2176  /* >>chng 03 jan 28, ignore lgBoundErr on lower boundary, low-T grains have negigible emission, PvH */
2177  Tlo = inv_ufunct(Ulo,nd,&lgBoundErr);
2178 
2179  fac2 = fwhm/Umax*sqrt(-log(PROB_CUTOFF_HI)/(4.*LN_TWO));
2180  /* >>chng 03 nov 10, safeguard against overflow, PvH */
2181  if( fac2 > log(DBL_MAX/10.) )
2182  {
2183  *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
2184  return;
2185  }
2186  Uhi = Umax*exp(fac2);
2187  Thi = inv_ufunct(Uhi,nd,&lgBoundErr);
2188 
2189  nbin = 0;
2190 
2191  T1 = Tlo;
2192  uu1 = ufunct(T1,nd,&lgErr);
2193  lgBoundErr = lgBoundErr || lgErr;
2194  help1 = log(uu1/Umax);
2195  p1 = c1*exp(-c2*pow2(help1));
2196  splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(T1),&y,&lgErr);
2197  lgBoundErr = lgBoundErr || lgErr;
2198  L1 = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
2199 
2200  /* >>chng 03 nov 10, safeguard against underflow, PvH */
2201  if( uu1*p1*L1 < 1.e5*DBL_MIN )
2202  {
2203  *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
2204  return;
2205  }
2206 
2207  /* >>chng 04 feb 01, change 2.*NQMIN -> 1.2*NQMIN, PvH */
2208  help1 = pow(Thi/Tlo,1./(1.2*NQMIN));
2209  mul_fac = MIN2(QT_RATIO,help1);
2210 
2211  sum = 0.;
2212  RadCooling = 0.;
2213 
2214  double z[8];
2215 
2216  do
2217  {
2218  double s0,s1,wid;
2219 
2220  T2 = T1*mul_fac;
2221  uu2 = ufunct(T2,nd,&lgErr);
2222  lgBoundErr = lgBoundErr || lgErr;
2223  help1 = log(uu2/Umax);
2224  p2 = c1*exp(-c2*pow2(help1));
2225  splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(T2),&y,&lgErr);
2226  lgBoundErr = lgBoundErr || lgErr;
2227  L2 = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
2228 
2229  wid = uu2 - uu1;
2230  vlog(z,uu1,uu2,p1,p2,p1*L1,p2*L2,1.,1.);
2231  s0 = log_integral(uu1,p1,uu2,p2,z[0],z[2],z[1],z[3]);
2232  s1 = log_integral(uu1,p1*L1,uu2,p2*L2,z[0],z[4],z[1],z[5]);
2233 
2234  qprob[nbin] = s0;
2235  Lambda[nbin] = s1/s0;
2236  xx = log(Lambda[nbin]*EN1RYD*gv.bin[nd]->cnv_GR_pH);
2237  splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,xx,&y,&lgErr);
2238  lgBoundErr = lgBoundErr || lgErr;
2239  qtemp[nbin] = exp(y);
2240  delu[nbin] = wid;
2241  p[nbin] = qprob[nbin]/delu[nbin];
2242  dPdlnT[nbin] = p[nbin]*qtemp[nbin]*uderiv(qtemp[nbin],nd);
2243 
2244  sum += qprob[nbin];
2245  RadCooling += qprob[nbin]*Lambda[nbin];
2246 
2247  T1 = T2;
2248  uu1 = uu2;
2249  p1 = p2;
2250  L1 = L2;
2251 
2252  ++nbin;
2253 
2254  } while( T2 < Thi && nbin < NQGRID );
2255 
2256  fac = RadCooling*EN1RYD*gv.bin[nd]->cnv_GR_pCM3/gv.bin[nd]->GrainHeat;
2257 
2258  for( i=0; i < nbin; ++i )
2259  {
2260  qprob[i] /= fac;
2261  dPdlnT[i] /= fac;
2262  }
2263 
2264  *tol = sum/fac;
2265  *qnbin = nbin;
2266  *new_tmin = qtemp[0];
2267  *ErrorCode = MAX2(*ErrorCode,QH_ANALYTIC);
2268 
2269  /* do quality control */
2270  if( TolFac > STRICT )
2271  *ErrorCode = MAX2(*ErrorCode,QH_ANALYTIC_RELAX);
2272 
2273  if( lgBoundErr )
2274  *ErrorCode = MAX2(*ErrorCode,QH_THIGH_FAIL);
2275 
2276  if( fabs(sum/fac-1.) > PROB_TOL )
2277  *ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);
2278 
2279  if( dPdlnT[0] > SAFETY*PROB_CUTOFF_LO || dPdlnT[nbin-1] > SAFETY*PROB_CUTOFF_HI )
2280  *ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
2281 
2282  if( trace.lgTrace && trace.lgDustBug )
2283  {
2284  fprintf( ioQQQ, " GetProbDistr_HighLimit found tol1 %.4e tol2 %.4e\n",
2285  fabs(sum-1.), fabs(sum/fac-1.) );
2286  fprintf( ioQQQ, " zone %ld %s nbin %ld total prob %.4e new_tmin %.4e\n",
2287  nzone,gv.bin[nd]->chDstLab,nbin,sum/fac,*new_tmin );
2288  }
2289  return;
2290 }
2291 
2292 
2293 /* calculate derivative of the enthalpy function dU/dT (aka the specific heat) at a given temperature, in Ryd/K */
2294 STATIC double uderiv(double temp,
2295  size_t nd)
2296 {
2297  enth_type ecase;
2298  long int i,
2299  j;
2300  double N_C,
2301  N_H;
2302  double deriv = 0.,
2303  hok[3] = {1275., 1670., 4359.},
2304  numer,
2305  dnumer,
2306  denom,
2307  ddenom,
2308  x;
2309 
2310 
2311  DEBUG_ENTRY( "uderiv()" );
2312 
2313  if( temp <= 0. )
2314  {
2315  fprintf( ioQQQ, " uderiv called with non-positive temperature: %.6e\n" , temp );
2316  cdEXIT(EXIT_FAILURE);
2317  }
2318  ASSERT( nd < gv.bin.size() );
2319 
2320  const double CALORIE = 4.184e7; /* this is the thermochemical calorie in erg */
2321 
2322  ecase = gv.which_enth[gv.bin[nd]->matType];
2323  switch( ecase )
2324  {
2325  case ENTH_CAR:
2326  numer = (4.15e-22/EN1RYD)*pow(temp,3.3);
2327  dnumer = (3.3*4.15e-22/EN1RYD)*pow(temp,2.3);
2328  denom = 1. + 6.51e-03*temp + 1.5e-06*temp*temp + 8.3e-07*pow(temp,2.3);
2329  ddenom = 6.51e-03 + 2.*1.5e-06*temp + 2.3*8.3e-07*pow(temp,1.3);
2330  /* dU/dT for pah/graphitic grains in Ryd/K, derived from:
2331  * >>refer grain physics Guhathakurta & Draine, 1989, ApJ, 345, 230 */
2332  deriv = (dnumer*denom - numer*ddenom)/pow2(denom);
2333  break;
2334  case ENTH_CAR2:
2335  /* dU/dT for graphite grains in Ryd/K, using eq 9 of */
2336  /* >>refer grain physics Draine B.T., and Li A., 2001, ApJ, 551, 807 */
2337  deriv = (DebyeDeriv(temp/863.,2) + 2.*DebyeDeriv(temp/2504.,2))*BOLTZMANN/EN1RYD;
2338  break;
2339  case ENTH_SIL:
2340  /* dU/dT for silicate grains (and grey grains) in Ryd/K */
2341  /* limits of tlim set above, 0 and DBL_MAX, so always OK */
2342  /* >>refer grain physics Guhathakurta & Draine, 1989, ApJ, 345, 230 */
2343  for( j = 0; j < 4; j++ )
2344  {
2345  if( temp > tlim[j] && temp <= tlim[j+1] )
2346  {
2347  deriv = cval[j]*pow(temp,ppower[j]);
2348  break;
2349  }
2350  }
2351  break;
2352  case ENTH_SIL2:
2353  /* dU/dT for silicate grains in Ryd/K, using eq 11 of */
2354  /* >>refer grain physics Draine B.T., and Li A., 2001, ApJ, 551, 807 */
2355  deriv = (2.*DebyeDeriv(temp/500.,2) + DebyeDeriv(temp/1500.,3))*BOLTZMANN/EN1RYD;
2356  break;
2357  case ENTH_PAH:
2358  /* dU/dT for PAH grains in Ryd/K, using eq A.4 of */
2359  /* >>refer grain physics Dwek E., Arendt R.G., Fixsen D.J. et al., 1997, ApJ, 475, 565 */
2360  /* this expression is only valid upto 2000K */
2361  x = log10(MIN2(temp,2000.));
2362  deriv = exp10(-21.26+3.1688*x-0.401894*pow2(x))/EN1RYD;
2363  break;
2364  case ENTH_PAH2:
2365  /* dU/dT for PAH grains in Ryd/K, approximately using eq 33 of */
2366  /* >>refer grain physics Draine B.T., and Li A., 2001, ApJ, 551, 807 */
2367  /* N_C and N_H should actually be nint(N_C) and nint(N_H),
2368  * but this can lead to FP overflow for very large grains */
2369  N_C = no_atoms(nd);
2370  if( N_C <= 25. )
2371  {
2372  N_H = 0.5*N_C;
2373  }
2374  else if( N_C <= 100. )
2375  {
2376  N_H = 2.5*sqrt(N_C);
2377  }
2378  else
2379  {
2380  N_H = 0.25*N_C;
2381  }
2382  deriv = 0.;
2383  for( i=0; i < 3; i++ )
2384  {
2385  double help1 = hok[i]/temp;
2386  if( help1 < 300. )
2387  {
2388  double help2 = exp(help1);
2389  double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
2390  deriv += N_H/(N_C-2.)*pow2(help1)*help2/pow2(help3)*BOLTZMANN/EN1RYD;
2391  }
2392  }
2393  deriv += (DebyeDeriv(temp/863.,2) + 2.*DebyeDeriv(temp/2504.,2))*BOLTZMANN/EN1RYD;
2394  break;
2395  case ENTH_SIC:
2396  /* c_p for alpha-SiC in cal/(mol*K), derived from Table 4 in */
2397  /* >>refer grain physics Chekhovskoy V. Ya., 1971, J. Chem. Thermodynamics, 3, 289 */
2398  if( temp < 300. )
2399  deriv = 13.25*(3.979984e-01*DebyeDeriv(temp/747.,3) + 6.020016e-01*DebyeDeriv(temp/1647.,3));
2400  else
2401  deriv = 13.250 - 2035./temp + 288.e5/pow2(temp)*exp(-5680./temp); /* Eq. 3 */
2402  // now convert to Ryd/K per atom
2403  deriv *= CALORIE/(EN1RYD*2.*AVOGADRO);
2404  break;
2405  default:
2406  fprintf( ioQQQ, " uderiv called with unknown type for enthalpy function: %d\n", ecase );
2407  cdEXIT(EXIT_FAILURE);
2408  }
2409 
2410  /* >>chng 00 mar 23, use formula 3.1 of Guhathakurtha & Draine, by PvH */
2411  /* >>chng 03 jan 17, use MAX2(..,1) to prevent crash for extremely small grains, PvH */
2412  deriv *= MAX2(no_atoms(nd)-2.,1.);
2413 
2414  if( deriv <= 0. )
2415  {
2416  fprintf( ioQQQ, " uderiv finds non-positive derivative: %.6e, what's up?\n" , deriv );
2417  cdEXIT(EXIT_FAILURE);
2418  }
2419  return deriv;
2420 }
2421 
2422 
2423 /* calculate the enthalpy of a grain at a given temperature, in Ryd */
2424 STATIC double ufunct(double temp,
2425  size_t nd,
2426  /*@out@*/ bool *lgBoundErr)
2427 {
2428  double enthalpy,
2429  y;
2430 
2431  DEBUG_ENTRY( "ufunct()" );
2432 
2433  if( temp <= 0. )
2434  {
2435  fprintf( ioQQQ, " ufunct called with non-positive temperature: %.6e\n" , temp );
2436  cdEXIT(EXIT_FAILURE);
2437  }
2438  ASSERT( nd < gv.bin.size() );
2439 
2440  /* >>chng 02 apr 22, interpolate in DustEnth array to get enthalpy, by PvH */
2441  splint_safe(gv.dsttmp,gv.bin[nd]->DustEnth,gv.bin[nd]->EnthSlp,NDEMS,log(temp),&y,lgBoundErr);
2442  enthalpy = exp(y);
2443 
2444  ASSERT( enthalpy > 0. );
2445  return enthalpy;
2446 }
2447 
2448 
2449 /* this is the inverse of ufunct: determine grain temperature as a function of enthalpy */
2450 STATIC double inv_ufunct(double enthalpy,
2451  size_t nd,
2452  /*@out@*/ bool *lgBoundErr)
2453 {
2454  double temp,
2455  y;
2456 
2457  DEBUG_ENTRY( "inv_ufunct()" );
2458 
2459  if( enthalpy <= 0. )
2460  {
2461  fprintf( ioQQQ, " inv_ufunct called with non-positive enthalpy: %.6e\n" , enthalpy );
2462  cdEXIT(EXIT_FAILURE);
2463  }
2464  ASSERT( nd < gv.bin.size() );
2465 
2466  /* >>chng 02 apr 22, interpolate in DustEnth array to get temperature, by PvH */
2467  splint_safe(gv.bin[nd]->DustEnth,gv.dsttmp,gv.bin[nd]->EnthSlp2,NDEMS,log(enthalpy),&y,lgBoundErr);
2468  temp = exp(y);
2469 
2470  ASSERT( temp > 0. );
2471  return temp;
2472 }
2473 
2474 
2475 /* initialize interpolation arrys for grain enthalpy */
2476 void InitEnthalpy()
2477 {
2478  DEBUG_ENTRY( "InitEnthalpy()" );
2479 
2480  double z[8];
2481  for( size_t nd=0; nd < gv.bin.size(); nd++ )
2482  {
2483  double tdust2 = GRAIN_TMIN;
2484  double C_V2 = uderiv(tdust2,nd);
2485  /* at low temps, C_V = C*T^3 -> U = C*T^4/4 = C_V*T/4 */
2486  gv.bin[nd]->DustEnth[0] = C_V2*tdust2/4.;
2487  double tdust1 = tdust2;
2488  double C_V1 = C_V2;
2489 
2490  for( long i=1; i < NDEMS; i++ )
2491  {
2492  double tmid,Cmid;
2493  tdust2 = exp(gv.dsttmp[i]);
2494  C_V2 = uderiv(tdust2,nd);
2495  tmid = sqrt(tdust1*tdust2);
2496  /* this ensures accuracy for silicate enthalpy */
2497  for( long j=1; j < 4; j++ )
2498  {
2499  if( tdust1 < tlim[j] && tlim[j] < tdust2 )
2500  {
2501  tmid = tlim[j];
2502  break;
2503  }
2504  }
2505  Cmid = uderiv(tmid,nd);
2506  vlog(z,tdust1,tmid,tdust2,C_V1,Cmid,C_V2,1.,1.);
2507  gv.bin[nd]->DustEnth[i] = gv.bin[nd]->DustEnth[i-1] +
2508  log_integral(tdust1,C_V1,tmid,Cmid,z[0],z[3],z[1],z[4]) +
2509  log_integral(tmid,Cmid,tdust2,C_V2,z[1],z[4],z[2],z[5]);
2510  tdust1 = tdust2;
2511  C_V1 = C_V2;
2512  }
2513  }
2514 
2515  /* conversion for logarithmic interpolation */
2516  for( size_t nd=0; nd < gv.bin.size(); nd++ )
2517  {
2518  vlog(gv.bin[nd]->DustEnth,gv.bin[nd]->DustEnth,0,NDEMS);
2519  /* set up coefficients for splint */
2520  spline(gv.dsttmp,gv.bin[nd]->DustEnth,NDEMS,2e31,2e31,gv.bin[nd]->EnthSlp);
2521  spline(gv.bin[nd]->DustEnth,gv.dsttmp,NDEMS,2e31,2e31,gv.bin[nd]->EnthSlp2);
2522  }
2523  return;
2524 }
2525 
2526 
2527 /* helper function for calculating specific heat, uses Debye theory */
2528 STATIC double DebyeDeriv(double x,
2529  long n)
2530 {
2531  long i,
2532  nn;
2533  double res;
2534 
2535  DEBUG_ENTRY( "DebyeDeriv()" );
2536 
2537  ASSERT( x > 0. );
2538  ASSERT( n == 2 || n == 3 );
2539 
2540  if( x < 0.001 )
2541  {
2542  /* for general n this is Gamma(n+2)*zeta(n+1)*powi(x,n) */
2543  if( n == 2 )
2544  {
2545  res = 6.*1.202056903159594*pow2(x);
2546  }
2547  else if( n == 3 )
2548  {
2549  res = 24.*1.082323233711138*pow3(x);
2550  }
2551  else
2552  /* added to keep lint happy - note that assert above confimred that n is 2 or 3,
2553  * but lint flagged possible flow without defn of res */
2554  TotalInsanity();
2555  }
2556  else
2557  {
2558  nn = 4*MAX2(4,2*(long)(0.05/x));
2559  vector<double> xx(nn);
2560  vector<double> rr(nn);
2561  vector<double> aa(nn);
2562  vector<double> ww(nn);
2563  gauss_legendre(nn,xx,aa);
2564  gauss_init(nn,0.,1.,xx,aa,rr,ww);
2565 
2566  res = 0.;
2567  for( i=0; i < nn; i++ )
2568  {
2569  double help1 = rr[i]/x;
2570  if( help1 < 300. )
2571  {
2572  double help2 = exp(help1);
2573  double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
2574  res += ww[i]*powi(rr[i],n+1)*help2/pow2(help3);
2575  }
2576  }
2577  res /= pow2(x);
2578  }
2579  return (double)n*res;
2580 }
t_mesh::emm
double emm() const
Definition: mesh.h:84
GrainVar::lgBakesPAH_heat
bool lgBakesPAH_heat
Definition: grainvar.h:481
LOOP_MAX
static const long LOOP_MAX
Definition: grains_qheat.cpp:65
ChargeBin::PotSurf
double PotSurf
Definition: grainvar.h:228
GetProbDistr_HighLimit
STATIC void GetProbDistr_HighLimit(long, double, double, double, vector< double > &, vector< double > &, vector< double > &, double *, long *, double *, QH_Code *)
Definition: grains_qheat.cpp:2117
QH_DELTA
Definition: grains_qheat.cpp:29
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:751
NSTARTUP
static const long NSTARTUP
Definition: grains_qheat.cpp:57
QH_DELTA_FAIL
Definition: grains_qheat.cpp:35
avx_ptr
Definition: vectorize.h:151
thermal
t_thermal thermal
Definition: thermal.cpp:6
gauss_legendre
void gauss_legendre(long int, vector< double > &, vector< double > &)
Definition: grains_mie.cpp:4604
InitEnthalpy
void InitEnthalpy()
Definition: grains_qheat.cpp:2476
t_trace::lgDustBug
bool lgDustBug
Definition: trace.h:76
ChargeBin::ehat
flex_arr< realnum > ehat
Definition: grainvar.h:242
avx_ptr::ptr0
T * ptr0()
Definition: vectorize.h:221
exp10
double exp10(double x)
Definition: cddefines.h:1383
GrainVar::QHSaveFile
FILE * QHSaveFile
Definition: grainvar.h:573
GrainVar::GrainHeatSum
double GrainHeatSum
Definition: grainvar.h:546
get_ptr
T * get_ptr(T *v)
Definition: cddefines.h:1113
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
RELAX
static const double RELAX
Definition: grains_qheat.cpp:95
t_mesh::widflx
double widflx(size_t i) const
Definition: mesh.h:147
QH_NBIN_FAIL
Definition: grains_qheat.cpp:42
no_atoms
double no_atoms(size_t nd)
Definition: grains_qheat.cpp:16
QH_ARRAY_FAIL
Definition: grains_qheat.cpp:32
spldrv_safe
void spldrv_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
Definition: thirdparty.h:239
uderiv
STATIC double uderiv(double, size_t)
Definition: grains_qheat.cpp:2294
GrainVar::GasCoolColl
double GasCoolColl
Definition: grainvar.h:546
DEF_FAC
static const double DEF_FAC
Definition: grains_qheat.cpp:68
strg_type
strg_type
Definition: grainvar.h:82
LIM1
static const realnum LIM1
Definition: grains_qheat.cpp:175
TorF
char TorF(bool l)
Definition: cddefines.h:753
GrainMakeDiffuse
void GrainMakeDiffuse()
Definition: grains_qheat.cpp:222
QH_ANALYTIC
Definition: grains_qheat.cpp:29
QH_BOUND_FAIL
Definition: grains_qheat.cpp:35
MAX_LOOP
static const long MAX_LOOP
Definition: grains_qheat.cpp:85
GrainVar::dsttmp
double dsttmp[NDEMS]
Definition: grainvar.h:566
sign
T sign(T x, T y)
Definition: cddefines.h:846
ChargeBin::HeatingRate2
double HeatingRate2
Definition: grainvar.h:279
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_rfield::SummedCon
double * SummedCon
Definition: rfield.h:163
ChargeBin::yhat_primary
flex_arr< realnum > yhat_primary
Definition: grainvar.h:241
LIM3
static const realnum LIM3
Definition: grains_qheat.cpp:177
tlim
static const double tlim[5]
Definition: grains_qheat.cpp:127
pow3
T pow3(T a)
Definition: cddefines.h:992
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
GrainVar::lgQHPunLast
bool lgQHPunLast
Definition: grainvar.h:572
QH_THIGH_FAIL
Definition: grains_qheat.cpp:32
nzone
long int nzone
Definition: cddefines.cpp:14
ufunct
STATIC double ufunct(double, size_t, bool *)
Definition: grains_qheat.cpp:2424
MIN2
#define MIN2(a, b)
Definition: cddefines.h:807
vlog
void vlog(const double x[], double y[], long nlo, long nhi)
Definition: vectorize_log.cpp:219
t_mesh::anu
double anu(size_t i) const
Definition: mesh.h:111
MAX_EVENTS
static const double MAX_EVENTS
Definition: grains_qheat.cpp:61
QH_STEP_FAIL
Definition: grains_qheat.cpp:39
ChargeBin::FracPop
double FracPop
Definition: grainvar.h:228
dense
t_dense dense
Definition: global.cpp:15
GrainVar::GraphiteEmission
vector< realnum > GraphiteEmission
Definition: grainvar.h:581
log_integral
STATIC double log_integral(double, double, double, double, double, double, double, double)
Definition: grains_qheat.cpp:1759
STRICT
static const double STRICT
Definition: grains_qheat.cpp:94
ChargeBin::yhat
flex_arr< realnum > yhat
Definition: grainvar.h:240
QH_REBIN_FAIL
Definition: grains_qheat.cpp:42
ppower
static const double ppower[4]
Definition: grains_qheat.cpp:128
DebyeDeriv
STATIC double DebyeDeriv(double, long)
Definition: grains_qheat.cpp:2528
QH_LOW_FAIL
Definition: grains_qheat.cpp:39
gauss_init
void gauss_init(long int, double, double, const vector< double > &, const vector< double > &, vector< double > &, vector< double > &)
Definition: grains_mie.cpp:4576
ChargeBin::DustZ
long DustZ
Definition: grainvar.h:224
t_rfield::nflux_with_check
long int nflux_with_check
Definition: rfield.h:51
FWHM_RATIO2
static const double FWHM_RATIO2
Definition: grains_qheat.cpp:82
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
trace
t_trace trace
Definition: trace.cpp:5
QHEAT_TOL
static const double QHEAT_TOL
Definition: grains_qheat.cpp:88
FWHM_RATIO
static const double FWHM_RATIO
Definition: grains_qheat.cpp:79
dprintf
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1198
t_mesh::ipointC
size_t ipointC(double anu) const
Definition: mesh.h:152
GrainVar::dstAbundThresholdNear
realnum dstAbundThresholdNear
Definition: grainvar.h:569
spline
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
Definition: thirdparty.h:150
QH_WIDE_FAIL
Definition: grains_qheat.cpp:42
DEN_SIL
static const double DEN_SIL
Definition: grains_qheat.cpp:121
QH_Code
QH_Code
Definition: grains_qheat.cpp:26
qheat
void qheat(vector< double > &, vector< double > &, long *, size_t)
GrainVar::GrainEmission
vector< realnum > GrainEmission
Definition: grainvar.h:580
GrainVar::which_strg
strg_type which_strg[MAT_TOP]
Definition: grainvar.h:522
ChargeBin::ipThresInf
long ipThresInf
Definition: grainvar.h:224
fast_expm1
realnum fast_expm1(realnum x)
Definition: grains_qheat.cpp:180
QH_RETRY
Definition: grains_qheat.cpp:35
STATIC
#define STATIC
Definition: cddefines.h:118
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
GrainVar::dstAbundThresholdFar
realnum dstAbundThresholdFar
Definition: grainvar.h:570
ChargeBin::ThresInfInc
double ThresInfInc
Definition: grainvar.h:228
rfield
t_rfield rfield
Definition: rfield.cpp:9
GrainVar::which_enth
enth_type which_enth[MAT_TOP]
Definition: grainvar.h:517
GetProbDistr_LowLimit
STATIC void GetProbDistr_LowLimit(size_t, double, double, double, const vector< double > &, const vector< double > &, vector< double > &, vector< double > &, vector< double > &, long *, double *, long *, QH_Code *)
realnum
float realnum
Definition: cddefines.h:124
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
avx_ptr::data
T * data()
Definition: vectorize.h:213
QH_CONV_FAIL
Definition: grains_qheat.cpp:35
max
long max(int a, long b)
Definition: cddefines.h:821
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:484
powi
double powi(double, long int)
Definition: service.cpp:690
min
long min(int a, long b)
Definition: cddefines.h:766
ScanProbDistr
STATIC void ScanProbDistr(const vector< double > &, const vector< double > &, long, double, long, double, long *, long *, long *, long *, QH_Code *)
Definition: grains_qheat.cpp:1784
enth_type
enth_type
Definition: grainvar.h:47
PROB_TOL
static const double PROB_TOL
Definition: grains_qheat.cpp:71
t_mesh::anu2
double anu2(size_t i) const
Definition: mesh.h:115
iterations
t_iterations iterations
Definition: iterations.cpp:6
ChargeBin::PotSurfInc
double PotSurfInc
Definition: grainvar.h:228
t_thermal::heating
double heating(long nelem, long ion)
Definition: thermal.h:186
t_mesh::anumin
double anumin(size_t i) const
Definition: mesh.h:139
WIDE_FAIL_MAX
static const long WIDE_FAIL_MAX
Definition: grains_qheat.cpp:91
QH_FATAL
Definition: grains_qheat.cpp:42
NQMIN
static const long NQMIN
Definition: grains_qheat.cpp:49
QH_OK
Definition: grains_qheat.cpp:29
SAFETY
static const double SAFETY
Definition: grains_qheat.cpp:54
GrainVar::lgWD01
bool lgWD01
Definition: grainvar.h:481
GRAIN_TMIN
const double GRAIN_TMIN
Definition: grainvar.h:19
inv_ufunct
STATIC double inv_ufunct(double, size_t, bool *)
Definition: grains_qheat.cpp:2450
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:617
t_iterations::lgLastIt
bool lgLastIt
Definition: iterations.h:47
ChargeBin::nfill
long nfill
Definition: grainvar.h:224
GrainVar::lgDHetOn
bool lgDHetOn
Definition: grainvar.h:491
ChargeBin::ipThresInfVal
long ipThresInfVal
Definition: grainvar.h:224
pow2
T pow2(T a)
Definition: cddefines.h:985
LIM2
static const realnum LIM2
Definition: grains_qheat.cpp:176
QH_LOOP_FAIL
Definition: grains_qheat.cpp:32
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
powpq
double powpq(double x, int p, int q)
Definition: service.cpp:726
QH_HIGH_FAIL
Definition: grains_qheat.cpp:39
NDEMS
const int NDEMS
Definition: grainvar.h:16
NQTEST
static const long NQTEST
Definition: grains_qheat.cpp:75
QH_NEGRATE_FAIL
Definition: grains_qheat.cpp:32
GrainVar::GrainHeatScaleFactor
realnum GrainHeatScaleFactor
Definition: grainvar.h:559
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:585
t_hmi::H2_total
double H2_total
Definition: hmi.h:25
t_thermal::ConstGrainTemp
realnum ConstGrainTemp
Definition: thermal.h:59
PROB_CUTOFF_LO
static const double PROB_CUTOFF_LO
Definition: grains_qheat.cpp:52
ChargeBin::fac1
flex_arr< double > fac1
Definition: grainvar.h:262
QH_ANALYTIC_RELAX
Definition: grains_qheat.cpp:29
splint_safe
void splint_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
Definition: thirdparty.h:199
qheat_init
STATIC void qheat_init(size_t, vector< double > &, double *)
Definition: grains_qheat.cpp:856
MAX2
#define MAX2(a, b)
Definition: cddefines.h:828
NQGRID
const int NQGRID
Definition: grainvar.h:34
CONSERV_TOL
const double CONSERV_TOL
Definition: grainvar.h:37
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
QH_NO_REBIN
Definition: grains_qheat.cpp:39
pow
double pow(double x, int i)
Definition: cddefines.h:782
MW_SIL
static const double MW_SIL
Definition: grains_qheat.cpp:124
QT_RATIO
static const double QT_RATIO
Definition: grains_qheat.cpp:106
MAX3
#define MAX3(a, b, c)
Definition: cddefines.h:833
GrainVar::SilicateEmission
vector< realnum > SilicateEmission
Definition: grainvar.h:582
t_mesh::anumax
double anumax(size_t i) const
Definition: mesh.h:143
PROB_CUTOFF_HI
static const double PROB_CUTOFF_HI
Definition: grains_qheat.cpp:53
gv
GrainVar gv
Definition: grainvar.cpp:5
hmi
t_hmi hmi
Definition: hmi.cpp:5
t_phycon::te
double te
Definition: phycon.h:21
frac
double frac(double d)
Definition: thirdparty.cpp:3597
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:348
vexpm1
void vexpm1(const double x[], double y[], long nlo, long nhi)
Definition: vectorize_exp.cpp:297
ChargeBin::ThresSurf
double ThresSurf
Definition: grainvar.h:228
t_rfield::nflux
long int nflux
Definition: rfield.h:48
GrainVar::GrainHeatChem
double GrainHeatChem
Definition: grainvar.h:546
t_rfield::nPositive
long int nPositive
Definition: rfield.h:55
RebinQHeatResults
STATIC long RebinQHeatResults(size_t, long, long, vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, QH_Code *)
Definition: grains_qheat.cpp:1890
GrainMakeDiffuseSingle
STATIC long GrainMakeDiffuseSingle(double Tgrain, double fracpop, avx_ptr< realnum > &flux, long nflux)
Definition: grains_qheat.cpp:194
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
TryDoubleStep
STATIC double TryDoubleStep(vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, const vector< double > &, const vector< double > &, double, double *, double, long, size_t, bool *)
Definition: grains_qheat.cpp:1586
cval
static const double cval[4]
Definition: grains_qheat.cpp:129
GrainVar::lgDColOn
bool lgDColOn
Definition: grainvar.h:496
ChargeBin::ThresSurfVal
double ThresSurfVal
Definition: grainvar.h:228
Generated by   doxygen 1.8.5