cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 */
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 doxygen 1.8.1.2