Changeset 789 for sans/XOP_Dev/MonteCarlo/DebyeSpheres.c
 Timestamp:
 Feb 10, 2011 1:33:02 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/MonteCarlo/DebyeSpheres.c
r758 r789 7 7 #include "DebyeSpheres.h" 8 8 9 #pragma XOP_SET_STRUCT_PACKING // All structures are 2bytealigned.9 //#pragma XOP_SET_STRUCT_PACKING // All structures are 2bytealigned. 10 10 11 11 // Prototypes … … 158 158 159 159 160 #pragma XOP_RESET_STRUCT_PACKING // All structures are 2bytealigned. 161 // All structures are 2bytealigned. 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 int i,j; 173 int npt; 174 175 176 177 // check for all of the required waves 178 if (p>zwavH == NIL) { 179 SetNaN64(&p>result); 180 return NON_EXISTENT_WAVE; 181 } 182 if (p>ywavH == NIL) { 183 SetNaN64(&p>result); 184 return NON_EXISTENT_WAVE; 185 } 186 if (p>xwavH == NIL) { 187 SetNaN64(&p>result); 188 return NON_EXISTENT_WAVE; 189 } 190 191 //check to see that all are double 192 if(WaveType(p>zwavH) != NT_FP64 ) { 193 SetNaN64(&p>result); 194 return kExpectedNT_FP64; 195 } 196 if(WaveType(p>ywavH) != NT_FP64 ) { 197 SetNaN64(&p>result); 198 return kExpectedNT_FP64; 199 } 200 if(WaveType(p>xwavH) != NT_FP64 ) { 201 SetNaN64(&p>result); 202 return kExpectedNT_FP64; 203 } 204 205 // 206 npt = (int) WavePoints(p>xwavH); //wavePoints returns long, number of XYZ points 207 xv = WaveData(p>xwavH); //xyz locations 208 yv = WaveData(p>ywavH); 209 zv = WaveData(p>zwavH); 210 211 dmax = 0; 212 //do the i!=j double loop, keeping the maximum distance 213 214 for(i=0;i<npt;i+=1) { 215 for(j=(i+1);j<npt;j+=1) { 216 // dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]); 217 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]); 218 if(dij > dmax) { 219 dmax = dij; 220 } 221 } 222 } 223 224 p>result = dmax; 225 226 return 0; 227 228 } 229 230 /* 231 232 given the distances XYZ as a triplet (on a unit grid) 233 return the binned histogram of distances 234 235 */ 236 int 237 binDistanceX(BinParamPtr p) 238 { 239 double *xv,*yv,*zv,*bv; //pointers to input xyz coordinates 240 int i,j; 241 int npt,numBins,binIndex; 242 double grid,binWidth,val; 243 244 245 246 // check for all of the required waves 247 if (p>bwavH == NIL) { 248 SetNaN64(&p>result); 249 return NON_EXISTENT_WAVE; 250 } 251 if (p>zwavH == NIL) { 252 SetNaN64(&p>result); 253 return NON_EXISTENT_WAVE; 254 } 255 if (p>ywavH == NIL) { 256 SetNaN64(&p>result); 257 return NON_EXISTENT_WAVE; 258 } 259 if (p>xwavH == NIL) { 260 SetNaN64(&p>result); 261 return NON_EXISTENT_WAVE; 262 } 263 264 //check to see that all are double 265 if(WaveType(p>bwavH) != NT_FP64 ) { 266 SetNaN64(&p>result); 267 return kExpectedNT_FP64; 268 } 269 if(WaveType(p>zwavH) != NT_FP64 ) { 270 SetNaN64(&p>result); 271 return kExpectedNT_FP64; 272 } 273 if(WaveType(p>ywavH) != NT_FP64 ) { 274 SetNaN64(&p>result); 275 return kExpectedNT_FP64; 276 } 277 if(WaveType(p>xwavH) != NT_FP64 ) { 278 SetNaN64(&p>result); 279 return kExpectedNT_FP64; 280 } 281 282 // 283 npt = (int) WavePoints(p>xwavH); //wavePoints returns long, number of XYZ points 284 numBins = (int) WavePoints(p>bwavH); //wavePoints returns long, number of points in bin wave 285 286 xv = WaveData(p>xwavH); //xyz locations 287 yv = WaveData(p>ywavH); 288 zv = WaveData(p>zwavH); 289 bv = WaveData(p>bwavH); 290 291 grid = p>grid; 292 binWidth = p>binWidth; 293 294 //do the i!=j double loop, 295 for(i=0;i<npt;i+=1) { 296 for(j=(i+1);j<npt;j+=1) { 297 val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid; 298 binIndex = (int)(val/binWidth0.5); 299 if(binIndex > numBins 1 ) { 300 //Print "bad index" 301 } else { 302 bv[binIndex] += 1; 303 } 304 305 } 306 } 307 308 p>result = 0; 309 310 return 0; 311 312 } 313 314 315 /* 316 317 given the distances XYZ as a triplet (on a unit grid) and the SLD at each point, 318 return the binned histogram of distances for each of the parwise interactions 319 320 The returned binning is a matrix, and has to be assigned as such 321 322 */ 323 int 324 binSLDDistanceX(BinSLDParamPtr p) 325 { 326 double *xv,*yv,*zv; //pointers to input xyz coordinates 327 double *rho,*SLDLook,*PSFid; // rho and the SLD lookup vector 328 int i,j; 329 int npt,numBins,binIndex; 330 double grid,binWidth,val,retVal; 331 332 // for accessing the 2D wave data to write the results 333 waveHndl wavH,PSFwavH; 334 // long numDimensions; 335 // long dimensionSizes[MAX_DIMENSIONS+1]; 336 double value[2]; // Pointers used for double data. 337 long indices[MAX_DIMENSIONS]; 338 // 339 long rhoi,rhoj,rii,rji,PSFIndex; 340 341 342 // check for all of the required waves 343 if (p>PSFidH == NIL) { 344 SetNaN64(&p>result); 345 return NON_EXISTENT_WAVE; 346 } 347 if (p>SLDLookH == NIL) { 348 SetNaN64(&p>result); 349 return NON_EXISTENT_WAVE; 350 } 351 if (p>bwavH == NIL) { 352 SetNaN64(&p>result); 353 return NON_EXISTENT_WAVE; 354 } 355 if (p>rhowavH == NIL) { 356 SetNaN64(&p>result); 357 return NON_EXISTENT_WAVE; 358 } 359 if (p>zwavH == NIL) { 360 SetNaN64(&p>result); 361 return NON_EXISTENT_WAVE; 362 } 363 if (p>ywavH == NIL) { 364 SetNaN64(&p>result); 365 return NON_EXISTENT_WAVE; 366 } 367 if (p>xwavH == NIL) { 368 SetNaN64(&p>result); 369 return NON_EXISTENT_WAVE; 370 } 371 372 //check to see that all are double 373 if(WaveType(p>PSFidH) != NT_FP64 ) { 374 SetNaN64(&p>result); 375 return kExpectedNT_FP64; 376 } 377 if(WaveType(p>SLDLookH) != NT_FP64 ) { 378 SetNaN64(&p>result); 379 return kExpectedNT_FP64; 380 } 381 if(WaveType(p>bwavH) != NT_FP64 ) { 382 SetNaN64(&p>result); 383 return kExpectedNT_FP64; 384 } 385 if(WaveType(p>rhowavH) != NT_FP64 ) { 386 SetNaN64(&p>result); 387 return kExpectedNT_FP64; 388 } 389 if(WaveType(p>zwavH) != NT_FP64 ) { 390 SetNaN64(&p>result); 391 return kExpectedNT_FP64; 392 } 393 if(WaveType(p>ywavH) != NT_FP64 ) { 394 SetNaN64(&p>result); 395 return kExpectedNT_FP64; 396 } 397 if(WaveType(p>xwavH) != NT_FP64 ) { 398 SetNaN64(&p>result); 399 return kExpectedNT_FP64; 400 } 401 402 403 // access the 2D wave data for writing using the direct method 404 wavH = p>bwavH; 405 if (wavH == NIL) 406 return NOWAV; 407 // 408 PSFwavH = p>PSFidH; 409 410 npt = (int) WavePoints(p>xwavH); //wavePoints returns long, number of XYZ points 411 numBins = (int) WavePoints(p>bwavH); //wavePoints returns long, number of points in bin wave 412 413 xv = WaveData(p>xwavH); //xyz locations 414 yv = WaveData(p>ywavH); 415 zv = WaveData(p>zwavH); 416 rho = WaveData(p>rhowavH); 417 SLDLook = WaveData(p>SLDLookH); 418 PSFid = WaveData(p>PSFidH); //this one is 2D 419 420 grid = p>grid; 421 binWidth = p>binWidth; 422 423 //do the i!=j double loop, 424 for(i=0;i<npt;i+=1) { 425 for(j=(i+1);j<npt;j+=1) { 426 val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid; 427 binIndex = (int)(val/binWidth0.5); 428 if(binIndex > numBins 1 ) { 429 //Print "bad index" 430 } else { 431 rhoi = (long) rho[i]; //get the rho value at i and j 432 rhoj = (long) rho[j]; 433 rii = (long) SLDLook[rhoi]; //rho i index 434 rji = (long) SLDLook[rhoj]; //rho j index 435 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 436 indices[0] = rii; 437 indices[1] = rji; 438 if (retVal = MDGetNumericWavePointValue(PSFwavH, indices, value)) 439 return retVal; 440 //PSFIndex = (long) PSFid[rii][rji]; //doesn't work 441 PSFIndex = (long) value[0]; 442 443 //now do the assignment to the 2D 444 // equivalent to binMatrix[binIndex][PSFIndex] 445 446 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 447 indices[0] = binIndex; 448 indices[1] = PSFIndex; 449 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 450 return retVal; 451 value[0] += 1; // Real part 452 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 453 return retVal; 454 455 } 456 457 } 458 } 459 460 p>result = 0; 461 462 return 0; 463 464 } 465 466 467 ///// this is directly from Numerical Recipes 468 //  I did change the float to double, since Igor treats all as double 469 // and n is an int, not a pointer (seemed unnecessary) 470 // 471 #define MAXBIT 30 472 #define MAXDIM 6 473 static int iminarg1,iminarg2; 474 #define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2)) 475 476 int 477 SobolX(SobolParamPtr p) 478 { 479 int j,k,l; 480 unsigned long i,im,ipp; 481 static double fac; 482 static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1]; 483 static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4}; 484 static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4}; 485 static unsigned long iv[MAXDIM*MAXBIT+1]={ 486 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}; 487 488 static int initDone=0; 489 char buf[256]; 490 491 int n=0; 492 double *x; //output x vector 493 494 // check for all of the required waves 495 if (p>bwavH == NIL) { 496 SetNaN64(&p>result); 497 return NON_EXISTENT_WAVE; 498 } 499 500 //check to see that all are double 501 if(WaveType(p>bwavH) != NT_FP64 ) { 502 SetNaN64(&p>result); 503 return kExpectedNT_FP64; 504 } 505 x = WaveData(p>bwavH); 506 n = (int)(p>nIn); // not sure that the negative input will be properly cast to int 507 508 // sprintf(buf, "input, recast n = %g %d\r",p>nIn, n); 509 // XOPNotice(buf); 510 511 if (n < 0) { 512 513 if(initDone) { 514 sprintf(buf, "Don't reinitialize\r"); 515 XOPNotice(buf); 516 return 0; 517 } 518 519 for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k]; 520 for (k=1;k<=MAXDIM;k++) { 521 for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBITj); 522 for (j=mdeg[k]+1;j<=MAXBIT;j++) { 523 ipp=ip[k]; 524 i=iu[jmdeg[k]][k]; 525 i ^= (i >> mdeg[k]); 526 for (l=mdeg[k]1;l>=1;l) { 527 if (ipp & 1) i ^= iu[jl][k]; 528 ipp >>= 1; 529 } 530 iu[j][k]=i; 531 } 532 } 533 fac=1.0/(1L << MAXBIT); 534 in=0; 535 536 initDone=1; 537 538 sprintf(buf, "Initialization loop done\r"); 539 XOPNotice(buf); 540 541 } else { 542 im=in; 543 for (j=1;j<=MAXBIT;j++) { 544 if (!(im & 1)) break; 545 im >>= 1; 546 } 547 if (j > MAXBIT) { 548 sprintf(buf, "MAXBIT too small in sobseq\r"); 549 XOPNotice(buf); 550 } 551 im=(j1)*MAXDIM; 552 for (k=1;k<=IMIN(n,MAXDIM);k++) { 553 ix[k] ^= iv[im+k]; 554 x[k1]=ix[k]*fac; /// this is a real array to send back, count this one from zero 555 //sprintf(buf, "calculate x[%d] = %g\r",k,ix[k]*fac); 556 //XOPNotice(buf); 557 } 558 in++; 559 } 560 561 // sprintf(buf, "x[0],x[1] = %g %g\r",x[0],x[1]); 562 // XOPNotice(buf); 563 564 565 p>result = 0; 566 567 return 0; 568 } 569 570 #undef MAXBIT 571 #undef MAXDIM 572 573 574 //#pragma XOP_RESET_STRUCT_PACKING // All structures are 2bytealigned.
Note: See TracChangeset
for help on using the changeset viewer.