Context Navigation

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

Last change on this file since 996 was 834, checked in by srkline, 11 years ago

Changes to the XOP code to upgrade to ToolKit? v6. Changes are the ones outlined in the Appendix A of the TK6 manual. SOme of the XOP support routines and the #pragma for 2-byte structures have changed. Per Howard Rodstein, there is no need to change the c files to cpp. c should work and compile just fine.

These changes work correctly on my mac. Next is to make sure that they work correctly on the two build machines.

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

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

Download in other formats: