Context Navigation

DebyeSpheres.c @ 1249

Last change on this file since 1249 was 997, checked in by srkline, 7 years ago
Updates for Xcode 7, Toolkit 7, and 64-bit version
Property svn:executable set to ``*
File size: 13.8 KB

Line
1	/* DebyeSpheres.c
2
3
4	*/
5
6	#include "XOPStandardHeaders.h" // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
7	#include "DebyeSpheres.h"
8
9	//#pragma XOP_SET_STRUCT_PACKING // All structures are 2-byte-aligned.
10
11	// Prototypes
12	HOST_IMPORT int XOPMain(IORecHandle ioRecHandle);
13
14	// Custom error codes
15	//#define REQUIRES_IGOR_200 1 + FIRST_XOP_ERR
16	#define NON_EXISTENT_WAVE 2 + FIRST_XOP_ERR
17	#define REQUIRES_SP_OR_DP_WAVE 3 + FIRST_XOP_ERR
18
19
20
21	/* Calculates the scattered intensity from a collection of spheres
22	using Debye's method.
23
24	everything was previously declared as float, some inputs were as double.
25
26	-- to avoid confusion, calculate everything in double, pass everything in as double
27	-- this is more precision than is necessary, but avoids strange results from type casting
28	behind the scenes
29	*/
30	int
31	DebyeSpheresX(AltiParamsPtr p)
32	{
33	double qv; // input q-value
34	double ival; //output intensity value
35	double xv,yv,zv,rv; //pointers to input xyz-rho coordinates
36	CountInt i,j;
37	CountInt npt;
38	double rval,grid,vol,fQR,dum,dij;
39
40
41
42	// check for all of the required waves
43	if (p->rhowavH == NIL) {
44	SetNaN64(&p->result);
45	return NON_EXISTENT_WAVE;
46	}
47	if (p->zwavH == NIL) {
48	SetNaN64(&p->result);
49	return NON_EXISTENT_WAVE;
50	}
51	if (p->ywavH == NIL) {
52	SetNaN64(&p->result);
53	return NON_EXISTENT_WAVE;
54	}
55	if (p->xwavH == NIL) {
56	SetNaN64(&p->result);
57	return NON_EXISTENT_WAVE;
58	}
59
60
61	//check to see that all are float, not double
62	/*
63	if(WaveType(p->rhowavH) != NT_FP32 ) {
64	SetNaN64(&p->result);
65	return kExpectedNT_FP32;
66	}
67	if(WaveType(p->zwavH) != NT_FP32 ) {
68	SetNaN64(&p->result);
69	return kExpectedNT_FP32;
70	}
71	if(WaveType(p->ywavH) != NT_FP32 ) {
72	SetNaN64(&p->result);
73	return kExpectedNT_FP32;
74	}
75	if(WaveType(p->xwavH) != NT_FP32 ) {
76	SetNaN64(&p->result);
77	return kExpectedNT_FP32;
78	}
79	*/
80
81	//check to see that all are double
82	if(WaveType(p->rhowavH) != NT_FP64 ) {
83	SetNaN64(&p->result);
84	return kExpectedNT_FP64;
85	}
86	if(WaveType(p->zwavH) != NT_FP64 ) {
87	SetNaN64(&p->result);
88	return kExpectedNT_FP64;
89	}
90	if(WaveType(p->ywavH) != NT_FP64 ) {
91	SetNaN64(&p->result);
92	return kExpectedNT_FP64;
93	}
94	if(WaveType(p->xwavH) != NT_FP64 ) {
95	SetNaN64(&p->result);
96	return kExpectedNT_FP64;
97	}
98
99
100	// (NO) -- convert the input to float. Do all calculations as float.
101	// do everything in double. no reason to do in float, except for the previous Altivec incarnation
102
103	rval = p->Rprimary; // primary sphere radius
104	grid = p->grid; // calling program should set this to 0.62*Rprimary
105	qv = p->qval;
106
107	//
108	npt = (CountInt) WavePoints(p->xwavH); //wavePoints returns long, number of XYZ points
109	xv = (double*)WaveData(p->xwavH); //xyz locations
110	yv = (double*)WaveData(p->ywavH);
111	zv = (double*)WaveData(p->zwavH);
112	rv = (double*)WaveData(p->rhowavH);
113
114
115	vol = 4.03.1415927/3.0rvalrvalrval;
116	ival = 0.0;
117	fQR = PhiQR(qv,rval);
118	//do the i=j sum
119	for(i=0;i<npt;i+=1) {
120	dum = rv[i]volfQR;
121	ival += dum*dum;
122	}
123	//do the i!=j double sum
124	dum = volvolfQR*fQR;
125	for(i=0;i<npt;i+=1) {
126	for(j=(i+1);j<npt;j+=1) {
127	dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]) * grid;
128	ival += 2.0rv[i]rv[j]dumsin(dij*qv)/dij/qv;
129	}
130	}
131
132	p->result = ival;
133
134	return 0;
135
136	}
137
138
139	double PhiQR(double qval, double rval)
140	{
141	double retval,qr;
142
143	qr = qval*rval;
144	retval = ( 3(sin(qr) - qrcos(qr))/qr/qr/qr );
145	return(retval);
146	}
147
148	double XYZDistance(double x1, double x2,double y1, double y2,double z1, double z2)
149	{
150	double retval,dx,dy,dz;
151
152	dx = (x1-x2);
153	dy = (y1-y2);
154	dz = (z1-z2);
155	retval = sqrt( dxdx + dydy + dz*dz );
156	return(retval);
157	}
158
159
160	/*
161
162	given the distances XYZ as a triplet (on a unit grid)
163	return the maximum distance. The calling program must multiply by
164	the grid dimension to get real distance
165
166	*/
167	int
168	maxDistanceX(DistParamPtr p)
169	{
170	double dmax,dij; //output dmax value, dij
171	double xv,yv,*zv; //pointers to input xyz coordinates
172	CountInt i,j;
173	CountInt npt;
174	CountInt p1,p2;
175
176	// check for all of the required waves
177	if (p->zwavH == NIL) {
178	SetNaN64(&p->result);
179	return NON_EXISTENT_WAVE;
180	}
181	if (p->ywavH == NIL) {
182	SetNaN64(&p->result);
183	return NON_EXISTENT_WAVE;
184	}
185	if (p->xwavH == NIL) {
186	SetNaN64(&p->result);
187	return NON_EXISTENT_WAVE;
188	}
189
190	//check to see that all are double
191	if(WaveType(p->zwavH) != NT_FP64 ) {
192	SetNaN64(&p->result);
193	return kExpectedNT_FP64;
194	}
195	if(WaveType(p->ywavH) != NT_FP64 ) {
196	SetNaN64(&p->result);
197	return kExpectedNT_FP64;
198	}
199	if(WaveType(p->xwavH) != NT_FP64 ) {
200	SetNaN64(&p->result);
201	return kExpectedNT_FP64;
202	}
203
204	//
205	npt = (CountInt) WavePoints(p->xwavH); //wavePoints returns long, number of XYZ points
206	xv = (double*)WaveData(p->xwavH); //xyz locations
207	yv = (double*)WaveData(p->ywavH);
208	zv = (double*)WaveData(p->zwavH);
209
210	p1 = (CountInt) p->p1;
211	p2 = (CountInt) p->p2;
212
213	dmax = 0;
214	//do the i!=j double loop, keeping the maximum distance
215
216	for(i=p1;i<p2;i+=1) {
217	for(j=(i+1);j<npt;j+=1) {
218	// dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]);
219	dij = (xv[i]-xv[j])(xv[i]-xv[j]) + (yv[i]-yv[j])(yv[i]-yv[j]) + (zv[i]-zv[j])*(zv[i]-zv[j]);
220	if(dij > dmax) {
221	dmax = dij;
222	}
223	}
224	}
225
226	p->result = dmax;
227
228	return 0;
229
230	}
231
232	/*
233
234	given the distances XYZ as a triplet (on a unit grid)
235	return the binned histogram of distances
236
237	*/
238	int
239	binDistanceX(BinParamPtr p)
240	{
241	double xv,yv,zv,bv; //pointers to input xyz coordinates
242	CountInt i,j;
243	CountInt npt,numBins,binIndex;
244	double grid,binWidth,val;
245	CountInt p1,p2;
246
247
248
249	// check for all of the required waves
250	if (p->bwavH == NIL) {
251	SetNaN64(&p->result);
252	return NON_EXISTENT_WAVE;
253	}
254	if (p->zwavH == NIL) {
255	SetNaN64(&p->result);
256	return NON_EXISTENT_WAVE;
257	}
258	if (p->ywavH == NIL) {
259	SetNaN64(&p->result);
260	return NON_EXISTENT_WAVE;
261	}
262	if (p->xwavH == NIL) {
263	SetNaN64(&p->result);
264	return NON_EXISTENT_WAVE;
265	}
266
267	//check to see that all are double
268	if(WaveType(p->bwavH) != NT_FP64 ) {
269	SetNaN64(&p->result);
270	return kExpectedNT_FP64;
271	}
272	if(WaveType(p->zwavH) != NT_FP64 ) {
273	SetNaN64(&p->result);
274	return kExpectedNT_FP64;
275	}
276	if(WaveType(p->ywavH) != NT_FP64 ) {
277	SetNaN64(&p->result);
278	return kExpectedNT_FP64;
279	}
280	if(WaveType(p->xwavH) != NT_FP64 ) {
281	SetNaN64(&p->result);
282	return kExpectedNT_FP64;
283	}
284
285	//
286	npt = (CountInt) WavePoints(p->xwavH); //wavePoints returns long, number of XYZ points
287	numBins = (CountInt) WavePoints(p->bwavH); //wavePoints returns long, number of points in bin wave
288	p1 = (CountInt) p->p1;
289	p2 = (CountInt) p->p2;
290
291
292	xv = (double*)WaveData(p->xwavH); //xyz locations
293	yv = (double*)WaveData(p->ywavH);
294	zv = (double*)WaveData(p->zwavH);
295	bv = (double*)WaveData(p->bwavH);
296
297	grid = p->grid;
298	binWidth = p->binWidth;
299
300	//do the i!=j double loop,
301	for(i=p1;i<p2;i+=1) {
302	for(j=(i+1);j<npt;j+=1) {
303	val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
304	binIndex = (CountInt)(val/binWidth-0.5);
305	if(binIndex > numBins -1 ) {
306	//Print "bad index"
307	} else {
308	bv[binIndex] += 1;
309	}
310
311	}
312	}
313
314	p->result = 0;
315
316	return 0;
317
318	}
319
320
321	/*
322
323	given the distances XYZ as a triplet (on a unit grid) and the SLD at each point,
324	return the binned histogram of distances for each of the parwise interactions
325
326	The returned binning is a matrix, and has to be assigned as such
327
328	*/
329	int
330	binSLDDistanceX(BinSLDParamPtr p)
331	{
332	double xv,yv,*zv; //pointers to input xyz coordinates
333	double rho,SLDLook,*PSFid; // rho and the SLD lookup vector
334	CountInt i,j;
335	CountInt npt,numBins,binIndex;
336	double grid,binWidth,val,retVal;
337	CountInt p1,p2;
338	CountInt intSLD;
339
340
341	// for accessing the 2D wave data to write the results
342	waveHndl wavH,PSFwavH;
343	// long numDimensions;
344	// long dimensionSizes[MAX_DIMENSIONS+1];
345	double value[2]; // Pointers used for double data.
346	IndexInt indices[MAX_DIMENSIONS];
347	//
348	CountInt rhoi,rhoj,rii,rji,PSFIndex;
349
350
351	// check for all of the required waves
352	if (p->PSFidH == NIL) {
353	SetNaN64(&p->result);
354	return NON_EXISTENT_WAVE;
355	}
356	if (p->SLDLookH == NIL) {
357	SetNaN64(&p->result);
358	return NON_EXISTENT_WAVE;
359	}
360	if (p->bwavH == NIL) {
361	SetNaN64(&p->result);
362	return NON_EXISTENT_WAVE;
363	}
364	if (p->rhowavH == NIL) {
365	SetNaN64(&p->result);
366	return NON_EXISTENT_WAVE;
367	}
368	if (p->zwavH == NIL) {
369	SetNaN64(&p->result);
370	return NON_EXISTENT_WAVE;
371	}
372	if (p->ywavH == NIL) {
373	SetNaN64(&p->result);
374	return NON_EXISTENT_WAVE;
375	}
376	if (p->xwavH == NIL) {
377	SetNaN64(&p->result);
378	return NON_EXISTENT_WAVE;
379	}
380
381	//check to see that all are double
382	if(WaveType(p->PSFidH) != NT_FP64 ) {
383	SetNaN64(&p->result);
384	return kExpectedNT_FP64;
385	}
386	if(WaveType(p->SLDLookH) != NT_FP64 ) {
387	SetNaN64(&p->result);
388	return kExpectedNT_FP64;
389	}
390	if(WaveType(p->bwavH) != NT_FP64 ) {
391	SetNaN64(&p->result);
392	return kExpectedNT_FP64;
393	}
394	if(WaveType(p->rhowavH) != NT_FP64 ) {
395	SetNaN64(&p->result);
396	return kExpectedNT_FP64;
397	}
398	if(WaveType(p->zwavH) != NT_FP64 ) {
399	SetNaN64(&p->result);
400	return kExpectedNT_FP64;
401	}
402	if(WaveType(p->ywavH) != NT_FP64 ) {
403	SetNaN64(&p->result);
404	return kExpectedNT_FP64;
405	}
406	if(WaveType(p->xwavH) != NT_FP64 ) {
407	SetNaN64(&p->result);
408	return kExpectedNT_FP64;
409	}
410
411
412	// access the 2D wave data for writing using the direct method
413	wavH = p->bwavH;
414	if (wavH == NIL)
415	return NOWAV;
416	//
417	PSFwavH = p->PSFidH;
418
419	npt = (CountInt) WavePoints(p->xwavH); //wavePoints returns long, number of XYZ points
420	numBins = (CountInt) WavePoints(p->bwavH); //wavePoints returns long, number of points in bin wave
421
422	xv = (double*)WaveData(p->xwavH); //xyz locations
423	yv = (double*)WaveData(p->ywavH);
424	zv = (double*)WaveData(p->zwavH);
425	rho = (double*)WaveData(p->rhowavH);
426	SLDLook = (double*)WaveData(p->SLDLookH);
427	PSFid = (double*)WaveData(p->PSFidH); //this one is 2D
428
429	p1 = (CountInt) p->p1;
430	p2 = (CountInt) p->p2;
431
432	intSLD = (CountInt) p->minSLD; //convert to int for use as index
433
434	grid = p->grid;
435	binWidth = p->binWidth;
436
437	//do the i!=j double loop,
438	for(i=p1;i<p2;i+=1) {
439	for(j=(i+1);j<npt;j+=1) {
440	val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
441	binIndex = (CountInt)(val/binWidth-0.5);
442	if(binIndex > numBins -1 ) {
443	//Print "bad index"
444	} else {
445	rhoi = (CountInt) rho[i]; //get the rho value at i and j
446	rhoj = (CountInt) rho[j];
447	rii = (CountInt) SLDLook[rhoi+intSLD]; //rho i index
448	rji = (CountInt) SLDLook[rhoj+intSLD]; //rho j index
449	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
450	indices[0] = (IndexInt)rii;
451	indices[1] = (IndexInt)rji;
452	if (retVal = MDGetNumericWavePointValue(PSFwavH, indices, value))
453	return retVal;
454	//PSFIndex = (long) PSFid[rii][rji]; //doesn't work
455	PSFIndex = (CountInt) value[0];
456
457	//now do the assignment to the 2D
458	// equivalent to binMatrix[binIndex][PSFIndex]
459
460	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
461	indices[0] = (IndexInt)binIndex;
462	indices[1] = (IndexInt)PSFIndex;
463	if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
464	return retVal;
465	value[0] += 1; // Real part
466	if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
467	return retVal;
468
469	}
470
471	}
472	}
473
474	p->result = 0;
475
476	return 0;
477
478	}
479
480
481	///// this is directly from Numerical Recipes
482	// -- I did change the float to double, since Igor treats all as double
483	// and n is an int, not a pointer (seemed unnecessary)
484	//
485	#define MAXBIT 30
486	#define MAXDIM 6
487	static int iminarg1,iminarg2;
488	#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))
489
490	int
491	SobolX(SobolParamPtr p)
492	{
493	int j,k,l;
494	unsigned long i,im,ipp;
495	static double fac;
496	static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
497	static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};
498	static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};
499	static unsigned long iv[MAXDIM*MAXBIT+1]={
500	0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};
501
502	static int initDone=0;
503	char buf[256];
504
505	int n=0;
506	double *x; //output x vector
507
508	// check for all of the required waves
509	if (p->bwavH == NIL) {
510	SetNaN64(&p->result);
511	return NON_EXISTENT_WAVE;
512	}
513
514	//check to see that all are double
515	if(WaveType(p->bwavH) != NT_FP64 ) {
516	SetNaN64(&p->result);
517	return kExpectedNT_FP64;
518	}
519	x = (double*)WaveData(p->bwavH);
520	n = (int)(p->nIn); // not sure that the negative input will be properly cast to int
521
522	// sprintf(buf, "input, recast n = %g %d\r",p->nIn, n);
523	// XOPNotice(buf);
524
525	if (n < 0) {
526
527	if(initDone) {
528	sprintf(buf, "Don't re-initialize\r");
529	XOPNotice(buf);
530	return 0;
531	}
532
533	for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k];
534	for (k=1;k<=MAXDIM;k++) {
535	for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j);
536	for (j=mdeg[k]+1;j<=MAXBIT;j++) {
537	ipp=ip[k];
538	i=iu[j-mdeg[k]][k];
539	i ^= (i >> mdeg[k]);
540	for (l=mdeg[k]-1;l>=1;l--) {
541	if (ipp & 1) i ^= iu[j-l][k];
542	ipp >>= 1;
543	}
544	iu[j][k]=i;
545	}
546	}
547	fac=1.0/(1L << MAXBIT);
548	in=0;
549
550	initDone=1;
551
552	sprintf(buf, "Initialization loop done\r");
553	XOPNotice(buf);
554
555	} else {
556	im=in;
557	for (j=1;j<=MAXBIT;j++) {
558	if (!(im & 1)) break;
559	im >>= 1;
560	}
561	if (j > MAXBIT) {
562	sprintf(buf, "MAXBIT too small in sobseq\r");
563	XOPNotice(buf);
564	}
565	im=(j-1)*MAXDIM;
566	for (k=1;k<=IMIN(n,MAXDIM);k++) {
567	ix[k] ^= iv[im+k];
568	x[k-1]=ix[k]*fac; /// this is a real array to send back, count this one from zero
569	//sprintf(buf, "calculate x[%d] = %g\r",k,ix[k]*fac);
570	//XOPNotice(buf);
571	}
572	in++;
573	}
574
575	// sprintf(buf, "x[0],x[1] = %g %g\r",x[0],x[1]);
576	// XOPNotice(buf);
577
578
579	p->result = 0;
580
581	return 0;
582	}
583
584	#undef MAXBIT
585	#undef MAXDIM
586
587
588	//#pragma XOP_RESET_STRUCT_PACKING // All structures are 2-byte-aligned.

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

Context Navigation

source: sans/XOP_Dev/MonteCarlo/DebyeSpheres.c @ 1249

Download in other formats: