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 "parser.h"
00010
00011 void ParseInterp(Parser &p)
00012 {
00013 bool lgDONE,
00014 lgHit;
00015 long int npairs;
00016 double cmax,
00017 cmin,
00018 fac;
00019
00020 DEBUG_ENTRY( "ParseInterp()" );
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 if( rfield.nShape >= LIMSPC-1 )
00032 {
00033 fprintf( ioQQQ, " Too many spectra entered. Increase LIMSPC\n" );
00034 cdEXIT(EXIT_FAILURE);
00035 }
00036
00037 if( !rfield.lgContMalloc[rfield.nShape] )
00038 {
00039 rfield.tNu[rfield.nShape].resize( NCELL );
00040 rfield.tslop[rfield.nShape].resize( NCELL );
00041 rfield.tFluxLog[rfield.nShape].resize( NCELL );
00042 rfield.lgContMalloc[rfield.nShape] = true;
00043 }
00044
00045 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
00046
00047
00048 npairs = 0;
00049
00050
00051 lgDONE = false;
00052
00053
00054 p.m_lgEOF = false;
00055 while( !lgDONE && !p.m_lgEOF )
00056 {
00057
00058
00059
00060
00061 while( !p.lgEOL() && npairs<NCELL )
00062 {
00063 rfield.tNu[rfield.nShape][npairs].set(
00064 p.FFmtRead());
00065 rfield.tFluxLog[rfield.nShape][npairs] =
00066 (realnum)p.FFmtRead();
00067 ++npairs;
00068 }
00069
00070 if( !p.lgEOL() )
00071 {
00072 fprintf( ioQQQ, "Too many continuum points were entered.\n" );
00073 fprintf( ioQQQ,
00074 "The current logic limits the number of possible points to the value of NCELL, which is %i.\n",NCELL );
00075 fprintf( ioQQQ,
00076 "Increase the value of NCELL in rfield.h.\nSorry.\n" );
00077 cdEXIT(EXIT_FAILURE);
00078 }
00079
00080
00081 --npairs;
00082
00083 do {
00084
00085 p.getline();
00086
00087
00088
00089 } while( !p.m_lgEOF && p.isComment());
00090
00091
00092 if( called.lgTalk && (p.strcmp("CONT")==0) )
00093 {
00094 p.echo();
00095 }
00096
00097
00098 if( p.strcmp("CONT") != 0 )
00099 {
00100
00101 lgDONE = true;
00102 }
00103
00104
00105 if( p.last() )
00106 p.m_lgEOF = true;
00107 }
00108
00109
00110 if( lgDONE )
00111 {
00112 input.nRead -= input.iReadWay;
00113 }
00114
00115
00116 --npairs;
00117 if( npairs < 1 )
00118 {
00119 fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
00120 cdEXIT(EXIT_FAILURE);
00121 }
00122
00123
00124
00125 else if( npairs >= NCELL - 2 )
00126 {
00127 fprintf( ioQQQ, " Too many continuum points entered.\n" );
00128 fprintf( ioQQQ,
00129 "The current logic limits the number of possible points to the value of NCELL, which is %i.\nSorry.\n",NCELL );
00130 fprintf( ioQQQ, " Increase NCELL (in cddefines.h) to more than the number of continuum points.\n" );
00131 cdEXIT(EXIT_FAILURE);
00132 }
00133
00134 if( rfield.tNu[rfield.nShape][0].Ryd() == 0. )
00135 {
00136
00137 if( rfield.tNu[rfield.nShape][1].Ryd() > 0. )
00138 {
00139
00140 rfield.tNu[rfield.nShape][0].set(rfield.emm);
00141 }
00142 else if( rfield.tNu[rfield.nShape][1].Ryd() < 0. )
00143 {
00144
00145 rfield.tNu[rfield.nShape][0].set(log10(rfield.emm));
00146 }
00147 else
00148 {
00149
00150 fprintf( ioQQQ,
00151 "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
00152 rfield.nShape );
00153 cdEXIT(EXIT_FAILURE);
00154 }
00155 }
00156
00157
00158 if( rfield.tNu[rfield.nShape][0].Ryd() >= 5. )
00159 {
00160 for( long i=0; i < (npairs + 1); i++ )
00161 {
00162 rfield.tNu[rfield.nShape][i].set(
00163 pow(10.,rfield.tNu[rfield.nShape][i].Ryd())/FR1RYD);
00164 }
00165 }
00166
00167 else if( rfield.tNu[rfield.nShape][0].Ryd() < 0. )
00168 {
00169
00170 for( long i=0; i < (npairs + 1); i++ )
00171 {
00172 rfield.tNu[rfield.nShape][i].set(
00173 pow(10.,(double)rfield.tNu[rfield.nShape][i].Ryd()));
00174 }
00175 }
00176
00177 else
00178 {
00179
00180 for( long i=0; i < (npairs + 1); i++ )
00181 {
00182 if( rfield.tNu[rfield.nShape][i].Ryd() == 0. )
00183 {
00184 fprintf( ioQQQ, "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
00185 i );
00186 cdEXIT(EXIT_FAILURE);
00187 }
00188 }
00189 }
00190
00191
00192 {
00193
00194 enum {DEBUG_LOC=false};
00195 if( DEBUG_LOC )
00196 {
00197 for( long i=0; i < npairs; i++ )
00198 {
00199 fprintf(ioQQQ,"%.4e\t%.3e\n",
00200 rfield.tNu[rfield.nShape][i].Ryd(),
00201 pow(10.,(double)rfield.tFluxLog[rfield.nShape][i]) * rfield.tNu[rfield.nShape][i].Ryd());
00202 }
00203 cdEXIT(EXIT_SUCCESS);
00204 }
00205 }
00206
00207 for( long i=0; i < npairs; i++ )
00208 {
00209
00210 if( rfield.tNu[rfield.nShape][i+1].Ryd() <= rfield.tNu[rfield.nShape][i].Ryd() )
00211 {
00212 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",
00213 i, rfield.tNu[rfield.nShape][i].Ryd() );
00214 cdEXIT(EXIT_FAILURE);
00215 }
00216
00217
00218 rfield.tslop[rfield.nShape][i] = (realnum)((rfield.tFluxLog[rfield.nShape][i+1] -
00219 rfield.tFluxLog[rfield.nShape][i])/log10(rfield.tNu[rfield.nShape][i+1].Ryd()/
00220 rfield.tNu[rfield.nShape][i].Ryd()));
00221 }
00222
00223
00224
00225
00226
00227 if( npairs + 2 < NCELL )
00228 {
00229
00230
00231
00232 for( long i=npairs + 1; i < NCELL; i++ )
00233 {
00234 rfield.tNu[rfield.nShape][i].set(0.);
00235 }
00236 }
00237
00238
00239 if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm )
00240 {
00241
00242 fprintf( ioQQQ,
00243 "\n NOTE The incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" );
00244 fprintf( ioQQQ,
00245 " NOTE You may be making a BIG mistake.\n\n" );
00246 }
00247
00248
00249 if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm )
00250 rfield.lgMMok = false;
00251
00252 if( rfield.tNu[rfield.nShape][0].Ryd() > 1/36. )
00253 rfield.lgHPhtOK = false;
00254
00255
00256 if( rfield.tNu[rfield.nShape][npairs].Ryd() < rfield.egamry )
00257 rfield.lgGamrOK = false;
00258
00259
00260 if( rfield.tNu[rfield.nShape][npairs].Ryd() < rfield.EnerGammaRay )
00261 rfield.lgXRayOK = false;
00262
00263
00264 cmax = -38.;
00265 cmin = 28;
00266
00267
00268
00269 for( long i=0; i <= npairs; i++ )
00270 {
00271 cmax = MAX2(cmax,rfield.tFluxLog[rfield.nShape][i]);
00272 cmin = MIN2(cmin,rfield.tFluxLog[rfield.nShape][i]);
00273 }
00274
00275
00276 if( cmax - cmin > 74. )
00277 {
00278 fprintf( ioQQQ, "The dynamic range of the specified continuum is too large.\nSorry.\n" );
00279 cdEXIT(EXIT_FAILURE);
00280 }
00281
00282 if( cmin < -37. )
00283 {
00284 fac = -cmin - 37.;
00285
00286 for( long i=0; i <= npairs; i++ )
00287 {
00288 rfield.tFluxLog[rfield.nShape][i] += (realnum)fac;
00289 }
00290 }
00291
00292 else if( cmax > 37. )
00293 {
00294 fac = cmax - 37.;
00295
00296 for( long i=0; i <= npairs; i++ )
00297 {
00298 rfield.tFluxLog[rfield.nShape][i] -= (realnum)fac;
00299 }
00300 }
00301
00302
00303 if( trace.lgConBug && trace.lgTrace )
00304 {
00305 fprintf( ioQQQ, " Table for this continuum;\ni\tTNU\tTFAC\tTSLOP, npairs=%li\n",
00306 npairs );
00307 for( long i=0; i < npairs; i++ )
00308 {
00309 fprintf( ioQQQ, "%li\t%.4e\t%.4e\t%.4e\n",
00310 i , rfield.tNu[rfield.nShape][i].Ryd(),
00311 rfield.tFluxLog[rfield.nShape][i], rfield.tslop[rfield.nShape][i] );
00312 }
00313 fprintf( ioQQQ, "%li\t%.4e\t%.4e\n",
00314 npairs , rfield.tNu[rfield.nShape][npairs].Ryd(),
00315 rfield.tFluxLog[rfield.nShape][npairs]);
00316 }
00317
00318
00319
00320
00321
00322 long n;
00323 for( n=0; n <= npairs; n++ )
00324 {
00325 if( rfield.tNu[rfield.nShape][n].Ryd() > 1. )
00326 break;
00327 }
00328
00329
00330 fac = rfield.tFluxLog[rfield.nShape][max(n-1,0)];
00331
00332 for( long i=0; i <= npairs; i++ )
00333 {
00334 rfield.tFluxLog[rfield.nShape][i] -= (realnum)fac;
00335
00336 }
00337
00338
00339 cmin = log10( FLT_MIN );
00340 cmax = log10( FLT_MAX );
00341 lgHit = false;
00342 for( long i=0; i <= npairs; i++ )
00343 {
00344 if( rfield.tFluxLog[rfield.nShape][i] <= cmin ||
00345 rfield.tFluxLog[rfield.nShape][i] >= cmax )
00346 {
00347 fprintf(ioQQQ,
00348 " The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i);
00349 fprintf(ioQQQ,
00350 " The frequency is %f and the log of the flux is %f.\n\n",
00351 rfield.tNu[rfield.nShape][i].Ryd() ,
00352 rfield.tFluxLog[rfield.nShape][i]);
00353 lgHit = true;
00354 }
00355 }
00356 if( lgHit )
00357 {
00358 fprintf(ioQQQ,"\n NOTE The log of the flux given in an interpolate command is outside the range of this cpu.\n");
00359 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");
00360 fprintf(ioQQQ," NOTE Note that the interpolate command only is used for the shape of the continuum.\n");
00361 fprintf(ioQQQ," NOTE The order of magnitude of the flux is not used in any way.\n");
00362 fprintf(ioQQQ," NOTE For safety this could be of order unity.\n\n");
00363 }
00364
00365 rfield.ncont[rfield.nShape] = npairs+1;
00366
00367
00368 for( long i=npairs+1; i < rfield.nupper; i++ )
00369 {
00370
00371
00372 rfield.tFluxLog[rfield.nShape][i] = 0.;
00373 rfield.tslop[rfield.nShape][i] = 0.;
00374 rfield.tNu[rfield.nShape][i].set(0.);
00375 }
00376
00377
00378 ++rfield.nShape;
00379
00380
00381
00382 ASSERT( rfield.nShape < LIMSPC );
00383
00384 return;
00385 }