Context Navigation

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

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

Download in other formats: