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 /*gridXspec handles all grid calculations, called by griddo */ 00004 /*gridFunc */ 00005 /*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */ 00006 #include "cddefines.h" 00007 #include "punch.h" 00008 #include "warnings.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 /*gridXspec handles all grid calculations, called by grid_do */ 00017 void gridXspec(realnum xc[], long int nInterpVars) 00018 { 00019 long int i, j, k; 00020 double averageChi2; 00021 00022 DEBUG_ENTRY( "gridXspec()" ); 00023 00024 if( nInterpVars > LIMPAR ) 00025 { 00026 fprintf( ioQQQ, "grid_do: too many parameters are varied, increase LIMPAR\n" ); 00027 cdEXIT(EXIT_FAILURE); 00028 } 00029 00030 optimize.nOptimiz = 0; 00031 grid.nintparm = nInterpVars; 00032 00033 /* if this is changed there must be some change made to actually 00034 * stuff the additional parameter information. */ 00035 grid.naddparm = 0; 00036 00037 ASSERT( grid.nintparm + grid.naddparm <= LIMPAR ); 00038 00039 grid.totNumModels = 1; 00040 /* >>chng 06 aug 21, allow the number of parameter values to be different for different parameters. */ 00041 for( i=0; i<nInterpVars; i++ ) 00042 { 00043 /* >>chng 06 sep 4, use grid variable instead of passing to routine. */ 00044 grid.totNumModels *= grid.numParamValues[i]; 00045 } 00046 /* grid.totNumModels = (long)pow((double)nParVals, (double)nInterpVars); */ 00047 ASSERT( grid.totNumModels > 1 ); 00048 00049 grid.paramNames = (char**)MALLOC(sizeof(char*)*(unsigned)(nInterpVars+grid.naddparm) ); 00050 grid.paramMethods = (long*)MALLOC(sizeof(long)*(unsigned)(nInterpVars+grid.naddparm) ); 00051 grid.paramRange = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) ); 00052 grid.paramData = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) ); 00053 grid.interpParameters = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) ); 00054 /* save abort flag */ 00055 grid.lgAbort = (bool*)MALLOC(sizeof(bool)*(unsigned)(grid.totNumModels) ); 00056 /* number of warnings */ 00057 grid.lgWarn = (bool*)MALLOC(sizeof(bool)*(unsigned)(grid.totNumModels) ); 00058 00059 for( i=0; i<nInterpVars+grid.naddparm; i++ ) 00060 { 00061 grid.paramNames[i] = (char*)MALLOC(sizeof(char)*(unsigned)(12) ); 00062 grid.paramRange[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(6) ); 00063 grid.paramData[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numParamValues[i]) ); 00064 00065 sprintf( grid.paramNames[i], "%s%ld", "PARAM", i+1 ); 00066 /* Method is 0 for linear, 1 for logarithmic, for now all linear. */ 00067 grid.paramMethods[i] = 0; 00068 /* Initial */ 00069 grid.paramRange[i][0] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)/2.f; 00070 /* Delta */ 00071 grid.paramRange[i][1] = grid.paramIncrements[i]/10.f; 00072 /* Minimum */ 00073 grid.paramRange[i][2] = xc[i]-grid.paramIncrements[i]/10.f; 00074 /* Bottom */ 00075 grid.paramRange[i][3] = xc[i]-grid.paramIncrements[i]/10.f; 00076 /* Top */ 00077 grid.paramRange[i][4] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f; 00078 /* Maximum */ 00079 grid.paramRange[i][5] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f; 00080 00081 for( j=0; j<grid.numParamValues[i]; j++ ) 00082 { 00083 grid.paramData[i][j] = xc[i]+grid.paramIncrements[i]*j; 00084 } 00085 } 00086 00087 for( i=0; i<grid.totNumModels; i++ ) 00088 { 00089 grid.interpParameters[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nInterpVars) ); 00090 } 00091 00092 /* >>chng 06 aug 23, this logic now allows non-square parameter spaces. */ 00093 for( i=0; i< grid.totNumModels; i++ ) 00094 { 00095 realnum variableVector[LIMPAR]; 00096 00097 for( j=0; j<nInterpVars; j++ ) 00098 { 00099 int index; 00100 long volumeOtherDimensions = 1; 00101 00102 /* by "volume", we simply mean the product of the parameter dimensions 00103 * AFTER the present index, i.e.: 00104 * first "volume" is product of grid.numParamValues[1]*grid.numParamValues[2]*....grid.numParamValues[n] 00105 * second "volume" is product of grid.numParamValues[2]*grid.numParamValues[3]*....grid.numParamValues[n] 00106 * last "volume" is unity. */ 00107 for( k=j+1; k<nInterpVars; k++ ) 00108 { 00109 volumeOtherDimensions *= grid.numParamValues[k]; 00110 } 00111 00112 /* For each successive parameter, the "volume" is less than the previous one. 00113 * So the left-hand side of this modulus operation increases for each parameter, 00114 * which causes the index of each parameter to be incremented more often than the 00115 * index of the previous parameter. Thus, the last dimension is incremented 00116 * every time, then the second to last dimension is incremented with each repeat 00117 * of the last dimension. This repeats until, finally, the first dimension is 00118 * incremented. For example, the indices of the parameter vectors for a 2x2x3 00119 * box would be ordered as such: 00120 * [0][0][0] 00121 * [0][0][1] 00122 * [0][0][2] 00123 * [0][1][0] 00124 * [0][1][1] 00125 * [0][1][2] 00126 * [1][0][0] 00127 * [1][0][1] 00128 * [1][0][2] 00129 * [1][1][0] 00130 * [1][1][1] 00131 * [1][1][2] 00132 */ 00133 index = (int)( (i/volumeOtherDimensions)%(grid.numParamValues[j]) ); 00134 00135 /* this prevents parameter incrementation for debugging purposes. */ 00136 if( grid.lgStrictRepeat ) 00137 variableVector[j] = xc[j]; 00138 else 00139 variableVector[j] = xc[j] + grid.paramIncrements[j]*index; 00140 00141 grid.interpParameters[i][j] = variableVector[j]; 00142 } 00143 00144 for( j=nInterpVars; j<LIMPAR; j++ ) 00145 { 00146 variableVector[j] = xc[j]; 00147 } 00148 00149 if( i == grid.totNumModels - 1 ) 00150 { 00151 called.lgTalk = true; 00152 called.lgTalkIsOK = true; 00153 prt.lgFaintOn = true; 00154 grid.lgGridDone = true; 00155 } 00156 00157 averageChi2 = optimize_func(variableVector); 00158 /* silly test to use the var averageChi2 for something - keeps lint and some 00159 * compilers happy */ 00160 if( averageChi2 < 0. && averageChi2 == 0 ) 00161 fprintf( ioQQQ , "DEBUG interesting impossible print has occurred.\n"); 00162 00163 } 00164 return; 00165 } 00166 00167 /* GridGatherAfterCloudy - gather information just before 00168 * punch output generated */ 00169 void GridGatherAfterCloudy( 00170 /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */ 00171 const char *chTime) 00172 { 00173 DEBUG_ENTRY( "GridGatherAfterCloudy()" ); 00174 00175 if( !grid.lgGrid ) 00176 { 00177 /* not grid - return */ 00178 return; 00179 } 00180 00181 /* this has to be either MIDL or LAST */ 00182 if( chTime[0] == 'L' ) 00183 { 00184 ASSERT( optimize.nOptimiz>= 0 && optimize.nOptimiz<grid.totNumModels ); 00185 grid.lgAbort[optimize.nOptimiz] = lgAbort; 00186 grid.lgWarn[optimize.nOptimiz] = warnings.lgWarngs; 00187 } 00188 else if( chTime[0] != 'M' ) 00189 TotalInsanity(); 00190 return; 00191 } 00192 00193 /*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */ 00194 void GridGatherInCloudy( void ) 00195 { 00196 long i; 00197 static bool lgFirstRun = true; 00198 00199 DEBUG_ENTRY( "GridGatherInCloudy()" ); 00200 00201 if( !grid.lgGrid ) 00202 { 00203 TotalInsanity(); 00204 /* not grid - return */ 00205 return; 00206 } 00207 00208 /* malloc some arrays if first call and save continuum energies. */ 00209 if( lgFirstRun ) 00210 { 00211 long i1, i2; 00212 grid.numEnergies = rfield.nupper-2; 00213 grid.Energies = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numEnergies) ); 00214 grid.Spectra = (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(NUM_OUTPUT_TYPES) ); 00215 00216 for( i1=0; i1< NUM_OUTPUT_TYPES; i1++ ) 00217 { 00218 if( grid.lgOutputTypeOn[i1] ) 00219 { 00220 grid.Spectra[i1] = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) ); 00221 for( i2=0; i2<grid.totNumModels; i2++ ) 00222 { 00223 grid.Spectra[i1][i2] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numEnergies) ); 00224 } 00225 } 00226 } 00227 00228 for( i1=0; i1<grid.numEnergies; i1++ ) 00229 { 00230 grid.Energies[i1] = rfield.AnuOrg[i1]; 00231 } 00232 lgFirstRun = false; 00233 } 00234 00235 ASSERT( lgFirstRun == false ); 00236 ASSERT( optimize.nOptimiz>=0 && optimize.nOptimiz<grid.totNumModels); 00237 00238 for( i=1; i< NUM_OUTPUT_TYPES; i++ ) 00239 { 00240 /* Grab spectrum for xspec printout 00241 * at this point nOptimiz has already been incremented for first model */ 00242 if( grid.lgOutputTypeOn[i] ) 00243 cdSPEC2( i, grid.numEnergies, grid.Spectra[i][optimize.nOptimiz]); 00244 } 00245 return; 00246 }