cloudy trunk
Loading...
Searching...
No Matches
atmdat.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#include "cddefines.h"
4#include "atmdat.h"
5#include "thirdparty.h"
7
8double ****HS_He1_Xsectn;
9double ****HS_He1_Energy;
10double *****OP_Helike_Xsectn;
11double *****OP_Helike_Energy;
13
14void ReadCollisionRateTable( CollRateCoeffArray& coll_rate_table, FILE* io, FunctPtr GetIndices, long nMolLevs, long nTemps, long nTrans )
15{
16 DEBUG_ENTRY( "ReadCollisionRateTable()" );
17
18 char chLine[INPUT_LINE_LENGTH] = "";
19
20 // negative nTrans and/or nTemps will be signal that the numbers are not already known
21 // They should never be zero.
22 ASSERT( nTemps != 0 && nTrans != 0 );
23
24 // Get the row of temperatures
25 while( read_whole_line( chLine, (int)sizeof(chLine), io ) != NULL )
26 {
27 if( chLine[0] == '!' || chLine[0] == '#' || chLine[0] == '\n' )
28 continue;
29 else
30 break;
31 }
32 ASSERT( strlen( chLine ) > 0 );
33
34 // fill the temperature array
35 {
36 char *chColltemp = strtok(chLine," \t\n");
37 while( chColltemp != NULL )
38 {
39 coll_rate_table.temps.push_back( atof(chColltemp) );
40 chColltemp = strtok(NULL," \t\n");
41 }
42
43 if( nTemps < 0 )
44 nTemps = coll_rate_table.temps.size();
45 else
46 ASSERT( (unsigned)nTemps == coll_rate_table.temps.size() );
47 }
48
49 ASSERT( coll_rate_table.collrates.size() == 0 );
50 coll_rate_table.collrates.reserve( nMolLevs );
51 for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
52 {
53 coll_rate_table.collrates.reserve( ipHi, nMolLevs );
54 for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
55 {
56 coll_rate_table.collrates.reserve( ipHi, ipLo, nTemps );
57 }
58 }
59 coll_rate_table.collrates.alloc();
60 // initialize to zero
61 coll_rate_table.collrates.zero();
62
63 long ipTrans = 0;
64 while( read_whole_line( chLine, (int)sizeof(chLine), io ) != NULL )
65 {
66 if( chLine[0] == '!' || chLine[0] == '#' || chLine[0] == '\n' )
67 continue;
68
69 long i = 1;
70 long ipHi = -1, ipLo = -1;
71 (*GetIndices)( ipHi, ipLo, chLine, i );
72 ipTrans++;
73
74 // sentinel that transition will not be used for whatever reason
75 if( ipHi == -1 && ipLo == -1 )
76 continue;
77
78 // skip these lines
79 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
80 {
81 if( nTrans > 0 && ipTrans == nTrans )
82 break;
83 else
84 continue;
85 }
86
87 /* Indices between the very highest levels seem to be reversed */
88 if( ipHi < ipLo )
89 {
90 ASSERT( ipLo == nMolLevs - 1);
91 long temp = ipHi;
92 ipHi = ipLo;
93 ipLo = temp;
94 }
95
96 ASSERT( ipHi >= 0 );
97 ASSERT( ipLo >= 0 );
98
99 bool lgEOL = false;
100 for( long j = 0; j < nTemps; ++j )
101 {
102 coll_rate_table.collrates[ipHi][ipLo][j] =
103 (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
104 ASSERT( !lgEOL );
105 }
106
107 // try to read one more and assert that it gets lgEOL
108 FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
109 ASSERT( lgEOL );
110
111 {
112 enum {DEBUG_LOC=false};
113 if( DEBUG_LOC )
114 {
115 printf("The values of up and lo are %ld & %ld \n",ipHi,ipLo);
116 printf("The collision rates are");
117 for( long i = 0; i < nTemps; ++i )
118 {
119 printf( "\n %e", coll_rate_table.collrates[ipHi][ipLo][i]);
120 }
121 printf("\n");
122 }
123 }
124
125 if( nTrans > 0 && ipTrans == nTrans )
126 break;
127 }
128 ASSERT( ipTrans > 0 );
129 if( nTrans > 0 )
130 ASSERT( ipTrans == nTrans );
131
132 return;
133}
134
135double InterpCollRate( const CollRateCoeffArray& rate_table, const long& ipHi, const long& ipLo, const double& ftemp)
136{
137 DEBUG_ENTRY("InterpCollRate()");
138 double ret_collrate = 0.;
139
140 if( rate_table.temps.size() == 0 )
141 {
142 return 0.;
143 }
144
145 if( ftemp < rate_table.temps[0] )
146 {
147 ret_collrate = rate_table.collrates[ipHi][ipLo][0];
148 }
149 else if( ftemp > rate_table.temps.back() )
150 {
151 ret_collrate = rate_table.collrates[ipHi][ipLo][ rate_table.temps.size()-1 ];
152 }
153 else if( rate_table.temps.size() == 1 )
154 {
155 // lamda data files can have only one temperature (see cn.dat as of Feb. 10, 2009)
156 ret_collrate = rate_table.collrates[ipHi][ipLo][0];
157 }
158 else
159 {
160 ret_collrate = linint(&rate_table.temps[0],
161 &rate_table.collrates[ipHi][ipLo][0],
162 rate_table.temps.size(),
163 ftemp);
164 }
165
166 ASSERT( !isnan( ret_collrate ) );
167 return(ret_collrate);
168}
169
double InterpCollRate(const CollRateCoeffArray &rate_table, const long &ipHi, const long &ipLo, const double &ftemp)
Definition atmdat.cpp:135
double ***** OP_Helike_Xsectn
Definition atmdat.cpp:10
double **** HS_He1_Xsectn
Definition atmdat.cpp:8
double **** HS_He1_Energy
Definition atmdat.cpp:9
void ReadCollisionRateTable(CollRateCoeffArray &coll_rate_table, FILE *io, FunctPtr GetIndices, long nMolLevs, long nTemps, long nTrans)
Definition atmdat.cpp:14
long **** OP_Helike_NumPts
Definition atmdat.cpp:12
t_atmdat atmdat
Definition atmdat.cpp:6
double ***** OP_Helike_Energy
Definition atmdat.cpp:11
Funct * FunctPtr
Definition atmdat.h:408
#define isnan
Definition cddefines.h:620
#define ASSERT(exp)
Definition cddefines.h:578
struct t_CollRatesArray CollRateCoeffArray
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void reserve(size_type i1)
size_type size() const
vector< double > temps
Definition cddefines.h:1267
multi_arr< double, 3 > collrates
Definition cddefines.h:1269
double linint(const double x[], const double y[], long n, double xval)