Context Navigation

MonteCarlo.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: 15.9 KB

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

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

Context Navigation

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

Download in other formats: