cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
temp_change.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 /*tfidle update some temperature dependent variables */
4 /*tauff compute optical depth where cloud is thin to free-free and plasma freq */
5 #include "cddefines.h"
6 #include "conv.h"
7 #include "opacity.h"
8 #include "dense.h"
9 #include "phycon.h"
10 #include "stopcalc.h"
11 #include "trace.h"
12 #include "rfield.h"
13 #include "doppvel.h"
14 #include "radius.h"
15 #include "wind.h"
16 #include "thermal.h"
17 #include "vectorize.h"
18 
19 /*tauff compute optical depth where cloud is thin to free-free and plasma freq */
20 STATIC void tauff();
21 /* On first run, fill GauntFF with gaunt factors */
22 
26 STATIC void tfidle(
27  bool lgForceUpdate);
28 
32  double TempNew ,
33  /* option to force update of all variables */
34  bool lgForceUpdate)
35 {
36 
37  DEBUG_ENTRY( "TempChange()" );
38 
39  /* set new temperature */
40  if( TempNew > phycon.TEMP_LIMIT_HIGH )
41  {
42  /* temp is too high */
43  fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
44  " is above the upper limit of the code, %.3eK.\n",
45  TempNew , phycon.TEMP_LIMIT_HIGH );
46  fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
47 
48  TempNew = phycon.TEMP_LIMIT_HIGH*0.99999;
49  lgAbort = true;
50  }
51  else if( TempNew < phycon.TEMP_LIMIT_LOW )
52  {
53  /* temp is too low */
54  fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
55  " is below the lower limit of the code, %.3eK.\n",
56  TempNew , phycon.TEMP_LIMIT_LOW );
57  fprintf(ioQQQ," Consider setting a lowest temperature with the SET TEMPERATURE FLOOR command.\n");
58  fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
59  TempNew = phycon.TEMP_LIMIT_LOW*1.00001;
60  lgAbort = true;
61  }
62  else if( TempNew < StopCalc.TeFloor )
63  {
64  if( trace.lgTrace || trace.nTrConvg>=2 )
65  fprintf(ioQQQ,"temp_change: temp change floor hit, TempNew=%.3e TeFloor=%.3e, "
66  "setting constant temperature, nTotalIoniz=%li\n",
67  TempNew , StopCalc.TeFloor , conv.nTotalIoniz);
68  /* temperature floor option -
69  * go to constant temperature calculation if temperature
70  * falls below floor */
74  /*fprintf(ioQQQ,"DEBUG TempChange hit temp floor, setting const temp to %.3e\n",
75  phycon.te );*/
76  }
77  else
78  {
79  /* temp is within range */
80  phycon.te = TempNew;
81  }
82 
83  /* now update related variables */
84  tfidle( lgForceUpdate );
85  return;
86 }
91  double TempNew )
92 {
93 
94  DEBUG_ENTRY( "TempChange()" );
95 
96  /* set new temperature */
97  if( TempNew > phycon.TEMP_LIMIT_HIGH )
98  {
99  /* temp is too high */
100  fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
101  " is above the upper limit of the code, %.3eK.\n",
102  TempNew , phycon.TEMP_LIMIT_HIGH );
103  fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
104 
105  TempNew = phycon.TEMP_LIMIT_HIGH*0.99999;
106  lgAbort = true;
107  }
108  else if( TempNew < phycon.TEMP_LIMIT_LOW )
109  {
110  /* temp is too low */
111  fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
112  " is below the lower limit of the code, %.3eK.\n",
113  TempNew , phycon.TEMP_LIMIT_LOW );
114  fprintf(ioQQQ," Consider setting a lowest temperature with the SET TEMPERATURE FLOOR command.\n");
115  fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
116  TempNew = phycon.TEMP_LIMIT_LOW*1.00001;
117  lgAbort = true;
118  }
119  else
120  {
121  /* temp is within range */
122  phycon.te = TempNew;
123  }
124 
125  /* now update related variables */
126  tfidle( false );
127  return;
128 }
129 
130 void tfidle(
131  /* option to force update of all variables */
132  bool lgForceUpdate)
133 {
134  static double tgffused=-1.;
135  static double ttused = 0.;
136  static bool lgZLogSet = false;
137 
138  DEBUG_ENTRY( "tfidle()" );
139 
140  /* called with lgForceUpdate true in zero.c, when we must update everything */
141  if( lgForceUpdate )
142  {
143  ttused = -1.;
144  tgffused = -1.;
145  }
146 
147  /* check that eden not negative */
148  if( dense.eden <= 0. )
149  {
150  fprintf( ioQQQ, "tfidle called with a zero or negative electron density,%10.2e\n",
151  dense.eden );
152  TotalInsanity();
153  }
154 
155  /* check that temperature not negative */
156  if( phycon.te <= 0. )
157  {
158  fprintf( ioQQQ, "tfidle called with a negative electron temperature,%10.2e\n",
159  phycon.te );
160  TotalInsanity();
161  }
162 
163  /* one time only, set up array of logs of charge squared */
164  if( !lgZLogSet )
165  {
166  for( long nelem=0; nelem<LIMELM; ++nelem )
167  {
168  /* this array is used to modify the log temperature array
169  * defined below, for hydrogenic species of charge nelem+1 */
170  phycon.sqlogz[nelem] = log10( POW2(nelem+1.) );
171  }
172  lgZLogSet = true;
173  }
174 
175  if( ! fp_equal( phycon.te, ttused ) )
176  {
177  ttused = phycon.te;
179  /* current temperature in various units */
180  phycon.te_eV = phycon.te/EVDEGK;
181  phycon.te_ryd = phycon.te/TE1RYD;
182  phycon.te_wn = phycon.te / T1CM;
183 
185  phycon.sqrte = sqrt(phycon.te);
186  thermal.halfte = 0.5/phycon.te;
187  thermal.tsq1 = 1./phycon.tesqrd;
189  phycon.teinv = 1./phycon.te;
190 
191  phycon.alogte = log10(phycon.te);
192  phycon.alnte = log(phycon.te);
193 
194  phycon.telogn[0] = phycon.alogte;
195  for( int i=1; i < 7; i++ )
196  {
197  phycon.telogn[i] = phycon.telogn[i-1]*phycon.telogn[0];
198  }
199 
200  phycon.te10 = pow(phycon.te,0.10);
206 
207  phycon.te01 = pow(phycon.te,0.01);
213 
214  phycon.te001 = pow(phycon.te,0.001);
220  /*>>>chng 06 June 30 -Humeshkar Nemala*/
221  phycon.te0001 = pow(phycon.te ,0.0001);
227 
228  }
229 
230  /* >>>chng 99 nov 23, removed this line, so back to old method of h coll */
231  /* used for hydrogenic collisions */
232  /*
233  * following electron density has approximate correction for neutrals
234  * corr of hi*1.7e-4 accounts for col ion by HI;
235  * >>refer H0 correction for collisional contribution Drawin, H.W. 1969, Zs Phys 225, 483.
236  * also quoted in Dalgarno & McCray 1972
237  * extensive discussion of this in
238  *>>refer H0 collisions Lambert, D.L.
239  * used EdenHCorr instead
240  * edhi = eden + hi * 1.7e-4
241  */
243  /* dense.HCorrFac is unity by default and changed with the set HCOR command */
244  dense.xIonDense[ipHYDROGEN][0]*1.7e-4 * dense.HCorrFac;
246 
247  /*>>chng 93 jun 04,
248  * term with hi added June 4, 93, to account for warm pdr */
249  /* >>chng 05 jan 05, Will Henney noticed that 1.e-4 used here is not same as
250  * 1.7e-4 used for EdenHCorr, which had rewritten the expression.
251  * change so that edensqte uses EdenHCorr rather than reevaluating */
252  /*dense.edensqte = ((dense.eden + dense.xIonDense[ipHYDROGEN][0]/1e4)/phycon.sqrte);*/
254  dense.cdsqte = dense.edensqte*COLL_CONST;
255  dense.SqrtEden = sqrt(dense.eden);
256 
257  /* rest have to do with radiation field and frequency mesh which may not be defined yet */
258  if( !lgRfieldMalloced || !rfield.lgMeshSetUp() )
259  return;
260 
261  /* correction factors for induced recombination,
262  * also used as Boltzmann factors
263  * check for zero is because ContBoltz is zeroed out in initialization
264  * of code, its possible this is a constant density grid of models
265  * in which the code is called as a subroutine */
266  /* >>chng 01 aug 21, must also test on size of continuum nflux because
267  * conintitemp can increase nflux then call this routine, although
268  * temp may not have changed */
269  if( ! fp_equal(tgffused, phycon.te) || rfield.ContBoltz[0] <= 0. )
270  {
271  tgffused = phycon.te;
272  for( long i=0; i < rfield.nflux_with_check; ++i )
274  /* atom_level2 uses ContBoltz to see whether the rates will be significant.
275  * If the numbers did not agree then this test would be flawed, resulting in
276  * div by zero */
278  for( long i=0; i < rfield.nflux_with_check; ++i )
280  /* this is Boltzmann factor averaged over the width of the cell */
282  for( long i=0; i < rfield.nflux_with_check; ++i )
285  for( long i=0; i < rfield.nflux_with_check; ++i )
286  {
289  }
291  /* ipMaxBolt is number of non-zero cells, so non-zero up through ipMaxBolt-1 */
293  for( long i=rfield.nflux_with_check-1; i >= 0; --i )
294  {
295  if( rfield.ContBoltz[i] > 0. )
296  break;
297  rfield.ipMaxBolt = i;
298  }
299  }
300 
301  /* find frequency where thin to bremsstrahlung or plasma frequency */
302  tauff();
303 }
304 
305 /*tauff compute optical depth where cloud is thin to free-free and plasma freq */
306 STATIC void tauff()
307 {
308  realnum fac;
309 
310  /* simply return if space not yet allocated */
311  if( !lgOpacMalloced )
312  return;
313 
314  DEBUG_ENTRY( "tauff()" );
315 
316  if( !conv.nTotalIoniz )
318 
319  /* routine sets variable ipEnergyBremsThin, index for lowest cont cell that is optically thin */
320  /* find frequency where continuum thin to free-free */
321  while( rfield.ipEnergyBremsThin < rfield.nflux &&
323  {
325  }
326 
327  /* TFF will be frequency where cloud becomes optically thin to bremsstrahlung
328  * >>chng 96 may 7, had been 2, change as per Kevin Volk bug report */
330  {
331  /* tau can be zero when plasma frequency is within energy grid, */
334  fac = MAX2(fac,0.f);
336  }
337  else
338  {
339  rfield.EnergyBremsThin = 0.f;
340  }
341 
342  /* did not include plasma freq before
343  * function returns larger of these two frequencies */
345 
346  /* now increment ipEnergyBremsThin still further, until above plasma frequency */
347  while( rfield.ipEnergyBremsThin < rfield.nflux &&
349  {
351  }
352  return;
353 }
354 
356 {
357  ASSERT( massAMU > 0. );
358  // force a fairly conservative upper limit
359  ASSERT( massAMU < 10000. );
360 
361  /* usually TurbVel =0, reset with turbulence command
362  * cm/s here, but was entered in km/s with command */
363  double turb2 = POW2(DoppVel.TurbVel);
364 
365  /* this is option to dissipate the turbulence. DispScale is entered with
366  * dissipate keyword on turbulence command. The velocity is reduced here,
367  * by an assumed exponential scale, and also adds heat */
368  if( DoppVel.DispScale > 0. )
369  {
370  /* square of exp depth dependence */
371  turb2 *= sexp( 2.*radius.depth / DoppVel.DispScale );
372  }
373 
374  /* in case of D-Critical flow include initial velocity as
375  * a component of turbulence */
376  if( ! ( wind.lgBallistic() || wind.lgStatic() ) )
377  {
378  turb2 += POW2(wind.windv0);
379  }
380 
381  realnum width = (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/massAMU+turb2);
382  ASSERT( width > 0.f );
383  return width;
384 }
385 
387 {
388 #if 0
389  /* usually TurbVel =0, reset with turbulence command
390  * cm/s here, but was entered in km/s with command */
391  double turb2 = POW2(DoppVel.TurbVel);
392 
393  /* this is option to dissipate the turbulence. DispScale is entered with
394  * dissipate keyword on turbulence command. The velocity is reduced here,
395  * by an assumed exponential scale, and also adds heat */
396  if( DoppVel.DispScale > 0. )
397  {
398  /* square of exp depth dependence */
399  turb2 *= sexp( 2.*radius.depth / DoppVel.DispScale );
400  }
401 
402  /* in case of D-Critical flow include initial velocity as
403  * a component of turbulence */
404  if( ! ( wind.lgBallistic() || wind.lgStatic() ) )
405  {
406  turb2 += POW2(wind.windv0);
407  }
408 #endif
409  DEBUG_ENTRY( "GetAveVelocity()" );
410 
411  /* this is average (NOT rms) particle speed for Maxwell distribution, Mihalas 70, 9-70 */
412  fixit("turbulence was included here for molecules but not ions. Now neither. Resolve.");
413  return (realnum)sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT*phycon.te/massAMU);
414 }
double depth
Definition: radius.h:31
double te007
Definition: phycon.h:58
t_thermal thermal
Definition: thermal.cpp:6
double te03
Definition: phycon.h:58
long int ipEnergyBremsThin
Definition: rfield.h:228
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double te_update
Definition: thermal.h:148
bool lgOpacMalloced
Definition: cdinit.cpp:40
double widflx(size_t i) const
Definition: mesh.h:147
double EdenHCorr
Definition: dense.h:227
t_opac opac
Definition: opacity.cpp:5
double * ContBoltzAvg
Definition: rfield.h:137
realnum EnergyBremsThin
Definition: rfield.h:229
STATIC void tauff()
realnum windv0
Definition: wind.h:11
double cdsqte
Definition: dense.h:246
realnum HCorrFac
Definition: dense.h:121
double te02
Definition: phycon.h:58
long int ipMaxBolt
Definition: rfield.h:232
t_StopCalc StopCalc
Definition: stopcalc.cpp:7
t_conv conv
Definition: conv.cpp:5
double te0002
Definition: phycon.h:58
t_phycon phycon
Definition: phycon.cpp:6
realnum TurbVel
Definition: doppvel.h:21
sys_float sexp(sys_float x)
Definition: service.cpp:1095
bool lgTemperatureConstant
Definition: thermal.h:44
FILE * ioQQQ
Definition: cddefines.cpp:7
double sqlogz[LIMELM]
Definition: phycon.h:86
t_DoppVel DoppVel
Definition: doppvel.cpp:5
double te_eV
Definition: phycon.h:24
double anu(size_t i) const
Definition: mesh.h:111
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:31
STATIC void tfidle(bool lgForceUpdate)
t_dense dense
Definition: global.cpp:15
double te005
Definition: phycon.h:58
double te003
Definition: phycon.h:58
void vexp(const double x[], double y[], long nlo, long nhi)
long int nflux_with_check
Definition: rfield.h:51
double te004
Definition: phycon.h:58
double edensqte
Definition: dense.h:241
const double TEMP_LIMIT_LOW
Definition: phycon.h:121
Wind wind
Definition: wind.cpp:5
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_trace trace
Definition: trace.cpp:5
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
double te07
Definition: phycon.h:58
bool lgBallistic(void) const
Definition: wind.h:31
double te01
Definition: phycon.h:58
#define POW2
Definition: cddefines.h:983
bool lgMeshSetUp() const
Definition: mesh.h:75
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
double tsq1
Definition: thermal.h:142
double te001
Definition: phycon.h:58
double teinv
Definition: phycon.h:33
double te0005
Definition: phycon.h:58
double te0004
Definition: phycon.h:58
const double TEMP_LIMIT_HIGH
Definition: phycon.h:123
t_rfield rfield
Definition: rfield.cpp:9
realnum DispScale
Definition: doppvel.h:51
float realnum
Definition: cddefines.h:124
double * ContBoltz
Definition: rfield.h:126
double alnte
Definition: phycon.h:95
realnum GetDopplerWidth(realnum massAMU)
realnum GetAveVelocity(realnum massAMU)
double te05
Definition: phycon.h:58
int nTrConvg
Definition: trace.h:27
double halfte
Definition: thermal.h:142
t_radius radius
Definition: radius.cpp:5
double anumin(size_t i) const
Definition: mesh.h:139
double te0007
Definition: phycon.h:58
long int nTotalIoniz
Definition: conv.h:159
#define ASSERT(exp)
Definition: cddefines.h:617
double * ContBoltzHelp1
Definition: rfield.h:129
double te0001
Definition: phycon.h:58
realnum ConstTemp
Definition: thermal.h:56
bool lgRfieldMalloced
Definition: cdinit.cpp:39
const int LIMELM
Definition: cddefines.h:307
double te0003
Definition: phycon.h:58
double te40
Definition: phycon.h:58
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double te_ryd
Definition: phycon.h:27
double te10
Definition: phycon.h:58
realnum EdenHCorr_f
Definition: dense.h:229
double te20
Definition: phycon.h:58
double * ContBoltzHelp2
Definition: rfield.h:130
double te002
Definition: phycon.h:58
double eden
Definition: dense.h:201
realnum ** TauAbsGeo
Definition: opacity.h:90
#define MAX2(a, b)
Definition: cddefines.h:828
double TeFloor
Definition: stopcalc.h:33
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
double te70
Definition: phycon.h:58
double te90
Definition: phycon.h:58
double alogte
Definition: phycon.h:92
double pow(double x, int i)
Definition: cddefines.h:782
double telogn[7]
Definition: phycon.h:86
double SqrtEden
Definition: dense.h:223
bool lgStatic(void) const
Definition: wind.h:24
double te_wn
Definition: phycon.h:30
double sqrte
Definition: phycon.h:58
#define fixit(a)
Definition: cddefines.h:416
double * vexp_arg
Definition: rfield.h:133
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:348
void vexpm1(const double x[], double y[], long nlo, long nhi)
double te04
Definition: phycon.h:58
realnum plsfrq
Definition: rfield.h:428
long int nflux
Definition: rfield.h:48
double te30
Definition: phycon.h:58
double te32
Definition: phycon.h:58
bool lgAbort
Definition: cddefines.cpp:10
double tesqrd
Definition: phycon.h:36