00001
00002
00003
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
00029
00030
00031
00032
00033
00034
00035
00036 if( rfield.nspec >= LIMSPC-1 )
00037 {
00038 fprintf( ioQQQ, " Too many spectra entered. Increase LIMSPC\n" );
00039 puts( "[Stop in ParseInterp]" );
00040 cdEXIT(EXIT_FAILURE);
00041 }
00042
00043 if( !rfield.lgContMalloc[rfield.nspec] )
00044 {
00045 rfield.tNuRyd[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) );
00046 rfield.tslop[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) );
00047 rfield.tFluxLog[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) );
00048 rfield.lgContMalloc[rfield.nspec] = true;
00049 }
00050
00051 strcpy( rfield.chSpType[rfield.nspec], "INTER" );
00052
00053
00054 npairs = 0;
00055
00056
00057 lgDONE = false;
00058
00059
00060 *lgEOF = false;
00061 while( !lgDONE && !*lgEOF )
00062 {
00063 i = 5;
00064 lgEOL = false;
00065
00066
00067
00068
00069
00070 while( !lgEOL && npairs<NCELL )
00071 {
00072 rfield.tNuRyd[rfield.nspec][npairs] =
00073 (float)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
00074 rfield.tFluxLog[rfield.nspec][npairs] =
00075 (float)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
00076 ++npairs;
00077 }
00078
00079 if( !lgEOL )
00080 {
00081 fprintf( ioQQQ, "Too many continuum points were entered.\n" );
00082 fprintf( ioQQQ,
00083 "The current logic limits the number of possible points to the value of NCELL, which is %i.\n",NCELL );
00084 fprintf( ioQQQ,
00085 "Increase the value of NCELL in rfield.h.\nSorry.\n" );
00086 puts( "[Stop in ParseInterp]" );
00087 cdEXIT(EXIT_FAILURE);
00088 }
00089
00090
00091 --npairs;
00092
00093
00094 input_readarray(chCard,lgEOF);
00095
00096
00097
00098 while( !*lgEOF && lgInputComment(chCard) )
00099 {
00100 input_readarray(chCard,lgEOF);
00101 }
00102
00103
00104 cap4( chLab4 , chCard );
00105
00106
00107 if( called.lgTalk && (strncmp(chLab4,"CONT",4)==0) )
00108 {
00109 fprintf( ioQQQ, " * ");
00110
00111 k=0;
00112 while( chCard[k]!='\0' )
00113 {
00114 fprintf(ioQQQ,"%c",chCard[k]);
00115 ++k;
00116 }
00117 while( k<80 )
00118 {
00119 fprintf(ioQQQ,"%c",' ');
00120 ++k;
00121 }
00122 fprintf( ioQQQ,"*\n");
00123 }
00124
00125
00126 caps(chCard);
00127
00128
00129 if( strncmp(chCard,"CONT",4) != 0 )
00130 {
00131
00132 lgDONE = true;
00133 }
00134
00135
00136 if( chCard[0] == ' ' )
00137 *lgEOF = true;
00138 }
00139
00140
00141 if( lgDONE )
00142 {
00143 input.nRead -= input.iReadWay;
00144 }
00145
00146
00147 --npairs;
00148 if( npairs < 1 )
00149 {
00150 fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
00151 puts( "[Stop in ParseInterp]" );
00152 cdEXIT(EXIT_FAILURE);
00153 }
00154
00155
00156
00157 else if( npairs >= NCELL - 2 )
00158 {
00159 fprintf( ioQQQ, " Too many continuum points entered.\n" );
00160 fprintf( ioQQQ,
00161 "The current logic limits the number of possible points to the value of NCELL, which is %i.\nSorry.\n",NCELL );
00162 fprintf( ioQQQ, " Increase NCELL (in cddefines.h) to more than the number of continuum points.\n" );
00163 puts( "[Stop in ParseInterp]" );
00164 cdEXIT(EXIT_FAILURE);
00165 }
00166
00167 if( rfield.tNuRyd[rfield.nspec][0] == 0. )
00168 {
00169
00170 if( rfield.tNuRyd[rfield.nspec][1] > 0. )
00171 {
00172
00173 rfield.tNuRyd[rfield.nspec][0] = (float)rfield.emm;
00174 }
00175 else if( rfield.tNuRyd[rfield.nspec][1] < 0. )
00176 {
00177
00178 rfield.tNuRyd[rfield.nspec][0] = (float)log10(rfield.emm);
00179 }
00180 else
00181 {
00182
00183 fprintf( ioQQQ,
00184 "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
00185 rfield.nspec );
00186 puts( "[Stop in ParseInterp]" );
00187 cdEXIT(EXIT_FAILURE);
00188 }
00189 }
00190
00191
00192 if( rfield.tNuRyd[rfield.nspec][0] >= 5. )
00193 {
00194 for( i=0; i < (npairs + 1); i++ )
00195 {
00196 rfield.tNuRyd[rfield.nspec][i] =
00197 (float)pow(10.,rfield.tNuRyd[rfield.nspec][i] - 15.517);
00198 }
00199 }
00200
00201 else if( rfield.tNuRyd[rfield.nspec][0] < 0. )
00202 {
00203
00204 for( i=0; i < (npairs + 1); i++ )
00205 {
00206 rfield.tNuRyd[rfield.nspec][i] = (float)pow(10.,(double)rfield.tNuRyd[rfield.nspec][i]);
00207 }
00208 }
00209
00210 else
00211 {
00212
00213 for( i=0; i < (npairs + 1); i++ )
00214 {
00215 if( rfield.tNuRyd[rfield.nspec][i] == 0. )
00216 {
00217 fprintf( ioQQQ, "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
00218 i );
00219 puts( "[Stop in ParseInterp]" );
00220 cdEXIT(EXIT_FAILURE);
00221 }
00222 }
00223 }
00224
00225
00226 {
00227
00228
00229 enum {DEBUG_LOC=false};
00230
00231 if( DEBUG_LOC )
00232 {
00233 for( i=0; i < npairs; i++ )
00234 {
00235 fprintf(ioQQQ,"%.4e\t%.3e\n",
00236 rfield.tNuRyd[rfield.nspec][i],
00237 pow(10.,(double)rfield.tFluxLog[rfield.nspec][i]) * rfield.tNuRyd[rfield.nspec][i]);
00238 }
00239 cdEXIT(EXIT_SUCCESS);
00240 }
00241 }
00242
00243 for( i=0; i < npairs; i++ )
00244 {
00245
00246 if( rfield.tNuRyd[rfield.nspec][i+1] <= rfield.tNuRyd[rfield.nspec][i] )
00247 {
00248 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",
00249 i, rfield.tNuRyd[rfield.nspec][i] );
00250 puts( "[Stop in ParseInterp]" );
00251 cdEXIT(EXIT_FAILURE);
00252 }
00253
00254
00255 rfield.tslop[rfield.nspec][i] = (float)((rfield.tFluxLog[rfield.nspec][i+1] -
00256 rfield.tFluxLog[rfield.nspec][i])/log10(rfield.tNuRyd[rfield.nspec][i+1]/
00257 rfield.tNuRyd[rfield.nspec][i]));
00258 }
00259
00260
00261
00262
00263
00264 if( npairs + 2 < NCELL )
00265 {
00266
00267
00268
00269 for( i=npairs + 1; i < NCELL; i++ )
00270 {
00271 rfield.tNuRyd[rfield.nspec][i] = 0.;
00272 }
00273 }
00274
00275
00276 if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm )
00277 {
00278
00279 fprintf( ioQQQ,
00280 " Warning: the incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" );
00281 fprintf( ioQQQ,
00282 " >>>>>>>>>>You are making a BIG mistake.\n" );
00283 }
00284
00285
00286 if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm )
00287 {
00288 rfield.lgMMok = false;
00289 }
00290
00291 if( rfield.tNuRyd[rfield.nspec][0] > 1/36. )
00292 {
00293 rfield.lgHPhtOK = false;
00294 }
00295
00296
00297 if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.egamry )
00298 {
00299 rfield.lgGamrOK = false;
00300 }
00301
00302
00303 if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.EnerGammaRay )
00304 {
00305 rfield.lgXRayOK = false;
00306 }
00307
00308
00309 cmax = -38.;
00310 cmin = 28;
00311
00312
00313
00314 for( i=0; i <= npairs; i++ )
00315 {
00316 cmax = MAX2(cmax,rfield.tFluxLog[rfield.nspec][i]);
00317 cmin = MIN2(cmin,rfield.tFluxLog[rfield.nspec][i]);
00318 }
00319
00320
00321 if( cmax - cmin > 74. )
00322 {
00323 fprintf( ioQQQ, "The dynamic range of the specified continuum is too large.\nSorry.\n" );
00324 puts( "[Stop in ParseInterp]" );
00325 cdEXIT(EXIT_FAILURE);
00326 }
00327
00328 if( cmin < -37. )
00329 {
00330 fac = -cmin - 37.;
00331
00332 for( i=0; i <= npairs; i++ )
00333 {
00334 rfield.tFluxLog[rfield.nspec][i] += (float)fac;
00335 }
00336 }
00337
00338 else if( cmax > 37. )
00339 {
00340 fac = cmax - 37.;
00341
00342 for( i=0; i <= npairs; i++ )
00343 {
00344 rfield.tFluxLog[rfield.nspec][i] -= (float)fac;
00345 }
00346 }
00347
00348
00349 if( trace.lgConBug && trace.lgTrace )
00350 {
00351 fprintf( ioQQQ, " Table for this continuum; TNU, TFAC, TSLOP\n" );
00352 for( i=0; i < (npairs + 1); i++ )
00353 {
00354 fprintf( ioQQQ, "%12.4e %12.4e %12.4e\n", rfield.tNuRyd[rfield.nspec][i],
00355 rfield.tFluxLog[rfield.nspec][i], rfield.tslop[rfield.nspec][i] );
00356 }
00357 }
00358
00359
00360
00361
00362 i = 0;
00363 while( rfield.tNuRyd[rfield.nspec][i] <= 1. && i < npairs-1 )
00364 {
00365 i += 1;
00366 }
00367
00368
00369
00370
00371 if( i < 1 )
00372 {
00373 fprintf( ioQQQ, "ParseInput finds insane i after rfield.tNuRyd loop\n");
00374 ShowMe();
00375 puts( "[Stop in ParseInterp]" );
00376 cdEXIT(EXIT_FAILURE);
00377 }
00378
00379
00380 fac = rfield.tFluxLog[rfield.nspec][i-1];
00381 for( i=0; i <= npairs; i++ )
00382 {
00383 rfield.tFluxLog[rfield.nspec][i] -= (float)fac;
00384
00385 }
00386
00387
00388 cmin = log10( FLT_MIN );
00389 cmax = log10( FLT_MAX );
00390 lgHit = false;
00391 for( i=0; i <= npairs; i++ )
00392 {
00393 if( rfield.tFluxLog[rfield.nspec][i] <= cmin ||
00394 rfield.tFluxLog[rfield.nspec][i] >= cmax )
00395 {
00396 fprintf(ioQQQ,
00397 " The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i);
00398 fprintf(ioQQQ,
00399 " The frequency is %f and the log of the flux is %f.\n\n",
00400 rfield.tNuRyd[rfield.nspec][i] ,
00401 rfield.tFluxLog[rfield.nspec][i]);
00402 lgHit = true;
00403 }
00404 }
00405 if( lgHit )
00406 {
00407 fprintf(ioQQQ,"\n The log of the flux given in an interpolate command is outside the range of this cpu.\n");
00408 fprintf(ioQQQ," I will try to renormalize it to be within the range of this cpu, but I crash, this is a likely reason.\n");
00409 fprintf(ioQQQ," Note that the interpolate command only is used for the shape of the continuum.\n");
00410 fprintf(ioQQQ," The order of magnitude of the flux is not used in any way.\n");
00411 fprintf(ioQQQ," For safety this could be of order unity.\n\n");
00412 }
00413
00414
00415 for( i=npairs; i < rfield.nupper; i++ )
00416 {
00417
00418
00419 rfield.tFluxLog[rfield.nspec][i] = 0.;
00420 rfield.tslop[rfield.nspec][i] = 0.;
00421 rfield.tNuRyd[rfield.nspec][i] = 0.;;
00422 }
00423
00424
00425 ++rfield.nspec;
00426 if( rfield.nspec >= LIMSPC )
00427 {
00428
00429 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
00430 puts( "[Stop in ParseCoronal]" );
00431 cdEXIT(EXIT_FAILURE);
00432 }
00433
00434 DEBUG_EXIT( "ParseInterp()" );
00435
00436 return;
00437 }