cloudy
trunk
Main Page
Related Pages
Data Structures
Files
File List
Globals
All
Data Structures
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
source
punch_average.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
/*punch_average routine to read in list of species to output as averages */
4
#include "
cddefines.h
"
5
#include "
cddrive.h
"
6
#include "
input.h
"
7
#include "
elementnames.h
"
8
#include "
punch.h
"
9
10
/* return value is number of lines, -1 if file could not be opened */
11
void
punch_average
(
12
/* the file we will write to */
13
long
int
ipPun,
14
/* this is the job we shall do - READ and PUNCH */
15
const
char
*chJob,
16
char
*chHeader)
17
{
18
long
int
i;
19
bool
lgEOF , lgEOL;
20
long
int
nLine;
21
22
char
chLine[
FILENAME_PATH_LENGTH_2
] ;
23
24
DEBUG_ENTRY
(
"punch_average()"
);
25
26
if
( strcmp(chJob,
"READ"
) == 0 )
27
{
28
char
chCap[
INPUT_LINE_LENGTH
];
29
char
chTemp[
INPUT_LINE_LENGTH
];
30
31
/* this is index for first line we will read. use this to count
32
* total number of species */
33
nLine =
input
.
nRead
;
34
/* very first time this routine is called, chJob is "READ" and we read
35
* in lines from the input stream. The species labels and other information
36
* are store in the punch structure. These are output in a later call
37
* to this routine with argument PUNCH */
38
39
/* number of lines we will save */
40
punch
.
nAverageList
[ipPun] = 0;
41
42
/* get the next line, and check for eof */
43
input_readarray
(chLine,&lgEOF);
44
if
( lgEOF )
45
{
46
fprintf(
ioQQQ
,
47
" Punch average hit EOF while reading list; use END to end list.\n"
);
48
cdEXIT
(EXIT_FAILURE);
49
}
50
51
/* convert line to caps */
52
strcpy( chCap, chLine );
53
caps
(chCap);
54
55
/* keep reading until we hit END */
56
while
( strncmp(chCap,
"END"
,3 ) != 0 )
57
{
58
/* count number of species we will save */
59
++
punch
.
nAverageList
[ipPun];
60
61
/* get next line and check for eof */
62
input_readarray
(chLine,&lgEOF);
63
if
( lgEOF )
64
{
65
fprintf(
ioQQQ
,
" punch averages hit EOF while reading species list; use END to end list.\n"
);
66
cdEXIT
(EXIT_FAILURE);
67
}
68
69
/* convert it to caps */
70
strcpy( chCap, chLine );
71
caps
(chCap);
72
}
73
/*# define PADEBUG*/
74
# ifdef PADEBUG
75
fprintf(
ioQQQ
,
"DEBUG punch_average %li species read in.\n"
,
76
punch
.
nAverageList
[ipPun] );
77
# endif
78
79
/* now create space that will be needed to hold these arrays */
80
81
if
(
punch
.
ipPnunit
[ipPun] == NULL )
82
{
83
/* nAverageIonList is set of ions for averages */
84
punch
.
nAverageIonList
[ipPun] = (
int
*)
MALLOC
(
85
(
size_t
)(
punch
.
nAverageList
[ipPun]*
sizeof
(int)) );
86
87
/* nAverage2ndPar is set of second parameters for averages */
88
punch
.
nAverage2ndPar
[ipPun] = (
int
*)
MALLOC
((
size_t
)
89
(
punch
.
nAverageList
[ipPun]*
sizeof
(int)) );
90
91
/* chAverageType is label for type of average */
92
punch
.
chAverageType
[ipPun] = (
char
**)
MALLOC
((
size_t
)
93
(
punch
.
nAverageList
[ipPun]*
sizeof
(
char
*)) );
94
95
/* chAverageSpeciesLabel is label for species */
96
punch
.
chAverageSpeciesLabel
[ipPun] = (
char
**)
MALLOC
((
size_t
)
97
(
punch
.
nAverageList
[ipPun]*
sizeof
(
char
*)) );
98
for
( i=0; i<
punch
.
nAverageList
[ipPun]; ++i )
99
{
100
/* create space for labels themselves */
101
punch
.
chAverageType
[ipPun][i] = (
char
*)
MALLOC
((
size_t
)
102
(5*
sizeof
(char)) );
103
104
/* chAverageSpeciesLabel is label for species */
105
punch
.
chAverageSpeciesLabel
[ipPun][i] = (
char
*)
MALLOC
((
size_t
)
106
(5*
sizeof
(char)) );
107
}
108
}
109
110
/* reset array input read to first of the species lines */
111
input
.
nRead
= nLine;
112
113
# ifdef PADEBUG
114
fprintf(
ioQQQ
,
"DEBUG punch_average %li species read in.\n"
,
115
punch
.
nAverageList
[ipPun] );
116
# endif
117
118
/* reread the input lines and save the data */
119
input_readarray
(chLine,&lgEOF);
120
if
( lgEOF )
121
{
122
fprintf(
ioQQQ
,
123
" Punch average hit EOF while reading list; use END to end list.\n"
);
124
cdEXIT
(EXIT_FAILURE);
125
}
126
127
/* convert line to caps */
128
strcpy( chCap, chLine );
129
caps
(chCap);
130
131
/* use this to count number of species, and will assert equal to
132
* total malloced above */
133
nLine = 0;
134
/* keep reading until we hit END */
135
while
( strncmp(chCap,
"END"
,3 ) != 0 )
136
{
137
/* count number of species we will save */
138
++nLine;
139
if
(
nMatch
(
"TEMP"
, chCap ))
140
{
141
/* temperature */
142
strcpy(
punch
.
chAverageType
[ipPun][nLine-1] ,
"TEMP"
);
143
}
144
else
if
(
nMatch
(
"COLU"
, chCap ))
145
{
146
/* column density */
147
strcpy(
punch
.
chAverageType
[ipPun][nLine-1] ,
"COLU"
);
148
}
149
else
if
(
nMatch
(
"IONI"
, chCap ))
150
{
151
/* ionization fraction */
152
strcpy(
punch
.
chAverageType
[ipPun][nLine-1] ,
"IONI"
);
153
}
154
else
155
{
156
fprintf(
ioQQQ
,
"PROBLEM one of the jobs TEMPerature, COLUmn density, or IONIzation, must appear.\n"
);
157
cdEXIT
(EXIT_FAILURE);
158
}
159
160
/* get element name, a string we shall pass on to the routine
161
* that computes final quantities */
162
if
( (i =
GetElem
( chCap )) < 0 )
163
{
164
/* no name found */
165
fprintf(
ioQQQ
,
"punch average did not an element on this line, sorry %s\n"
,
166
chCap );
167
cdEXIT
(EXIT_FAILURE);
168
}
169
strcpy(
punch
.
chAverageSpeciesLabel
[ipPun][nLine-1] ,
elementnames
.
chElementNameShort
[i]);
170
171
/* now get ionization stage */
172
i = 5;
173
punch
.
nAverageIonList
[ipPun][nLine-1] = (int)
174
FFmtRead
( chCap , &i,
INPUT_LINE_LENGTH
,&lgEOL );
175
if
( lgEOF )
176
{
177
/* error - needed that ionization stage */
178
NoNumb
( chCap );
179
}
180
181
/* look for volume keyword, otherwise will be radius
182
* only used for some options */
183
if
(
nMatch
(
"VOLU"
, chCap ) )
184
{
185
/* volume */
186
punch
.
nAverage2ndPar
[ipPun][nLine-1] = 1;
187
}
188
else
189
{
190
/* radius */
191
punch
.
nAverage2ndPar
[ipPun][nLine-1] = 0;
192
}
193
194
/* get next line and check for eof */
195
input_readarray
(chLine,&lgEOF);
196
if
( lgEOF )
197
{
198
fprintf(
ioQQQ
,
" punch averages hit EOF while reading species list; use END to end list.\n"
);
199
cdEXIT
(EXIT_FAILURE);
200
}
201
202
/* convert it to caps */
203
strcpy( chCap, chLine );
204
caps
(chCap);
205
}
206
207
/* these must equal or we have a major logical error */
208
ASSERT
( nLine ==
punch
.
nAverageList
[ipPun]);
209
210
# ifdef PADEBUG
211
for
( i=0; i<nLine ; ++i )
212
{
213
fprintf(
ioQQQ
,
"PDDEBUG %s %s %i %i\n"
,
214
punch
.
chAverageType
[ipPun][i],
215
punch
.
chAverageSpeciesLabel
[ipPun][i] ,
216
punch
.
nAverageIonList
[ipPun][i] ,
217
punch
.
nAverage2ndPar
[ipPun][i] );
218
}
219
# endif
220
221
/* punch headers */
222
sprintf(chHeader,
"#averages"
);
223
for
( i=0; i<nLine ; ++i )
224
{
225
sprintf(chTemp,
"\t %s %s %i %i"
,
226
punch
.
chAverageType
[ipPun][i],
227
punch
.
chAverageSpeciesLabel
[ipPun][i] ,
228
punch
.
nAverageIonList
[ipPun][i] ,
229
punch
.
nAverage2ndPar
[ipPun][i] );
230
strcat( chHeader, chTemp );
231
}
232
strcat( chHeader,
"\n"
);
233
}
234
else
if
( strcmp(chJob,
"PUNCH"
) == 0 )
235
{
236
/* do the output */
237
for
( i=0; i<
punch
.
nAverageList
[ipPun] ; ++i )
238
{
239
double
result;
240
char
chWeight[7];
241
if
(
punch
.
nAverage2ndPar
[ipPun][i] == 0 )
242
strcpy( chWeight ,
"RADIUS"
);
243
else
244
strcpy( chWeight ,
"VOLUME"
);
245
246
if
( strncmp(
punch
.
chAverageType
[ipPun][i] ,
"TEMP"
, 4 ) == 0)
247
{
248
/* temperature */
249
if
(
cdTemp
(
250
punch
.
chAverageSpeciesLabel
[ipPun][i] ,
251
punch
.
nAverageIonList
[ipPun][i] ,
252
&result ,
253
chWeight ) )
254
{
255
fprintf(
ioQQQ
,
" punch average temperature could not identify the species.\n"
);
256
cdEXIT
(EXIT_FAILURE);
257
}
258
/* will report log of temperature */
259
result = log10( result );
260
}
261
else
if
( strncmp(
punch
.
chAverageType
[ipPun][i] ,
"IONI"
, 4 ) == 0)
262
{
263
/* ionization fraction
264
* H2 is a special case, HYDRO 0 requests
265
* the H2 fraction, n(H2)/n(H) */
266
if
( strncmp(
"HYDR"
,
267
punch
.
chAverageSpeciesLabel
[ipPun][i] , 4)==0 &&
268
punch
.
nAverageIonList
[ipPun][i]== 0 )
269
strncpy(
punch
.
chAverageSpeciesLabel
[ipPun][i],
270
"H2 "
, 4 );
271
if
(
cdIonFrac
(
272
punch
.
chAverageSpeciesLabel
[ipPun][i] ,
273
punch
.
nAverageIonList
[ipPun][i] ,
274
&result ,
275
chWeight ,
276
false
277
) )
278
{
279
fprintf(
ioQQQ
,
" punch average ionization fraction could not identify the species.\n"
);
280
cdEXIT
(EXIT_FAILURE);
281
}
282
/* will report log of ionization fraction */
283
result = log10( result );
284
}
285
else
if
( strncmp(
punch
.
chAverageType
[ipPun][i] ,
"COLU"
, 4 ) == 0)
286
{
287
/* column density */
288
if
(
cdColm
(
289
punch
.
chAverageSpeciesLabel
[ipPun][i] ,
290
punch
.
nAverageIonList
[ipPun][i] ,
291
&result ) )
292
{
293
fprintf(
ioQQQ
,
" punch average column density fraction could not identify the species.\n"
);
294
cdEXIT
(EXIT_FAILURE);
295
}
296
/* will report log of column density */
297
result = log10( result );
298
}
299
else
300
TotalInsanity
();
301
302
fprintf(
punch
.
ipPnunit
[ipPun],
"\t %.3f"
, result );
303
}
304
fprintf(
punch
.
ipPnunit
[ipPun],
"\n"
);
305
}
306
else
307
{
308
/* total insanity */
309
TotalInsanity
();
310
}
311
/* return */
312
return
;
313
}
Generated for cloudy by
1.8.1.1