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

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

a number of changes, mostly to allow everything to compile.

added conditional compile to ensure that XML code would not be compiled if VSANS was present, since it's not XML-aware.

modified V_MainPanel to avoid conflicts with the SANS version. There still may be some functions hidden in procedures that do not have the V_ prefix yet, but these are either for functions that should point to a common file, or procedures that have been hidden from the VSANS panel

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