cloudy trunk
Loading...
Searching...
No Matches
temp_change.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*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 "physconst.h"
7#include "conv.h"
8#include "opacity.h"
9#include "iso.h"
10#include "dense.h"
11#include "phycon.h"
12#include "stopcalc.h"
13#include "continuum.h"
14#include "trace.h"
15#include "rfield.h"
16#include "doppvel.h"
17#include "radius.h"
18#include "wind.h"
19#include "thermal.h"
20#include "conv.h"
21
22/*tauff compute optical depth where cloud is thin to free-free and plasma freq */
23STATIC void tauff(void);
24/* On first run, fill GauntFF with gaunt factors */
25STATIC void FillGFF(void);
26/* Interpolate on GauntFF to calc gaunt at current temp, phycon.te */
27STATIC realnum InterpolateGff( long charge , double ERyd );
28STATIC int LinterpTable(realnum **t, realnum *v, long int lta, long int ltb, realnum x, realnum *a, long int *pipx);
29STATIC int LinterpVector(realnum **t, realnum *v, long lta , long ltb, realnum *yy , long ny, realnum **a);
30STATIC void fhunt(realnum *xx, long int n, realnum x, long int *j);
31
35STATIC void tfidle(
36 bool lgForceUpdate);
37
38static long lgGffNotFilled = true;
39
40const long N_TE_GFF = 41;
41static long N_PHOTON_GFF; /* will cover full energy range full range in one-tenth dec steps */
42static realnum ***GauntFF;
44/* the array of logs of temperatures at which GauntFF is defined */
46/* the array of logs of u at which GauntFF is defined */
48
52 double TempNew ,
53 /* option to force update of all variables */
54 bool lgForceUpdate)
55{
56
57 DEBUG_ENTRY( "TempChange()" );
58
59 /* set new temperature */
60 if( TempNew > phycon.TEMP_LIMIT_HIGH )
61 {
62 /* temp is too high */
63 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
64 " is above the upper limit of the code, %.3eK.\n",
65 TempNew , phycon.TEMP_LIMIT_HIGH );
66 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
67
68 TempNew = phycon.TEMP_LIMIT_HIGH*0.99999;
69 lgAbort = true;
70 }
71 else if( TempNew < phycon.TEMP_LIMIT_LOW )
72 {
73 /* temp is too low */
74 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
75 " is below the lower limit of the code, %.3eK.\n",
76 TempNew , phycon.TEMP_LIMIT_LOW );
77 fprintf(ioQQQ," Consider setting a lowest temperature with the SET TEMPERATURE FLOOR command.\n");
78 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
79 TempNew = phycon.TEMP_LIMIT_LOW*1.00001;
80 lgAbort = true;
81 }
82 else if( TempNew < StopCalc.TeFloor )
83 {
84 if( trace.lgTrace || trace.nTrConvg>=2 )
85 fprintf(ioQQQ,"temp_change: temp change floor hit, TempNew=%.3e TeFloor=%.3e, "
86 "setting constant temperature, nTotalIoniz=%li\n",
87 TempNew , StopCalc.TeFloor , conv.nTotalIoniz);
88 /* temperature floor option -
89 * go to constant temperature calculation if temperature
90 * falls below floor */
91 thermal.lgTemperatureConstant = true;
92 thermal.ConstTemp = (realnum)StopCalc.TeFloor;
93 phycon.te = thermal.ConstTemp;
94 /*fprintf(ioQQQ,"DEBUG TempChange hit temp floor, setting const temp to %.3e\n",
95 phycon.te );*/
96 }
97 else
98 {
99 /* temp is within range */
100 phycon.te = TempNew;
101 }
102
103 /* now update related variables */
104 tfidle( lgForceUpdate );
105 return;
106}
107
111 double TempNew )
112{
113
114 DEBUG_ENTRY( "TempChange()" );
115
116 /* set new temperature */
117 if( TempNew > phycon.TEMP_LIMIT_HIGH )
118 {
119 /* temp is too high */
120 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
121 " is above the upper limit of the code, %.3eK.\n",
122 TempNew , phycon.TEMP_LIMIT_HIGH );
123 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
124
125 TempNew = phycon.TEMP_LIMIT_HIGH*0.99999;
126 lgAbort = true;
127 }
128 else if( TempNew < phycon.TEMP_LIMIT_LOW )
129 {
130 /* temp is too low */
131 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
132 " is below the lower limit of the code, %.3eK.\n",
133 TempNew , phycon.TEMP_LIMIT_LOW );
134 fprintf(ioQQQ," Consider setting a lowest temperature with the SET TEMPERATURE FLOOR command.\n");
135 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
136 TempNew = phycon.TEMP_LIMIT_LOW*1.00001;
137 lgAbort = true;
138 }
139 else
140 {
141 /* temp is within range */
142 phycon.te = TempNew;
143 }
144
145 /* now update related variables */
146 tfidle( false );
147 return;
148}
149
151 /* option to force update of all variables */
152 bool lgForceUpdate)
153{
154 static double tgffused=-1.,
155 tgffsued2=-1.;
156 static long int nff_defined=-1;
157 static long maxion = 0, oldmaxion = 0;
158 static double ttused = 0.;
159 static bool lgZLogSet = false;
160 bool lgGauntF;
161 long int ion;
162 long int i,
163 nelem,
164 if1,
165 ipTe,
166 ret;
167 double fac,
168 fanu;
169
170 DEBUG_ENTRY( "tfidle()" );
171
172 /* called with lgForceUpdate true in zero.c, when we must update everything */
173 if( lgForceUpdate )
174 {
175 ttused = -1.;
176 tgffused = -1.;
177 tgffsued2 = -1.;
178 }
179
180 /* check that eden not negative */
181 if( dense.eden <= 0. )
182 {
183 fprintf( ioQQQ, "tfidle called with a zero or negative electron density,%10.2e\n",
184 dense.eden );
186 }
187
188 /* check that temperature not negative */
189 if( phycon.te <= 0. )
190 {
191 fprintf( ioQQQ, "tfidle called with a negative electron temperature,%10.2e\n",
192 phycon.te );
194 }
195
196 /* one time only, set up array of logs of charge squared */
197 if( !lgZLogSet )
198 {
199 for( nelem=0; nelem<LIMELM; ++nelem )
200 {
201 /* this array is used to modify the log temperature array
202 * defined below, for hydrogenic species of charge nelem+1 */
203 phycon.sqlogz[nelem] = log10( POW2(nelem+1.) );
204 }
205 lgZLogSet = true;
206 }
207
208 if( ! fp_equal( phycon.te, ttused ) )
209 {
210 ttused = phycon.te;
211 thermal.te_update = phycon.te;
212 /* current temperature in various units */
213 phycon.te_eV = phycon.te/EVDEGK;
214 phycon.te_ryd = phycon.te/TE1RYD;
215 phycon.te_wn = phycon.te / T1CM;
216
217 phycon.tesqrd = POW2(phycon.te);
218 phycon.sqrte = sqrt(phycon.te);
219 thermal.halfte = 0.5/phycon.te;
220 thermal.tsq1 = 1./phycon.tesqrd;
221 phycon.te32 = phycon.te*phycon.sqrte;
222 phycon.teinv = 1./phycon.te;
223
224 phycon.alogte = log10(phycon.te);
225 phycon.alnte = log(phycon.te);
226
227 phycon.telogn[0] = phycon.alogte;
228 for( i=1; i < 7; i++ )
229 {
230 phycon.telogn[i] = phycon.telogn[i-1]*phycon.telogn[0];
231 }
232
233 phycon.te10 = pow(phycon.te,0.10);
234 phycon.te20 = phycon.te10 * phycon.te10;
235 phycon.te30 = phycon.te20 * phycon.te10;
236 phycon.te40 = phycon.te30 * phycon.te10;
237 phycon.te70 = phycon.sqrte * phycon.te20;
238 phycon.te90 = phycon.te70 * phycon.te20;
239
240 phycon.te01 = pow(phycon.te,0.01);
241 phycon.te02 = phycon.te01 * phycon.te01;
242 phycon.te03 = phycon.te02 * phycon.te01;
243 phycon.te04 = phycon.te02 * phycon.te02;
244 phycon.te05 = phycon.te03 * phycon.te02;
245 phycon.te07 = phycon.te05 * phycon.te02;
246
247 phycon.te001 = pow(phycon.te,0.001);
248 phycon.te002 = phycon.te001 * phycon.te001;
249 phycon.te003 = phycon.te002 * phycon.te001;
250 phycon.te004 = phycon.te002 * phycon.te002;
251 phycon.te005 = phycon.te003 * phycon.te002;
252 phycon.te007 = phycon.te005 * phycon.te002;
253 /*>>>chng 06 June 30 -Humeshkar Nemala*/
254 phycon.te0001 = pow(phycon.te ,0.0001);
255 phycon.te0002 = phycon.te0001 * phycon.te0001;
256 phycon.te0003 = phycon.te0002 * phycon.te0001;
257 phycon.te0004 = phycon.te0002 * phycon.te0002;
258 phycon.te0005 = phycon.te0003 * phycon.te0002;
259 phycon.te0007 = phycon.te0005 * phycon.te0002;
260
261 }
262
263 /* >>>chng 99 nov 23, removed this line, so back to old method of h coll */
264 /* used for hydrogenic collisions */
265 /*
266 * following electron density has approximate correction for neutrals
267 * corr of hi*1.7e-4 accounts for col ion by HI;
268 * >>refer H0 correction for collisional contribution Drawin, H.W. 1969, Zs Phys 225, 483.
269 * also quoted in Dalgarno & McCray 1972
270 * extensive discussion of this in
271 *>>refer H0 collisions Lambert, D.L.
272 * used EdenHCorr instead
273 * edhi = eden + hi * 1.7e-4
274 */
275 dense.EdenHCorr = dense.eden +
276 /* dense.HCorrFac is unity by default and changed with the set HCOR command */
277 dense.xIonDense[ipHYDROGEN][0]*1.7e-4 * dense.HCorrFac;
278 dense.EdenHCorr_f = (realnum)dense.EdenHCorr;
279
280 /*>>chng 93 jun 04,
281 * term with hi added June 4, 93, to account for warm pdr */
282 /* >>chng 05 jan 05, Will Henney noticed that 1.e-4 used here is not same as
283 * 1.7e-4 used for EdenHCorr, which had rewritten the expression.
284 * change so that edensqte uses EdenHCorr rather than reevaluating */
285 /*dense.edensqte = ((dense.eden + dense.xIonDense[ipHYDROGEN][0]/1e4)/phycon.sqrte);*/
286 dense.edensqte = dense.EdenHCorr/phycon.sqrte;
287 dense.cdsqte = dense.edensqte*COLL_CONST;
288 dense.SqrtEden = sqrt(dense.eden);
289
290 /* rest have to do with radiation field and frequency mesh which may not be defined yet */
291 if( !lgRfieldMalloced || !rfield.lgMeshSetUp )
292 {
293 return;
294 }
295
296 /* correction factors for induced recombination,
297 * also used as Boltzmann factors
298 * check for zero is because ContBoltz is zeroed out in initialization
299 * of code, its possible this is a constant density grid of models
300 * in which the code is called as a subroutine */
301 /* >>chng 01 aug 21, must also test on size of continuum nflux because
302 * conintitemp can increase nflux then call this routine, although
303 * temp may not have changed */
304 if( ! fp_equal(tgffused, phycon.te)
305 || rfield.ContBoltz[0] <= 0. || nff_defined<rfield.nflux )
306 {
307 tgffused = phycon.te;
308 fac = TE1RYD/phycon.te;
309 i = 0;
310 fanu = fac*rfield.anu[i];
311 /* SEXP_LIMIT is 84 in cddefines.h - it is the -ln of smallest number
312 * that sexp can handle, and is used elsewhere in the code.
313 * atom_level2 uses ContBoltz to see whether the rates will be significant.
314 * If the numbers did not agree then this test would be flawed, resulting in
315 * div by zero */
316 while( i < rfield.nupper && fanu < SEXP_LIMIT )
317 {
318 rfield.ContBoltz[i] = exp(-fanu);
319 ++i;
320 /* this is Boltzmann factor for NEXT cell */
321 fanu = fac*rfield.anu[i];
322 }
323 /* ipMaxBolt is number of cells defined, so defined up through ipMaxBolt-1 */
324 rfield.ipMaxBolt = i;
325
326 /* zero out remainder */
327 /* >>chng 01 apr 14, upper limit has been ipMaxBolt+1, off by one */
328 for( i=rfield.ipMaxBolt; i < rfield.nupper; i++ )
329 {
330 rfield.ContBoltz[i] = 0.;
331 }
332 }
333
334 /* find frequency where thin to bremsstrahlung or plasma frequency */
335 tauff();
336
337 oldmaxion = maxion;
338 maxion = 0;
339
340 /* Find highest maximum stage of ionization */
341 for( nelem = 0; nelem < LIMELM; nelem++ )
342 {
343 maxion = MAX2( maxion, dense.IonHigh[nelem] );
344 }
345
346 /* reevaluate if temperature or number of cells has changed */
347 if( !fp_equal(phycon.te,tgffsued2) ||
348 /* this test - reevaluate if upper bound of defined values is
349 * above nflux, the highest point. This only triggers in
350 * large grids when continuum source gets harder */
351 nff_defined < rfield.nflux ||
352 /* this occurs when now have more stages of ionization than in previous defn */
353 maxion > oldmaxion )
354 {
355 static bool lgFirstRunDone = false;
356 long lowion;
357 /* >>chng 02 jan 10, only reevaluate needed states of ionization */
358 if( fp_equal(phycon.te,tgffsued2) && nff_defined >= rfield.nflux &&
359 maxion > oldmaxion )
360 {
361 /* this case temperature did not change by much, but upper
362 * stage of ionization increased. only need to evaluate
363 * stages that have not been already done */
364 lowion = oldmaxion;
365 }
366 else
367 {
368 /* temperature changed - do them all */
369 lowion = 1;
370 }
371
372 /* if1 will certainly be set to some positive value in gffsub */
373 if1 = 1;
374
375 /* chng 02 may 16, by Ryan...one gaunt factor array for all charges */
376 /* First index is EFFECTIVE CHARGE! */
377 /* highest stage of ionization is LIMELM,
378 * index goes from 1 to LIMELM */
379
380 nff_defined = rfield.nflux;
381
382 ASSERT( if1 >= 0 );
383
384 tgffsued2 = phycon.te;
385 lgGauntF = true;
386
387 /* new gaunt factors */
388 if( lgGffNotFilled )
389 {
390 FillGFF();
391 }
392
393 if( lgFirstRunDone == false )
394 {
395 maxion = LIMELM;
396 lgFirstRunDone = true;
397 }
398
399 /* >> chng 03 jan 23, rjrw -- move a couple of loops down into
400 * subroutines, and make those subroutines generic
401 * (i.e. dependences only on arguments, may be useful elsewhere...) */
402
403 /* Make Gaunt table for new temperature */
404 ipTe = -1;
405 for( ion=1; ion<=LIMELM; ion++ )
406 {
407 /* Interpolate for all tabulated photon energies at this temperature */
408 ret = LinterpTable(GauntFF[ion], TeGFF, N_PHOTON_GFF, N_TE_GFF, (realnum)phycon.alogte, GauntFF_T[ion], &ipTe);
409 if(ret == -1)
410 {
411 fprintf(ioQQQ," LinterpTable for GffTable called with te out of bounds \n");
413 }
414 }
415
416 /* Interpolate for all ions at required photon energies -- similar
417 * to LinterpTable, but not quite similar enough... */
418 /* >>chng 04 jun 30, change nflux to nupper */
419 ret = LinterpVector(GauntFF_T+lowion, PhoGFF, maxion-lowion+1, N_PHOTON_GFF,
420 rfield.anulog, rfield.nupper,/*rfield.nflux,*/ rfield.gff + lowion);
421 if(ret == -1)
422 {
423 fprintf(ioQQQ," LinterpVector for GffTable called with photon energy out of bounds \n");
425 }
426 }
427 else
428 {
429 /* this is flag that would have been set in gffsub, and
430 * printed in debug statement below. We are not evaluating
431 * so set to -1 */
432 if1 = -1;
433 lgGauntF = false;
434 }
435
436 if( trace.lgTrace && trace.lgTrGant )
437 {
438 fprintf( ioQQQ, " tfidle; gaunt factors?" );
439 fprintf( ioQQQ, "%2c", TorF(lgGauntF) );
440
441 fprintf( ioQQQ, "%2f g 1 2=%10.2e%9.2ld flag%2f guv(1,n)%10.2e\n",
442 rfield.gff[1][0], rfield.gff[1][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1],
443 if1, rfield.gff[1][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon],
444 rfield.gff[1][rfield.nflux-1] );
445 }
446 return;
447}
448
449/*tauff compute optical depth where cloud is thin to free-free and plasma freq */
450STATIC void tauff(void)
451{
452 realnum fac;
453
454 /* simply return if space not yet allocated */
455 if( !lgOpacMalloced )
456 return;
457
458 DEBUG_ENTRY( "tauff()" );
459
460 if( !conv.nTotalIoniz )
461 rfield.ipEnergyBremsThin = 0;
462
463 /* routine sets variable ipEnergyBremsThin, index for lowest cont cell that is optically thin */
464 /* find frequency where continuum thin to free-free */
465 while( rfield.ipEnergyBremsThin < rfield.nflux &&
466 opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] >= 1. )
467 {
468 ++rfield.ipEnergyBremsThin;
469 }
470
471 /* TFF will be frequency where cloud becomes optically thin to bremsstrahlung
472 * >>chng 96 may 7, had been 2, change as per Kevin Volk bug report */
473 if( rfield.ipEnergyBremsThin > 1 && opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] > 0.001 )
474 {
475 /* tau can be zero when plasma frequency is within energy grid, */
476 fac = (1.f - opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1])/(opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] -
477 opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1]);
478 fac = MAX2(fac,0.f);
479 rfield.EnergyBremsThin = rfield.anu[rfield.ipEnergyBremsThin-1] + rfield.widflx[rfield.ipEnergyBremsThin-1]*fac;
480 }
481 else
482 {
483 rfield.EnergyBremsThin = 0.f;
484 }
485
486 /* did not include plasma freq before
487 * function returns larger of these two frequencies */
488 rfield.EnergyBremsThin = MAX2(rfield.plsfrq,rfield.EnergyBremsThin);
489
490 /* now increment ipEnergyBremsThin still further, until above plasma frequency */
491 while( rfield.ipEnergyBremsThin < rfield.nflux &&
492 rfield.anu[rfield.ipEnergyBremsThin] <= rfield.EnergyBremsThin )
493 {
494 ++rfield.ipEnergyBremsThin;
495 }
496 return;
497}
498
500{
501 ASSERT( massAMU > 0. );
502 // force a fairly conservative upper limit
503 ASSERT( massAMU < 10000. );
504
505 /* usually TurbVel =0, reset with turbulence command
506 * cm/s here, but was entered in km/s with command */
507 double turb2 = POW2(DoppVel.TurbVel);
508
509 /* this is option to dissipate the turbulence. DispScale is entered with
510 * dissipate keyword on turbulence command. The velocity is reduced here,
511 * by an assumed exponential scale, and also adds heat */
512 if( DoppVel.DispScale > 0. )
513 {
514 /* square of exp depth dependence */
515 turb2 *= sexp( 2.*radius.depth / DoppVel.DispScale );
516 }
517
518 /* in case of D-Critical flow include initial velocity as
519 * a component of turbulence */
520 if( ! ( wind.lgBallistic() || wind.lgStatic() ) )
521 {
522 turb2 += POW2(wind.windv0);
523 }
524
525 realnum width = (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/massAMU+turb2);
526 ASSERT( width > 0.f );
527 return width;
528}
529
531{
532#if 0
533 /* usually TurbVel =0, reset with turbulence command
534 * cm/s here, but was entered in km/s with command */
535 double turb2 = POW2(DoppVel.TurbVel);
536
537 /* this is option to dissipate the turbulence. DispScale is entered with
538 * dissipate keyword on turbulence command. The velocity is reduced here,
539 * by an assumed exponential scale, and also adds heat */
540 if( DoppVel.DispScale > 0. )
541 {
542 /* square of exp depth dependence */
543 turb2 *= sexp( 2.*radius.depth / DoppVel.DispScale );
544 }
545
546 /* in case of D-Critical flow include initial velocity as
547 * a component of turbulence */
548 if( ! ( wind.lgBallistic() || wind.lgStatic() ) )
549 {
550 turb2 += POW2(wind.windv0);
551 }
552#endif
553
554 /* this is average (NOT rms) particle speed for Maxwell distribution, Mihalas 70, 9-70 */
555 fixit(); // turbulence was included here for molecules but not ions. Now neither. Resolve.
556 return (realnum)sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT*phycon.te/massAMU);
557}
558
559STATIC void FillGFF( void )
560{
561
562 long i1,i2,i3,j,charge,GffMAGIC = 100804;
563 double Temp, photon;
564 bool lgEOL;
565
566 DEBUG_ENTRY( "FillGFF()" );
567
568 const int chLine_LENGTH = 1000;
569 char chLine[chLine_LENGTH];
570
571 FILE *ioDATA;
572
573 for( long i = 0; i < N_TE_GFF; i++ )
574 {
575 TeGFF[i] = 0.25f*i;
576 }
577 /* >>chng 06 feb 14, assert thrown at T == 1e10 K, Ryan Porter proposes this fix */
578 TeGFF[N_TE_GFF-1] += 0.01f;
579
580 /* number of photon energies */
581 N_PHOTON_GFF = (int)( 20. * ( log10(rfield.egamry) - log10(rfield.emm) ) ) + 20;
582
583 PhoGFF = (realnum*)MALLOC( sizeof(realnum)*(unsigned)N_PHOTON_GFF );
584
585 for( long i = 0; i < N_PHOTON_GFF; i++ )
586 {
587 PhoGFF[i] = 0.05f*i + log10(rfield.emm) - 0.5;
588 }
589
590 GauntFF = (realnum***)MALLOC( sizeof(realnum**)*(unsigned)(LIMELM+2) );
591 for( long i = 1; i <= LIMELM; i++ )
592 {
593 GauntFF[i] = (realnum**)MALLOC( sizeof(realnum*)*(unsigned)N_PHOTON_GFF );
594 for( j = 0; j < N_PHOTON_GFF; j++ )
595 {
596 GauntFF[i][j] = (realnum*)MALLOC( sizeof(realnum)*(unsigned)N_TE_GFF );
597 }
598 }
599
600 GauntFF_T = (realnum**)MALLOC( sizeof(realnum*)*(unsigned)(LIMELM+2) );
601 for( long i = 1; i <= LIMELM; i++ )
602 {
603 GauntFF_T[i] = (realnum*)MALLOC( sizeof(realnum)*(unsigned)N_PHOTON_GFF );
604 }
605
606 if( !rfield.lgCompileGauntFF )
607 {
608 if( trace.lgTrace )
609 fprintf( ioQQQ," FillGFF opening gauntff.dat:");
610
611 try
612 {
613 ioDATA = open_data( "gauntff.dat", "r" );
614 }
615 catch( cloudy_exit )
616 {
617 fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
618 fprintf( ioQQQ, "but the code runs much faster if you compile gauntff.dat!\n");
619 ioDATA = NULL;
620 }
621
622 if( ioDATA == NULL )
623 {
624 /* Do on the fly computation of Gff's instead. */
625 for( charge=1; charge<=LIMELM; charge++ )
626 {
627 for( long i=0; i < N_PHOTON_GFF; i++ )
628 {
629 photon = pow((realnum)10.f,PhoGFF[i]);
630 for(j=0; j<N_TE_GFF; j++)
631 {
632 Temp = pow((realnum)10.f,TeGFF[j]);
633 GauntFF[charge][i][j] = (realnum)cont_gaunt_calc( Temp, (double)charge, photon );
634 }
635 }
636 }
637 }
638 else
639 {
640 /* check that magic number is ok */
641 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
642 {
643 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
645 }
646 long i = 1;
647 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
648 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
649 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
650
651 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
652 {
653 fprintf( ioQQQ,
654 " FillGFF: the version of gauntff.dat is not the current version.\n" );
655 fprintf( ioQQQ,
656 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
657 GffMAGIC ,
659 N_TE_GFF,
660 i1 , i2 , i3 );
661 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
662 fprintf( ioQQQ,
663 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
665 }
666
667 /* now read in the data */
668 for( charge = 1; charge <= LIMELM; charge++ )
669 {
670 for( i = 0; i < N_PHOTON_GFF; i++ )
671 {
672 /* get next line image */
673 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
674 {
675 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
677 }
678 /* each line starts with charge and energy level ( index in rfield ) */
679 i3 = 1;
680 i1 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
681 i2 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
682 /* check that these numbers are correct */
683 if( i1!=charge || i2!=i )
684 {
685 fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n");
686 fprintf( ioQQQ,
687 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
689 }
690
691 /* loop over temperatures to produce array of free free gaunt factors */
692 for(j = 0; j < N_TE_GFF; j++)
693 {
694 GauntFF[charge][i][j] = (realnum)FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
695
696 if( lgEOL )
697 {
698 fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n");
699 fprintf( ioQQQ,
700 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
702 }
703 }
704 }
705
706 }
707
708 /* check that magic number is ok */
709 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
710 {
711 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
713 }
714 i = 1;
715 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
716 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
717 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
718
719 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
720 {
721 fprintf( ioQQQ,
722 " FillGFF: the version of gauntff.dat is not the current version.\n" );
723 fprintf( ioQQQ,
724 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
725 GffMAGIC ,
727 N_TE_GFF,
728 i1 , i2 , i3 );
729 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
730 fprintf( ioQQQ,
731 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
733 }
734
735 /* close the data file */
736 fclose( ioDATA );
737 }
738 if( trace.lgTrace )
739 fprintf( ioQQQ," - opened and read ok.\n");
740 }
741 else
742 {
743 /* option to create table of gaunt factors,
744 * executed with the compile gaunt command */
745 FILE *ioGFF;
746
747 /*GffMAGIC the magic number for the table of recombination coefficients */
748 ioGFF = open_data( "gauntff.dat" , "w", AS_LOCAL_ONLY );
749 fprintf(ioGFF,"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
750 GffMAGIC ,
752 N_TE_GFF,
754 N_TE_GFF );
755
756 for( charge = 1; charge <= LIMELM; charge++ )
757 {
758 for( long i=0; i < N_PHOTON_GFF; i++ )
759 {
760 fprintf(ioGFF, "%li\t%li", charge, i );
761 /* loop over temperatures to produce array of gaunt factors */
762 for(j = 0; j < N_TE_GFF; j++)
763 {
764 /* Store gaunt factor in N_TE_GFF half dec steps */
765 Temp = pow((realnum)10.f,TeGFF[j]);
766 photon = pow((realnum)10.f,PhoGFF[i]);
767 GauntFF[charge][i][j] = (realnum)cont_gaunt_calc( Temp, (double)charge, photon );
768 fprintf(ioGFF, "\t%f", GauntFF[charge][i][j] );
769 }
770 fprintf(ioGFF, "\n" );
771 }
772 }
773
774 /* end the file with the same information */
775 fprintf(ioGFF,"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
776 GffMAGIC ,
778 N_TE_GFF,
780 N_TE_GFF );
781
782 fclose( ioGFF );
783
784 fprintf( ioQQQ, "FillGFF: compilation complete, gauntff.dat created.\n" );
785 }
786
787 lgGffNotFilled = false;
788
789 /* We have already checked the validity of cont_gaunt_calc in sanitycheck.c.
790 * Now we check to see if the InterpolateGff routine also works correctly. */
791 {
792 double gaunt, error;
793 double tempsave = phycon.te;
794 long logu, loggamma2;
795
796 for( logu=-4; logu<=4; logu++)
797 {
798 for(loggamma2=-4; loggamma2<=4; loggamma2++)
799 {
800 double SutherlandGff[9][9]=
801 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
802 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
803 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
804 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
805 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
806 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
807 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
808 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
809 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};
810
811 phycon.te = (TE1RYD/pow(10.,(double)loggamma2));
812 phycon.alogte = log10(phycon.te);
813 double ERyd = pow(10.,(double)(logu-loggamma2));
814 if( ERyd > rfield.emm && ERyd < rfield.egamry )
815 {
816 gaunt = InterpolateGff( 1, ERyd );
817 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
818 if( error>0.05 )
819 {
820 fprintf(ioQQQ," PROBLEM DISASTER tfidle found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
821 logu, loggamma2, error );
823 }
824 }
825 }
826 }
827 phycon.te = tempsave;
828 phycon.alogte = log10(phycon.te);
829 }
830
831 return;
832}
833
834/* Interpolate Gff at some temperature */
835STATIC realnum InterpolateGff( long charge , double ERyd )
836{
837 double GauntAtLowerPho, GauntAtUpperPho;
838 static long int ipTe=-1, ipPho=-1;
839 double gaunt = 0.;
840 double xmin , xmax;
841 long i;
842
843 DEBUG_ENTRY( "InterpolateGff()" );
844
845 if( ipTe < 0 )
846 {
847 /* te totally unknown */
848 if( phycon.alogte < TeGFF[0] || phycon.alogte > TeGFF[N_TE_GFF-1] )
849 {
850 fprintf(ioQQQ," InterpolateGff called with te out of bounds \n");
852 }
853 /* now search for temperature */
854 for( i=0; i<N_TE_GFF-1; ++i )
855 {
856 if( phycon.alogte > TeGFF[i] && phycon.alogte <= TeGFF[i+1] )
857 {
858 /* found the temperature - use it */
859 ipTe = i;
860 break;
861 }
862 }
863 ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1) );
864
865 }
866 else if( phycon.alogte < TeGFF[ipTe] )
867 {
868 /* temp is too low, must also lower ipTe */
869 ASSERT( phycon.alogte > TeGFF[0] );
870 /* decrement ipTe until it is correct */
871 while( phycon.alogte < TeGFF[ipTe] && ipTe > 0)
872 --ipTe;
873 }
874 else if( phycon.alogte > TeGFF[ipTe+1] )
875 {
876 /* temp is too high */
877 ASSERT( phycon.alogte <= TeGFF[N_TE_GFF-1] );
878 /* increment ipTe until it is correct */
879 while( phycon.alogte > TeGFF[ipTe+1] && ipTe < N_TE_GFF-1)
880 ++ipTe;
881 }
882
883 ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1) );
884
885 /* ipTe should now be valid */
886 ASSERT( phycon.alogte >= TeGFF[ipTe] && phycon.alogte <= TeGFF[ipTe+1] && ipTe < N_TE_GFF-1 );
887
888 /***************/
889 /* This bit is completely analogous to the above, but for the photon vector instead of temp. */
890 if( ipPho < 0 )
891 {
892 if( log10(ERyd) < PhoGFF[0] || log10(ERyd) > PhoGFF[N_PHOTON_GFF-1] )
893 {
894 fprintf(ioQQQ," InterpolateGff called with photon energy out of bounds \n");
896 }
897 for( i=0; i<N_PHOTON_GFF-1; ++i )
898 {
899 if( log10(ERyd) > PhoGFF[i] && log10(ERyd) <= PhoGFF[i+1] )
900 {
901 ipPho = i;
902 break;
903 }
904 }
905 ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1) );
906
907 }
908 else if( log10(ERyd) < PhoGFF[ipPho] )
909 {
910 ASSERT( log10(ERyd) >= PhoGFF[0] );
911 while( log10(ERyd) < PhoGFF[ipPho] && ipPho > 0)
912 --ipPho;
913 }
914 else if( log10(ERyd) > PhoGFF[ipPho+1] )
915 {
916 ASSERT( log10(ERyd) <= PhoGFF[N_PHOTON_GFF-1] );
917 while( log10(ERyd) > PhoGFF[ipPho+1] && ipPho < N_PHOTON_GFF-1)
918 ++ipPho;
919 }
920 ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1) );
921 ASSERT( log10(ERyd)>=PhoGFF[ipPho]
922 && log10(ERyd)<=PhoGFF[ipPho+1] && ipPho<N_PHOTON_GFF-1 );
923
924 /* Calculate the answer...must interpolate on two variables.
925 * First interpolate on T, at both the lower and upper photon energies.
926 * Then interpolate between these results for the right photon energy. */
927
928 GauntAtLowerPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) *
929 (GauntFF[charge][ipPho][ipTe+1] - GauntFF[charge][ipPho][ipTe]) + GauntFF[charge][ipPho][ipTe];
930
931 GauntAtUpperPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) *
932 (GauntFF[charge][ipPho+1][ipTe+1] - GauntFF[charge][ipPho+1][ipTe]) + GauntFF[charge][ipPho+1][ipTe];
933
934 gaunt = (log10(ERyd) - PhoGFF[ipPho]) / (PhoGFF[ipPho+1] - PhoGFF[ipPho]) *
935 (GauntAtUpperPho - GauntAtLowerPho) + GauntAtLowerPho;
936
937 xmax = MAX4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1],
938 GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] );
939 ASSERT( gaunt <= xmax );
940
941 xmin = MIN4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1],
942 GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] );
943 ASSERT( gaunt >= xmin );
944
945 ASSERT( gaunt > 0. );
946
947 return (realnum)gaunt;
948}
949
950/* Interpolate in table t[lta][ltb], with physical values for the
951 second index given by v[ltb], for values x, and put results in
952 a[lta]; store the index found if that's useful; assumes v[] is
953 sorted */
954STATIC int LinterpTable(realnum **t, realnum *v, long int lta, long int ltb, realnum x, realnum *a, long int *pipx)
955{
956 long int ipx=-1;
957 realnum frac;
958 long i;
959
960 DEBUG_ENTRY( "LinterpTable()" );
961
962 if(pipx != NULL)
963 ipx = *pipx;
964
965 fhunt (v,ltb,x,&ipx); /* search for index */
966 if(pipx != NULL)
967 *pipx = ipx;
968
969 if( ipx == -1 || ipx == ltb )
970 {
971 return -1;
972 }
973
974 ASSERT( (ipx >=0) && (ipx < ltb-1) );
975 ASSERT( x >= v[ipx] && x <= v[ipx+1]);
976
977 frac = (x - v[ipx]) / (v[ipx+1] - v[ipx]);
978 for( i=0; i<lta; i++ )
979 {
980 a[i] = frac*t[i][ipx+1]+(1.f-frac)*t[i][ipx];
981 }
982
983 return 0;
984}
985
986/* Interpolate in table t[lta][ltb], with physical values for the second index given by v[ltb],
987 for values yy[ny], and put results in a[lta][ly] */
988STATIC int LinterpVector(realnum **t, realnum *v, long lta , long ltb, realnum *yy, long ly, realnum **a)
989{
990 realnum yl, frac;
991 long i, j, n;
992
993 DEBUG_ENTRY( "LinterpVector()" );
994
995 if( yy[0] < v[0] || yy[ly-1] > v[ltb-1] )
996 {
997 return -1;
998 }
999
1000 n = 0;
1001 yl = yy[n];
1002 for(j = 1; j < ltb && n < ly; j++) {
1003 while (yl < v[j] && n < ly) {
1004 frac = (yl-v[j-1])/(v[j]-v[j-1]);
1005 for(i = 0; i < lta; i++)
1006 a[i][n] = frac*t[i][j]+(1.f-frac)*t[i][j-1];
1007 n++;
1008 if(n == ly)
1009 break;
1010 ASSERT( yy[n] >= yy[n-1] );
1011 yl = yy[n];
1012 }
1013 }
1014
1015 return 0;
1016}
1017STATIC void fhunt(realnum *xx, long int n, realnum x, long int *j)
1018{
1019 /*lint -e731 boolean argument to equal / not equal */
1020 long int jl, jm, jh, in;
1021 int up;
1022
1023 jl = *j;
1024 up = (xx[n-1] > xx[0]);
1025 if(jl < 0 || jl >= n)
1026 {
1027 jl = -1;
1028 jh = n;
1029 }
1030 else
1031 {
1032 in = 1;
1033 if((x >= xx[jl]) == up)
1034 {
1035 if(jl == n-1)
1036 {
1037 *j = jl;
1038 return;
1039 }
1040 jh = jl + 1;
1041 while ((x >= xx[jh]) == up)
1042 {
1043 jl = jh;
1044 in += in;
1045 jh += in;
1046 if(jh >= n)
1047 {
1048 jh = n;
1049 break;
1050 }
1051 }
1052 }
1053 else
1054 {
1055 if(jl == 0)
1056 {
1057 *j = -1;
1058 return;
1059 }
1060 jh = jl--;
1061 while ((x < xx[jl]) == up)
1062 {
1063 jh = jl;
1064 in += in;
1065 jl -= in;
1066 if(jl <= 0)
1067 {
1068 jl = 0;
1069 break;
1070 }
1071 }
1072 }
1073 }
1074 while (jh-jl != 1)
1075 {
1076 jm = (jh+jl)/2;
1077 if((x > xx[jm]) == up)
1078 jl = jm;
1079 else
1080 jh = jm;
1081 }
1082 *j = jl;
1083 return;
1084}
1085 /*lint +e731 boolean argument to equal / not equal */
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
#define MIN4(a, b, c, d)
Definition cddefines.h:771
const double SEXP_LIMIT
Definition cddefines.h:1476
#define MALLOC(exp)
Definition cddefines.h:501
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
char TorF(bool l)
Definition cddefines.h:710
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define MAX4(a, b, c, d)
Definition cddefines.h:792
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
bool lgRfieldMalloced
Definition cdinit.cpp:98
bool lgOpacMalloced
Definition cdinit.cpp:100
double cont_gaunt_calc(double temp, double z, double photon)
t_conv conv
Definition conv.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
@ AS_LOCAL_ONLY
Definition cpu.h:208
t_dense dense
Definition dense.cpp:24
t_DoppVel DoppVel
Definition doppvel.cpp:5
#define chLine_LENGTH
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double PI
Definition physconst.h:29
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double T1CM
Definition physconst.h:167
UNUSED const double EVDEGK
Definition physconst.h:186
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
UNUSED const double COLL_CONST
Definition physconst.h:229
UNUSED const double TE1RYD
Definition physconst.h:183
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_StopCalc StopCalc
Definition stopcalc.cpp:5
STATIC void tauff(void)
void TempChange(double TempNew, bool lgForceUpdate)
STATIC void tfidle(bool lgForceUpdate)
STATIC int LinterpTable(realnum **t, realnum *v, long int lta, long int ltb, realnum x, realnum *a, long int *pipx)
realnum GetDopplerWidth(realnum massAMU)
STATIC void FillGFF(void)
static long N_PHOTON_GFF
realnum GetAveVelocity(realnum massAMU)
static long lgGffNotFilled
static realnum TeGFF[N_TE_GFF]
STATIC realnum InterpolateGff(long charge, double ERyd)
static realnum * PhoGFF
static realnum *** GauntFF
const long N_TE_GFF
static realnum ** GauntFF_T
STATIC void fhunt(realnum *xx, long int n, realnum x, long int *j)
STATIC int LinterpVector(realnum **t, realnum *v, long lta, long ltb, realnum *yy, long ny, realnum **a)
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
Wind wind
Definition wind.cpp:5