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