source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_SuperFormula.ipf @ 1089

Last change on this file since 1089 was 1089, checked in by srkline, 5 years ago

adding superformula calculation to the FFT method. SuperEllipsoid?, SuperToroid?, and SuperQuartic? are also included. From Panel input, draws the shape and places in mat ready for FFT or Debye.

File size: 26.3 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3/// http://en.wikipedia.org/wiki/Superformula
4//
5// and the SuperEllipsoid
6// and the superquadric
7// and the supertoroid
8//
9//
10// The superEllipsoid is in an implicit form, so it can easily be converted
11// to voxels, since the implicit form defines inside/outside of the surface.
12//
13// -- need to clean up the superEllipsoid version so that it can be incorporated into the real-space
14//  modeling functions - since this would be a rather unique, if not totally useless thing to do.
15//
16//
17// the general superformula is given in polar coordinates, so it is going to require a
18// bit of math to get the voxel representation.
19//
20// With an implicit equation for the surface - it's a snap to generate the voxels.
21//
22//
23
24
25
26// you can also use:
27//
28// mat_as_3dCloud()
29//
30// Gizmo_superSurface()
31//
32
33
34
35Macro Setup_SuperFormulas()
36
37        DoWindow/F SuperFormulaPanel
38        if(V_flag == 0)
39                Setup_super3D_waves(128)
40       
41                Variable/G gSuperRadioVal= 4            //start with superFormula checked
42
43                SuperFormulaPanel()
44        endif
45               
46End
47
48Proc SuperFormulaPanel()
49        PauseUpdate; Silent 1           // building window...
50        NewPanel /W=(634,45,934,394) /K=1
51       
52       
53        DoWindow/C SuperFormulaPanel
54//      ShowTools/A
55        SetDrawLayer UserBack
56        SetVariable setvar0,pos={37.00,92.00},size={80.00,18.00},title="r",fSize=12
57        SetVariable setvar0,limits={-20,20,1},value= _NUM:3,disable=2
58        SetVariable setvar1,pos={37.00,118.00},size={80.00,18.00},title="t",fSize=12
59        SetVariable setvar1,limits={-20,20,1},value= _NUM:2,disable=2
60        SetVariable setvar2,pos={37.00,145.00},size={80.00,18.00},title="s",fSize=12
61        SetVariable setvar2,limits={-20,20,1},value= _NUM:4,disable=2
62        SetVariable setvar3,pos={37.00,172.00},size={80.00,18.00},title="rad",fSize=12
63        SetVariable setvar3,limits={0,50,1},value= _NUM:15,disable=2
64       
65        SetVariable setvar4,pos={162.00,91.00},size={80.00,18.00},title="m",fSize=12
66        SetVariable setvar4,limits={-20,20,1},value= _NUM:15
67        SetVariable setvar5,pos={162.00,117.00},size={80.00,18.00},title="n1",fSize=12
68        SetVariable setvar5,limits={-20,20,1},value= _NUM:1     
69        SetVariable setvar6,pos={162.00,143.00},size={80.00,18.00},title="n2",fSize=12
70        SetVariable setvar6,limits={-20,20,1},value= _NUM:2
71        SetVariable setvar7,pos={162.00,170.00},size={80.00,18.00},title="n3",fSize=12
72        SetVariable setvar7,limits={-20,20,1},value= _NUM:6
73
74        SetVariable setvar8,pos={35,226.00},size={120.00,18.00},title="x-scaling"
75        SetVariable setvar8,fSize=12,limits={0,50,1},value= _NUM:20
76        SetVariable setvar9,pos={35,253.00},size={120.00,18.00},title="y-scaling"
77        SetVariable setvar9,fSize=12,limits={0,50,1},value= _NUM:20
78        SetVariable setvar10,pos={35,281.00},size={120.00,18.00},title="z-scaling"
79        SetVariable setvar10,fSize=12,limits={0,50,1},value= _NUM:20
80       
81        CheckBox check0,pos={43.00,15.00},size={64.00,15.00},title="Ellipsoid",fSize=12
82        CheckBox check0,value= 0,mode=1,proc=SuperCheckProc
83        CheckBox check1,pos={44.00,38.00},size={54.00,15.00},title="Toroid",fSize=12
84        CheckBox check1,value= 0,mode=1,proc=SuperCheckProc
85        CheckBox check2,pos={149.00,14.00},size={61.00,15.00},title="Quadric",fSize=12
86        CheckBox check2,value= 0,mode=1,proc=SuperCheckProc
87        CheckBox check3,pos={149.00,38.00},size={96.00,15.00},title="SuperFormula"
88        CheckBox check3,fSize=12,value= 1,mode=1,proc=SuperCheckProc
89
90        Button button0,pos={208.00,215.00},size={70.00,20.00},proc=SuperCalcButtonProc,title="Calculate"
91
92End
93
94Function SuperCheckProc(cba) : CheckBoxControl
95        STRUCT WMCheckboxAction &cba
96
97        switch( cba.eventCode )
98                case 2: // mouse up
99                        Variable checked = cba.checked
100                       
101                        // which radio button?
102                        // be sure that the others are "off"
103                       
104                        // disable the parameters not needed
105                        // enable the parameters that are needed
106                       
107                        NVAR gRadioVal= root:gSuperRadioVal
108       
109                        strswitch (cba.ctrlName)
110                                case "check0":          // Ellipsoid
111                                        gRadioVal= 1
112                                       
113                                        SetVariable setvar0 disable=0
114                                        SetVariable setvar1 disable=0
115                                        SetVariable setvar2 disable=2
116                                        SetVariable setvar3 disable=2
117                                        SetVariable setvar4 disable=2
118                                        SetVariable setvar5 disable=2
119                                        SetVariable setvar6 disable=2
120                                        SetVariable setvar7 disable=2
121                                       
122                                        break
123                                case "check1":          // Toroid
124                                        gRadioVal= 2
125                                        SetVariable setvar0 disable=0
126                                        SetVariable setvar1 disable=0
127                                        SetVariable setvar2 disable=2
128                                        SetVariable setvar3 disable=0
129                                        SetVariable setvar4 disable=2
130                                        SetVariable setvar5 disable=2
131                                        SetVariable setvar6 disable=2
132                                        SetVariable setvar7 disable=2
133                                        break
134                                case "check2":          // Quadric
135                                        gRadioVal= 3
136                                        SetVariable setvar0 disable=0
137                                        SetVariable setvar1 disable=0
138                                        SetVariable setvar2 disable=0
139                                        SetVariable setvar3 disable=2
140                                        SetVariable setvar4 disable=2
141                                        SetVariable setvar5 disable=2
142                                        SetVariable setvar6 disable=2
143                                        SetVariable setvar7 disable=2
144                                        break
145                                case "check3":          // SuperFormula
146                                        gRadioVal= 4
147                                        SetVariable setvar0 disable=2
148                                        SetVariable setvar1 disable=2
149                                        SetVariable setvar2 disable=2
150                                        SetVariable setvar3 disable=2
151                                        SetVariable setvar4 disable=0
152                                        SetVariable setvar5 disable=0
153                                        SetVariable setvar6 disable=0
154                                        SetVariable setvar7 disable=0
155                                        break
156                        endswitch
157                        CheckBox check0,value= gRadioVal==1
158                        CheckBox check1,value= gRadioVal==2
159                        CheckBox check2,value= gRadioVal==3
160                        CheckBox check3,value= gRadioVal==4
161                       
162
163                       
164                        break
165                case -1: // control being killed
166                        break
167        endswitch
168
169        return 0
170End
171
172Function SuperCalcButtonProc(ba) : ButtonControl
173        STRUCT WMButtonAction &ba
174
175        switch( ba.eventCode )
176                case 2: // mouse up
177                        // click code here
178                       
179                        // switch to the proper function call
180                        // read in the values from the panel
181                        // calculate the shape
182                        NVAR gRadioVal= root:gSuperRadioVal
183
184                        Variable aa,bb,cc
185                        Variable mm,n1,n2,n3
186                        Variable rr,ss,tt,rad
187                       
188                        ControlInfo setvar8
189                        aa = V_Value
190                        ControlInfo setvar9
191                        bb = V_Value
192                        ControlInfo setvar10
193                        cc = V_Value
194
195
196                        switch (gRadioVal)
197                                case 1:         // Ellipsoid
198                                        ControlInfo setvar0
199                                        rr = V_Value
200                                        ControlInfo setvar1
201                                        tt = V_Value                           
202                                       
203                                        isInsideSuperEllipsoid(rr,tt,aa,bb,cc)
204                                       
205                                        break
206                                case 2:         // Toroid
207                                        ControlInfo setvar0
208                                        rr = V_Value
209                                        ControlInfo setvar1
210                                        tt = V_Value                           
211                                        ControlInfo setvar3
212                                        rad = V_Value
213                       
214                                        isInsideSuperToroid(rr,tt,aa,bb,cc,rad)
215
216                                        break
217                                case 3:         // Quadric
218                                        ControlInfo setvar0
219                                        rr = V_Value
220                                        ControlInfo setvar1
221                                        tt = V_Value                           
222                                        ControlInfo setvar2
223                                        ss = V_Value
224
225                                        isInsideSuperQuadric(rr,ss,tt,aa,bb,cc)
226                                       
227                                        break
228                                case 4:         // SuperFormula
229                                        ControlInfo setvar4
230                                        mm = V_Value
231                                        ControlInfo setvar5
232                                        n1 = V_Value                           
233                                        ControlInfo setvar6
234                                        n2 = V_Value
235                                        ControlInfo setvar7
236                                        n3 = V_Value
237
238                                        isInsideSuperFormula(mm,n1,n2,n3,aa,bb,cc)
239
240                                        break
241                        endswitch
242
243                       
244                       
245                        break
246                case -1: // control being killed
247                        break
248        endswitch
249
250        return 0
251End
252
253
254
255
256
257
258
259
260
261
262
263// the 3D version
264// partially converted
265//
266//              superFormula_3d(1.5,1.5,7,7,1,1)
267//      superFormula_3d(1.1,1.1,.2,.2,1,1)
268//
269//
270Function superFormula_3d(n1,n2,n3,n4,a1,a2)
271        Variable n1,n2,n3,n4,a1,a2
272
273        variable np,ii,jj,nu,nv,raux1,raux2,sig,r1,r2
274       
275        np=1000
276        Make/O/D/N=(np) u3,v3
277       
278        Make/O/D/N=(np+1,np+1,3) M_Parametric//,r1,r2
279
280//      u3 = (x/np)*2*pi                // range from (0,2Pi)           //wrong definition
281//      v3 = (x/np)*pi          // range from (0,pi)   
282               
283        u3 = (x/np)*2*pi - pi           // range from (-pi, pi)
284        v3 = (x/np)*pi - pi/2           // range from (-pi/2, pi/2)
285
286        nu = np
287        nv = np
288        sig = 1
289       
290        for(ii=0;ii<nu;ii+=1)
291                for(jj=0;jj<nv;jj+=1)
292                        raux1 = abs(1/a1*abs(cos(n1*u3[ii]/4)))^n3+abs(1/a2*abs(sin(n1*u3[ii]/4)))^n4
293                        r1 = abs(raux1)^(-1/n2)
294                        raux2 = abs(1/a1*abs(cos(n1*v3[jj]/4)))^n3+abs(1/a2*abs(sin(n1*v3[jj]/4)))^n4
295                        r2 = abs(raux2)^(-1/n2)
296
297// the three values here are XYZ calculated from the polar representation of the superformula
298//
299                        M_Parametric[ii][jj][0] = r1*cos(u3[ii])*r2*cos(v3[jj])
300                        M_Parametric[ii][jj][1] = r1*sin(u3[ii])*r2*cos(v3[jj])
301                        M_Parametric[ii][jj][2] = sig*r2*sin(v3[jj])                   
302                endfor
303        endfor
304                       
305        return(0)
306end
307
308
309// the 3D version
310// partially converted
311//
312//              superFormula_3d_v2(3.5,3.5,1,1,1,1,1)
313//      superFormula_3d_v2(1.1,1.1,.5,.5,.5,1,1)
314//              superFormula_3d_v2(3.5,3.5,2,7,15,1,1)
315//
316// a more general version - note that there is a 7th parameter
317//
318Function superFormula_3d_v2(m1,m2,n1,n2,n3,aa,bb)
319        Variable m1,m2,n1,n2,n3,aa,bb
320
321        variable np,ii,jj,nu,nv,raux1,raux2,sig,r1,r2
322       
323        np=1000
324        Make/O/D/N=(np) u3,v3
325       
326        Make/O/D/N=(np+1,np+1,3) M_Parametric//,r1,r2
327
328//      u3 = (x/np)*2*pi                // range from (0,2Pi)           //this is the standard definition of spherical coordinates
329//      v3 = (x/np)*pi          // range from (0,pi)            //but not what was used in this case
330               
331        u3 = (x/np)*2*pi - pi           // range from (-pi, pi)
332        v3 = (x/np)*pi - pi/2           // range from (-pi/2, pi/2)
333
334        nu = np
335        nv = np
336        sig = 1
337       
338        for(ii=0;ii<nu;ii+=1)
339                for(jj=0;jj<nv;jj+=1)
340                        raux1 = abs(1/aa*abs(cos(m1*u3[ii]/4)))^n2+abs(1/bb*abs(sin(m2*u3[ii]/4)))^n3
341                        r1 = abs(raux1)^(-1/n1)
342                        raux2 = abs(1/aa*abs(cos(m1*v3[jj]/4)))^n2+abs(1/bb*abs(sin(m2*v3[jj]/4)))^n3
343                        r2 = abs(raux2)^(-1/n1)
344
345// the three values here are XYZ calculated from the polar representation of the superformula
346//
347                        M_Parametric[ii][jj][0] = r1*cos(u3[ii])*r2*cos(v3[jj])
348                        M_Parametric[ii][jj][1] = r1*sin(u3[ii])*r2*cos(v3[jj])
349                        M_Parametric[ii][jj][2] = sig*r2*sin(v3[jj])                   
350                endfor
351        endfor
352                       
353        return(0)
354end
355
356
357
358 
359Function SuperEllipsoid_3d(rr,tt,aa,bb,cc)
360        Variable rr,tt,aa,bb,cc
361
362        variable np,ii,jj,nu,nv,r1,r2
363       
364        np=1000
365        Make/O/D/N=(np) uu,vv
366       
367        Make/O/D/N=(np+1,np+1,3) M_Parametric
368
369//      u3 = (x/np)*2*pi                // range from (0,2Pi)           //wrong definition
370//      v3 = (x/np)*pi          // range from (0,pi)   
371               
372        uu = (x/np)*2*pi - pi           // range from (-pi, pi)
373        vv = (x/np)*pi - pi/2           // range from (-pi/2, pi/2)
374
375        nu = np
376        nv = np
377       
378        for(ii=0;ii<nu;ii+=1)
379                for(jj=0;jj<nv;jj+=1)
380
381                        M_Parametric[ii][jj][0] = aa*SE_C(vv[jj],2/tt)*SE_C(uu[ii],2/rr)
382                        M_Parametric[ii][jj][1] = bb*SE_C(vv[jj],2/tt)*SE_S(uu[ii],2/rr)
383                        M_Parametric[ii][jj][2] = cc*SE_S(vv[jj],2/tt)
384                endfor
385        endfor
386                       
387        return(0)
388end
389
390Function SE_C(ww,mm)
391        Variable ww,mm
392       
393        Variable ans
394       
395        ans = sign(cos(ww)) * (abs(cos(ww)))^mm
396       
397        return(ans)
398End
399
400
401Function SE_S(ww,mm)
402        Variable ww,mm
403       
404        Variable ans
405       
406        ans = sign(sin(ww)) * (abs(sin(ww)))^mm
407       
408        return(ans)
409End
410
411
412// will plot either the 2d or 3d version, whichever was most recently
413// calculated -- the M_parametric wave is plotted
414Window Gizmo_superSurface() : GizmoPlot
415        PauseUpdate; Silent 1           // building window...
416        // Building Gizmo 7 window...
417        NewGizmo/T="Gizmo0"/W=(286,371,959,941)
418        ModifyGizmo startRecMacro=700
419        ModifyGizmo scalingOption=63
420        AppendToGizmo Surface=root:M_Parametric,name=surface0
421        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ surfaceColorType,1}
422        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ srcMode,4}
423        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ frontColor,1,0,0,1}
424        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ backColor,0,0,1,1}
425        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ SurfaceCTABScaling,16}
426        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ textureType,1}
427        ModifyGizmo modifyObject=surface0,objectType=Surface,property={calcNormals,1}
428        AppendToGizmo light=Directional,name=light0
429        ModifyGizmo modifyObject=light0,objectType=light,property={ position,-0.832778,0.305999,0.461353,0.000000}
430        ModifyGizmo modifyObject=light0,objectType=light,property={ direction,-0.832778,0.305999,0.461353}
431        ModifyGizmo modifyObject=light0,objectType=light,property={ ambient,0.533333,0.533333,0.533333,1.000000}
432        ModifyGizmo modifyObject=light0,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
433        AppendToGizmo freeAxesCue={0,0,0,1.3},name=freeAxesCue0
434        AppendToGizmo Axes=boxAxes,name=axes0
435        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
436        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
437        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
438        ModifyGizmo setDisplayList=0, opName=pushAttribute0, operation=pushAttribute, data=1
439        ModifyGizmo setDisplayList=1, object=light0
440        ModifyGizmo setDisplayList=2, object=surface0
441        ModifyGizmo setDisplayList=3, opName=popAttribute0, operation=popAttribute
442        ModifyGizmo setDisplayList=4, object=freeAxesCue0
443        ModifyGizmo setDisplayList=5, object=axes0
444        ModifyGizmo autoscaling=1
445        ModifyGizmo currentGroupObject=""
446        ModifyGizmo showInfo
447        ModifyGizmo infoWindow={1129,602,1946,899}
448        ModifyGizmo endRecMacro
449        ModifyGizmo SETQUATERNION={-0.111733,-0.610876,-0.767227,-0.160412}
450EndMacro
451
452
453
454///////
455
456//
457Proc Setup_super3D_waves(num)
458        Variable num=128
459        Make/O/D/N=(num,num,num) superW
460        SetScale/P x -(num/2),1,"",superW
461        SetScale/P y -(num/2),1,"",superW
462        SetScale/P z -(num/2),1,"",superW
463       
464End
465
466////
467//
468Function isInsideSuperEllipsoid(rr,tt,aa,bb,cc)
469        Variable rr,tt,aa,bb,cc
470
471        if(exists("root:mat") == 0)
472                FFT_MakeMatrixButtonProc("")
473        endif
474               
475        Duplicate/O root:mat inW, voxW
476        Redimension/D inW
477        Variable num=DimSize(inW, 0)
478        SetScale/P x -(num/2),1,"",inW
479        SetScale/P y -(num/2),1,"",inW
480        SetScale/P z -(num/2),1,"",inW
481       
482        // using the wave scaling
483        inW = ( abs(x/aa)^rr + abs(y/bb)^rr )^(tt/rr) + abs(z/cc)^tt
484        voxW = inW[p][q][r] <= 1 ? 1 : 0
485
486        Wave w = root:mat
487        w = voxW
488       
489        return(0)
490end
491
492// pass superW to this function, really just for the scaling...
493// so I need to fix this in the future
494Function isInsideSuperToroid(rr,tt,aa,bb,cc,rad)
495        Variable rr,tt,aa,bb,cc,rad
496       
497        Variable dd
498        if(exists("root:mat") == 0)
499                FFT_MakeMatrixButtonProc("")
500        endif
501               
502        Duplicate/O root:mat inW, voxW
503        Redimension/D inW
504        Variable num=DimSize(inW, 0)
505        SetScale/P x -(num/2),1,"",inW
506        SetScale/P y -(num/2),1,"",inW
507        SetScale/P z -(num/2),1,"",inW
508       
509        // using the wave scaling
510        dd = rad/sqrt(aa^2 + bb^2)
511        inW = ( ( abs(x/aa)^rr + abs(y/bb)^rr )^(1/rr) - dd )^tt + abs(z/cc)^tt
512        voxW = inW[p][q][r] <= 1 ? 1 : 0
513
514        Wave w = root:mat
515        w = voxW
516       
517        return(0)
518end
519
520// pass superW to this function, really just for the scaling...
521// so I need to fix this in the future
522//
523// See https://en.wikipedia.org/w/index.php?title=Superquadrics&oldid=770878683
524//
525Function isInsideSuperQuadric(rr,ss,tt,aa,bb,cc)
526        Variable rr,ss,tt,aa,bb,cc
527       
528        if(exists("root:mat") == 0)
529                FFT_MakeMatrixButtonProc("")
530        endif
531               
532        Duplicate/O root:mat inW, voxW
533        Redimension/D inW
534        Variable num=DimSize(inW, 0)
535        SetScale/P x -(num/2),1,"",inW
536        SetScale/P y -(num/2),1,"",inW
537        SetScale/P z -(num/2),1,"",inW
538       
539        // using the wave scaling
540        inW = ( abs(x/aa)^rr + abs(y/bb)^ss + abs(z/cc)^tt )
541        voxW = inW[p][q][r] <= 1 ? 1 : 0
542
543        Wave w = root:mat
544        w = voxW
545       
546        return(0)
547end
548
549
550// pass superW to this function, really just for the scaling...
551// so I need to fix this in the future
552//
553// this is a modified version of the superformula
554//
555// I have fixed a == b == 1. almost all examples have this, and the
556// equation goes haywire if you stray far from one. For ease in use with voxels,
557// I've added aa,bb,cc as scaling values for x,y,z, in the same way it is done for the
558// ellipsoid and qudric equations
559//
560//
561// This could be further extended by using a separate set of parameters (m,n1,n2,n3) for
562// r_phi and r_theta.
563//
564// The "m" values could also be different (m1, m2) for each r_ calculation.
565//
566Function isInsideSuperFormula(m,n1,n2,n3,aa,bb,cc)
567        Variable m,n1,n2,n3,aa,bb,cc
568       
569        if(exists("root:mat") == 0)
570                FFT_MakeMatrixButtonProc("")
571        endif
572               
573        Duplicate/O root:mat inW, voxW
574        Redimension/D inW
575        Variable num=DimSize(inW, 0)
576        SetScale/P x -(num/2),1,"",inW
577        SetScale/P y -(num/2),1,"",inW
578        SetScale/P z -(num/2),1,"",inW
579       
580        // using the wave scaling
581       
582        // calculate r
583        // (not needed - calculate as part of asin()
584       
585        // calculate phi
586        Duplicate/O inW phi
587        phi = atan2(y,x)
588       
589        // calculate theta
590        Duplicate/O inW theta
591        theta = asin(z/sqrt(x^2+y^2+z^2))
592       
593        // calculate r(phi)
594        Duplicate/O inW r_phi
595//      r_phi = ( abs(1/aa*cos(m*phi/4))^n2+abs(1/bb*sin(m*phi/4))^n3 )^(-1/n1)
596        r_phi = ( abs(cos(m*phi/4))^n2+abs(sin(m*phi/4))^n3 )^(-1/n1)
597       
598        // calculate r(theta)
599        Duplicate/O inW r_theta
600//      r_theta = ( abs(1/aa*cos(m*theta/4))^n2+abs(1/bb*sin(m*theta/4))^n3 )^(-1/n1)
601        r_theta = ( abs(cos(m*theta/4))^n2+abs(sin(m*theta/4))^n3 )^(-1/n1)
602               
603        // evaluate for inside/outside
604//      inW = (x/(r_theta*r_phi))^2 + (y/(r_theta*r_phi))^2 + (z/(r_theta))^2           // ?? add cc here to scale z
605        inW = (x/(r_theta*r_phi*aa))^2 + (y/(r_theta*r_phi*bb))^2 + (z/(r_theta*cc))^2          // ?? add cc here to scale z
606       
607        voxW = inW[p][q][r] <= 1 ? 1 : 0
608
609        Wave w = root:mat
610        Redimension/B voxW
611        w = voxW
612
613        KillWaves/Z phi,theta,r_phi,r_theta     
614        return(0)
615end
616
617
618//
619// plots the voxelgram generated by isInsideSuperEllipse()
620//
621// execute mat = voxW, then the voxelgram can be FFT'd, or Debye's method
622//
623Window Gizmo_superVox() : GizmoPlot
624        PauseUpdate; Silent 1           // building window...
625        // Building Gizmo 7 window...
626        NewGizmo/W=(1810,345,2232,750)
627        ModifyGizmo startRecMacro=700
628        ModifyGizmo scalingOption=63
629        AppendToGizmo voxelgram=root:voxW,name=voxelgram0
630        ModifyGizmo ModifyObject=voxelgram0,objectType=voxelgram,property={ valueRGBA,0,1,0.000015,0.195544,0.800000,0.100008}
631        ModifyGizmo ModifyObject=voxelgram0,objectType=voxelgram,property={ mode,0}
632        ModifyGizmo ModifyObject=voxelgram0,objectType=voxelgram,property={ pointSize,2}
633        AppendToGizmo freeAxesCue={0,0,0,1},name=freeAxesCue0
634        AppendToGizmo Axes=boxAxes,name=axes0
635        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
636        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
637        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,3}
638        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,3}
639        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,3}
640        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
641        AppendToGizmo light=Directional,name=light0
642        ModifyGizmo modifyObject=light0,objectType=light,property={ position,0.514700,0.447400,-0.731400,0.000000}
643        ModifyGizmo modifyObject=light0,objectType=light,property={ direction,0.514700,0.447400,-0.731400}
644        ModifyGizmo modifyObject=light0,objectType=light,property={ ambient,0.933300,0.933300,0.933300,1.000000}
645        ModifyGizmo modifyObject=light0,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
646        AppendToGizmo attribute blendFunc={770,771},name=blendFunc0
647        ModifyGizmo setDisplayList=0, object=freeAxesCue0
648        ModifyGizmo setDisplayList=1, object=light0
649        ModifyGizmo setDisplayList=2, attribute=blendFunc0
650        ModifyGizmo setDisplayList=3, object=voxelgram0
651        ModifyGizmo setDisplayList=4, object=axes0
652        ModifyGizmo autoscaling=1
653        ModifyGizmo currentGroupObject=""
654        ModifyGizmo showInfo
655        ModifyGizmo infoWindow={1476,23,2293,321}
656        ModifyGizmo endRecMacro
657        ModifyGizmo SETQUATERNION={0.140038,0.591077,0.776874,0.165764}
658EndMacro
659
660
661
662/////////////////////////////////////////////////////
663/// below is from spherical harmonic demo -- and needs to be adapted to the SuperFormula
664
665Function calcParametric(L,M,size)
666        Variable L,M,size
667       
668        Make/O/N=(size+1,size+1,3) M_Parametric
669       
670        Variable i,j,dt,df,rr,theta,phi
671        dt=pi/size
672        df=2*dt
673       
674        SetScale/P x 0,1,"", M_Parametric
675        SetScale/P y 0,1,"", M_Parametric
676       
677        for(i=0;i<=size;i+=1)
678                theta=i*dt
679                for(j=0;j<=size;j+=1)
680                        phi=j*df
681                        rr=sqrt(magsqr(sphericalHarmonics(L,m,theta,phi)))
682                        M_Parametric[i][j][0]=rr*sin(theta)*cos(phi)
683                        M_Parametric[i][j][1]=rr*sin(theta)*sin(phi)
684                        M_Parametric[i][j][2]=rr*cos(theta)
685                endfor
686        endfor
687End
688
689Function updateParamSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
690        String ctrlName
691        Variable varNum
692        String varStr
693        String varName
694       
695        NVAR LL,MM,resolution
696       
697        if(abs(MM)>LL)
698                if(MM<0)
699                        MM=-LL
700                else
701                        MM=LL
702                endif
703        endif
704       
705        calcParametric(LL,MM,resolution)
706End
707
708
709//// from another WM Demo
710Function makeSphere(pointsx,pointsy)
711        Variable pointsx,pointsy
712       
713        Variable i,j,rad
714        Make/O/n=(pointsx,pointsy,3) sphereData
715        Variable anglePhi,angleTheta
716        Variable dPhi,dTheta
717       
718       
719        dPhi=2*pi/(pointsx-1)
720        dTheta=pi/(pointsy-1)
721        Variable xx,yy,zz
722        Variable sig
723       
724        for(j=0;j<pointsy;j+=1)
725                angleTheta=j*dTheta
726                zz=sin(angleTheta)
727                if(angleTheta>pi/2)
728                        sig=-1
729                else
730                        sig=1
731                endif
732                for(i=0;i<pointsx;i+=1)
733                        anglePhi=i*dPhi
734                        xx=zz*cos(anglePhi)
735                        yy=zz*sin(anglePhi)
736                        sphereData[i][j][0]=xx
737                        sphereData[i][j][1]=yy
738                        sphereData[i][j][2]=sig*sqrt(1-xx*xx-yy*yy)
739//                      sphereData[i][j][2]=sqrt(1-xx*xx-yy*yy)
740                endfor
741        endfor
742End
743
744
745//////////////////////////////////////////////////////////////////////////////////////////////
746Function makeDoubleCone(pointsx,pointsy)
747        Variable pointsx,pointsy
748       
749        Variable i,j,rad
750        Make/O/n=(pointsx,pointsy,3) sphereData
751        Variable anglePhi,angleTheta
752        Variable dPhi,dTheta
753       
754       
755        dPhi=2*pi/pointsx
756        dTheta=pi/pointsy
757        Variable xx,yy,zz
758       
759        for(j=0;j<pointsy;j+=1)
760                angleTheta=j*dTheta
761                zz=cos(angleTheta)
762                for(i=0;i<pointsx;i+=1)
763                        anglePhi=i*dPhi
764                        xx=zz*cos(anglePhi)
765                        yy=zz*sin(anglePhi)
766                        sphereData[i][j][0]=xx
767                        sphereData[i][j][1]=yy
768                        sphereData[i][j][2]=zz
769                endfor
770        endfor
771End
772
773/////////////////////////////
774
775
776
777////////////////////////////////////////////////
778//
779// another way to plot data as a cloud of points of an isosurface
780// change the isoValue to plot a different surface
781//
782Window Gizmo_Isosurface() : GizmoPlot
783        PauseUpdate; Silent 1           // building window...
784        // Building Gizmo 7 window...
785        NewGizmo/W=(35,45,550,505)
786        ModifyGizmo startRecMacro=700
787        ModifyGizmo scalingOption=63
788        AppendToGizmo isoSurface=root:superW,name=isoSurface0
789        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ surfaceColorType,1}
790        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineColorType,0}
791        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidthType,0}
792        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ fillMode,4}
793        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidth,1}
794        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ isoValue,30}
795        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ frontColor,1,0,0,1}
796        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ backColor,0,0,1,1}
797        AppendToGizmo Axes=boxAxes,name=axes0
798        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
799        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
800        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,2}
801        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,2}
802        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,2}
803        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
804        AppendToGizmo freeAxesCue={0,0,0,1},name=freeAxesCue0
805        ModifyGizmo setDisplayList=0, object=freeAxesCue0
806        ModifyGizmo setDisplayList=1, object=axes0
807        ModifyGizmo setDisplayList=2, object=isoSurface0
808        ModifyGizmo autoscaling=1
809        ModifyGizmo currentGroupObject=""
810        ModifyGizmo showInfo
811        ModifyGizmo infoWindow={551,23,1368,322}
812        ModifyGizmo endRecMacro
813        ModifyGizmo SETQUATERNION={0.055218,0.485731,0.843569,0.222280}
814EndMacro
815
816//
817// another cloud. would need a way to set the isoValue here as well, to pick the layer that is
818// viewed
819//
820Window mat_as_3dCloud() : GizmoPlot
821        PauseUpdate; Silent 1           // building window...
822        // Building Gizmo 7 window...
823        NewGizmo/W=(35,45,550,505)
824        ModifyGizmo startRecMacro=700
825        ModifyGizmo scalingOption=63
826        AppendToGizmo isoSurface=root:mat,name=isoSurface0
827        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ surfaceColorType,1}
828        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineColorType,0}
829        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidthType,0}
830        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ fillMode,4}
831        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidth,1}
832        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ isoValue,1}
833        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ frontColor,1,0,0,1}
834        ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ backColor,0,0,1,1}
835        AppendToGizmo Axes=boxAxes,name=axes0
836        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
837        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}
838        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
839        AppendToGizmo freeAxesCue={0,0,0,1.5},name=freeAxesCue0
840        AppendToGizmo voxelgram=root:mat,name=voxelgram0
841        ModifyGizmo ModifyObject=voxelgram0,objectType=voxelgram,property={ valueRGBA,0,1,1.000000,0.000000,0.000000,1.000000}
842        ModifyGizmo ModifyObject=voxelgram0,objectType=voxelgram,property={ mode,0}
843        ModifyGizmo setDisplayList=0, object=axes0
844        ModifyGizmo setDisplayList=1, object=freeAxesCue0
845        ModifyGizmo setDisplayList=2, object=isoSurface0
846        ModifyGizmo autoscaling=1
847        ModifyGizmo currentGroupObject=""
848        ModifyGizmo showInfo
849        ModifyGizmo infoWindow={551,23,1368,322}
850        ModifyGizmo endRecMacro
851        ModifyGizmo SETQUATERNION={0.205660,0.395226,0.765282,0.464591}
852EndMacro
853
854
855////////////////////////////////////////
856//
857// I have some of the 2D representations here for testing
858//
859Function superFormula_2d(n1,n2,n3,n4,a1,a2)
860        Variable n1,n2,n3,n4,a1,a2
861       
862        Variable np
863        np=1000
864       
865        Make/O/D/N=(np) u,xx,yy,raux,rr
866        u = (x/np)*(2*pi)
867       
868        raux=abs(1/a1*abs(cos(n1*u/4)))^n3+abs(1/a2*abs(sin(n1*u/4)))^n4
869        rr=abs(raux)^(-1/n2)
870        xx=rr*cos(u)
871        yy=rr*sin(u)
872 
873        return(0)
874end
875
876
877Function is_2DPointInside(xPt,yPt)
878        Variable xPt,yPt
879       
880        Variable uVal, rVal
881        uVal = atan2(yPt,xPt) + pi
882        rVal = xPt/cos(uVal)
883       
884//      Print uVal,rVal
885       
886        WAVE rr = root:rr
887        WAVE u = root:u
888       
889        Variable val
890        val = interp(uVal, u, rr )
891//      Print val
892       
893        if( rval^2 > val^2)
894//              Print "outside"
895                return(0)
896        else
897//              Print "inside"
898                return(1)
899        endif
900       
901        return(0)
902end
903
904Function testInsideOutside(range)
905        variable range
906        Make/O/D/N=10000 testX, testY, testZ
907        testX = enoise(-range,range)
908        testY = enoise(-range,range)
909        testZ = is_2DPointInside(-testx,-testy)         //not sure why, but it makes the graph correct
910End
911
912Window SF_2DGraph() : Graph
913        PauseUpdate; Silent 1           // building window...
914        Display /W=(1002,44,1607,594) yy vs xx
915        ModifyGraph mode=3
916        ModifyGraph marker=19
917        ModifyGraph msize=1
918EndMacro
919
920Window SF_2D_InsideOutside() : Graph
921        PauseUpdate; Silent 1           // building window...
922        Display /W=(1616,44,2170,600) testY vs testX
923        ModifyGraph mode=3
924        ModifyGraph marker=19
925        ModifyGraph msize=2
926        ModifyGraph zColor(testY)={testz,*,*,Rainbow}
927EndMacro
928
929
930Window SF_2DGraph_overlay() : Graph
931        PauseUpdate; Silent 1           // building window...
932        Display /W=(1002,44,1607,594) testY vs testX
933        AppendToGraph yy vs xx
934        ModifyGraph mode=3
935        ModifyGraph marker(testY)=16,marker(yy)=19
936        ModifyGraph msize=2
937        ModifyGraph zColor(testY)={testZ,*,*,BlueRedGreen256}
938        Cursor/P A yy 1
939        ShowInfo
940EndMacro
941
942//// end 2D versions
Note: See TracBrowser for help on using the repository browser.