Context Navigation

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

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

Download in other formats: