cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_gaunt.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 
4 #include "cddefines.h"
5 #include "phycon.h"
6 #include "rfield.h"
7 #include "opacity.h"
8 #include "iso.h"
9 #include "dense.h"
10 #include "mole.h"
11 #include "vectorize.h"
12 #include "atmdat_gaunt.h"
13 
14 // use 3rd-order interpolation
15 static const size_t N_GFF_INTERP = 4;
16 
17 static const char* gauntff_file = "gauntff_merged_Zxx.dat";
18 
19 static const long gaunt_magic = 20140510L;
20 
21 void t_gaunt::p_read_table(const char* fnam, long nelem)
22 {
23  DEBUG_ENTRY( "t_gaunt::p_read_table()" );
24 
25  if( p_gff.size() > 0 && !isnan(p_gff[nelem][0][0][0]) )
26  return;
27 
28  // generate file name
29  string fnam1 = fnam;
30  ostringstream oss;
31  oss << setw(2) << setfill('0') << nelem+1;
32  string::size_type ptr = fnam1.find( "xx" );
33  ASSERT( ptr != string::npos );
34  fnam1.replace( ptr, 2, oss.str() );
35 
36  fstream io;
37  open_data( io, fnam1.c_str(), mode_r );
38  string line;
39  // read magic number
40  getline( io, line );
41  istringstream iss( line );
42  long magic;
43  iss >> magic;
44  if( magic != gaunt_magic )
45  {
46  fprintf( ioQQQ, " t_gaunt::p_read_table() found wrong magic number in file %s.\n", fnam );
47  fprintf( ioQQQ, " I expected to find version: %ld\n", gaunt_magic );
49  }
50  // read dimensions of the table
51  getline( io, line );
52  iss.str(line);
53  iss >> p_np_gam2 >> p_np_u;
54  ASSERT( p_np_gam2 >= N_GFF_INTERP && p_np_u >= N_GFF_INTERP );
55  // read start value for log(gamma^2)
56  getline( io, line );
57  iss.str(line);
58  iss >> p_lg_gam2_min;
59  // read start value for log(u)
60  getline( io, line );
61  iss.str(line);
62  iss >> p_lg_u_min;
63  // read step size in dex
64  getline( io, line );
65  iss.str(line);
66  iss >> p_step;
67  // read atomic number
68  getline( io, line );
69  iss.str(line);
70  int Z;
71  iss >> Z;
72  ASSERT( Z == nelem+1 );
73 
74  if( p_gff.size() == 0 )
75  {
77  // set to NaN to avoid uninitialized use
78  p_gff.invalidate();
79 
80  p_lg_gam2.resize( p_np_gam2 );
81  for( size_t ipg2=0; ipg2 < p_np_gam2; ++ipg2 )
82  p_lg_gam2[ipg2] = p_lg_gam2_min + double(ipg2)*p_step;
83  p_lg_gam2_max = p_lg_gam2.back();
84 
85  p_lg_u.resize( p_np_u );
86  for( size_t ipu=0; ipu < p_np_u; ++ipu )
87  p_lg_u[ipu] = p_lg_u_min + double(ipu)*p_step;
88  p_lg_u_max = p_lg_u.back();
89  }
90  else
91  {
92  // check that table parameters are the same
93  ASSERT( p_np_gam2 == p_lg_gam2.size() );
94  ASSERT( p_np_u == p_lg_u.size() );
95  ASSERT( fp_equal( p_lg_gam2_min, p_lg_gam2[0] ) );
96  ASSERT( fp_equal( p_lg_u_min, p_lg_u[0] ) );
97  ASSERT( fp_equal( p_step, p_lg_u[1]-p_lg_u[0], 2*int(abs(p_lg_u[0]/p_step)) ) );
98  }
99 
100  // next two lines are comments
101  getline( io, line );
102  getline( io, line );
103 
104  for( size_t ipu=0; ipu < p_np_u; ++ipu )
105  {
106  getline( io, line );
107  iss.str(line);
108  for( size_t ipg2=0; ipg2 < p_np_gam2; ++ipg2 )
109  {
110  double val;
111  iss >> val;
112  p_gff[nelem][0][ipu][ipg2] = log(val);
113  }
114  // reset eof flag, otherwise subsequent reads would fail
115  iss.clear( iss.rdstate() & ~iss.eofbit );
116 
117  // calculate higher derivatives for Newton interpolation
118  for( size_t i=1; i < N_GFF_INTERP; ++i )
119  for( size_t ipg2=0; ipg2 < p_np_gam2-i; ++ipg2 )
120  p_gff[nelem][i][ipu][ipg2] =
121  (p_gff[nelem][i-1][ipu][ipg2+1]-p_gff[nelem][i-1][ipu][ipg2])/
122  (double(i)*p_step);
123  }
124 
125  // the remainder of the file contains the uncertainties of the gff values
126  // we will not read those here...
127 
128  if( io.fail() || iss.fail() )
129  {
130  fprintf( ioQQQ, " An error occurred while reading the file %s. Bailing out.\n", fnam );
132  }
133 }
134 
135 double t_gaunt::gauntff(long Z, double Te, double anu)
136 {
137  DEBUG_ENTRY( "t_gaunt::gauntff()" );
138 
139  ASSERT( Z > 0 && Z <= LIMELM );
140 
142 
143  double gam2 = pow2(double(Z))*TE1RYD/Te;
144  double u = anu*TE1RYD/Te;
145 
146  ASSERT( gam2 > 0. && u > 0. );
147 
148  double lg_gam2 = log10(gam2);
149  double lg_u = log10(u);
150 
151  ASSERT( p_lg_gam2_min <= lg_gam2 && lg_gam2 <= p_lg_gam2_max );
152  ASSERT( p_lg_u_min <= lg_u && lg_u <= p_lg_u_max );
153 
154  long ipg2 = min(max(size_t((lg_gam2-p_lg_gam2_min)/p_step),1)-1,p_np_gam2-N_GFF_INTERP);
155  long ipu = min(max(size_t((lg_u-p_lg_u_min)/p_step),1)-1,p_np_u-N_GFF_INTERP);
156 
157  double interp[N_GFF_INTERP];
158 
159  for( size_t i=0; i < N_GFF_INTERP; ++i )
160  interp[i] = lagrange(get_ptr(&p_lg_gam2[ipg2]), get_ptr(&p_gff[Z-1][0][ipu+i][ipg2]),
161  N_GFF_INTERP, lg_gam2);
162 
163  return exp(lagrange(get_ptr(&p_lg_u[ipu]), interp, N_GFF_INTERP, lg_u));
164 }
165 
166 void t_gaunt::p_gauntff_vec(long Z, double Te, const double anulog10[], long nflux)
167 {
168  DEBUG_ENTRY( "t_gaunt::p_gauntff_vec()" );
169 
171 
172  if( fp_equal( Te, p_gff_ion[Z].Te_used )
173  && p_gff_ion[Z].nflux_used >= nflux )
174  return;
175 
176  if( p_vexp_arg.size() == 0 )
178 
179  if( p_gff_ion[Z].gff.size() == 0 )
180  {
181  p_gff_ion[Z].gff.resize( rfield.nflux_with_check );
183  }
184 
185  if( !fp_equal( Te, p_gff_ion[Z].Te_used ) )
186  p_gauntff_vec_sub( Z, Te, anulog10, 0, nflux );
187  else if( p_gff_ion[Z].nflux_used < nflux )
188  p_gauntff_vec_sub( Z, Te, anulog10, p_gff_ion[Z].nflux_used, nflux );
189 
190  p_gff_ion[Z].Te_used = Te;
191  p_gff_ion[Z].nflux_used = nflux;
192 }
193 
194 void t_gaunt::p_gauntff_vec_sub(long Z, double Te, const double anulog10[], long nmin, long nmax)
195 {
196  DEBUG_ENTRY( "t_gaunt::p_gauntff_vec_sub()" );
197 
198  ASSERT( Z > 0 && Z <= LIMELM );
199 
200  vector< realnum, allocator_avx<realnum> >& gff = p_gff_ion[Z].gff;
201 
202  double lg_gam2 = log10(pow2(double(Z))*TE1RYD/Te);
203  ASSERT( p_lg_gam2_min <= lg_gam2 && lg_gam2 <= p_lg_gam2_max );
204  long ipg2 = min(max(size_t((lg_gam2-p_lg_gam2_min)/p_step),1)-1,p_np_gam2-N_GFF_INTERP);
205 
206  double lg_u0 = log10(TE1RYD/Te);
207  ASSERT( p_lg_u_min <= lg_u0+anulog10[nmin] && lg_u0+anulog10[nmax-1] <= p_lg_u_max );
208  size_t ipumin = min(max(size_t((lg_u0+anulog10[nmin]-p_lg_u_min)/p_step),1)-1,p_np_u-N_GFF_INTERP);
209  size_t ipumax = min(max(size_t((lg_u0+anulog10[nmax-1]-p_lg_u_min)/p_step),1)-1,p_np_u-N_GFF_INTERP)+N_GFF_INTERP;
210 
212  for( size_t ipu=ipumin; ipu < ipumax; ++ipu )
213  {
214  // Newton interpolation in the gamma^2 direction
215  size_t i = 1;
216  const double* plu = get_ptr(&p_lg_gam2[ipg2]);
217  double val = p_gff[Z-1][0][ipu][ipg2];
218  double fac = lg_gam2 - *plu++;
219  while( true )
220  {
221  val += fac*p_gff[Z-1][i][ipu][ipg2];
222  if( ++i == N_GFF_INTERP )
223  break;
224  fac *= lg_gam2 - *plu++;
225  }
226  interp[0][ipu] = val;
227  }
228 
229  // fill in higher derivatives for Newton interpolation
230  for( size_t i=1; i < N_GFF_INTERP; ++i )
231  for( size_t ipu=ipumin; ipu < ipumax-i; ++ipu )
232  interp[i][ipu] = (interp[i-1][ipu+1]-interp[i-1][ipu])/(double(i)*p_step);
233 
234  size_t ipu = ipumin;
235 
236  // this loop burns lots of CPU time, so it should be well optimized
237  for( long j=nmin; j < nmax; ++j )
238  {
239  double lg_u = lg_u0 + anulog10[j];
240  while( lg_u >= p_lg_u[ipu+2] )
241  ipu++;
242  // Newton interpolation in the u direction
243  size_t i = 1;
244  const double* plu = get_ptr(&p_lg_u[ipu]);
245  const double* pin = get_ptr(&interp[0][ipu]);
246  double gval = *pin;
247  double fac = lg_u - *plu++;
248  while( true )
249  {
250  pin += p_np_u;
251  gval += fac*(*pin);
252  if( ++i == N_GFF_INTERP )
253  break;
254  fac *= lg_u - *plu++;
255  }
256  p_vexp_arg[j] = realnum(gval);
257  }
258  vexp( get_ptr(p_vexp_arg), get_ptr(gff), nmin, nmax );
259 }
260 
261 void t_gaunt::p_setup_brems(long ion, double Te)
262 {
263  DEBUG_ENTRY( "t_gaunt::p_setup_brems()" );
264 
265  long limit = min( rfield.ipMaxBolt, rfield.nflux );
266 
267  if( ion > 0 && ion < LIMELM+1 )
268  {
269  if( fp_equal( Te, p_cache_vec[ion].Te_used )
270  && p_cache_vec[ion].nhi >= limit )
271  return;
272 
273  if( p_cache_vec[ion].brems_vec.size() == 0 )
274  p_cache_vec[ion].brems_vec.resize( rfield.nflux_with_check );
275 
276  p_gauntff_vec( ion, Te, rfield.anulog10ptr(), limit );
277 
278  p_cache_vec[ion].Te_used = Te;
279  p_cache_vec[ion].nhi = limit;
280  double ion2 = pow2(double(ion));
281 
282  for( long i=0; i < p_cache_vec[ion].nhi; ++i )
283  p_cache_vec[ion].brems_vec[i] =
284  ion2*p_gff_ion[ion].gff[i]*rfield.widflx(i)*rfield.ContBoltz[i];
285  }
286  else if( ion == -1 )
287  {
288  // special case for H-
289  if( fp_equal( Te, p_hminus_vec.Te_used )
290  && p_hminus_vec.nhi >= limit )
291  return;
292 
293  if( p_hminus_vec.brems_vec.size() == 0 )
295 
296  p_gauntff_vec( 1, Te, rfield.anulog10ptr(), limit );
297 
298  p_hminus_vec.Te_used = Te;
299  p_hminus_vec.nhi = limit;
300 
301  for( long i=0; i < p_hminus_vec.nhi; ++i )
303  p_gff_ion[1].gff[i]*opac.OpacStack[i-1+opac.iphmra]*
305  }
306  else
307  {
308  TotalInsanity();
309  }
310 }
311 
312 double t_gaunt::brems_cool(long ion, double Te)
313 {
314  DEBUG_ENTRY( "t_gaunt::brems_cool()" );
315 
316  ASSERT( Te > 0. );
317 
318  long limit = min( rfield.ipMaxBolt, rfield.nflux );
319 
320  if( ion > 0 && ion < LIMELM+1 )
321  {
322  if( fp_equal( Te, p_cache[ion].Te_used )
323  && p_cache[ion].nlo == rfield.ipEnergyBremsThin
324  && p_cache[ion].nhi == limit )
325  return p_cache[ion].brems_sum;
326 
327  p_setup_brems( ion, Te );
328 
329  p_cache[ion].Te_used = Te;
331  p_cache[ion].nhi = limit;
332 
333  p_cache[ion].brems_sum = 0.;
334  for( long i=p_cache[ion].nlo; i < p_cache[ion].nhi; ++i )
335  p_cache[ion].brems_sum += p_cache_vec[ion].brems_vec[i];
336 
337  return p_cache[ion].brems_sum;
338  }
339  else if( ion == -1 )
340  {
341  // special case for H-
342  if( fp_equal( Te, p_hminus.Te_used )
344  && p_hminus.nhi == limit )
345  return p_hminus.brems_sum;
346 
347  p_setup_brems( ion, Te );
348 
349  p_hminus.Te_used = Te;
351  p_hminus.nhi = limit;
352 
353  p_hminus.brems_sum = 0.;
354  for( long i=p_hminus.nlo; i < p_hminus.nhi; ++i )
355  // OpacStack contains the ratio of the H- to H brems cross section
357 
358  return p_hminus.brems_sum;
359  }
360  else
361  {
362  TotalInsanity();
363  }
364 }
365 
366 void t_gaunt::brems_rt(long ion, double Te, double abun, vector<double>& arr)
367 {
368  DEBUG_ENTRY( "t_gaunt::brems_rt()" );
369 
370  ASSERT( Te > 0. );
371 
372  size_t limit = min( rfield.ipMaxBolt, rfield.nflux );
373 
374  if( ion > 0 && ion < LIMELM+1 )
375  {
376  ASSERT( arr.size() >= limit );
377  p_setup_brems( ion, Te );
378  for( size_t i=0; i < limit; ++i )
379  arr[i] += abun*p_cache_vec[ion].brems_vec[i];
380  }
381  else if( ion == -1 )
382  {
383  ASSERT( arr.size() >= limit );
384  p_setup_brems( ion, Te );
385  for( size_t i=0; i < limit; ++i )
386  arr[i] += abun*p_hminus_vec.brems_vec[i];
387  }
388  else
389  {
390  TotalInsanity();
391  }
392 }
393 
394 void t_gaunt::brems_opac(long ion, double Te, double abun, double arr[])
395 {
396  DEBUG_ENTRY( "t_gaunt::brems_opac()" );
397 
398  ASSERT( Te > 0. );
399 
400  if( ion > 0 && ion < LIMELM+1 )
401  {
403  for( long i=0; i < rfield.nflux; ++i )
404  arr[i] += abun*pow2(double(ion))*p_gff_ion[ion].gff[i];
405  }
406  else if( ion == -1 )
407  {
409  for( long i=0; i < rfield.nflux; ++i )
410  arr[i] += abun*p_gff_ion[1].gff[i]*opac.OpacStack[i-1+opac.iphmra];
411  }
412  else
413  {
414  TotalInsanity();
415  }
416 }
417 
419 {
420  DEBUG_ENTRY( "t_gaunt::brems_sum_ions()" );
421 
422  sum.den_Hm = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
423  sum.den_Hp = dense.xIonDense[ipHYDROGEN][1];
424  sum.den_Hp += findspecieslocal("H2+")->den;
425  sum.den_Hp += findspecieslocal("H3+")->den;
426  sum.den_Hep = dense.xIonDense[ipHELIUM][1];
427  sum.den_Hepp = dense.xIonDense[ipHELIUM][2];
428  sum.den_ion[0] = 0.;
429  for( long ion=1; ion < LIMELM+1; ++ion )
430  {
431  sum.den_ion[ion] = 0.;
432  for( long nelem=ipLITHIUM; nelem < LIMELM; ++nelem )
433  if( dense.lgElmtOn[nelem] && ion <= nelem+1 )
434  sum.den_ion[ion] += dense.xIonDense[nelem][ion];
435  }
436  /* add molecular ions */
437  for( long ipMol = 0; ipMol < mole_global.num_calc; ipMol++ )
438  {
439  if( !mole_global.list[ipMol]->isMonatomic() &&
440  mole_global.list[ipMol]->isIsotopicTotalSpecies() &&
441  mole_global.list[ipMol]->charge > 0 &&
442  mole_global.list[ipMol]->label != "H2+" &&
443  mole_global.list[ipMol]->label != "H3+" )
444  {
445  ASSERT( mole_global.list[ipMol]->charge < LIMELM+1 );
446  sum.den_ion[mole_global.list[ipMol]->charge] += mole.species[ipMol].den;
447  }
448  }
449 }
450 
451 void dgaunt_check(double Z)
452 {
453  DEBUG_ENTRY( "dgaunt_check()" );
454 
455  exit_type exit_status = ES_SUCCESS;
456 
457  double toler = 0.002;
458 
459  for( long logu=-13; logu <= 12; logu++ )
460  {
461  int i = 0;
462  double gam2[3], gaunt[3];
463  double u = exp10((double)(logu));
464  for( long loggamma2=-500; loggamma2 <= 900; loggamma2++ )
465  {
466  gam2[i] = exp10(double(loggamma2)/100.);
467  double Te = pow2(Z)*(TE1RYD/gam2[i]);
468  double ERyd = pow2(Z)*u/gam2[i];
469  if( Te > phycon.TEMP_LIMIT_LOW && Te < phycon.TEMP_LIMIT_HIGH &&
470  ERyd > rfield.emm() && ERyd < rfield.egamry() )
471  {
472  gaunt[i] = t_gaunt::Inst().gauntff( long(Z), Te, ERyd );
473  if( i < 2 )
474  ++i;
475  else {
476  if( abs((gaunt[2]+gaunt[0])/(2.*gaunt[1]) - 1.) > toler ) {
477  fprintf( ioQQQ, " PROBLEM DISASTER cont_gaunt_calc discontinuity in"
478  " gam2 direction at gam2 = %e u = %e\n", gam2[i], u );
479  fprintf( ioQQQ, " gam2 = %e gaunt = %e\n", gam2[0], gaunt[0] );
480  fprintf( ioQQQ, " gam2 = %e gaunt = %e\n", gam2[1], gaunt[1] );
481  fprintf( ioQQQ, " gam2 = %e gaunt = %e\n", gam2[2], gaunt[2] );
482  exit_status = ES_FAILURE;
483  i = 0;
484  }
485  else {
486  gam2[0] = gam2[1];
487  gam2[1] = gam2[2];
488  gaunt[0] = gaunt[1];
489  gaunt[1] = gaunt[2];
490  }
491  }
492  }
493  }
494  }
495 
496  for( long loggamma2=-5; loggamma2 <= 9; loggamma2++ )
497  {
498  int i = 0;
499  double u[3], gaunt[3];
500  double gam2 = exp10(double(loggamma2));
501  double Te = pow2(Z)*(TE1RYD/gam2);
502  for( long logu=-1300; logu <= 1200; logu++ )
503  {
504  u[i] = exp10((double)(logu)/100.);
505  double ERyd = pow2(Z)*u[i]/gam2;
506  if( Te > phycon.TEMP_LIMIT_LOW && Te < phycon.TEMP_LIMIT_HIGH &&
507  ERyd > rfield.emm() && ERyd < rfield.egamry() )
508  {
509  gaunt[i] = t_gaunt::Inst().gauntff( long(Z), Te, ERyd );
510  if( i < 2 )
511  ++i;
512  else {
513  if( abs((gaunt[2]+gaunt[0])/(2.*gaunt[1]) - 1.) > toler ) {
514  fprintf( ioQQQ, " PROBLEM DISASTER cont_gaunt_calc discontinuity in"
515  " u direction at gam2 = %e u = %e\n", gam2, u[2] );
516  fprintf( ioQQQ, " u = %e gaunt = %e\n", u[0], gaunt[0] );
517  fprintf( ioQQQ, " u = %e gaunt = %e\n", u[1], gaunt[1] );
518  fprintf( ioQQQ, " u = %e gaunt = %e\n", u[2], gaunt[2] );
519  exit_status = ES_FAILURE;
520  i = 0;
521  }
522  else {
523  u[0] = u[1];
524  u[1] = u[2];
525  gaunt[0] = gaunt[1];
526  gaunt[1] = gaunt[2];
527  }
528  }
529  }
530  }
531  }
532 
533  if( exit_status != ES_SUCCESS )
534  cdEXIT(exit_status);
535 }
double emm() const
Definition: mesh.h:93
t_mole_global mole_global
Definition: mole.cpp:7
vector< double, allocator_avx< double > > ContBoltz
Definition: rfield.h:124
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
double brems_cool(long ion, double Te)
double * OpacStack
Definition: opacity.h:164
qList st
Definition: iso.h:482
double exp10(double x)
Definition: cddefines.h:1368
long int ipEnergyBremsThin
Definition: rfield.h:226
T * get_ptr(T *v)
Definition: cddefines.h:1109
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
double widflx(size_t i) const
Definition: mesh.h:156
t_opac opac
Definition: opacity.cpp:5
int num_calc
Definition: mole.h:362
size_t p_np_gam2
Definition: atmdat_gaunt.h:49
void set_NaN(sys_float &x)
Definition: cpu.cpp:906
t_brems_vec p_cache_vec[LIMELM+1]
Definition: atmdat_gaunt.h:66
void p_gauntff_vec(long Z, double Te, const double anulog10[], long nflux)
double den_Hepp
Definition: atmdat_gaunt.h:40
double p_lg_gam2_max
Definition: atmdat_gaunt.h:53
double den_Hep
Definition: atmdat_gaunt.h:39
long int ipMaxBolt
Definition: rfield.h:230
size_type size() const
t_phycon phycon
Definition: phycon.cpp:6
void dgaunt_check(double Z)
FILE * ioQQQ
Definition: cddefines.cpp:7
molezone * findspecieslocal(const char buf[])
exit_type
Definition: cddefines.h:142
long nflux_used
Definition: atmdat_gaunt.h:13
size_t p_np_u
Definition: atmdat_gaunt.h:50
double brems_sum
Definition: atmdat_gaunt.h:23
long int iphmra
Definition: opacity.h:223
t_dense dense
Definition: global.cpp:15
static t_gaunt & Inst()
Definition: cddefines.h:209
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
void vexp(const double x[], double y[], long nlo, long nhi)
long int nflux_with_check
Definition: rfield.h:49
const double TEMP_LIMIT_LOW
Definition: phycon.h:121
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
t_brems_vec p_hminus_vec
Definition: atmdat_gaunt.h:67
const ios_base::openmode mode_r
Definition: cpu.h:267
void p_setup_brems(long ion, double Te)
const int ipH1s
Definition: iso.h:29
double Te_used
Definition: atmdat_gaunt.h:29
t_mole_local mole
Definition: mole.cpp:8
const double TEMP_LIMIT_HIGH
Definition: phycon.h:123
t_rfield rfield
Definition: rfield.cpp:9
double p_step
Definition: atmdat_gaunt.h:56
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
#define EXIT_FAILURE
Definition: cddefines.h:168
double Te_used
Definition: atmdat_gaunt.h:20
static const long gaunt_magic
void p_read_table(const char *fnam, long nelem)
long max(int a, long b)
Definition: cddefines.h:817
#define cdEXIT(FAIL)
Definition: cddefines.h:482
long min(int a, long b)
Definition: cddefines.h:762
vector< realnum, allocator_avx< realnum > > gff
Definition: atmdat_gaunt.h:14
double den_ion[LIMELM+1]
Definition: atmdat_gaunt.h:42
vector< double > brems_vec
Definition: atmdat_gaunt.h:31
double p_lg_u_min
Definition: atmdat_gaunt.h:54
bool lgElmtOn[LIMELM]
Definition: dense.h:160
double gauntff(long Z, double Te, double anu)
double lagrange(const double x[], const double y[], long n, double xval)
double p_lg_gam2_min
Definition: atmdat_gaunt.h:52
t_brems_sum p_cache[LIMELM+1]
Definition: atmdat_gaunt.h:64
vector< double > p_lg_u
Definition: atmdat_gaunt.h:59
static const char * gauntff_file
vector< realnum, allocator_avx< realnum > > p_vexp_arg
Definition: atmdat_gaunt.h:61
#define ASSERT(exp)
Definition: cddefines.h:613
double den_Hm
Definition: atmdat_gaunt.h:37
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
T pow2(T a)
Definition: cddefines.h:981
double den
Definition: mole.h:421
void p_gauntff_vec_sub(long Z, double Te, const double anulog10[], long nmin, long nmax)
double p_lg_u_max
Definition: atmdat_gaunt.h:55
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
#define isnan
Definition: cddefines.h:659
const int ipHELIUM
Definition: cddefines.h:350
double egamry() const
Definition: mesh.h:97
double Te_used
Definition: atmdat_gaunt.h:12
t_gff p_gff_ion[LIMELM+1]
Definition: atmdat_gaunt.h:63
multi_arr< double, 4 > p_gff
Definition: atmdat_gaunt.h:60
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
MoleculeList list
Definition: mole.h:365
vector< double > p_lg_gam2
Definition: atmdat_gaunt.h:58
void brems_sum_ions(t_brems_den &sum) const
void brems_opac(long ion, double Te, double abun, double arr[])
const int ipHYDROGEN
Definition: cddefines.h:349
long int nflux
Definition: rfield.h:46
const int ipLITHIUM
Definition: cddefines.h:351
void brems_rt(long ion, double Te, double abun, vector< double > &arr)
const double * anulog10ptr() const
Definition: mesh.h:136
double den_Hp
Definition: atmdat_gaunt.h:38
static const size_t N_GFF_INTERP
t_brems_sum p_hminus
Definition: atmdat_gaunt.h:65