Context Navigation

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

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

Download in other formats: