source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_ConnectedRods.ipf @ 798

Last change on this file since 798 was 798, checked in by srkline, 11 years ago

I'm not sure this is a great idea, but I'm putting the FFT / Debye sphere work that I have completed into SVN. It's all really rough - the math, I believe is correct, but the interface if really, really rough. But it's not going to develop without help.
File size: 10.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3// I should be able to calculate the average number of connections from a search of the
4// connXYZ table. I might have other counting routines elsewhere for this too.
5//
6//
7//              Set the Euler angles 55,-65,-150 for a consistent view
8//
9//
10
11Proc Overnight()
12        mTestConnRods()
13End
14
15
16Proc mTestConnRods()
17
18        String suffix = ""
19        Variable ranType = 1
20        Variable t1 = ticks
21       
22// constant for all steps
23        root:FFT_SolventSLD = 0
24
25// always start fresh when changing dimensions
26//      root:FFT_T = 40
27//      root:FFT_N = 256
28//      FFT_MakeMatrixButtonProc("")
29//      FFTEraseMatrixButtonProc("")
30       
31//      DoWindow/F Gizmo_VoxelMat
32       
33//      suffix = "a"
34//      CalculateTriangulated(40,40,20,"_CR"+suffix,ranType)
35//      DoUpdate/W=Gizmo_VoxelMat
36//      ExportGizmo wave as "p_connRod_"+suffix
37//     
38//      suffix = "b"
39//      CalculateTriangulated(100,40,20,"_CR"+suffix,ranType)
40//      DoUpdate/W=Gizmo_VoxelMat
41//      ExportGizmo wave as "p_connRod_"+suffix
42//     
43//      suffix = "c"
44//      CalculateTriangulated(200,40,20,"_CR"+suffix,ranType)
45//      DoUpdate/W=Gizmo_VoxelMat
46//      ExportGizmo wave as "p_connRod_"+suffix
47//     
48//      suffix = "d"
49//      CalculateTriangulated(300,40,20,"_CR"+suffix,ranType)
50//      DoUpdate/W=Gizmo_VoxelMat
51//      ExportGizmo wave as "p_connRod_"+suffix
52//     
53//      suffix = "e"
54//      CalculateTriangulated(400,40,20,"_CR"+suffix,ranType)
55//      DoUpdate/W=Gizmo_VoxelMat
56//      ExportGizmo wave as "p_connRod_"+suffix
57//     
58//      suffix = "f"
59//      CalculateTriangulated(500,40,20,"_CR"+suffix,ranType)
60//      DoUpdate/W=Gizmo_VoxelMat
61//      ExportGizmo wave as "p_connRod_"+suffix
62//     
63//     
64       
65// always start fresh when changing dimensions
66        root:FFT_T = 5
67        root:FFT_N = 256
68        FFT_MakeMatrixButtonProc("")
69        FFTEraseMatrixButtonProc("")
70       
71        suffix = "gg"
72        CalculateTriangulated(6,40,20,"_CR"+suffix,ranType)
73        DoUpdate/W=Gizmo_VoxelMat
74        ExportGizmo wave as "p_connRod_"+suffix
75       
76        suffix = "hh"
77        CalculateTriangulated(7,40,20,"_CR"+suffix,ranType)
78        DoUpdate/W=Gizmo_VoxelMat
79        ExportGizmo wave as "p_connRod_"+suffix
80       
81        suffix = "ii"
82        CalculateTriangulated(9,40,20,"_CR"+suffix,ranType)
83        DoUpdate/W=Gizmo_VoxelMat
84        ExportGizmo wave as "p_connRod_"+suffix
85       
86        suffix = "jj"
87        CalculateTriangulated(10,40,20,"_CR"+suffix,ranType)
88        DoUpdate/W=Gizmo_VoxelMat
89        ExportGizmo wave as "p_connRod_"+suffix
90       
91        suffix = "kk"
92        CalculateTriangulated(20,40,20,"_CR"+suffix,ranType)
93        DoUpdate/W=Gizmo_VoxelMat
94        ExportGizmo wave as "p_connRod_"+suffix
95       
96        suffix = "ll"
97        CalculateTriangulated(30,40,20,"_CR"+suffix,ranType)
98        DoUpdate/W=Gizmo_VoxelMat
99        ExportGizmo wave as "p_connRod_"+suffix
100       
101       
102        // then repeat everything with the random sequence
103        ranType = 0
104       
105        // always start fresh when changing dimensions
106//      root:FFT_T = 40
107//      root:FFT_N = 256
108//      FFT_MakeMatrixButtonProc("")
109//      FFTEraseMatrixButtonProc("")
110       
111//      DoWindow/F Gizmo_VoxelMat
112       
113//      suffix = "m"
114//      CalculateTriangulated(40,40,20,"_CR"+suffix,ranType)
115//      DoUpdate/W=Gizmo_VoxelMat
116//      ExportGizmo wave as "p_connRod_"+suffix
117//     
118//      suffix = "n"
119//      CalculateTriangulated(100,40,20,"_CR"+suffix,ranType)
120//      DoUpdate/W=Gizmo_VoxelMat
121//      ExportGizmo wave as "p_connRod_"+suffix
122//     
123//      suffix = "o"
124//      CalculateTriangulated(200,40,20,"_CR"+suffix,ranType)
125//      DoUpdate/W=Gizmo_VoxelMat
126//      ExportGizmo wave as "p_connRod_"+suffix
127//     
128//      suffix = "p"
129//      CalculateTriangulated(300,40,20,"_CR"+suffix,ranType)
130//      DoUpdate/W=Gizmo_VoxelMat
131//      ExportGizmo wave as "p_connRod_"+suffix
132//     
133//      suffix = "q"
134//      CalculateTriangulated(400,40,20,"_CR"+suffix,ranType)
135//      DoUpdate/W=Gizmo_VoxelMat
136//      ExportGizmo wave as "p_connRod_"+suffix
137//     
138//      suffix = "r"
139//      CalculateTriangulated(500,40,20,"_CR"+suffix,ranType)
140//      DoUpdate/W=Gizmo_VoxelMat
141//      ExportGizmo wave as "p_connRod_"+suffix
142       
143       
144       
145// always start fresh when changing dimensions
146        root:FFT_T = 5
147        root:FFT_N = 256
148        FFT_MakeMatrixButtonProc("")
149        FFTEraseMatrixButtonProc("")
150       
151        suffix = "ss"
152        CalculateTriangulated(6,40,20,"_CR"+suffix,ranType)
153        DoUpdate/W=Gizmo_VoxelMat
154        ExportGizmo wave as "p_connRod_"+suffix
155       
156        suffix = "tt"
157        CalculateTriangulated(7,40,20,"_CR"+suffix,ranType)
158        DoUpdate/W=Gizmo_VoxelMat
159        ExportGizmo wave as "p_connRod_"+suffix
160       
161        suffix = "uu"
162        CalculateTriangulated(9,40,20,"_CR"+suffix,ranType)
163        DoUpdate/W=Gizmo_VoxelMat
164        ExportGizmo wave as "p_connRod_"+suffix
165       
166        suffix = "vv"
167        CalculateTriangulated(10,40,20,"_CR"+suffix,ranType)
168        DoUpdate/W=Gizmo_VoxelMat
169        ExportGizmo wave as "p_connRod_"+suffix
170       
171        suffix = "ww"
172        CalculateTriangulated(20,40,20,"_CR"+suffix,ranType)
173        DoUpdate/W=Gizmo_VoxelMat
174        ExportGizmo wave as "p_connRod_"+suffix
175       
176        suffix = "xx"
177        CalculateTriangulated(30,40,20,"_CR"+suffix,ranType)
178        DoUpdate/W=Gizmo_VoxelMat
179        ExportGizmo wave as "p_connRod_"+suffix
180       
181       
182        Print "Total elapsed time (s) = ",(ticks-t1)/60.15
183        Print "Total elapsed time (min) = ",(ticks-t1)/60.15/60
184        Print "Total elapsed time (hr) = ",(ticks-t1)/60.15/60/60
185       
186
187end
188
189Function PrintWaveNote(kwStr)
190        String kwStr
191       
192        String alph="abcdefghijklmnopqrstuvwxyz"
193        String str=""
194        Variable ii
195       
196        for(ii=1;ii<12;ii+=1)
197                Wave w = $("root:iBin_"+alph[ii]+"p")
198                str = note(w)
199                if(strlen(kwStr) > 0)
200                        Print NumberByKey(kwStr, Str ,"=",";")
201                else
202                        Print str
203                endif
204        endfor
205       
206        return(0)
207End
208
209// npts is the number of nodes
210// diam is the diameter of the rods
211// nPass is the number of averaging passes
212// tagStr is the identifying tag for the files
213// ranType selects the RNG => 1 = Sobol, other = enoise
214Function CalculateTriangulated(nPts,diam,nPass,tagStr,ranType)
215        Variable nPts,diam,nPass
216        String tagStr
217        Variable ranType
218       
219        Variable ii,np,frac,nocc
220        Wave m=root:mat
221        NVAR grid=root:FFT_T
222        NVAR Nedge = root:FFT_N
223
224        np = 0
225        frac = 0
226        for(ii=0;ii<nPass;ii+=1)                //number of averaging passes
227                m=0
228               
229                if(ranType == 1)
230                        SobolFill3DMat(m,nPts)
231                else
232                        RandomFill3DMat(m,nPts)
233                endif
234                ParseMatrix3D_rho(m)                            // get the triplets of points to connect
235               
236                Print "connecting"
237                ConnectTriangulated(diam)
238               
239                Execute "DoFFT()"
240                Print "step ii=",ii
241                if(ii==0)
242                        Duplicate/O iBin, $("iBin"+tagStr),$("qBin"+tagStr),$("i2"+tagStr),$("eBin"+tagStr)
243                        wave ib=$("iBin"+tagStr)
244                        wave qb=$("qBin"+tagStr)
245                        wave i2=$("i2"+tagStr)
246                        wave eb=$("eBin"+tagStr)
247                        Wave iBin=iBin
248                        Wave qbin=qBin
249                        qb = qBin
250                        ib = iBin
251                        i2 = iBin*iBin
252                        eb = 0
253                else
254                        ib += iBin
255                        i2 += iBin*iBin
256                endif
257               
258                // get an average of the fill parameters, since it's somewhat random
259                nocc = NonZeroValues(m)
260                np += nocc
261                frac += nocc/DimSize(m,0)^3
262                Print "np,frac = ",np,frac
263        endfor
264
265        ib /= npass
266        i2 /= npass
267        if(nPass > 1)
268                eb = sqrt((i2-ib^2)/(npass-1))
269        else
270                eb = 1e-10
271        endif
272        np /= npass
273        frac /= npass
274       
275
276        String nStr
277        sprintf nstr,"T=%d;N=%d;Npass=%d;NNodes=%d;VolFrac=%g;",grid,Nedge,nPass,nPts,frac
278        Note ib, nStr
279
280        return(0)
281End
282
283
284
285
286
287// at this point, the matrix has been filled with points
288// and parsed to 3 xyz waves with ParseMatrix3D_rho()
289//
290Function ConnectTriangulated(diam)
291        Variable diam
292
293        Variable ii,num,thick,fill
294        Variable x1,x2,y1,y2,z1,z2
295        String testStr=""
296       
297        NVAR FFT_T = root:FFT_T
298        thick = diam/FFT_T / 2
299//      Print "diam, grid, thick = ",diam,FFT_T,thick
300       
301        Wave m=root:mat
302        Wave x3d=root:x3d
303        Wave y3d=root:y3d
304        Wave z3d=root:z3d
305       
306        ConvertXYZto3N(x3d,y3d,z3d)
307        Wave xyzTrip=root:matGiz
308       
309        Triangulate3d/OUT=2 xyzTrip
310        // M_TetraPath is generated
311        Wave M_TetraPath=M_TetraPath
312        M_TetraPath = round(M_TetraPath)
313       
314
315        num = DimSize(M_TetraPath,0)
316        fill = 1
317       
318        // make the list of connected points (to save time)
319        Make/O/T/N=0 connXYZ
320       
321        Print "num = ",num
322        for(ii=0;ii<num-1;ii+=1)
323                if(numtype(M_tetraPath[ii][0]) == 0 && numtype(M_tetraPath[ii+1][0]) == 0)              // 0 = a normal number
324                        x1 = M_TetraPath[ii][0]
325                        x2 = M_TetraPath[ii+1][0]
326                        y1 = M_TetraPath[ii][1]
327                        y2 = M_TetraPath[ii+1][1]
328                        z1 = M_TetraPath[ii][2]
329                        z2 = M_TetraPath[ii+1][2]
330                       
331                        testStr = XYZstr(x1,y1,z1,x2,y2,z2)
332                        if(!XYZ_are_Connected(connXYZ,testStr))         // are points already connected?
333                                // if not, connect them
334                                ConnectPoints3D(m, x1,y1,z1,x2,y2,z2,thick,fill)
335//                              print "connected"
336                                // and list them as connected
337                                AddXYZConnList(connXYZ,x1,y1,z1,x2,y2,z2)
338//                              Print "fraction done = ",ii/num
339                        endif
340                else
341                        //Print "ii = ",ii
342                endif
343        endfor
344
345        Print "done connecting"
346//      Execute "NumberOfPoints(\"mat\")"
347       
348        return(0)
349End
350
351// adds two entries to the list
352// pt1->pt2
353// pt2->pt1
354//
355Function AddXYZConnList(w,x1,y1,z1,x2,y2,z2)
356        Wave/T w
357        Variable x1,y1,z1,x2,y2,z2
358
359        String str1=XYZstr(x1,y1,z1,x2,y2,z2)
360        String str2=XYZstr(x2,y2,z2,x1,y1,z1)
361
362        Variable num=numpnts(w)
363        InsertPoints num,2,w
364        w[num] = str1
365        w[num+1] = str2
366       
367        return(0)
368End
369
370// wave w is a text wave
371//
372Function XYZ_are_Connected(w,str)
373        Wave/T w
374        String str
375       
376        Variable ans=0
377       
378        // do the check
379        FindValue/TEXT=str/TXOP=4 w
380       
381        //return the answer
382        if(V_Value == -1)
383                return(0)               //not connected
384        else
385                return(1)               //yes, connected
386        endif
387End
388
389Function/S XYZstr(x1,y1,z1,x2,y2,z2)
390        Variable x1,y1,z1,x2,y2,z2
391
392        String str=""
393        str += num_to_3char(x1)
394        str += num_to_3char(y1)
395        str += num_to_3char(z1)
396        str += num_to_3char(x2)
397        str += num_to_3char(y2)
398        str += num_to_3char(z2)
399       
400        return(str)
401End
402
403Function/S num_to_3char(num)
404        Variable num
405       
406        String numStr=""
407        if(num<10)
408                numStr = "00"+num2str(num)
409        else
410                if(num<100)
411                        numStr = "0"+num2str(num)
412                else
413                        numStr = num2str(num)
414                Endif
415        Endif
416        return(numstr)
417End
418
419
420/// this does not work ---
421//
422// the fCOnnectDots3D function does not work - and I
423// don't know what's wronk with it - so I switched to the
424// Triangulate3D operation
425// Jan 2011
426//
427//Function TestConnRods(nPts,nConn,nPass,tagStr)
428//      Variable nPts,nConn,nPass
429//      String tagStr
430//     
431//      Variable ii,np,frac,num
432//      Variable row,col,lay,xt,yt,zt,err
433//     
434//      Wave m=mat
435//      NVAR grid=root:FFT_T
436//      NVAR Nedge = root:FFT_N
437//      WAVE/Z x3d=root:x3d
438//     
439//      for(ii=0;ii<nPass;ii+=1)                //number of averaging passes
440//              //fill the spheres into the box
441//              m=0
442//              RandomFill3DMat(m,nPts)
443//              ParseMatrix3D_rho(m)                            // get the triplets of points to connect
444//             
445//              num=numpnts(x3d)
446//              Make/O/N=(num) numConnection3D=0
447//              Make/O/T/N=(num) connectedTo3D=""
448//              fConnectDots3D(nConn)
449//             
450//              Execute "DoFFT()"
451//              Print "step ii=",ii
452//              if(ii==0)
453//                      Duplicate/O iBin, $("iBin"+tagStr),$("qBin"+tagStr),$("i2"+tagStr),$("eBin"+tagStr)
454//                      wave ib=$("iBin"+tagStr)
455//                      wave qb=$("qBin"+tagStr)
456//                      wave i2=$("i2"+tagStr)
457//                      wave eb=$("eBin"+tagStr)
458//                      Wave iBin=iBin
459//                      Wave qbin=qBin
460//                      qb = qBin
461//                      ib = iBin
462//                      i2 = iBin*iBin
463//                      eb = 0
464//              else
465//                      ib += iBin
466//                      i2 += iBin*iBin
467//              endif
468//      endfor
469//
470//      ib /= npass
471//      i2 /= npass
472//      if(nPass > 1)
473//              eb = sqrt((i2-ib^2)/(npass-1))
474//      else
475//              eb = 1e-10
476//      endif
477//     
478//      np = NonZeroValues(m)
479//      frac = np/DimSize(m,0)^3
480//      String nStr
481//      sprintf nstr,"T=%d;N=%d;Npass=%d;NNodes=%d;NConn=%d;VolFrac=%g;",grid,Nedge,nPass,nPts,nConn,frac
482//      Note ib, nStr
483//
484//     
485//End
486
Note: See TracBrowser for help on using the repository browser.