Context Navigation

MonteCarlo3.c @ 787

Last change on this file since 787 was 711, checked in by srkline, 13 years ago

Added wavelength distribution as part of the simulation

Fixed bug where pixel locations of (-1,0) row or column were incorrectly rounded towards zero, and ended up on the detector in row 0 or column 0, effectively doubling the count there.

Separated the 4 functions into 4 separate files, for ease in checking that they are all the same. Now the only difference in each file should be the RNG used. And they had better be different - I'm almost positive now that cross-calling of the same RNG by different threads is responsible for the crashing.

File size: 16.0 KB

Line
1	/*
2	* MonteCarlo3.c
3	* SANSMonteCarlo
4	*
5	* Created by Steve Kline on 7/1/10.
6	* Copyright 2010 NCNR. All rights reserved.
7	*
8	*/
9
10	#include "XOPStandardHeaders.h" // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
11	#include "MonteCarlo.h"
12
13	//static int gCallSpinProcess = 1; // Set to 1 to all user abort (cmd dot) and background processing.
14
15	// these versions are DIRECT COPIES of the main version in MonteCarlo.c
16	// make changes there and copy them here. All that changes here is that the random
17	// number calls are different.
18	//
19	// version X uses ran3
20	// version X2 uses ran1
21	// version X3 uses ran3a
22	// version X4 usus ran1a
23
24	int
25	Monte_SANSX3(MC_ParamsPtr p) {
26	double inputWave; / pointer to double precision wave data */
27	double ran_dev; / pointer to double precision wave data */
28	double nt; / pointer to double precision wave data */
29	double j1; / pointer to double precision wave data */
30	double j2; / pointer to double precision wave data */
31	double nn; / pointer to double precision wave data */
32	// double MC_linear_data; / pointer to double precision wave data */
33	double results; / pointer to double precision wave data */
34	double retVal; //return value
35
36	long imon;
37	double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas;
38	long ind,index,n_index;
39	double qmax,theta_max,q0,zpow;
40	long n1,n2,n3;
41	double dth,zz,xx,yy,phi;
42	double theta,ran,ll,rr;
43	long done,find_theta,err; //used as logicals
44	long xPixel,yPixel;
45	double vx,vy,vz,theta_z;
46	double sig_abs,ratio,sig_total;
47	double testQ,testPhi,left,delta,dummy,pi;
48	double sigabs_0,num_bins;
49	long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
50	long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
51	long NMultipleCoherent,NCoherentEvents;
52	double deltaLam,v1,v2,currWavelength,rsq,fac; //for simulating wavelength distribution
53
54
55	// for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
56	waveHndl wavH;
57	// int waveType,hState;
58	long numDimensions;
59	long dimensionSizes[MAX_DIMENSIONS+1];
60	// char* dataStartPtr;
61	// long dataOffset;
62	// long numRows, numColumns;
63	long numRows_ran_dev;
64	// double dp0, dp;
65	double value[2]; // Pointers used for double data.
66	long seed;
67	long indices[MAX_DIMENSIONS];
68
69	// char buf[256];
70
71	/* check that wave handles are all valid */
72	if (p->inputWaveH == NIL) {
73	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
74	return(NON_EXISTENT_WAVE);
75	}
76	if (p->ran_devH == NIL) {
77	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
78	return(NON_EXISTENT_WAVE);
79	}
80	if (p->ntH == NIL) {
81	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
82	return(NON_EXISTENT_WAVE);
83	}
84	if (p->j1H == NIL) {
85	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
86	return(NON_EXISTENT_WAVE);
87	}
88	if (p->j2H == NIL) {
89	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
90	return(NON_EXISTENT_WAVE);
91	}
92	if (p->nnH == NIL) {
93	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
94	return(NON_EXISTENT_WAVE);
95	}
96	if (p->MC_linear_dataH == NIL) {
97	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
98	return(NON_EXISTENT_WAVE);
99	}
100	if (p->resultsH == NIL) {
101	SetNaN64(&p->retVal); /* return NaN if wave is not valid */
102	return(NON_EXISTENT_WAVE);
103	}
104
105	p->retVal = 0;
106
107	// trusting that all inputs are DOUBLE PRECISION WAVES!!!
108	inputWave = WaveData(p->inputWaveH);
109	ran_dev = WaveData(p->ran_devH);
110	nt = WaveData(p->ntH);
111	j1 = WaveData(p->j1H);
112	j2 = WaveData(p->j2H);
113	nn = WaveData(p->nnH);
114	// MC_linear_data = WaveData(p->MC_linear_dataH);
115	results = WaveData(p->resultsH);
116
117	seed = (long)results[0];
118
119	// sprintf(buf, "input seed = %ld\r", seed);
120	// XOPNotice(buf);
121
122	if(seed >= 0) {
123	seed = -1234509876;
124	}
125
126	dummy = ran3a(&seed); //initialize the random sequence by passing in a negative value
127	seed = 12348765; //non-negative after that does nothing
128
129	imon = (int)inputWave[0];
130	r1 = inputWave[1];
131	r2 = inputWave[2];
132	xCtr = inputWave[3];
133	yCtr = inputWave[4];
134	sdd = inputWave[5];
135	pixSize = inputWave[6];
136	thick = inputWave[7];
137	wavelength = inputWave[8];
138	sig_incoh = inputWave[9];
139	sig_sas = inputWave[10];
140	deltaLam = inputWave[11];
141	xCtr_long = (long)(xCtr+0.5);
142	yCtr_long = (long)(yCtr+0.5);
143
144	dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left); //0 is the rows
145	if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
146	return retVal;
147	numRows_ran_dev = dimensionSizes[0];
148
149	pi = 4.0*atan(1.0);
150
151	// access the 2D wave data for writing using the direct method
152	wavH = p->MC_linear_dataH;
153	if (wavH == NIL)
154	return NOWAV;
155
156	// waveType = WaveType(wavH);
157	// if (waveType & NT_CMPLX)
158	// return NO_COMPLEX_WAVE;
159	// if (waveType==TEXT_WAVE_TYPE)
160	// return NUMERIC_ACCESS_ON_TEXT_WAVE;
161	// if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes))
162	// return retVal;
163	// numRows = dimensionSizes[0];
164	// numColumns = dimensionSizes[1];
165
166	// if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset))
167	// return retVal;
168
169	// hState = MoveLockHandle(wavH); // So wave data can't move. Remember to call HSetState when done.
170	// dataStartPtr = (char)(wavH) + dataOffset;
171	// dp0 = (double*)dataStartPtr; // Pointer to the start of the 2D wave data.
172
173	//scattering power and maximum qvalue to bin
174	// zpow = .1 //scattering power, calculated below
175	qmax = 4.0*pi/wavelength; //maximum Q to bin 1D data. (A-1) (not really used)
176	sigabs_0 = 0.0; // ignore absorption cross section/wavelength [1/(cm A)]
177	n_index = 50; // maximum number of scattering events per neutron
178	num_bins = 200; //number of 1-D bins (not really used)
179
180	//c total SAS cross-section
181	//
182	zpow = sig_sas*thick; //since I now calculate the sig_sas from the model
183	sig_abs = sigabs_0 * wavelength;
184	sig_total = sig_abs + sig_sas + sig_incoh;
185	// Print "The TOTAL XSECTION. (CM-1) is ",sig_total
186	// Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
187	// results[0] = sig_total;
188	// results[1] = sig_sas;
189	// RATIO = SIG_ABS / SIG_TOTAL
190	ratio = sig_incoh / sig_total;
191
192	theta_max = wavelengthqmax/(2pi);
193	//C SET Theta-STEP SIZE.
194	dth = theta_max/num_bins;
195	// Print "theta bin size = dth = ",dth
196
197	//C INITIALIZE COUNTERS.
198	n1 = 0;
199	n2 = 0;
200	n3 = 0;
201	NSingleIncoherent = 0;
202	NSingleCoherent = 0;
203	NDoubleCoherent = 0;
204	NMultipleScatter = 0;
205	NScatterEvents = 0;
206	NMultipleCoherent = 0;
207	NCoherentEvents = 0;
208
209	isOn = 0;
210
211	//C MONITOR LOOP - looping over the number of incedent neutrons
212	//note that zz, is the z-position in the sample - NOT the scattering power
213	// NOW, start the loop, throwing neutrons at the sample.
214	do {
215	////SpinProcess() IS A CALLBACK, and not good for Threading!
216	// if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) { // Spins cursor and allows background processing.
217	// retVal = -1; // User aborted.
218	// break;
219	// }
220
221	vx = 0.0; // Initialize direction vector.
222	vy = 0.0;
223	vz = 1.0;
224
225	theta = 0.0; // Initialize scattering angle.
226	phi = 0.0; // Intialize azimuthal angle.
227	n1 += 1; // Increment total number neutrons counter.
228	done = 0; // True when neutron is absorbed or when scattered out of the sample.
229	index = 0; // Set counter for number of scattering events.
230	zz = 0.0; // Set entering dimension of sample.
231	incoherentEvent = 0;
232	coherentEvent = 0;
233
234	do { // Makes sure position is within circle.
235	ran = ran3a(&seed); //[0,1]
236	xx = 2.0r1(ran-0.5); //X beam position of neutron entering sample.
237	ran = ran3a(&seed); //[0,1]
238	yy = 2.0r1(ran-0.5); //Y beam position ...
239	rr = sqrt(xxxx+yyyy); //Radial position of neutron in incident beam.
240	} while(rr>r1);
241
242	//pick the wavelength out of the wavelength spread, approximate as a gaussian
243	// from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the
244	// 2.35 converts to a gaussian std dev.
245	do {
246	v1=2.0*ran3a(&seed)-1.0;
247	v2=2.0*ran3a(&seed)-1.0;
248	rsq=v1v1+v2v2;
249	} while (rsq >= 1.0 \|\| rsq == 0.0);
250	fac=sqrt(-2.0*log(rsq)/rsq);
251
252	// gset=v1*fac //technically, I'm throwing away one of the two values
253	currWavelength = (v2fac)deltaLam*wavelength/2.35 + wavelength;
254
255	do { //Scattering Loop, will exit when "done" == 1
256	// keep scattering multiple times until the neutron exits the sample
257	ran = ran3a(&seed); //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH
258	ll = path_len(ran,sig_total);
259	//Determine new scattering direction vector.
260	err = NewDirection(&vx,&vy,&vz,theta,phi); //vx,vy,vz updated, theta, phi unchanged by function
261
262	//X,Y,Z-POSITION OF SCATTERING EVENT.
263	xx += ll*vx;
264	yy += ll*vy;
265	zz += ll*vz;
266	rr = sqrt(xxxx+yyyy); //radial position of scattering event.
267
268	//sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
269	//XOPNotice(buf);
270
271	//Check whether interaction occurred within sample volume.
272	if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
273	//NEUTRON INTERACTED.
274	//sprintf(buf,"neutron interacted\r");
275	//XOPNotice(buf);
276
277	index += 1; //Increment counter of scattering events.
278	if (index == 1) {
279	n2 += 1; //Increment # of scat. neutrons
280	}
281	ran = ran3a(&seed); //[0,1]
282	//Split neutron interactions into scattering and absorption events
283	if (ran > ratio ) { //C NEUTRON SCATTERED coherently
284	//sprintf(buf,"neutron scatters coherently\r");
285	//XOPNotice(buf);
286	coherentEvent += 1;
287	find_theta = 0; //false
288	do {
289	// pick a q-value from the deviate function
290	// pnt2x truncates the point to an integer before returning the x
291	// so get it from the wave scaling instead
292	// q0 =left + binarysearchinterp(ran_dev,ran3a(seed))*delta;
293
294	q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta;
295	theta = q0/2/pi*currWavelength; //SAS approximation
296
297	find_theta = 1; //always accept
298
299	//sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
300	//XOPNotice(buf);
301
302	} while(!find_theta);
303
304	ran = ran3a(&seed); //[0,1]
305	phi = 2.0piran; //Chooses azimuthal scattering angle.
306	} else {
307	//NEUTRON scattered incoherently
308	//sprintf(buf,"neutron scatters incoherent\r");
309	//XOPNotice(buf);
310	incoherentEvent += 1;
311	// phi and theta are random over the entire sphere of scattering
312	// !can't just choose random theta and phi, won't be random over sphere solid angle
313
314	ran = ran3a(&seed); //[0,1]
315	theta = acos(2.0*ran-1);
316
317	ran = ran3a(&seed); //[0,1]
318	phi = 2.0piran; //Chooses azimuthal scattering angle.
319	} //(ran > ratio)
320	} else {
321	//NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
322	done = 1; //done = true, will exit from loop
323	//Increment #scattering events array
324	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
325	indices[0] =index; //this sets access to nn[index]
326	if (index <= n_index) {
327	if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value))
328	return retVal;
329	value[0] += 1; // add one to the value
330	if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value))
331	return retVal;
332	// nn[index] += 1;
333	}
334
335	if( index != 0) { //neutron was scattered, figure out where it went
336	theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis.
337	testQ = 2pisin(theta_z)/currWavelength;
338
339	// pick a random phi angle, and see if it lands on the detector
340	// since the scattering is isotropic, I can safely pick a new, random value
341	// this would not be true if simulating anisotropic scattering.
342	testPhi = ran3a(&seed)2pi;
343
344	// is it on the detector?
345	FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
346
347	if(xPixel != -1 && yPixel != -1) {
348	isOn += 1;
349	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
350	indices[0] = xPixel;
351	indices[1] = yPixel;
352	if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
353	return retVal;
354	value[0] += 1; // Real part
355	if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
356	return retVal;
357	//if(index==1) // only the single scattering events
358	//dp = dp0 + xPixel + yPixel*numColumns; //offset the pointer to the exact memory location
359	//*dp += 1; //increment the value there
360	//endif
361	}
362
363
364	/* is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
365	if(theta_z < theta_max) {
366	//Choose index for scattering angle array.
367	//IND = NINT(THETA_z/DTH + 0.4999999)
368	ind = (long)(theta_z/dth + 0.4999999); //round is eqivalent to nint()
369	nt[ind] += 1; //Increment bin for angle.
370	//Increment angle array for single scattering events.
371	if (index == 1) {
372	j1[ind] += 1;
373	}
374	//Increment angle array for double scattering events.
375	if (index == 2) {
376	j2[ind] += 1;
377	}
378	}
379	/**/
380
381	// increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
382	NScatterEvents += index; //total number of scattering events
383	if(index == 1 && incoherentEvent == 1) {
384	NSingleIncoherent += 1;
385	}
386	if(index == 1 && coherentEvent == 1) {
387	NSingleCoherent += 1;
388	}
389	if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
390	NDoubleCoherent += 1;
391	}
392	if(index > 1) {
393	NMultipleScatter += 1;
394	}
395	if(coherentEvent >= 1 && incoherentEvent == 0) {
396	NCoherentEvents += 1;
397	}
398	if(coherentEvent > 1 && incoherentEvent == 0) {
399	NMultipleCoherent += 1;
400	}
401
402	} else { // index was zero, neutron must be transmitted, so just increment the proper counters and data
403	isOn += 1;
404	nt[0] += 1;
405	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
406	//indices[0] = xCtr_long; //don't put everything in one pixel
407	//indices[1] = yCtr_long;
408	indices[0] = (long)(xCtr+xx/pixSize+0.5);
409	indices[1] = (long)(yCtr+yy/pixSize+0.5);
410	// check for valid indices - got an XOP error, probably from here
411	if(indices[0] > 127) indices[0] = 127;
412	if(indices[0] < 0) indices[0] = 0;
413	if(indices[1] > 127) indices[1] = 127;
414	if(indices[1] < 0) indices[1] = 0;
415
416	if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
417	return retVal;
418	value[0] += 1; // Real part
419	if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
420	return retVal;
421	}
422	}
423	} while (!done);
424	} while(n1 < imon);
425
426	// assign the results to the wave
427
428	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
429	value[0] = (double)n1;
430	indices[0] = 0;
431	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
432	return retVal;
433	value[0] = (double)n2;
434	indices[0] = 1;
435	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
436	return retVal;
437	value[0] = (double)isOn;
438	indices[0] = 2;
439	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
440	return retVal;
441	value[0] = (double)NScatterEvents;
442	indices[0] = 3;
443	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
444	return retVal;
445	value[0] = (double)NSingleCoherent;
446	indices[0] = 4;
447	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
448	return retVal;
449	value[0] = (double)NMultipleCoherent;
450	indices[0] = 5;
451	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
452	return retVal;
453	value[0] = (double)NMultipleScatter;
454	indices[0] = 6;
455	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
456	return retVal;
457	value[0] = (double)NCoherentEvents;
458	indices[0] = 7;
459	if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
460	return retVal;
461
462	// HSetState((Handle)wavH, hState); //release the handle of the 2D data wave
463	// WaveHandleModified(wavH); // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
464
465	return(0);
466	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: sans/XOP_Dev/MonteCarlo/MonteCarlo3.c @ 787

Download in other formats: