cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*punch_average routine to read in list of species to output as averages */ 00004 #include "cddefines.h" 00005 #include "cddrive.h" 00006 #include "input.h" 00007 #include "elementnames.h" 00008 #include "punch.h" 00009 00010 /* return value is number of lines, -1 if file could not be opened */ 00011 void punch_average( 00012 /* the file we will write to */ 00013 long int ipPun, 00014 /* this is the job we shall do - READ and PUNCH */ 00015 const char *chJob, 00016 char *chHeader) 00017 { 00018 long int i; 00019 bool lgEOF , lgEOL; 00020 long int nLine; 00021 00022 char chLine[FILENAME_PATH_LENGTH_2] ; 00023 00024 DEBUG_ENTRY( "punch_average()" ); 00025 00026 if( strcmp(chJob,"READ") == 0 ) 00027 { 00028 char chCap[INPUT_LINE_LENGTH]; 00029 char chTemp[INPUT_LINE_LENGTH]; 00030 00031 /* this is index for first line we will read. use this to count 00032 * total number of species */ 00033 nLine = input.nRead; 00034 /* very first time this routine is called, chJob is "READ" and we read 00035 * in lines from the input stream. The species labels and other information 00036 * are store in the punch structure. These are output in a later call 00037 * to this routine with argument PUNCH */ 00038 00039 /* number of lines we will save */ 00040 punch.nAverageList[ipPun] = 0; 00041 00042 /* get the next line, and check for eof */ 00043 input_readarray(chLine,&lgEOF); 00044 if( lgEOF ) 00045 { 00046 fprintf( ioQQQ, 00047 " Punch average hit EOF while reading list; use END to end list.\n" ); 00048 cdEXIT(EXIT_FAILURE); 00049 } 00050 00051 /* convert line to caps */ 00052 strcpy( chCap, chLine ); 00053 caps(chCap); 00054 00055 /* keep reading until we hit END */ 00056 while( strncmp(chCap, "END" ,3 ) != 0 ) 00057 { 00058 /* count number of species we will save */ 00059 ++punch.nAverageList[ipPun]; 00060 00061 /* get next line and check for eof */ 00062 input_readarray(chLine,&lgEOF); 00063 if( lgEOF ) 00064 { 00065 fprintf( ioQQQ, " punch averages hit EOF while reading species list; use END to end list.\n" ); 00066 cdEXIT(EXIT_FAILURE); 00067 } 00068 00069 /* convert it to caps */ 00070 strcpy( chCap, chLine ); 00071 caps(chCap); 00072 } 00073 /*# define PADEBUG*/ 00074 # ifdef PADEBUG 00075 fprintf(ioQQQ , "DEBUG punch_average %li species read in.\n", 00076 punch.nAverageList[ipPun] ); 00077 # endif 00078 00079 /* now create space that will be needed to hold these arrays */ 00080 00081 if( punch.ipPnunit[ipPun] == NULL ) 00082 { 00083 /* nAverageIonList is set of ions for averages */ 00084 punch.nAverageIonList[ipPun] = (int *)MALLOC( 00085 (size_t)(punch.nAverageList[ipPun]*sizeof(int)) ); 00086 00087 /* nAverage2ndPar is set of second parameters for averages */ 00088 punch.nAverage2ndPar[ipPun] = (int *)MALLOC((size_t) 00089 (punch.nAverageList[ipPun]*sizeof(int)) ); 00090 00091 /* chAverageType is label for type of average */ 00092 punch.chAverageType[ipPun] = (char **)MALLOC((size_t) 00093 (punch.nAverageList[ipPun]*sizeof(char *)) ); 00094 00095 /* chAverageSpeciesLabel is label for species */ 00096 punch.chAverageSpeciesLabel[ipPun] = (char **)MALLOC((size_t) 00097 (punch.nAverageList[ipPun]*sizeof(char *)) ); 00098 for( i=0; i<punch.nAverageList[ipPun]; ++i ) 00099 { 00100 /* create space for labels themselves */ 00101 punch.chAverageType[ipPun][i] = (char *)MALLOC((size_t) 00102 (5*sizeof(char)) ); 00103 00104 /* chAverageSpeciesLabel is label for species */ 00105 punch.chAverageSpeciesLabel[ipPun][i] = (char *)MALLOC((size_t) 00106 (5*sizeof(char)) ); 00107 } 00108 } 00109 00110 /* reset array input read to first of the species lines */ 00111 input.nRead = nLine; 00112 00113 # ifdef PADEBUG 00114 fprintf(ioQQQ , "DEBUG punch_average %li species read in.\n", 00115 punch.nAverageList[ipPun] ); 00116 # endif 00117 00118 /* reread the input lines and save the data */ 00119 input_readarray(chLine,&lgEOF); 00120 if( lgEOF ) 00121 { 00122 fprintf( ioQQQ, 00123 " Punch average hit EOF while reading list; use END to end list.\n" ); 00124 cdEXIT(EXIT_FAILURE); 00125 } 00126 00127 /* convert line to caps */ 00128 strcpy( chCap, chLine ); 00129 caps(chCap); 00130 00131 /* use this to count number of species, and will assert equal to 00132 * total malloced above */ 00133 nLine = 0; 00134 /* keep reading until we hit END */ 00135 while( strncmp(chCap, "END" ,3 ) != 0 ) 00136 { 00137 /* count number of species we will save */ 00138 ++nLine; 00139 if( nMatch("TEMP" , chCap )) 00140 { 00141 /* temperature */ 00142 strcpy( punch.chAverageType[ipPun][nLine-1] , "TEMP" ); 00143 } 00144 else if( nMatch("COLU" , chCap )) 00145 { 00146 /* column density */ 00147 strcpy( punch.chAverageType[ipPun][nLine-1] , "COLU" ); 00148 } 00149 else if( nMatch("IONI" , chCap )) 00150 { 00151 /* ionization fraction */ 00152 strcpy( punch.chAverageType[ipPun][nLine-1] , "IONI" ); 00153 } 00154 else 00155 { 00156 fprintf(ioQQQ,"PROBLEM one of the jobs TEMPerature, COLUmn density, or IONIzation, must appear.\n"); 00157 cdEXIT(EXIT_FAILURE); 00158 } 00159 00160 /* get element name, a string we shall pass on to the routine 00161 * that computes final quantities */ 00162 if( (i = GetElem( chCap )) < 0 ) 00163 { 00164 /* no name found */ 00165 fprintf(ioQQQ, "punch average did not an element on this line, sorry %s\n", 00166 chCap ); 00167 cdEXIT(EXIT_FAILURE); 00168 } 00169 strcpy( punch.chAverageSpeciesLabel[ipPun][nLine-1] , elementnames.chElementNameShort[i]); 00170 00171 /* now get ionization stage */ 00172 i = 5; 00173 punch.nAverageIonList[ipPun][nLine-1] = (int) 00174 FFmtRead( chCap , &i,INPUT_LINE_LENGTH,&lgEOL ); 00175 if( lgEOF ) 00176 { 00177 /* error - needed that ionization stage */ 00178 NoNumb( chCap ); 00179 } 00180 00181 /* look for volume keyword, otherwise will be radius 00182 * only used for some options */ 00183 if( nMatch( "VOLU" , chCap ) ) 00184 { 00185 /* volume */ 00186 punch.nAverage2ndPar[ipPun][nLine-1] = 1; 00187 } 00188 else 00189 { 00190 /* radius */ 00191 punch.nAverage2ndPar[ipPun][nLine-1] = 0; 00192 } 00193 00194 /* get next line and check for eof */ 00195 input_readarray(chLine,&lgEOF); 00196 if( lgEOF ) 00197 { 00198 fprintf( ioQQQ, " punch averages hit EOF while reading species list; use END to end list.\n" ); 00199 cdEXIT(EXIT_FAILURE); 00200 } 00201 00202 /* convert it to caps */ 00203 strcpy( chCap, chLine ); 00204 caps(chCap); 00205 } 00206 00207 /* these must equal or we have a major logical error */ 00208 ASSERT( nLine == punch.nAverageList[ipPun]); 00209 00210 # ifdef PADEBUG 00211 for( i=0; i<nLine ; ++i ) 00212 { 00213 fprintf(ioQQQ, "PDDEBUG %s %s %i %i\n", 00214 punch.chAverageType[ipPun][i], 00215 punch.chAverageSpeciesLabel[ipPun][i] , 00216 punch.nAverageIonList[ipPun][i] , 00217 punch.nAverage2ndPar[ipPun][i] ); 00218 } 00219 # endif 00220 00221 /* punch headers */ 00222 sprintf(chHeader, "#averages"); 00223 for( i=0; i<nLine ; ++i ) 00224 { 00225 sprintf(chTemp, "\t %s %s %i %i", 00226 punch.chAverageType[ipPun][i], 00227 punch.chAverageSpeciesLabel[ipPun][i] , 00228 punch.nAverageIonList[ipPun][i] , 00229 punch.nAverage2ndPar[ipPun][i] ); 00230 strcat( chHeader, chTemp ); 00231 } 00232 strcat( chHeader, "\n"); 00233 } 00234 else if( strcmp(chJob,"PUNCH") == 0 ) 00235 { 00236 /* do the output */ 00237 for( i=0; i<punch.nAverageList[ipPun] ; ++i ) 00238 { 00239 double result; 00240 char chWeight[7]; 00241 if( punch.nAverage2ndPar[ipPun][i] == 0 ) 00242 strcpy( chWeight , "RADIUS"); 00243 else 00244 strcpy( chWeight , "VOLUME"); 00245 00246 if( strncmp( punch.chAverageType[ipPun][i] , "TEMP" , 4 ) == 0) 00247 { 00248 /* temperature */ 00249 if( cdTemp( 00250 punch.chAverageSpeciesLabel[ipPun][i] , 00251 punch.nAverageIonList[ipPun][i] , 00252 &result , 00253 chWeight ) ) 00254 { 00255 fprintf( ioQQQ, " punch average temperature could not identify the species.\n" ); 00256 cdEXIT(EXIT_FAILURE); 00257 } 00258 /* will report log of temperature */ 00259 result = log10( result ); 00260 } 00261 else if( strncmp( punch.chAverageType[ipPun][i] , "IONI" , 4 ) == 0) 00262 { 00263 /* ionization fraction 00264 * H2 is a special case, HYDRO 0 requests 00265 * the H2 fraction, n(H2)/n(H) */ 00266 if( strncmp( "HYDR" , 00267 punch.chAverageSpeciesLabel[ipPun][i] , 4)==0 && 00268 punch.nAverageIonList[ipPun][i]== 0 ) 00269 strncpy( punch.chAverageSpeciesLabel[ipPun][i], 00270 "H2 " , 4 ); 00271 if( cdIonFrac( 00272 punch.chAverageSpeciesLabel[ipPun][i] , 00273 punch.nAverageIonList[ipPun][i] , 00274 &result , 00275 chWeight , 00276 false 00277 ) ) 00278 { 00279 fprintf( ioQQQ, " punch average ionization fraction could not identify the species.\n" ); 00280 cdEXIT(EXIT_FAILURE); 00281 } 00282 /* will report log of ionization fraction */ 00283 result = log10( result ); 00284 } 00285 else if( strncmp( punch.chAverageType[ipPun][i] , "COLU" , 4 ) == 0) 00286 { 00287 /* column density */ 00288 if( cdColm( 00289 punch.chAverageSpeciesLabel[ipPun][i] , 00290 punch.nAverageIonList[ipPun][i] , 00291 &result ) ) 00292 { 00293 fprintf( ioQQQ, " punch average column density fraction could not identify the species.\n" ); 00294 cdEXIT(EXIT_FAILURE); 00295 } 00296 /* will report log of column density */ 00297 result = log10( result ); 00298 } 00299 else 00300 TotalInsanity(); 00301 00302 fprintf(punch.ipPnunit[ipPun], "\t %.3f", result ); 00303 } 00304 fprintf(punch.ipPnunit[ipPun], "\n"); 00305 } 00306 else 00307 { 00308 /* total insanity */ 00309 TotalInsanity(); 00310 } 00311 /* return */ 00312 return; 00313 }