Context Navigation

source: sans/XOP_Dev/MonteCarlo/MonteCarlo2.c @ 553

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

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

Download in other formats: