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

Generated for cloudy by doxygen 1.8.3.1