1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | |
---|
3 | // VERY incomplete methods for filling volumes with connected cylinders |
---|
4 | // |
---|
5 | // most is still being worked out and verified in 2D |
---|
6 | // |
---|
7 | /// search for "crude core-shell cylinder" for an example of how this can be done |
---|
8 | // -- maybe the same filling, but no mat=0 (ask) |
---|
9 | |
---|
10 | |
---|
11 | |
---|
12 | //These macros may or may not work - designed to help automate the |
---|
13 | //process of calculating a series of configurations. |
---|
14 | |
---|
15 | |
---|
16 | Proc ConnectDots3D(w,np,maxNumConn,sobol,fill) |
---|
17 | String w="mat" |
---|
18 | Variable np=100,maxNumConn=2,sobol=1,fill=10 |
---|
19 | Prompt w,"matrix" |
---|
20 | Prompt np,"number of points" |
---|
21 | prompt maxNumConn,"number of connections per point" |
---|
22 | Prompt sobol,"1=Sobol, 0=random" |
---|
23 | Prompt fill,"fill SLD value" |
---|
24 | |
---|
25 | // $w=0 |
---|
26 | |
---|
27 | if(sobol) |
---|
28 | SobolFill3DMat($w,np,fill) |
---|
29 | else |
---|
30 | RandomFill3DMat($w,np,fill) |
---|
31 | endif |
---|
32 | |
---|
33 | ParseMatrix3DToXYZ(w) |
---|
34 | |
---|
35 | Variable num=numpnts(x3d) |
---|
36 | Make/O/N=(num) numConnection3D=0 |
---|
37 | Make/O/T/N=(num) connectedTo3D="" |
---|
38 | fConnectDots3D(maxNumConn,fill) |
---|
39 | end |
---|
40 | |
---|
41 | // connect the dots, no periodic boundary conditions |
---|
42 | //thickness is set to 1 |
---|
43 | Function fConnectDots3D(maxNumConn,fill) |
---|
44 | Variable maxNumConn,fill |
---|
45 | |
---|
46 | |
---|
47 | WAVE x3d=x3d |
---|
48 | WAVE y3d=y3d |
---|
49 | WAVE z3d=z3d |
---|
50 | WAVE numConnection3D = numConnection3D |
---|
51 | WAVE/T connectedTo3D = connectedTo3D |
---|
52 | WAVE mat=mat |
---|
53 | Variable num=numpnts(x3d),ii=0 |
---|
54 | Variable matSize=DimSize(mat,0),jj |
---|
55 | |
---|
56 | Variable nnX,nnY,nnZ,nnInd,nnDist //return values of the nearest neighbor |
---|
57 | Variable thick = 1 |
---|
58 | do |
---|
59 | for(ii=0;ii<num;ii+=1) |
---|
60 | if(numConnection3D[ii] < maxNumConn) |
---|
61 | FindNearestNeighbor3D(x3d,y3d,z3d,ii,maxNumConn,nnX,nnY,nnZ,nnInd,nnDist) |
---|
62 | |
---|
63 | numConnection3D[nnInd] += 1 |
---|
64 | numConnection3D[ii] += 1 |
---|
65 | connectedTo3D[nnInd] += num2str(ii)+"," |
---|
66 | connectedTo3D[ii] += num2str(nnInd)+"," |
---|
67 | |
---|
68 | ConnectPoints3D(mat, x3d[ii],y3d[ii],z3d[ii],x3d[nnInd],y3d[nnInd],z3d[nnInd],thick,fill) //always fill |
---|
69 | endif |
---|
70 | endfor |
---|
71 | //without rotation, a very closed network structure tends to form (bias?) |
---|
72 | // rotate (10?) points so that we're not starting from the same place every time |
---|
73 | // -- but the "connectedTo" information is LOST... |
---|
74 | //Rotate 10, x2d,y2d, numConnection,connectedTo |
---|
75 | //Print sum(numConnection,-inf,inf), num*maxNumConn |
---|
76 | while(sum(numConnection3D,-inf,inf) < num*maxNumConn) |
---|
77 | |
---|
78 | End |
---|
79 | |
---|
80 | Function FindNearestNeighbor3D(xw,yw,zw,startPt,maxConn,nnX,nnY,nnZ,nnInd,nnDist) |
---|
81 | Wave xw,yw,zw |
---|
82 | Variable startPt,maxConn,&nnX,&nnY,&nnZ,&nnInd,&nnDist |
---|
83 | |
---|
84 | Variable num,matSize,numToFind |
---|
85 | |
---|
86 | //make waves to hold the answers |
---|
87 | num=numpnts(xw) |
---|
88 | Make/O/N=(num) distWave3D,distIndex3D |
---|
89 | |
---|
90 | //now calculate the distances |
---|
91 | Variable ii,dist,testPt |
---|
92 | |
---|
93 | //this is probably the slowest step by far in the process...(tiny XOP?) |
---|
94 | for(ii=0;ii<num;ii+=1) |
---|
95 | distWave3D[ii] = (xw[ii]-xw[startPt])^2 + (yw[ii]-yw[startPt])^2 + (zw[ii]-zw[startPt])^2 //d^2 |
---|
96 | endfor |
---|
97 | |
---|
98 | MakeIndex distWave3D distIndex3D //then distWave[distIndex[ii]] will loop through in order of distance |
---|
99 | WAVE numConnection3D = numConnection3D |
---|
100 | WAVE/T connectedTo3D = connectedTo3D |
---|
101 | WAVE mat = mat |
---|
102 | |
---|
103 | for(ii=1;ii<num;ii+=1) //[0] point is the test point, dist == 0 |
---|
104 | testPt = distIndex3D[ii] //index |
---|
105 | if(numConnection3D[testPt] < maxConn ) //can test pt accept another connection? |
---|
106 | if(WhichListItem(num2str(testPt),connectedTo3D[startPt],",",0) == -1) // not already connected to the starting point? |
---|
107 | // Print ii,testPt |
---|
108 | // Printf "nearest point is (%d,%d), dist^2 = %g\r",xw[testPt],yw[testPt],distWave[testPt] |
---|
109 | nnX = xw[testPt] //return values as pass-by-reference |
---|
110 | nnY = yw[testPt] |
---|
111 | nnZ = zw[testPt] |
---|
112 | nnInd = testPt |
---|
113 | nnDist = distWave3D[testPt] |
---|
114 | |
---|
115 | return(0) //found a point, return |
---|
116 | endif |
---|
117 | endif |
---|
118 | endfor |
---|
119 | |
---|
120 | return (1) //did not find a point, return error |
---|
121 | End |
---|
122 | |
---|
123 | ///!!! fails if the points are along an axis (or even close) - then the box becomes infinitely thin |
---|
124 | Function ConnectPoints3D_old(w, x1,y1,z1,x2,y2,z2,thick,fill) |
---|
125 | Wave w |
---|
126 | Variable x1,y1,z1,x2,y2,z2,thick,fill |
---|
127 | |
---|
128 | //find all of the points within dist=thick of the line defined by p1 and p2 |
---|
129 | // and give these points (in w) the value 1 |
---|
130 | |
---|
131 | // search through the box defined by p1 and p2 |
---|
132 | Variable ii,jj,kk,dist |
---|
133 | |
---|
134 | for(ii=min(x1,x2);ii<=max(x1,x2);ii+=1) |
---|
135 | for(jj=min(y1,y2);jj<=max(y1,y2);jj+=1) |
---|
136 | for(kk=min(z1,z2);kk<=max(z1,z2);kk+=1) |
---|
137 | //find distance between current point and P1P2 |
---|
138 | dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2) |
---|
139 | //print dist |
---|
140 | if(dist<=thick) |
---|
141 | w[ii][jj][kk] = fill //add the point |
---|
142 | endif |
---|
143 | endfor |
---|
144 | endfor |
---|
145 | endfor |
---|
146 | |
---|
147 | End |
---|
148 | |
---|
149 | // |
---|
150 | // thickness is given in voxels |
---|
151 | // |
---|
152 | // the search box is expanded by the thickness in each direction (within the box) |
---|
153 | // to account for connections that are along an axis direction |
---|
154 | // |
---|
155 | Function ConnectPoints3D(w, x1,y1,z1,x2,y2,z2,thick,fill) |
---|
156 | Wave w |
---|
157 | Variable x1,y1,z1,x2,y2,z2,thick,fill |
---|
158 | |
---|
159 | //find all of the points within dist=thick of the line defined by p1 and p2 |
---|
160 | // and give these points (in w) the value 1 |
---|
161 | |
---|
162 | // search through the box defined by p1 and p2 |
---|
163 | Variable ii,jj,kk,dist,minX,minY,minZ,maxX,maxY,maxZ |
---|
164 | Variable boxMax,slop |
---|
165 | boxMax = DimSize(w,0) - 1 //max index of box |
---|
166 | slop = thick + 0 // to make the bounding box large enough the get it all |
---|
167 | |
---|
168 | minX=min(x1,x2) //give me the smaller x |
---|
169 | minX=max(minX-slop,0) //give me zero if it's negative |
---|
170 | minY=min(y1,y2) |
---|
171 | minY=max(minY-slop,0) |
---|
172 | minZ=min(z1,z2) |
---|
173 | minZ=max(minZ-slop,0) |
---|
174 | |
---|
175 | maxX=max(x1,x2) //give me the larger x |
---|
176 | maxX=min(maxX+slop,boxMax) //give the the box dimension if it's too big |
---|
177 | maxY=max(y1,y2) |
---|
178 | maxY=min(maxY+slop,boxMax) |
---|
179 | maxZ=max(z1,z2) |
---|
180 | maxZ=min(maxZ+slop,boxMax) |
---|
181 | |
---|
182 | for(ii=minX;ii<=maxX;ii+=1) |
---|
183 | for(jj=minY;jj<=maxY;jj+=1) |
---|
184 | for(kk=minZ;kk<=maxZ;kk+=1) |
---|
185 | //find distance between current point and P1P2 |
---|
186 | dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2) |
---|
187 | //print dist |
---|
188 | if(dist<=thick) |
---|
189 | w[ii][jj][kk] = fill //add the point |
---|
190 | endif |
---|
191 | endfor |
---|
192 | endfor |
---|
193 | endfor |
---|
194 | |
---|
195 | End |
---|
196 | |
---|
197 | //solution courtesy of the math forum @ Drexel "Ask Dr. Math" |
---|
198 | Function Pt_Line_Distance3D(xx,yy,zz, x1,y1,z1,x2,y2,z2) |
---|
199 | Variable xx,yy,zz, x1,y1,z1,x2,y2,z2 |
---|
200 | |
---|
201 | Variable dist,x3,y3,z3,x4,y4,z4,aa,ab |
---|
202 | |
---|
203 | //AB |
---|
204 | x3 = x2-x1 |
---|
205 | y3 = y2-y1 |
---|
206 | z3 = z2-z1 |
---|
207 | //AP |
---|
208 | x4 = xx-x1 |
---|
209 | y4 = yy-y1 |
---|
210 | z4 = zz-z1 |
---|
211 | |
---|
212 | //aa = length of ABxAP |
---|
213 | aa = sqrt( (y3*z4-z3*y4)^2 + (x3*z4-x4*z3)^2 + (x3*y4-x4*y3)^2 ) |
---|
214 | |
---|
215 | //ab = length of AB |
---|
216 | ab = sqrt(x3^2+y3^2+z3^2) |
---|
217 | |
---|
218 | dist = aa/ab |
---|
219 | |
---|
220 | return(dist) |
---|
221 | End |
---|
222 | |
---|
223 | ////////////////////////// |
---|
224 | //// random rod fills - just a voxel at each point, so use T=diameter |
---|
225 | |
---|
226 | // fill a volume with "rods" that are made up of a numer of spheres |
---|
227 | // lined up in a random direction |
---|
228 | // |
---|
229 | // all of the spheres/rods that are added must be connected = "adjacent" to |
---|
230 | // an occupied voxel. |
---|
231 | // |
---|
232 | // this implementation is a random walk, end-to-end connection |
---|
233 | // |
---|
234 | Function ConnectedRodFill(mat,len,num,periodic,fill) |
---|
235 | Wave mat |
---|
236 | variable len,num //length in direction, number of spheres to add |
---|
237 | Variable periodic //==1 if periodic boundaries |
---|
238 | Variable fill |
---|
239 | |
---|
240 | Variable nptx,npty,nptz,ii,fail=0,nAdd |
---|
241 | Make/O/I/N=3 stVec,dirVec,endVec |
---|
242 | |
---|
243 | nptx=DimSize(mat,0) |
---|
244 | npty=DimSize(mat,1) |
---|
245 | nptz=DimSize(mat,2) |
---|
246 | |
---|
247 | //place first rod randomly |
---|
248 | stVec[0]=trunc(abs(enoise(nptx))) //random integer distr betw (0,npt) |
---|
249 | stVec[1]=trunc(abs(enoise(npty))) |
---|
250 | stVec[2]=trunc(abs(enoise(nptz))) |
---|
251 | dirVec[0] = PickDirection() |
---|
252 | dirVec[1] = PickDirection() |
---|
253 | dirVec[2] = PickDirection() |
---|
254 | |
---|
255 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic,fill) |
---|
256 | |
---|
257 | stVec = endVec //end point is the new starting point |
---|
258 | |
---|
259 | //loop over remaining |
---|
260 | ii=nAdd |
---|
261 | do |
---|
262 | |
---|
263 | dirVec[0] = PickDirection() |
---|
264 | dirVec[1] = PickDirection() |
---|
265 | dirVec[2] = PickDirection() |
---|
266 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic,fill) |
---|
267 | stVec = endVec //end point is the new starting point |
---|
268 | |
---|
269 | ii+=nAdd |
---|
270 | |
---|
271 | if(mod(ii,100000) == 0) |
---|
272 | print "Point = ",ii |
---|
273 | endif |
---|
274 | |
---|
275 | while(ii<num) |
---|
276 | |
---|
277 | |
---|
278 | return(0) |
---|
279 | End |
---|
280 | |
---|
281 | // fill a volume with "rods" that are made up of a numer of spheres |
---|
282 | // lined up in a random direction |
---|
283 | // |
---|
284 | // |
---|
285 | Function UnConnectedRodFill(mat,len,num,periodic,fill) |
---|
286 | Wave mat |
---|
287 | variable len,num //length in direction, number of spheres to add |
---|
288 | Variable periodic //==1 if periodic boundaries |
---|
289 | Variable fill |
---|
290 | |
---|
291 | Variable nptx,npty,nptz,ii,fail=0,nAdd |
---|
292 | Make/O/I/N=3 stVec,dirVec,endVec |
---|
293 | |
---|
294 | nptx=DimSize(mat,0) |
---|
295 | npty=DimSize(mat,1) |
---|
296 | nptz=DimSize(mat,2) |
---|
297 | |
---|
298 | //place first rod randomly |
---|
299 | stVec[0]=trunc(abs(enoise(nptx))) //random integer distr betw (0,npt) |
---|
300 | stVec[1]=trunc(abs(enoise(npty))) |
---|
301 | stVec[2]=trunc(abs(enoise(nptz))) |
---|
302 | dirVec[0] = PickDirection() |
---|
303 | dirVec[1] = PickDirection() |
---|
304 | dirVec[2] = PickDirection() |
---|
305 | |
---|
306 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic,fill) |
---|
307 | //loop over remaining |
---|
308 | ii=nAdd |
---|
309 | do |
---|
310 | //pick a new point |
---|
311 | stVec[0]=trunc(abs(enoise(nptx))) //random integer distr betw (0,npt) |
---|
312 | stVec[1]=trunc(abs(enoise(npty))) |
---|
313 | stVec[2]=trunc(abs(enoise(nptz))) |
---|
314 | dirVec[0] = PickDirection() |
---|
315 | dirVec[1] = PickDirection() |
---|
316 | dirVec[2] = PickDirection() |
---|
317 | //add as many spheres as possible in that direction |
---|
318 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic,fill) |
---|
319 | ii+=nAdd |
---|
320 | //Print "point ",ii |
---|
321 | |
---|
322 | while(ii<num) |
---|
323 | Print "failures = ",fail |
---|
324 | |
---|
325 | // ParseMatrix3D_rho(mat) //convert to XYZ |
---|
326 | |
---|
327 | return(0) |
---|
328 | End |
---|
329 | |
---|
330 | |
---|
331 | Function PickDirection() |
---|
332 | return( trunc(abs(enoise(3))) - 1 ) //integer -1,0,1 |
---|
333 | End |
---|
334 | |
---|
335 | //need to add periodic boundary conditions (done) |
---|
336 | // |
---|
337 | // need to insist on connectivity |
---|
338 | // |
---|
339 | // |
---|
340 | Function AddSpheresInDirection(m,stVec,dirVec,endVec,nTry,periodic,fill) |
---|
341 | WAVE m,stVec,dirVec,endVec |
---|
342 | Variable nTry,periodic,fill |
---|
343 | |
---|
344 | Variable added,nx,ny,nz,nDimx,nDimy,nDimz |
---|
345 | |
---|
346 | NVAR solventSLD = root:FFT_SolventSLD |
---|
347 | |
---|
348 | nDimx=DimSize(m,0)-1 |
---|
349 | nDimy=DimSize(m,1)-1 |
---|
350 | nDimz=DimSize(m,2)-1 |
---|
351 | added=0 |
---|
352 | nx=stVec[0] |
---|
353 | ny=stVec[1] |
---|
354 | nz=stVec[2] |
---|
355 | |
---|
356 | do |
---|
357 | nx+=dirVec[0] //coordinates on next test point to try to add in d(xyz) direction |
---|
358 | ny+=dirVec[1] |
---|
359 | nz+=dirVec[2] |
---|
360 | //check for out-of bounds - if periodic ==1, reflect to opposite face of cube, else return |
---|
361 | if( (nx>nDimx) || (nx<0) ) |
---|
362 | if(periodic) |
---|
363 | nx = abs(nDimx+1 - abs(nx)) |
---|
364 | else |
---|
365 | break |
---|
366 | endif |
---|
367 | Endif |
---|
368 | if( (ny>nDimy) || (ny<0) ) |
---|
369 | if(periodic) |
---|
370 | ny = abs(nDimy+1 - abs(ny)) |
---|
371 | else |
---|
372 | break |
---|
373 | endif |
---|
374 | Endif |
---|
375 | if( (nz>nDimz) || (nz<0) ) |
---|
376 | if(periodic) |
---|
377 | nz = abs(nDimz+1 - abs(nz)) |
---|
378 | else |
---|
379 | break |
---|
380 | endif |
---|
381 | Endif |
---|
382 | //add the point or exit |
---|
383 | if(m[nx][ny][nz] != solventSLD) |
---|
384 | break //point is already occupied |
---|
385 | else |
---|
386 | m[nx][ny][nz] = fill |
---|
387 | // nx +=1 //we'll back this off before returning |
---|
388 | // ny +=1 |
---|
389 | // nz +=1 |
---|
390 | added += 1 |
---|
391 | endif |
---|
392 | while(added<nTry) |
---|
393 | |
---|
394 | endVec[0] = nx - dirVec[0] //back up to the last point actually added |
---|
395 | endVec[1] = ny - dirVec[1] |
---|
396 | endVec[2] = nz - dirVec[2] |
---|
397 | return added |
---|
398 | End |
---|
399 | |
---|
400 | Proc PutRandomCylindersAtPoints(w,num,len,periodic,fill) |
---|
401 | String w="mat" |
---|
402 | Variable num=1000,len=10,periodic=1,fill = 10 |
---|
403 | Prompt w,"matrix" |
---|
404 | Prompt num,"number of starting points" |
---|
405 | prompt len,"length of cylinders" |
---|
406 | prompt periodic"1=periodic, 0=non-periodic boundaries" |
---|
407 | Prompt fill,"fill SLD value" |
---|
408 | |
---|
409 | |
---|
410 | // $w=0 |
---|
411 | RandomFill3DMat($w,num,fill) |
---|
412 | CylindersAtPoints($w,len,periodic,fill) |
---|
413 | |
---|
414 | NumberOfPoints() |
---|
415 | end |
---|
416 | |
---|
417 | Proc PutXAxisCylindersAtPoints(w,num,rad,len,periodic,sobol,fill) |
---|
418 | String w="mat" |
---|
419 | Variable num=100,rad=20,len=300,periodic=1,sobol=1,fill=10 |
---|
420 | Prompt w,"matrix" |
---|
421 | Prompt num,"number of starting points" |
---|
422 | prompt rad,"radius of cylinders" |
---|
423 | prompt len,"length of cylinders" |
---|
424 | prompt periodic,"1=periodic, 0=non-periodic boundaries" |
---|
425 | Prompt sobol,"1=Sobol, 0=random" |
---|
426 | Prompt fill,"fill SLD value" |
---|
427 | |
---|
428 | |
---|
429 | // $w=0 |
---|
430 | X_CylindersAtPoints($w,num,rad,len,sobol,periodic,fill) |
---|
431 | |
---|
432 | NumberOfPoints() |
---|
433 | end |
---|
434 | |
---|
435 | Proc PutXAxisCylindersSquare(w,rad,len,sep,fill) |
---|
436 | String w="mat" |
---|
437 | Variable rad=20,len=300,sep=80,fill=10 |
---|
438 | Prompt w,"matrix" |
---|
439 | prompt rad,"radius of cylinders" |
---|
440 | prompt len,"length of cylinders" |
---|
441 | prompt sep,"center-to-center separation of cylinders" |
---|
442 | Prompt fill,"fill SLD value" |
---|
443 | |
---|
444 | |
---|
445 | // $w=0 |
---|
446 | X_CylindersSquareGrid($w,rad,len,sep,fill) |
---|
447 | |
---|
448 | NumberOfPoints() |
---|
449 | end |
---|
450 | |
---|
451 | Proc PutXAxisCylindersHexagonal(w,rad,len,sep,fill) |
---|
452 | String w="mat" |
---|
453 | Variable rad=20,len=300,sep=80,fill=10 |
---|
454 | Prompt w,"matrix" |
---|
455 | prompt rad,"radius of cylinders" |
---|
456 | prompt len,"length of cylinders" |
---|
457 | prompt sep,"center-to-center separation of cylinders" |
---|
458 | Prompt fill,"fill SLD value" |
---|
459 | |
---|
460 | // $w=0 |
---|
461 | X_CylindersHexagonalGrid($w,rad,len,sep,fill) |
---|
462 | |
---|
463 | NumberOfPoints() |
---|
464 | end |
---|
465 | |
---|
466 | Function CylindersAtPoints(mat,len,periodic,fill) |
---|
467 | Wave mat |
---|
468 | variable len //length in direction |
---|
469 | Variable periodic //==1 if periodic boundaries |
---|
470 | Variable fill |
---|
471 | |
---|
472 | Make/O/B/N=3 stVec,dirVec,endVec |
---|
473 | |
---|
474 | ParseMatrix3D_rho(mat) |
---|
475 | Wave x3d=x3d |
---|
476 | Wave y3d=y3d |
---|
477 | Wave z3d=z3d |
---|
478 | |
---|
479 | Variable ii=0,num=numpnts(x3d) |
---|
480 | NVAR grid=root:FFT_T |
---|
481 | |
---|
482 | Variable nAdd |
---|
483 | |
---|
484 | for(ii=0;ii<num;ii+=1) |
---|
485 | dirVec[0] = PickDirection() |
---|
486 | dirVec[1] = PickDirection() |
---|
487 | dirVec[2] = PickDirection() |
---|
488 | stVec[0] = x3d[ii] |
---|
489 | stVec[1] = y3d[ii] |
---|
490 | stVec[2] = z3d[ii] |
---|
491 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic,fill) |
---|
492 | //FillSphereRadius(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],1) |
---|
493 | endfor |
---|
494 | |
---|
495 | return(0) |
---|
496 | End |
---|
497 | |
---|
498 | /////////////// |
---|
499 | |
---|
500 | Function X_CylindersAtPoints(mat,num,rad,len,sobol,periodic,fill) |
---|
501 | Wave mat |
---|
502 | variable num,rad,len //length in direction |
---|
503 | Variable periodic //==1 if periodic boundaries |
---|
504 | Variable sobol //==1 if sobol selection of points (2D) |
---|
505 | Variable fill |
---|
506 | |
---|
507 | NVAR solventSLD = root:FFT_SolventSLD |
---|
508 | |
---|
509 | Variable np |
---|
510 | np = DimSize(mat,0) // assumes that all dimensions are the same |
---|
511 | |
---|
512 | // fill a 2D plane with points |
---|
513 | Make/O/B/N=(np,np) plane |
---|
514 | plane = solventSLD |
---|
515 | |
---|
516 | if(sobol) |
---|
517 | if(WaveExists($"SobolInit") == 0) |
---|
518 | Make/O/D/N=2 SobolInit |
---|
519 | SobolX(-1,SobolInit) |
---|
520 | else |
---|
521 | SobolPoints2D(plane,num,fill) |
---|
522 | endif |
---|
523 | else |
---|
524 | RandomPoints2D(plane,num,fill) |
---|
525 | endif |
---|
526 | // put it in the proper plane of the matrix |
---|
527 | mat[np/2][][] = plane[q][r] // in the YZ plane |
---|
528 | |
---|
529 | ParseMatrix3D_rho(mat) |
---|
530 | Wave x3d=x3d |
---|
531 | Wave y3d=y3d |
---|
532 | Wave z3d=z3d |
---|
533 | |
---|
534 | Variable ii=0 |
---|
535 | NVAR grid=root:FFT_T |
---|
536 | |
---|
537 | for(ii=0;ii<num;ii+=1) |
---|
538 | FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,fill) //cylinder 1 |
---|
539 | endfor |
---|
540 | |
---|
541 | return(0) |
---|
542 | End |
---|
543 | |
---|
544 | |
---|
545 | Function X_CylindersSquareGrid(mat,rad,len,sep,fill) |
---|
546 | Wave mat |
---|
547 | variable rad,len //length of cylinders |
---|
548 | Variable sep // EDGE separation, in same units as cylinder |
---|
549 | Variable fill |
---|
550 | |
---|
551 | NVAR solventSLD = root:FFT_SolventSLD |
---|
552 | |
---|
553 | NVAR grid=root:FFT_T |
---|
554 | Variable np,spacing |
---|
555 | np = DimSize(mat,0) // assumes that all dimensions are the same |
---|
556 | |
---|
557 | // fill a 2D plane with points |
---|
558 | Make/O/B/N=(np,np) plane |
---|
559 | plane = solventSLD |
---|
560 | |
---|
561 | |
---|
562 | spacing = round(sep/grid) // so it's an integer |
---|
563 | FillPlaneSquareGrid(plane,spacing,fill) |
---|
564 | |
---|
565 | |
---|
566 | |
---|
567 | // put it in the proper plane of the matrix |
---|
568 | mat[np/2][][] = plane[q][r] // in the YZ plane |
---|
569 | |
---|
570 | |
---|
571 | ParseMatrix3D_rho(mat) |
---|
572 | Wave x3d=x3d |
---|
573 | Wave y3d=y3d |
---|
574 | Wave z3d=z3d |
---|
575 | |
---|
576 | Variable ii=0,num |
---|
577 | num = numpnts(x3d) |
---|
578 | |
---|
579 | for(ii=0;ii<num;ii+=1) |
---|
580 | FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,fill) //cylinder 1 |
---|
581 | endfor |
---|
582 | |
---|
583 | return(0) |
---|
584 | End |
---|
585 | |
---|
586 | Function X_CylindersHexagonalGrid(mat,rad,len,sep,fill) |
---|
587 | Wave mat |
---|
588 | variable rad,len //length of cylinders |
---|
589 | Variable sep // EDGE separation, in same units as cylinder |
---|
590 | Variable fill |
---|
591 | |
---|
592 | NVAR solventSLD = root:FFT_SolventSLD |
---|
593 | |
---|
594 | NVAR grid=root:FFT_T |
---|
595 | Variable np,spacing |
---|
596 | np = DimSize(mat,0) // assumes that all dimensions are the same |
---|
597 | |
---|
598 | // fill a 2D plane with points |
---|
599 | Make/O/B/N=(np,np) plane |
---|
600 | plane = solventSLD |
---|
601 | |
---|
602 | spacing = round(sep/grid) // so it's an integer |
---|
603 | FillPlaneHexagonal(plane,spacing,fill) |
---|
604 | |
---|
605 | // put it in the proper plane of the matrix |
---|
606 | mat[np/2][][] = plane[q][r] // in the YZ plane |
---|
607 | |
---|
608 | ParseMatrix3D_rho(mat) |
---|
609 | Wave x3d=x3d |
---|
610 | Wave y3d=y3d |
---|
611 | Wave z3d=z3d |
---|
612 | |
---|
613 | Variable ii=0,num |
---|
614 | num = numpnts(x3d) |
---|
615 | |
---|
616 | for(ii=0;ii<num;ii+=1) |
---|
617 | FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,fill) //cylinder 1 |
---|
618 | Print "cyl = ",ii,num |
---|
619 | endfor |
---|
620 | |
---|
621 | //// makes a crude core-shell cylinder |
---|
622 | // for(ii=0;ii<num;ii+=1) |
---|
623 | // FillXCylinder(mat,grid,rad-10,x3d[ii],y3d[ii],z3d[ii],len,3*fill) //cylinder 1 |
---|
624 | // endfor |
---|
625 | // |
---|
626 | return(0) |
---|
627 | End |
---|
628 | |
---|
629 | |
---|
630 | |
---|
631 | Function FillPlaneSquareGrid(w,sp,fill) |
---|
632 | Wave w |
---|
633 | Variable sp //spacing in pixels |
---|
634 | Variable fill |
---|
635 | |
---|
636 | Variable np,ii,jj |
---|
637 | |
---|
638 | np = DimSize(w,0) |
---|
639 | |
---|
640 | for(ii=sp;ii<(np-sp/2);ii+=sp) |
---|
641 | for(jj=sp;jj<(np-sp/2);jj+=sp) |
---|
642 | w[ii][jj] = fill |
---|
643 | endfor |
---|
644 | endfor |
---|
645 | |
---|
646 | |
---|
647 | return(0) |
---|
648 | End |
---|
649 | |
---|
650 | // sp is the |
---|
651 | Function FillPlaneHexagonal(w,sp,fill) |
---|
652 | Wave w |
---|
653 | Variable sp //spacing in pixels |
---|
654 | Variable fill |
---|
655 | |
---|
656 | Variable np,ii,jj,xx,yy |
---|
657 | |
---|
658 | np = DimSize(w,0) |
---|
659 | |
---|
660 | |
---|
661 | // do even columns, then odds |
---|
662 | for(ii=sp;ii<(np-sp/2);ii+=(sp*sqrt(3)) ) |
---|
663 | for(jj=sp;jj<(np-sp/2);jj+=sp) |
---|
664 | w[ii][jj] = fill |
---|
665 | endfor |
---|
666 | endfor |
---|
667 | |
---|
668 | // now odds |
---|
669 | for(ii=(sp+sp*sqrt(3)/2);ii<(np-sp/2);ii+=(sp*sqrt(3))) |
---|
670 | for(jj=sp/2;jj<(np-sp/2);jj+=sp) |
---|
671 | w[ii][jj] = fill |
---|
672 | endfor |
---|
673 | endfor |
---|
674 | |
---|
675 | return(0) |
---|
676 | End |
---|
677 | |
---|
678 | Proc HexagonalCylinders() |
---|
679 | |
---|
680 | Variable xc,yc,zc,len,rad,aa,cosa,sina,meanRad,pd,meanLen |
---|
681 | |
---|
682 | Variable grid=root:FFT_T |
---|
683 | Variable npt=root:FFT_N |
---|
684 | |
---|
685 | xc=npt/2 |
---|
686 | yc=npt/2 |
---|
687 | zc=npt/2 |
---|
688 | |
---|
689 | meanRad=20 |
---|
690 | pd=0.01 |
---|
691 | meanLen=400 |
---|
692 | len=meanLen |
---|
693 | |
---|
694 | // aa=meanRad/grid*2 +40 //grid units, +2 gives 2 units of separation |
---|
695 | aa=meanRad/grid*2 +2 //grid units, +2 gives 2 units of separation |
---|
696 | cosa = round(0.5*aa) |
---|
697 | sina = round(0.866*aa) |
---|
698 | |
---|
699 | |
---|
700 | Variable ii,nPass,x1,x2,x3 |
---|
701 | ii=0 |
---|
702 | nPass=1 |
---|
703 | |
---|
704 | |
---|
705 | x1 = Gauss(meanRad-pd*meanRad,meanRad,pd*meanRad) |
---|
706 | x2 = Gauss(meanRad,meanRad,pd*meanRad) |
---|
707 | x3 = Gauss(meanRad+pd*meanRad,meanRad,pd*meanRad) |
---|
708 | |
---|
709 | Print x1,x2,x3 |
---|
710 | // |
---|
711 | /////////// average over 3 point distribution -- |
---|
712 | // //calculate and sum the three contributions |
---|
713 | // FastOp mat=0 |
---|
714 | // rad = meanRad-pd*meanRad |
---|
715 | //// len = meanLen + gnoise(pd*meanLen) |
---|
716 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
717 | // Execute "DoFFT()" |
---|
718 | // Duplicate/O iBin iBin_1g_128 |
---|
719 | // iBin_1g_128 *= x1 |
---|
720 | // |
---|
721 | // FastOp mat=0 |
---|
722 | // rad = meanRad |
---|
723 | //// len = meanLen + gnoise(pd*meanLen) |
---|
724 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
725 | // Execute "DoFFT()" |
---|
726 | // iBin_1g_128 += x2*iBin |
---|
727 | // |
---|
728 | // FastOp mat=0 |
---|
729 | // rad = meanRad+pd*meanRad |
---|
730 | //// len = meanLen + gnoise(pd*meanLen) |
---|
731 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
732 | // Execute "DoFFT()" |
---|
733 | // iBin_1g_128 += x3*iBin |
---|
734 | // |
---|
735 | // iBin_1g_128 /= (x1+x2+x3) //renormalize |
---|
736 | // |
---|
737 | ////// end 3 pt average |
---|
738 | |
---|
739 | ////////////////////// |
---|
740 | // paired cylinders |
---|
741 | ///////////////////////// |
---|
742 | // //calculate and sum the three contributions |
---|
743 | // FastOp mat=0 |
---|
744 | // rad = meanRad-pd*meanRad |
---|
745 | //// len = meanLen + gnoise(pd*meanLen) |
---|
746 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
747 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
748 | // Execute "DoFFT()" |
---|
749 | // Duplicate/O iBin iBin_2g_256_sep40 |
---|
750 | // iBin_2g_256_sep40 *= x1 |
---|
751 | // |
---|
752 | // FastOp mat=0 |
---|
753 | // rad = meanRad |
---|
754 | //// len = meanLen + gnoise(pd*meanLen) |
---|
755 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
756 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
757 | // Execute "DoFFT()" |
---|
758 | // iBin_2g_256_sep40 += x2*iBin |
---|
759 | // |
---|
760 | // FastOp mat=0 |
---|
761 | // rad = meanRad+pd*meanRad |
---|
762 | //// len = meanLen + gnoise(pd*meanLen) |
---|
763 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
764 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
765 | // Execute "DoFFT()" |
---|
766 | // iBin_2g_256_sep40 += x3*iBin |
---|
767 | // |
---|
768 | // iBin_2g_256_sep40 /= (x1+x2+x3) //renormalize |
---|
769 | |
---|
770 | ///////////////////////// |
---|
771 | // |
---|
772 | // |
---|
773 | //// |
---|
774 | // ii = 0 |
---|
775 | // do |
---|
776 | // FastOp mat=0 |
---|
777 | // rad = meanRad + gnoise(pd*meanRad) |
---|
778 | // print ii,rad |
---|
779 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
780 | // Execute "DoFFT()" |
---|
781 | // if(ii==0) |
---|
782 | // Duplicate/O iBin iBin_1g_256 |
---|
783 | // else |
---|
784 | // iBin_1g_256 += iBin |
---|
785 | // endif |
---|
786 | // ii+=1 |
---|
787 | // while(ii<nPass) |
---|
788 | // iBin_1g_256 /= nPass |
---|
789 | // |
---|
790 | ///// |
---|
791 | // ii=0 |
---|
792 | // nPass=1 |
---|
793 | // do |
---|
794 | // FastOp mat=0 |
---|
795 | // rad = meanRad + gnoise(pd*meanRad) |
---|
796 | // print ii,rad |
---|
797 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) //cylinder 1 |
---|
798 | // rad = meanRad + gnoise(pd*meanRad) |
---|
799 | // print ii,rad |
---|
800 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
801 | // Execute "DoFFT()" |
---|
802 | // if(ii==0) |
---|
803 | // Duplicate/O iBin iBin_2g_256 |
---|
804 | // else |
---|
805 | // iBin_2g_256 += iBin |
---|
806 | // endif |
---|
807 | // ii+=1 |
---|
808 | // while(ii<nPass) |
---|
809 | // iBin_2g_256 /= nPass |
---|
810 | // |
---|
811 | ii=0 |
---|
812 | npass=1 |
---|
813 | do |
---|
814 | FastOp mat=0 |
---|
815 | rad = meanRad + gnoise(pd*meanRad) |
---|
816 | print ii,rad |
---|
817 | FillZCylinder(mat,grid,rad,xc,yc,zc,len,10) //cylinder 1 |
---|
818 | rad = meanRad + gnoise(pd*meanRad) |
---|
819 | print ii,rad |
---|
820 | FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,10) //cylinder 2 |
---|
821 | rad = meanRad + gnoise(pd*meanRad) |
---|
822 | print ii,rad |
---|
823 | FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,10) //cylinder 3 |
---|
824 | Execute "DoFFT()" |
---|
825 | if(ii==0) |
---|
826 | Duplicate/O iBin iBin_3g_256 |
---|
827 | else |
---|
828 | iBin_3g_256 += iBin |
---|
829 | endif |
---|
830 | ii+=1 |
---|
831 | while(ii<nPass) |
---|
832 | iBin_3g_256 /= nPass |
---|
833 | // |
---|
834 | rad = meanRad + gnoise(pd*meanRad) |
---|
835 | print rad |
---|
836 | FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,10) |
---|
837 | Execute "DoFFT()" |
---|
838 | Duplicate/O iBin iBin_2 |
---|
839 | /// |
---|
840 | rad = meanRad + gnoise(pd*meanRad) |
---|
841 | print rad |
---|
842 | FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,10) |
---|
843 | Execute "DoFFT()" |
---|
844 | Duplicate/O iBin iBin_3 |
---|
845 | /// |
---|
846 | rad = meanRad + gnoise(pd*meanRad) |
---|
847 | print rad |
---|
848 | FillZCylinder(mat,grid,rad,xc+sina,yc-cosa,zc,len,10) |
---|
849 | Execute "DoFFT()" |
---|
850 | Duplicate/O iBin iBin_4 |
---|
851 | /// |
---|
852 | rad = meanRad + gnoise(pd*meanRad) |
---|
853 | print rad |
---|
854 | FillZCylinder(mat,grid,rad,xc,yc-aa,zc,len,10) |
---|
855 | Execute "DoFFT()" |
---|
856 | Duplicate/O iBin iBin_5 |
---|
857 | /// |
---|
858 | rad = meanRad + gnoise(pd*meanRad) |
---|
859 | print rad |
---|
860 | FillZCylinder(mat,grid,rad,xc-sina,yc-cosa,zc,len,10) |
---|
861 | Execute "DoFFT()" |
---|
862 | Duplicate/O iBin iBin_6 |
---|
863 | /// |
---|
864 | rad = meanRad + gnoise(pd*meanRad) |
---|
865 | print rad |
---|
866 | FillZCylinder(mat,grid,rad,xc-sina,yc+cosa,zc,len,10) |
---|
867 | Execute "DoFFT()" |
---|
868 | Duplicate/O iBin iBin_7 |
---|
869 | |
---|
870 | end |
---|
871 | |
---|
872 | |
---|
873 | Proc Vary_N_In_Direction() |
---|
874 | |
---|
875 | Variable fill=10 |
---|
876 | |
---|
877 | FFTEraseMatrixButtonProc("") |
---|
878 | ConnectedRodFill(mat,2,30000*4,1,fill) |
---|
879 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
880 | Duplicate/O ival_XOP ival_XOP_2 |
---|
881 | Duplicate/O qval_XOP qval_XOP_2 |
---|
882 | |
---|
883 | FFTEraseMatrixButtonProc("") |
---|
884 | ConnectedRodFill(mat,3,30000*4,1,fill) |
---|
885 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
886 | Duplicate/O ival_XOP ival_XOP_3 |
---|
887 | Duplicate/O qval_XOP qval_XOP_3 |
---|
888 | |
---|
889 | FFTEraseMatrixButtonProc("") |
---|
890 | ConnectedRodFill(mat,4,30000*4,1,fill) |
---|
891 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
892 | Duplicate/O ival_XOP ival_XOP_4 |
---|
893 | Duplicate/O qval_XOP qval_XOP_4 |
---|
894 | |
---|
895 | FFTEraseMatrixButtonProc("") |
---|
896 | ConnectedRodFill(mat,5,30000*4,1,fill) |
---|
897 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
898 | Duplicate/O ival_XOP ival_XOP_5 |
---|
899 | Duplicate/O qval_XOP qval_XOP_5 |
---|
900 | |
---|
901 | FFTEraseMatrixButtonProc("") |
---|
902 | ConnectedRodFill(mat,6,30000*4,1,fill) |
---|
903 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
904 | Duplicate/O ival_XOP ival_XOP_6 |
---|
905 | Duplicate/O qval_XOP qval_XOP_6 |
---|
906 | |
---|
907 | FFTEraseMatrixButtonProc("") |
---|
908 | ConnectedRodFill(mat,7,30000*4,1,fill) |
---|
909 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
910 | Duplicate/O ival_XOP ival_XOP_7 |
---|
911 | Duplicate/O qval_XOP qval_XOP_7 |
---|
912 | |
---|
913 | FFTEraseMatrixButtonProc("") |
---|
914 | ConnectedRodFill(mat,8,30000*4,1,fill) |
---|
915 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
916 | Duplicate/O ival_XOP ival_XOP_8 |
---|
917 | Duplicate/O qval_XOP qval_XOP_8 |
---|
918 | |
---|
919 | FFTEraseMatrixButtonProc("") |
---|
920 | ConnectedRodFill(mat,9,30000*4,1,fill) |
---|
921 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
922 | Duplicate/O ival_XOP ival_XOP_9 |
---|
923 | Duplicate/O qval_XOP qval_XOP_9 |
---|
924 | |
---|
925 | FFTEraseMatrixButtonProc("") |
---|
926 | ConnectedRodFill(mat,10,30000*4,1,fill) |
---|
927 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
928 | Duplicate/O ival_XOP ival_XOP_10 |
---|
929 | Duplicate/O qval_XOP qval_XOP_10 |
---|
930 | |
---|
931 | End |
---|