cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat.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 #include "cddefines.h"
4 #include "atmdat.h"
5 #include "thirdparty.h"
7 
9 {
10  DEBUG_ENTRY( "t_atmdat::zero()" );
11  nsbig = 0;
12 
13  /***************************************************
14  * charge transfer ionization and recombination
15  ***************************************************/
16  /* HCharHeatMax, HCharCoolMax are largest fractions of heating in cur zone
17  * or cooling due to ct */
18  HCharHeatMax = 0.;
19  HCharCoolMax = 0.;
20 
21  for ( long nelem=0; nelem < t_atmdat::NCX; ++nelem)
22  {
23  CharExcIonTotal[nelem] = 0.;
24  CharExcRecTotal[nelem] = 0.;
25  }
26  HIonFrac = 0.;
27  HIonFracMax = 0.;
28  /* option to turn off all charge transfer, turned off with no charge transfer command */
29  lgCTOn = true;
30 
31  /* flag saying that charge transfer heating should be included,
32  * turned off with no CTHeat commmand */
33  HCharHeatOn = 1.;
34  for( long nelem1=0; nelem1 < t_atmdat::NCX; ++nelem1)
35  {
36  for( long nelem=0; nelem< LIMELM; ++nelem )
37  {
38  for( long ion=0; ion<LIMELM; ++ion )
39  {
40  CharExcIonOf[nelem1][nelem][ion] = 0.;
41  CharExcRecTo[nelem1][nelem][ion] = 0.;
42  }
43  }
44  }
45 
46  /* >>chng 97 jan 6, from 0 to 8.5e-10*q as per Alex Dalgarno e-mail
47  * >>chng 97 feb 6, from 8.5e-10*q 1.92e-9x as per Alex Dalgarno e-mail */
48  HCTAlex = 1.92e-9;
49 
50  lgInnerShellLine_on = true;
52  lgUTAprint = false;
53 }
54 
58 const double t_atmdat::aulThreshold = 1e-29;
59 
60 double ****HS_He1_Xsectn;
61 double ****HS_He1_Energy;
62 double *****OP_Helike_Xsectn;
63 double *****OP_Helike_Energy;
65 
66 void ReadCollisionRateTable( CollRateCoeffArray& coll_rate_table, FILE* io, FunctPtr GetIndices, long nMolLevs, long nTemps, long nTrans )
67 {
68  DEBUG_ENTRY( "ReadCollisionRateTable()" );
69 
70  char chLine[INPUT_LINE_LENGTH] = "";
71 
72  // negative nTrans and/or nTemps will be signal that the numbers are not already known
73  // They should never be zero.
74  ASSERT( nTemps != 0 && nTrans != 0 );
75 
76  // Get the row of temperatures
77  while( read_whole_line( chLine, (int)sizeof(chLine), io ) != NULL )
78  {
79  if( chLine[0] == '!' || chLine[0] == '#' || chLine[0] == '\n' )
80  continue;
81  else
82  break;
83  }
84  ASSERT( strlen( chLine ) > 0 );
85 
86  // fill the temperature array
87  {
88  char *chColltemp = strtok(chLine," \t\n");
89  while( chColltemp != NULL )
90  {
91  coll_rate_table.temps.push_back( atof(chColltemp) );
92  chColltemp = strtok(NULL," \t\n");
93  }
94 
95  if( nTemps < 0 )
96  nTemps = coll_rate_table.temps.size();
97  else
98  ASSERT( (unsigned)nTemps == coll_rate_table.temps.size() );
99  }
100 
101  ASSERT( coll_rate_table.collrates.size() == 0 );
102  coll_rate_table.collrates.reserve( nMolLevs );
103  for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
104  {
105  coll_rate_table.collrates.reserve( ipHi, nMolLevs );
106  for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
107  {
108  coll_rate_table.collrates.reserve( ipHi, ipLo, nTemps );
109  }
110  }
111  coll_rate_table.collrates.alloc();
112  // initialize to zero
113  coll_rate_table.collrates.zero();
114 
115  long ipTrans = 0;
116  while( read_whole_line( chLine, (int)sizeof(chLine), io ) != NULL )
117  {
118  if( chLine[0] == '!' || chLine[0] == '#' || chLine[0] == '\n' )
119  continue;
120 
121  long i = 1;
122  long ipHi = -1, ipLo = -1;
123  (*GetIndices)( ipHi, ipLo, chLine, i );
124  ipTrans++;
125 
126  // sentinel that transition will not be used for whatever reason
127  if( ipHi == -1 && ipLo == -1 )
128  continue;
129 
130  // skip these lines
131  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
132  {
133  if( nTrans > 0 && ipTrans == nTrans )
134  break;
135  else
136  continue;
137  }
138 
139  /* Indices between the very highest levels seem to be reversed */
140  if( ipHi < ipLo )
141  {
142  ASSERT( ipLo == nMolLevs - 1);
143  long temp = ipHi;
144  ipHi = ipLo;
145  ipLo = temp;
146  }
147 
148  ASSERT( ipHi >= 0 );
149  ASSERT( ipLo >= 0 );
150 
151  bool lgEOL = false;
152  for( long j = 0; j < nTemps; ++j )
153  {
154  coll_rate_table.collrates[ipHi][ipLo][j] =
155  (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
156  ASSERT( !lgEOL );
157  }
158 
159  // try to read one more and assert that it gets lgEOL
160  FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
161  ASSERT( lgEOL );
162 
163  {
164  enum {DEBUG_LOC=false};
165  if( DEBUG_LOC )
166  {
167  printf("The values of up and lo are %ld & %ld \n",ipHi,ipLo);
168  printf("The collision rates are");
169  for( long i = 0; i < nTemps; ++i )
170  {
171  printf( "\n %e", coll_rate_table.collrates[ipHi][ipLo][i]);
172  }
173  printf("\n");
174  }
175  }
176 
177  if( nTrans > 0 && ipTrans == nTrans )
178  break;
179  }
180  ASSERT( ipTrans > 0 );
181  if( nTrans > 0 )
182  ASSERT( ipTrans == nTrans );
183 
184  return;
185 }
186 
187 double InterpCollRate( const CollRateCoeffArray& rate_table, const long& ipHi, const long& ipLo, const double& ftemp)
188 {
189  DEBUG_ENTRY( "InterpCollRate()" );
190  double ret_collrate = 0.;
191 
192  if( rate_table.temps.size() == 0 )
193  {
194  return 0.;
195  }
196 
197  if( ftemp < rate_table.temps[0] )
198  {
199  ret_collrate = rate_table.collrates[ipHi][ipLo][0];
200  }
201  else if( ftemp > rate_table.temps.back() )
202  {
203  ret_collrate = rate_table.collrates[ipHi][ipLo][ rate_table.temps.size()-1 ];
204  }
205  else if( rate_table.temps.size() == 1 )
206  {
207  // lamda data files can have only one temperature (see cn.dat as of Feb. 10, 2009)
208  ret_collrate = rate_table.collrates[ipHi][ipLo][0];
209  }
210  else
211  {
212  ret_collrate = linint(&rate_table.temps[0],
213  &rate_table.collrates[ipHi][ipLo][0],
214  rate_table.temps.size(),
215  ftemp);
216  }
217 
218  ASSERT( !isnan( ret_collrate ) );
219  return(ret_collrate);
220 }
221 
t_atmdat atmdat
Definition: atmdat.cpp:6
multi_arr< double, 3 > collrates
Definition: atmdat.h:18
double HCharHeatMax
Definition: atmdat.h:300
double HCharHeatOn
Definition: atmdat.h:300
size_type size() const
double HCharCoolMax
Definition: atmdat.h:300
long int nsbig
Definition: atmdat.h:345
double HIonFracMax
Definition: atmdat.h:317
double CharExcRecTotal[NCX]
Definition: atmdat.h:310
double CharExcIonTotal[NCX]
Definition: atmdat.h:310
Definition: atmdat.h:461
static const double aulThreshold
Definition: atmdat.h:437
double **** HS_He1_Xsectn
Definition: atmdat.cpp:60
double InterpCollRate(const CollRateCoeffArray &rate_table, const long &ipHi, const long &ipLo, const double &ftemp)
Definition: atmdat.cpp:187
double ***** OP_Helike_Xsectn
Definition: atmdat.cpp:62
void zero()
Definition: atmdat.cpp:8
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
bool lgUTAprint
Definition: atmdat.h:433
#define ASSERT(exp)
Definition: cddefines.h:613
bool lgInnerShell_Kisielius
Definition: atmdat.h:431
void reserve(size_type i1)
const int LIMELM
Definition: cddefines.h:308
vector< double > temps
Definition: atmdat.h:16
double **** HS_He1_Energy
Definition: atmdat.cpp:61
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
#define isnan
Definition: cddefines.h:659
double ***** OP_Helike_Energy
Definition: atmdat.cpp:63
double linint(const double x[], const double y[], long n, double xval)
long **** OP_Helike_NumPts
Definition: atmdat.cpp:64
bool lgInnerShellLine_on
Definition: atmdat.h:429
double HIonFrac
Definition: atmdat.h:314
bool lgCTOn
Definition: atmdat.h:325
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
double HCTAlex
Definition: atmdat.h:321
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:300
void ReadCollisionRateTable(CollRateCoeffArray &coll_rate_table, FILE *io, FunctPtr GetIndices, long nMolLevs, long nTemps, long nTrans)
Definition: atmdat.cpp:66
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:300
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381