cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
grid_xspec.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*gridXspec handles all grid calculations, called by griddo */
4 /*gridFunc */
5 /*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */
6 #include "cddefines.h"
7 #include "punch.h"
8 #include "warnings.h"
9 #include "optimize.h"
10 #include "cddrive.h"
11 #include "rfield.h"
12 #include "grid.h"
13 #include "called.h"
14 #include "prt.h"
15 
16 /*gridXspec handles all grid calculations, called by grid_do */
17 void gridXspec(realnum xc[], long int nInterpVars)
18 {
19  long int i, j, k;
20  double averageChi2;
21 
22  DEBUG_ENTRY( "gridXspec()" );
23 
24  if( nInterpVars > LIMPAR )
25  {
26  fprintf( ioQQQ, "grid_do: too many parameters are varied, increase LIMPAR\n" );
27  cdEXIT(EXIT_FAILURE);
28  }
29 
30  optimize.nOptimiz = 0;
31  grid.nintparm = nInterpVars;
32 
33  /* if this is changed there must be some change made to actually
34  * stuff the additional parameter information. */
35  grid.naddparm = 0;
36 
38 
39  grid.totNumModels = 1;
40  /* >>chng 06 aug 21, allow the number of parameter values to be different for different parameters. */
41  for( i=0; i<nInterpVars; i++ )
42  {
43  /* >>chng 06 sep 4, use grid variable instead of passing to routine. */
45  }
46  /* grid.totNumModels = (long)pow((double)nParVals, (double)nInterpVars); */
47  ASSERT( grid.totNumModels > 1 );
48 
49  grid.paramNames = (char**)MALLOC(sizeof(char*)*(unsigned)(nInterpVars+grid.naddparm) );
50  grid.paramMethods = (long*)MALLOC(sizeof(long)*(unsigned)(nInterpVars+grid.naddparm) );
51  grid.paramRange = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
52  grid.paramData = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
53  grid.interpParameters = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) );
54  /* save abort flag */
55  grid.lgAbort = (bool*)MALLOC(sizeof(bool)*(unsigned)(grid.totNumModels) );
56  /* number of warnings */
57  grid.lgWarn = (bool*)MALLOC(sizeof(bool)*(unsigned)(grid.totNumModels) );
58 
59  for( i=0; i<nInterpVars+grid.naddparm; i++ )
60  {
61  grid.paramNames[i] = (char*)MALLOC(sizeof(char)*(unsigned)(12) );
62  grid.paramRange[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(6) );
63  grid.paramData[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numParamValues[i]) );
64 
65  sprintf( grid.paramNames[i], "%s%ld", "PARAM", i+1 );
66  /* Method is 0 for linear, 1 for logarithmic, for now all linear. */
67  grid.paramMethods[i] = 0;
68  /* Initial */
69  grid.paramRange[i][0] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)/2.f;
70  /* Delta */
71  grid.paramRange[i][1] = grid.paramIncrements[i]/10.f;
72  /* Minimum */
73  grid.paramRange[i][2] = xc[i]-grid.paramIncrements[i]/10.f;
74  /* Bottom */
75  grid.paramRange[i][3] = xc[i]-grid.paramIncrements[i]/10.f;
76  /* Top */
77  grid.paramRange[i][4] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f;
78  /* Maximum */
79  grid.paramRange[i][5] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f;
80 
81  for( j=0; j<grid.numParamValues[i]; j++ )
82  {
83  grid.paramData[i][j] = xc[i]+grid.paramIncrements[i]*j;
84  }
85  }
86 
87  for( i=0; i<grid.totNumModels; i++ )
88  {
89  grid.interpParameters[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nInterpVars) );
90  }
91 
92  /* >>chng 06 aug 23, this logic now allows non-square parameter spaces. */
93  for( i=0; i< grid.totNumModels; i++ )
94  {
95  realnum variableVector[LIMPAR];
96 
97  for( j=0; j<nInterpVars; j++ )
98  {
99  int index;
100  long volumeOtherDimensions = 1;
101 
102  /* by "volume", we simply mean the product of the parameter dimensions
103  * AFTER the present index, i.e.:
104  * first "volume" is product of grid.numParamValues[1]*grid.numParamValues[2]*....grid.numParamValues[n]
105  * second "volume" is product of grid.numParamValues[2]*grid.numParamValues[3]*....grid.numParamValues[n]
106  * last "volume" is unity. */
107  for( k=j+1; k<nInterpVars; k++ )
108  {
109  volumeOtherDimensions *= grid.numParamValues[k];
110  }
111 
112  /* For each successive parameter, the "volume" is less than the previous one.
113  * So the left-hand side of this modulus operation increases for each parameter,
114  * which causes the index of each parameter to be incremented more often than the
115  * index of the previous parameter. Thus, the last dimension is incremented
116  * every time, then the second to last dimension is incremented with each repeat
117  * of the last dimension. This repeats until, finally, the first dimension is
118  * incremented. For example, the indices of the parameter vectors for a 2x2x3
119  * box would be ordered as such:
120  * [0][0][0]
121  * [0][0][1]
122  * [0][0][2]
123  * [0][1][0]
124  * [0][1][1]
125  * [0][1][2]
126  * [1][0][0]
127  * [1][0][1]
128  * [1][0][2]
129  * [1][1][0]
130  * [1][1][1]
131  * [1][1][2]
132  */
133  index = (int)( (i/volumeOtherDimensions)%(grid.numParamValues[j]) );
134 
135  /* this prevents parameter incrementation for debugging purposes. */
136  if( grid.lgStrictRepeat )
137  variableVector[j] = xc[j];
138  else
139  variableVector[j] = xc[j] + grid.paramIncrements[j]*index;
140 
141  grid.interpParameters[i][j] = variableVector[j];
142  }
143 
144  for( j=nInterpVars; j<LIMPAR; j++ )
145  {
146  variableVector[j] = xc[j];
147  }
148 
149  if( i == grid.totNumModels - 1 )
150  {
151  called.lgTalk = true;
152  called.lgTalkIsOK = true;
153  prt.lgFaintOn = true;
154  grid.lgGridDone = true;
155  }
156 
157  averageChi2 = optimize_func(variableVector);
158  /* silly test to use the var averageChi2 for something - keeps lint and some
159  * compilers happy */
160  if( averageChi2 < 0. && averageChi2 == 0 )
161  fprintf( ioQQQ , "DEBUG interesting impossible print has occurred.\n");
162 
163  }
164  return;
165 }
166 
167 /* GridGatherAfterCloudy - gather information just before
168  * punch output generated */
170  /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
171  const char *chTime)
172 {
173  DEBUG_ENTRY( "GridGatherAfterCloudy()" );
174 
175  if( !grid.lgGrid )
176  {
177  /* not grid - return */
178  return;
179  }
180 
181  /* this has to be either MIDL or LAST */
182  if( chTime[0] == 'L' )
183  {
187  }
188  else if( chTime[0] != 'M' )
189  TotalInsanity();
190  return;
191 }
192 
193 /*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */
194 void GridGatherInCloudy( void )
195 {
196  long i;
197  static bool lgFirstRun = true;
198 
199  DEBUG_ENTRY( "GridGatherInCloudy()" );
200 
201  if( !grid.lgGrid )
202  {
203  TotalInsanity();
204  /* not grid - return */
205  return;
206  }
207 
208  /* malloc some arrays if first call and save continuum energies. */
209  if( lgFirstRun )
210  {
211  long i1, i2;
213  grid.Energies = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numEnergies) );
214  grid.Spectra = (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(NUM_OUTPUT_TYPES) );
215 
216  for( i1=0; i1< NUM_OUTPUT_TYPES; i1++ )
217  {
218  if( grid.lgOutputTypeOn[i1] )
219  {
220  grid.Spectra[i1] = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) );
221  for( i2=0; i2<grid.totNumModels; i2++ )
222  {
223  grid.Spectra[i1][i2] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numEnergies) );
224  }
225  }
226  }
227 
228  for( i1=0; i1<grid.numEnergies; i1++ )
229  {
230  grid.Energies[i1] = rfield.AnuOrg[i1];
231  }
232  lgFirstRun = false;
233  }
234 
235  ASSERT( lgFirstRun == false );
237 
238  for( i=1; i< NUM_OUTPUT_TYPES; i++ )
239  {
240  /* Grab spectrum for xspec printout
241  * at this point nOptimiz has already been incremented for first model */
242  if( grid.lgOutputTypeOn[i] )
244  }
245  return;
246 }

Generated for cloudy by doxygen 1.8.4