00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "punch.h"
00009 #include "optimize.h"
00010 #include "cddrive.h"
00011 #include "rfield.h"
00012 #include "grid.h"
00013 #include "called.h"
00014 #include "prt.h"
00015
00016 static int FirstRun = true;
00017
00018
00019
00020
00021 void gridXspec(float xc[], long int nInterpVars)
00022 {
00023 long int i, j, k;
00024 double averageChi2;
00025
00026 DEBUG_ENTRY( "grid_xspec()" );
00027
00028 if( nInterpVars > LIMPAR )
00029 {
00030 fprintf( ioQQQ, "grid_do: too many parameters are varied, increase LIMPAR\n" );
00031 puts( "[Stop]" );
00032 cdEXIT(EXIT_FAILURE);
00033 }
00034
00035 optimize.nOptimiz = 0;
00036
00037 grid.nintparm = nInterpVars;
00038
00039
00040
00041 grid.naddparm = 0;
00042
00043 ASSERT( grid.nintparm + grid.naddparm <= LIMPAR );
00044
00045 grid.totNumModels = 1;
00046
00047 for( i=0; i<nInterpVars; i++ )
00048 {
00049
00050 grid.totNumModels *= grid.numParamValues[i];
00051 }
00052
00053 ASSERT( grid.totNumModels > 1 );
00054
00055 grid.paramNames = (char**)MALLOC(sizeof(char*)*(unsigned)(nInterpVars+grid.naddparm) );
00056
00057 grid.paramMethods = (long*)MALLOC(sizeof(long)*(unsigned)(nInterpVars+grid.naddparm) );
00058
00059 grid.paramRange = (float**)MALLOC(sizeof(float*)*(unsigned)(nInterpVars+grid.naddparm) );
00060
00061 grid.paramData = (float**)MALLOC(sizeof(float*)*(unsigned)(nInterpVars+grid.naddparm) );
00062
00063 grid.interpParameters = (float**)MALLOC(sizeof(float*)*(unsigned)(grid.totNumModels) );
00064
00065 for( i=0; i<nInterpVars+grid.naddparm; i++ )
00066 {
00067 grid.paramNames[i] = (char*)MALLOC(sizeof(char)*(unsigned)(12) );
00068
00069 grid.paramRange[i] = (float*)MALLOC(sizeof(float)*(unsigned)(6) );
00070
00071 grid.paramData[i] = (float*)MALLOC(sizeof(float)*(unsigned)(grid.numParamValues[i]) );
00072
00073 sprintf( grid.paramNames[i], "%s%ld", "PARAM", i+1 );
00074
00075 grid.paramMethods[i] = 0;
00076
00077 grid.paramRange[i][0] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)/2.f;
00078
00079 grid.paramRange[i][1] = grid.paramIncrements[i]/10.f;
00080
00081 grid.paramRange[i][2] = xc[i]-grid.paramIncrements[i]/10.f;
00082
00083 grid.paramRange[i][3] = xc[i]-grid.paramIncrements[i]/10.f;
00084
00085 grid.paramRange[i][4] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f;
00086
00087 grid.paramRange[i][5] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f;
00088
00089 for( j=0; j<grid.numParamValues[i]; j++ )
00090 {
00091 grid.paramData[i][j] = xc[i]+grid.paramIncrements[i]*j;
00092 }
00093 }
00094
00095 for( i=0; i<grid.totNumModels; i++ )
00096 {
00097 grid.interpParameters[i] = (float*)MALLOC(sizeof(float)*(unsigned)(nInterpVars) );
00098 }
00099
00100
00101 for( i=0; i< grid.totNumModels; i++ )
00102 {
00103 float variableVector[LIMPAR];
00104
00105 for( j=0; j<nInterpVars; j++ )
00106 {
00107 int index;
00108 long volumeOtherDimensions = 1;
00109
00110
00111
00112
00113
00114
00115 for( k=j+1; k<nInterpVars; k++ )
00116 {
00117 volumeOtherDimensions *= grid.numParamValues[k];
00118 }
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 index = (int)( (i/volumeOtherDimensions)%(grid.numParamValues[j]) );
00142
00143 variableVector[j] = xc[j] + grid.paramIncrements[j]*index;
00144
00145 grid.interpParameters[i][j] = variableVector[j];
00146 }
00147 for( j=nInterpVars; j<LIMPAR; j++ )
00148 {
00149 variableVector[j] = xc[j];
00150 }
00151
00152 if( i == grid.totNumModels - 1 )
00153 {
00154 called.lgTalk = true;
00155 called.lgTalkIsOK = true;
00156 prt.lgFaintOn = true;
00157 punch.lgOpenUnits = true;
00158 grid.lgGridDone = true;
00159 }
00160
00161 averageChi2 = optimize_func(variableVector);
00162
00163
00164 if( averageChi2 < 0. && averageChi2 == 0 )
00165 fprintf( ioQQQ , "DEBUG interesting impossible print has occurred.\n");
00166
00167
00168 ASSERT( FirstRun == false );
00169
00170 }
00171
00172 DEBUG_EXIT( "grid_xspec()" );
00173 return;
00174 }
00175
00176
00177
00178 #if 0
00179 static void gridFunc(float param[] )
00180 {
00181 bool lgBAD, lgLimOK;
00182
00183 long int i;
00184
00185 double snorm;
00186
00187 DEBUG_ENTRY( "grid_func()" );
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204 zero();
00205
00206 if( optimize.lgOptimFlow )
00207 {
00208 fprintf( ioQQQ, " trace, grid_func variables" );
00209 for( i=0; i < optimize.nvary; i++ )
00210 {
00211 fprintf( ioQQQ, "%10.2e", param[i] );
00212 }
00213 fprintf( ioQQQ, "\n" );
00214 }
00215
00216 for( i=0; i < optimize.nvary; i++ )
00217 {
00218 optimize.vparm[0][i] = param[i];
00219 }
00220
00221
00222
00223 ++optimize.nOptimiz;
00224
00225
00226
00227 vary_input(&lgLimOK);
00228
00229 if( !lgLimOK )
00230 {
00231
00232
00233 fprintf( ioQQQ, " Iteration %ld not within range.\n",
00234 optimize.nOptimiz );
00235
00236 cdEXIT(EXIT_FAILURE);
00237 }
00238
00239 for( i=0; i < optimize.nvary; i++ )
00240 {
00241 optimize.varmax[i] = (float)MAX2(optimize.varmax[i],optimize.vpused[i]);
00242 optimize.varmin[i] = (float)MIN2(optimize.varmin[i],optimize.vpused[i]);
00243 }
00244
00245 lgBAD = cloudy();
00246
00247 if( lgBAD )
00248 {
00249 fprintf( ioQQQ, " Cloudy returned error condition - what happened?\n" );
00250 }
00251
00252 GridInitialize();
00253
00254 if( LineSave.ipNormWavL < 0 )
00255 {
00256 fprintf( ioQQQ,
00257 " Normalization line array index is bad. What has gone wrong?\n" );
00258 puts( "[Stop in grid_func]" );
00259 cdEXIT(EXIT_FAILURE);
00260 }
00261 snorm = LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent];
00262
00263 if( snorm == 0. )
00264 {
00265 fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity. What has gone wrong?\n" );
00266 fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
00267 puts( "[Stop in grid_func]" );
00268 cdEXIT(EXIT_FAILURE);
00269 }
00270
00271
00272 cdWarnings(ioQQQ);
00273
00274 DEBUG_EXIT( "grid_func()" );
00275 return;
00276 }
00277 #endif
00278
00279
00280 void GridInitialize( void )
00281 {
00282 long i,j;
00283
00284 DEBUG_ENTRY( "GridInitialize()" );
00285
00286
00287 if( FirstRun )
00288 {
00289 grid.numEnergies = rfield.nupper-2;
00290
00291 grid.Energies = (float*)MALLOC(sizeof(float)*(unsigned)(grid.numEnergies) );
00292
00293 grid.Spectra = (float***)MALLOC(sizeof(float**)*(unsigned)(NUM_OUTPUT_TYPES) );
00294
00295 for( i=0; i< NUM_OUTPUT_TYPES; i++ )
00296 {
00297
00298 {
00299 grid.Spectra[i] = (float**)MALLOC(sizeof(float*)*(unsigned)(grid.totNumModels) );
00300
00301 for( j=0; j<grid.totNumModels; j++ )
00302 {
00303 grid.Spectra[i][j] = (float*)MALLOC(sizeof(float)*(unsigned)(grid.numEnergies) );
00304 }
00305 }
00306 }
00307
00308 for( i=0; i<grid.numEnergies; i++ )
00309 {
00310 grid.Energies[i] = rfield.AnuOrg[i];
00311 grid.Energies[i] = rfield.AnuOrg[i];
00312 }
00313
00314 FirstRun = false;
00315 }
00316
00317 for( i=1; i< NUM_OUTPUT_TYPES; i++ )
00318 {
00319
00320 {
00321
00322 cdSPEC2( i, grid.numEnergies, grid.Spectra[i][optimize.nOptimiz-1]);
00323 }
00324 }
00325
00326 DEBUG_EXIT( "grid_func()" );
00327
00328 return;
00329 }
00330
00331
00332 void GridGather( void )
00333 {
00334 DEBUG_ENTRY( "GridGather()" );
00335
00336 long i;
00337
00338 if( optimize.nOptimiz==grid.totNumModels && grid.lgGrid )
00339 {
00340 ASSERT( optimize.nOptimiz >= 1 );
00341 for( i=1; i< NUM_OUTPUT_TYPES; i++ )
00342 {
00343
00344 {
00345
00346 cdSPEC2( i, grid.numEnergies, grid.Spectra[i][optimize.nOptimiz-1]);
00347 }
00348 }
00349 }
00350
00351 DEBUG_EXIT( "GridGather()" );
00352
00353 return;
00354 }