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 /*ParseInterp parse parameters on interpolate command */ 00004 #include "cddefines.h" 00005 #include "called.h" 00006 #include "rfield.h" 00007 #include "trace.h" 00008 #include "input.h" 00009 #include "parse.h" 00010 00011 void ParseInterp(char *chCard, 00012 bool *lgEOF) 00013 { 00014 char chLab4[5]; 00015 bool lgDONE, 00016 lgEOL, 00017 lgHit; 00018 long int i, 00019 k, 00020 npairs; 00021 double cmax, 00022 cmin, 00023 fac; 00024 00025 DEBUG_ENTRY( "ParseInterp()" ); 00026 00027 /* 00028 * this sub reads in the "interpolate" command, and has 00029 * logic for the "continue" lines as well. 00030 * OUTPUT: 00031 * TNU is vector of energies where the grid is defined 00032 * TSLOP initially is vector of log fnu at each freq 00033 * converted into slopes here too 00034 */ 00035 00036 if( rfield.nspec >= LIMSPC-1 ) 00037 { 00038 fprintf( ioQQQ, " Too many spectra entered. Increase LIMSPC\n" ); 00039 cdEXIT(EXIT_FAILURE); 00040 } 00041 /* >>chng 06 jul 16, fix memory leak when optimizing, PvH */ 00042 if( !rfield.lgContMalloc[rfield.nspec] ) 00043 { 00044 rfield.tNuRyd[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) ); 00045 rfield.tslop[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) ); 00046 rfield.tFluxLog[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) ); 00047 rfield.lgContMalloc[rfield.nspec] = true; 00048 } 00049 00050 strcpy( rfield.chSpType[rfield.nspec], "INTER" ); 00051 00052 /* read all of the numbers on a line */ 00053 npairs = 0; 00054 00055 /* this is flag saying that all numbers are in */ 00056 lgDONE = false; 00057 00058 /* this flag says we hit end of command stream */ 00059 *lgEOF = false; 00060 while( !lgDONE && !*lgEOF ) 00061 { 00062 i = 5; 00063 lgEOL = false; 00064 00065 /* keep scanning numbers until we hit eol for current line image */ 00066 /* >>chng 01 aug 10, add check on npairs not exceeding NC_ELL as per 00067 * Ilse van Bemmel bug report */ 00068 /*while( !lgEOL && npairs<rfield.nupper )*/ 00069 while( !lgEOL && npairs<NCELL ) 00070 { 00071 rfield.tNuRyd[rfield.nspec][npairs] = 00072 (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL); 00073 rfield.tFluxLog[rfield.nspec][npairs] = 00074 (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL); 00075 ++npairs; 00076 } 00077 /* check if we fell through due to too many points, did not hit EOL */ 00078 if( !lgEOL ) 00079 { 00080 fprintf( ioQQQ, "Too many continuum points were entered.\n" ); 00081 fprintf( ioQQQ, 00082 "The current logic limits the number of possible points to the value of NCELL, which is %i.\n",NCELL ); 00083 fprintf( ioQQQ, 00084 "Increase the value of NCELL in rfield.h.\nSorry.\n" ); 00085 cdEXIT(EXIT_FAILURE); 00086 } 00087 00088 /* added too many above, so back off */ 00089 --npairs; 00090 00091 /* read a new line, checking for EOF */ 00092 input_readarray(chCard,lgEOF); 00093 00094 /* option to ignore all lines starting with #, *, or %. 00095 * >>chng 06 sep 04 use routine to check for comments */ 00096 while( !*lgEOF && lgInputComment(chCard)/*((chCard[0] == '#' || chCard[0] == '*') || chCard[0] == '%')*/ ) 00097 { 00098 input_readarray(chCard,lgEOF); 00099 } 00100 00101 /* get capital first four char of chCard */ 00102 cap4( chLab4 , chCard ); 00103 00104 /* print the line, but only if it is a continue line */ 00105 if( called.lgTalk && (strncmp(chLab4,"CONT",4)==0) ) 00106 { 00107 fprintf( ioQQQ, " * "); 00108 /* next logic will print command plus spaces to right justified * */ 00109 k=0; 00110 while( chCard[k]!='\0' ) 00111 { 00112 fprintf(ioQQQ,"%c",chCard[k]); 00113 ++k; 00114 } 00115 while( k<80 ) 00116 { 00117 fprintf(ioQQQ,"%c",' '); 00118 ++k; 00119 } 00120 fprintf( ioQQQ,"*\n"); 00121 } 00122 00123 /* now convert entire line to caps */ 00124 caps(chCard); 00125 00126 /* is this a continue line? */ 00127 if( strncmp(chCard,"CONT",4) != 0 ) 00128 { 00129 /* we have a line but next command, not continue */ 00130 lgDONE = true; 00131 } 00132 00133 /* this is another way to hit end of input stream - blank lines */ 00134 if( chCard[0] == ' ' ) 00135 *lgEOF = true; 00136 } 00137 00138 /* if valid next line, backup one line */ 00139 if( lgDONE ) 00140 { 00141 input.nRead -= input.iReadWay; 00142 } 00143 00144 /* done reading all of the possible lines */ 00145 --npairs; 00146 if( npairs < 1 ) 00147 { 00148 fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" ); 00149 cdEXIT(EXIT_FAILURE); 00150 } 00151 00152 /* >>chng 02 apr 19, from > to >=, test did not trigger due to overrun protection 00153 * in main loop */ 00154 else if( npairs >= NCELL - 2 ) 00155 { 00156 fprintf( ioQQQ, " Too many continuum points entered.\n" ); 00157 fprintf( ioQQQ, 00158 "The current logic limits the number of possible points to the value of NCELL, which is %i.\nSorry.\n",NCELL ); 00159 fprintf( ioQQQ, " Increase NCELL (in cddefines.h) to more than the number of continuum points.\n" ); 00160 cdEXIT(EXIT_FAILURE); 00161 } 00162 00163 if( rfield.tNuRyd[rfield.nspec][0] == 0. ) 00164 { 00165 /* special case - if first energy is zero then it is low energy */ 00166 if( rfield.tNuRyd[rfield.nspec][1] > 0. ) 00167 { 00168 /* second energy positive, assume linear Ryd */ 00169 rfield.tNuRyd[rfield.nspec][0] = (realnum)rfield.emm; 00170 } 00171 else if( rfield.tNuRyd[rfield.nspec][1] < 0. ) 00172 { 00173 /* second energy negative, assume log of Ryd */ 00174 rfield.tNuRyd[rfield.nspec][0] = (realnum)log10(rfield.emm); 00175 } 00176 else 00177 { 00178 /* second energy zero, not allowed */ 00179 fprintf( ioQQQ, 00180 "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n", 00181 rfield.nspec ); 00182 cdEXIT(EXIT_FAILURE); 00183 } 00184 } 00185 00186 /* convert from log(Hz) to Ryd if first nu>5 */ 00187 if( rfield.tNuRyd[rfield.nspec][0] >= 5. ) 00188 { 00189 for( i=0; i < (npairs + 1); i++ ) 00190 { 00191 rfield.tNuRyd[rfield.nspec][i] = 00192 (realnum)pow(10.,rfield.tNuRyd[rfield.nspec][i] - 15.517); 00193 } 00194 } 00195 00196 else if( rfield.tNuRyd[rfield.nspec][0] < 0. ) 00197 { 00198 /* energies entered as logs */ 00199 for( i=0; i < (npairs + 1); i++ ) 00200 { 00201 rfield.tNuRyd[rfield.nspec][i] = (realnum)pow(10.,(double)rfield.tNuRyd[rfield.nspec][i]); 00202 } 00203 } 00204 00205 else 00206 { 00207 /* numbers are linear Rydbergs */ 00208 for( i=0; i < (npairs + 1); i++ ) 00209 { 00210 if( rfield.tNuRyd[rfield.nspec][i] == 0. ) 00211 { 00212 fprintf( ioQQQ, "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n", 00213 i ); 00214 cdEXIT(EXIT_FAILURE); 00215 } 00216 } 00217 } 00218 00219 /* option to print debugging file then exit */ 00220 { 00221 /* following should be set true to print information */ 00222 enum {DEBUG_LOC=false}; 00223 if( DEBUG_LOC ) 00224 { 00225 for( i=0; i < npairs; i++ ) 00226 { 00227 fprintf(ioQQQ,"%.4e\t%.3e\n", 00228 rfield.tNuRyd[rfield.nspec][i], 00229 pow(10.,(double)rfield.tFluxLog[rfield.nspec][i]) * rfield.tNuRyd[rfield.nspec][i]); 00230 } 00231 cdEXIT(EXIT_SUCCESS); 00232 } 00233 } 00234 00235 for( i=0; i < npairs; i++ ) 00236 { 00237 /* check that frequencies are monotonically increasing */ 00238 if( rfield.tNuRyd[rfield.nspec][i+1] <= rfield.tNuRyd[rfield.nspec][i] ) 00239 { 00240 fprintf( ioQQQ, "The energies MUST be in increasing order. Energy #%3ld=%10.2e Ryd was greater than or equal to the next one.\nSorry.\n", 00241 i, rfield.tNuRyd[rfield.nspec][i] ); 00242 cdEXIT(EXIT_FAILURE); 00243 } 00244 00245 /* TFAC is energy, and TSLOP is slope in f_nu; not photons */ 00246 rfield.tslop[rfield.nspec][i] = (realnum)((rfield.tFluxLog[rfield.nspec][i+1] - 00247 rfield.tFluxLog[rfield.nspec][i])/log10(rfield.tNuRyd[rfield.nspec][i+1]/ 00248 rfield.tNuRyd[rfield.nspec][i])); 00249 } 00250 00251 /*>>chng 06 may 17, must insure that rNuRyd is sane through NCELL values, not just 00252 * nupper values. This bit when stellar continuum had more cells than cloudy grid, 00253 * so npairs is much greater than nupper */ 00254 /*if( npairs + 2 < rfield.nupper )*/ 00255 if( npairs + 2 < NCELL ) 00256 { 00257 /* zero out remainder of array */ 00258 /*>>chng 06 may 17, same as above */ 00259 /*for( i=npairs + 1; i < rfield.nupper; i++ )*/ 00260 for( i=npairs + 1; i < NCELL; i++ ) 00261 { 00262 rfield.tNuRyd[rfield.nspec][i] = 0.; 00263 } 00264 } 00265 00266 /* now check that array is defined over all energies */ 00267 if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm ) 00268 { 00269 /* not defined over low energy part of array */ 00270 fprintf( ioQQQ, 00271 "\n NOTE The incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" ); 00272 fprintf( ioQQQ, 00273 " NOTE You may be making a BIG mistake.\n\n" ); 00274 } 00275 00276 /* check on IR */ 00277 if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm ) 00278 rfield.lgMMok = false; 00279 00280 if( rfield.tNuRyd[rfield.nspec][0] > 1/36. ) 00281 rfield.lgHPhtOK = false; 00282 00283 /* gamma ray, EGAMRY is roughly 100MeV */ 00284 if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.egamry ) 00285 rfield.lgGamrOK = false; 00286 00287 /* EnerGammaRay is roughly 100keV; high is gamma ray */ 00288 if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.EnerGammaRay ) 00289 rfield.lgXRayOK = false; 00290 00291 /* find min and max of continuum */ 00292 cmax = -38.; 00293 cmin = 28; 00294 00295 /* rfield.tFluxLog can be very large or small */ 00296 /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */ 00297 for( i=0; i <= npairs; i++ ) 00298 { 00299 cmax = MAX2(cmax,rfield.tFluxLog[rfield.nspec][i]); 00300 cmin = MIN2(cmin,rfield.tFluxLog[rfield.nspec][i]); 00301 } 00302 00303 /* check on dynamic range of input continuum */ 00304 if( cmax - cmin > 74. ) 00305 { 00306 fprintf( ioQQQ, "The dynamic range of the specified continuum is too large.\nSorry.\n" ); 00307 cdEXIT(EXIT_FAILURE); 00308 } 00309 00310 if( cmin < -37. ) 00311 { 00312 fac = -cmin - 37.; 00313 /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */ 00314 for( i=0; i <= npairs; i++ ) 00315 { 00316 rfield.tFluxLog[rfield.nspec][i] += (realnum)fac; 00317 } 00318 } 00319 00320 else if( cmax > 37. ) 00321 { 00322 fac = cmax - 37.; 00323 /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */ 00324 for( i=0; i <= npairs; i++ ) 00325 { 00326 rfield.tFluxLog[rfield.nspec][i] -= (realnum)fac; 00327 } 00328 } 00329 00330 /* option to print out results at this stage - "trace continuum" */ 00331 if( trace.lgConBug && trace.lgTrace ) 00332 { 00333 fprintf( ioQQQ, " Table for this continuum;\ni\tTNU\tTFAC\tTSLOP, npairs=%li\n", 00334 npairs ); 00335 for( i=0; i < npairs; i++ ) 00336 { 00337 fprintf( ioQQQ, "%li\t%.4e\t%.4e\t%.4e\n", 00338 i , rfield.tNuRyd[rfield.nspec][i], 00339 rfield.tFluxLog[rfield.nspec][i], rfield.tslop[rfield.nspec][i] ); 00340 } 00341 i = npairs; 00342 fprintf( ioQQQ, "%li\t%.4e\t%.4e\n", 00343 i , rfield.tNuRyd[rfield.nspec][i], 00344 rfield.tFluxLog[rfield.nspec][i]); 00345 } 00346 00347 /* renormalize continuum so that flux we will interpolate upon passes through unity 00348 * at near 1 Ryd. but first we must find 1 Ryd in the array. 00349 * find 1 Ryd, npairs is one less than number of continuum pairs */ 00350 i = 0; 00351 while( rfield.tNuRyd[rfield.nspec][i] <= 1. && i < npairs-1 ) 00352 { 00353 i += 1; 00354 } 00355 /* i is now the table point where rfield.tNuRyd[i-1] is <= 1 ryd, 00356 * and rfield.tNuRyd[i] > 1 ryd */ 00357 00358 /* following is impossible but sanity check */ 00359 if( i < 1 ) 00360 { 00361 fprintf( ioQQQ, "ParseInput finds insane i after rfield.tNuRyd loop\n"); 00362 ShowMe(); 00363 cdEXIT(EXIT_FAILURE); 00364 } 00365 00366 /* do the renormalization so 1 ryd is about unity */ 00367 fac = rfield.tFluxLog[rfield.nspec][i-1]; 00368 for( i=0; i <= npairs; i++ ) 00369 { 00370 rfield.tFluxLog[rfield.nspec][i] -= (realnum)fac; 00371 /*fprintf(ioQQQ,"DEBUG parse %li %e \n", i , rfield.tFluxLog[rfield.nspec][i] );*/ 00372 } 00373 00374 /* finally check that we are within dynamic range of this cpu */ 00375 cmin = log10( FLT_MIN ); 00376 cmax = log10( FLT_MAX ); 00377 lgHit = false; 00378 for( i=0; i <= npairs; i++ ) 00379 { 00380 if( rfield.tFluxLog[rfield.nspec][i] <= cmin || 00381 rfield.tFluxLog[rfield.nspec][i] >= cmax ) 00382 { 00383 fprintf(ioQQQ, 00384 " The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i); 00385 fprintf(ioQQQ, 00386 " The frequency is %f and the log of the flux is %f.\n\n", 00387 rfield.tNuRyd[rfield.nspec][i] , 00388 rfield.tFluxLog[rfield.nspec][i]); 00389 lgHit = true; 00390 } 00391 } 00392 if( lgHit ) 00393 { 00394 fprintf(ioQQQ,"\n NOTE The log of the flux given in an interpolate command is outside the range of this cpu.\n"); 00395 fprintf(ioQQQ," NOTE I will try to renormalize it to be within the range of this cpu, but if I crash, this is a likely reason.\n"); 00396 fprintf(ioQQQ," NOTE Note that the interpolate command only is used for the shape of the continuum.\n"); 00397 fprintf(ioQQQ," NOTE The order of magnitude of the flux is not used in any way.\n"); 00398 fprintf(ioQQQ," NOTE For safety this could be of order unity.\n\n"); 00399 } 00400 00401 /* zero out remainder of array */ 00402 for( i=npairs+1; i < rfield.nupper; i++ ) 00403 { 00404 /* >>chng 05 jul 27, next three indices had been ipspec, changed to nspec 00405 * bug caught by Ralf Quast */ 00406 rfield.tFluxLog[rfield.nspec][i] = 0.; 00407 rfield.tslop[rfield.nspec][i] = 0.; 00408 rfield.tNuRyd[rfield.nspec][i] = 0.; 00409 } 00410 00411 /* now increment number of continua */ 00412 ++rfield.nspec; 00413 if( rfield.nspec >= LIMSPC ) 00414 { 00415 /* too many continua were entered */ 00416 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00417 cdEXIT(EXIT_FAILURE); 00418 } 00419 00420 return; 00421 }