source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Panel.ipf

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

added a few corrections to the reduction:

Added sample apertures to the patch panel so that they can be corrected

A flag is now written to the data files if the "flip" has been done, and it will refuse to flip again. This flag can be reset if something goes wrong.

Multiple reduce now allows run numbers to be entered as is for SANS,

Filtering of files for the protocol panel should be better behaved now.

File size: 35.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3//
4// utility functions and procedures for displaying information
5// setting the matrix, and doing the calculations
6//
7
8
9//
10// TO DO:
11// x- I need to change a lot of routines (most notably Gizmo) to take "10" as the default SLD
12//               rather than "1"
13//
14//
15// -- incorporate utility function that reads the wave note from a binary file
16//                      Function/S LoadNoteFunc(PName,FName[,FileRef]) - use this as an "inspector?"
17////////////
18
19//
20
21//***********
22//For a 3D slicer view of dmat:
23//
24//Duplicate/O dmat dmatView
25//
26//Open a New 3D voxelgram. Don't set any of the levels
27//Open the 3D Slicer under the Gizmo menu
28//Make the X, Y, and Z slices
29//(Yellow Hot seems to be the best)
30//
31//// for a better view
32//dmatView = log(dmat)
33//
34//// if it all goes blank after the log transform get rid of the INF
35//dmatView = numtype(dmatView)== 0 ? dmatView : 0
36//
37//*************
38
39
40///PANEL
41///// panel procedures
42
43Proc Init_FFT()
44        DoWindow/F FFT_Panel
45        if(V_flag==0)
46                SetDataFolder root:
47                //init the globals
48                Variable/G root:FFT_N=128
49                Variable/G root:FFT_T=5
50                Variable/G root:FFT_Qmax = 0
51                Variable/G root:FFT_QmaxReal = 0
52                Variable/G root:FFT_DQ=0
53                Variable/G root:FFT_Qmin=0
54                Variable/G root:FFT_estTime = 0
55               
56                Variable/G root:FFT_SolventSLD = 0
57                Variable/G root:FFT_delRho = 1e-7                       //multiplier for SLD (other value is 1e-7)
58               
59                FFT_Qmax :=2*pi/FFT_T
60                FFT_QmaxReal := FFT_Qmax/2
61                FFT_DQ := pi/(FFT_N*FFT_T)
62                FFT_Qmin := 2*pi/(FFT_N*FFT_T)
63                //empirical fit (cubic) of time vs N=50 to N=256
64                FFT_estTime := 0.56 - 0.0156*FFT_N + 0.000116*FFT_N^2 + 8e-7*FFT_N^3
65//              FFT_estTime := FFT_N/128
66               
67                FFT_Panel()
68        endif
69End
70
71
72Proc FFT_Panel()
73        PauseUpdate; Silent 1           // building window...
74        NewPanel /W=(1452,44,1768,531)/K=1 as "FFT_Panel"
75        DoWindow/C FFT_Panel
76        SetDrawLayer UserBack
77        DrawLine 5,68,311,68
78        DrawLine 5,142,311,142
79        DrawLine 5,250,311,250
80        SetVariable FFTSetVar_0,pos={7,7},size={150,15},title="Cells per edge (N)"
81        SetVariable FFTSetVar_0,limits={50,512,2},value= FFT_N,live= 1
82        SetVariable FFTSetVar_1,pos={7,27},size={150,15},title="Length per Cell (T)"
83        SetVariable FFTSetVar_1,limits={1,5000,0.2},value= FFT_T,live= 1
84        SetVariable FFTSetVar_2,pos={183,7},size={120,15},title="Real Qmax"
85        SetVariable FFTSetVar_2,limits={0,0,0},value= FFT_QmaxReal,noedit= 1,live= 1
86        SetVariable FFTSetVar_3,pos={183,47},size={120,15},title="delta Q (A)"
87        SetVariable FFTSetVar_3,limits={0,0,0},value= FFT_DQ,noedit= 1,live= 1
88        SetVariable FFTSetVar_6,pos={183,27},size={120,15},title="Real Qmin (A)"
89        SetVariable FFTSetVar_6,limits={0,0,0},value= FFT_Qmin,noedit= 1,live= 1
90        Button FFTButton_0,pos={15,79},size={90,20},proc=FFT_MakeMatrixButtonProc,title="Make Matrix"
91
92        Button FFTButton_3,pos={14,265},size={70,20},proc=DoTheFFT_ButtonProc,title="Do FFT"
93        Button FFTButton_4,pos={180,264},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot FFT Results"
94       
95        Button FFTButton_1,pos={14,150},size={90,20},proc=FFTMakeGizmoButtonProc,title="Make Gizmo"
96        Button FFTButton_2,pos={14,175},size={100,20},proc=FFTDrawSphereButtonProc,title="Draw Sphere"
97        Button FFTButton_5,pos={14,200},size={130,20},proc=FFTDrawZCylinderButtonProc,title="Draw XYZ Cylinder"
98        Button FFTButton_20,pos={14,225},size={130,20},proc=FFTDrawRotCylinderButton,title="Draw Rot Cylinder"
99       
100//      Button FFTButton_6,pos={134,79},size={90,20},proc=FFTEraseMatrixButtonProc,title="Erase Matrix"
101        Button FFTButton_6a,pos={140,79},size={50,20},proc=FFTSaveMatrixButtonProc,title="Save"
102        Button FFTButton_6b,pos={200,79},size={50,20},proc=FFTLoadMatrixButtonProc,title="Load"
103        Button FFTButton_6c,pos={260,79},size={50,20},proc=FFT_AddMatrixButtonProc,title="Add"
104        Button FFTButton_7,pos={13,329},size={130,20},proc=FFT_BinnedSpheresButtonProc,title="Do Binned Debye"
105        Button FFTButton_7a,pos={180,329},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Binned Results"
106
107        Button FFTButton_8,pos={13,297},size={130,20},proc=FFT_AltiSpheresButtonProc,title="Do Debye Spheres"
108        Button FFTButton_8a,pos={180,297},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Debye Results"
109
110        Button FFTButton_14,pos={13,360},size={130,20},proc=FFT_BinnedSLDButtonProc,title="Do Binned SLD"
111        Button FFTButton_14a,pos={180,360},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot SLD Results"
112
113        SetVariable FFTSetVar_4,pos={7,47},size={100,15},title="FFT time(s)"
114        SetVariable FFTSetVar_4,limits={0,0,0},value= FFT_estTime,noedit= 1,live= 1,format="%g"
115        Button FFTButton_9,pos={200,400},size={100,20},proc=FFT_Get2DSlice,title="Get 2D Slice"
116
117        Button FFTButton_19,pos={168,150},size={130,20},proc=FFT_ChangeMatrixValuesButton,title="Replace Voxels"
118        Button FFTButton_12,pos={168,175},size={130,20},proc=FFT_ReplaceSolventButton,title="Replace Solvent"
119        Button FFTButton_11,pos={169,200},size={130,20},proc=FFT_RotateMat,title="Rotate Matrix"
120        Button FFTButton_10,pos={169,225},size={130,20},proc=FFT_TransposeMat,title="Transpose Matrix"
121
122        Button FFTButton_13,pos={14,109},size={120,20},proc=FFTFillSolventMatrixProc,title="Solvent Matrix"
123        SetVariable FFTSetVar_5,pos={155,111},size={150,15},title="Solvent SLD (10^-7)"
124        SetVariable FFTSetVar_5,limits={-99,99,1},value= FFT_SolventSLD,live= 1
125        Button FFTButton_15,pos={209,430},size={90,20},proc=Interp2DSliceButton,title="Interp 2D"
126        Button FFTButton_16,pos={14,460},size={70,20},proc=FFTHelpButton,title="Help"
127       
128        Button FFTButton_17,pos={13,400},size={120,20},proc=FFT_Iso2USANS,title="Iso to USANS"
129        Button FFTButton_18,pos={13,430},size={120,20},proc=FFT_Aniso2USANS,title="Aniso to USANS"
130
131EndMacro
132
133// Save a matrix wave, plus the N, T, and solvent values in the wave note for reloading
134Function FFTDrawRotCylinderButton(ba) : ButtonControl
135        STRUCT WMButtonAction &ba
136       
137        String win = ba.win
138
139        switch (ba.eventCode)
140                case 2:
141                        // click code here
142                        Execute "FFTDrawRotCylinder()"
143                       
144                        break
145        endswitch
146
147        return 0
148End
149
150// Save a matrix wave, plus the N, T, and solvent values in the wave note for reloading
151Function FFTSaveMatrixButtonProc(ba) : ButtonControl
152        STRUCT WMButtonAction &ba
153       
154        String win = ba.win
155
156        switch (ba.eventCode)
157                case 2:
158                        // click code here
159                        String fileStr=""
160                        SaveMyMatrix(fileStr)
161                       
162                        break
163        endswitch
164
165        return 0
166End
167
168// this will wave as Igor Binary, so be sure to use the ".ibw extension.
169// - this could possibly be enforced, but that's maybe not necessary at this stage.
170//
171Function SaveMyMatrix(fileStr)
172        String fileStr
173       
174        WAVE mat=root:mat
175        NVAR FFT_T = root:FFT_T
176        NVAR FFT_N = root:FFT_N
177        NVAR FFT_SolventSLD = root:FFT_SolventSLD
178        String str=""
179        sprintf str,"FFT_T=%g;FFT_N=%d;FFT_SolventSLD=%d;",FFT_T,FFT_N,FFT_SolventSLD
180        Note mat,str
181        Save/C/P=home mat as fileStr    //will ask for a file name if fileStr="" save as Igor Binary
182        Note/K mat                      //kill wave note on exiting since I don't properly update this anywhere else
183                       
184        return(0)
185end
186
187
188// load in a previously saved matrix, and reset FFT_N, FFT_T and solvent
189// from the wave note when saved
190Function FFTLoadMatrixButtonProc(ba) : ButtonControl
191        STRUCT WMButtonAction &ba
192       
193        String win = ba.win
194
195        switch (ba.eventCode)
196                case 2:
197                        // click code here
198                        String fileStr=""
199                        ReloadMatrix(fileStr)
200                       
201                        break
202        endswitch
203
204        return 0
205End
206
207// /H flag on the LoadWave command severs the connection with the binary file
208// -- this seems to be important - otherwise I get odd results and the wave (on disk) can change!
209//
210Function ReloadMatrix(fileStr)
211        String fileStr
212       
213        LoadWave/H/M/O/W/P=home         fileStr         //will ask for a file, Igor Binary format is assumed here
214        if(V_flag == 0)
215                return(0)               //user cancel
216        endif
217       
218        String str
219        str=note(mat)
220        NVAR FFT_T = root:FFT_T
221        NVAR FFT_N = root:FFT_N
222        NVAR FFT_SolventSLD = root:FFT_SolventSLD
223       
224        FFT_T = NumberByKey("FFT_T", str, "=" ,";")
225        FFT_N = NumberByKey("FFT_N", str, "=" ,";")
226        FFT_SolventSLD = NumberByKey("FFT_SolventSLD", str, "=" ,";")
227
228// if I got bad values, put in default values                   
229        if(numtype(FFT_T) != 0 )
230                FFT_T = 5
231        endif
232        if(numtype(FFT_N) != 0 )
233                FFT_N = DimSize(mat,0)
234        endif
235        if(numtype(FFT_SolventSLD) != 0 )
236                FFT_SolventSLD = 0
237        endif                   
238
239        ColorizeGizmo()
240       
241        Print "Loaded matrix parameters = ",str
242        Execute "NumberOfPoints()"
243                       
244        return(0)
245end
246
247// load in a previously saved matrix, and reset FFT_N, FFT_T and solvent
248// from the wave note when saved
249Function FFT_AddMatrixButtonProc(ba) : ButtonControl
250        STRUCT WMButtonAction &ba
251       
252        String win = ba.win
253
254        switch (ba.eventCode)
255                case 2:
256                        // click code here
257                        String fileStr=""
258                        LoadAndAddMatrix(fileStr)
259                       
260                        break
261        endswitch
262
263        return 0
264End
265
266// This function will load and add the matrix to whatever is currently present
267// - it is checked that the N, T, and SolventSLD are the same
268//
269// - it currently SUMS the voxels - so if there is overlap, you get a new value
270//
271Function LoadAndAddMatrix(fileStr)
272        String fileStr
273
274        String toLoadStr=""
275
276// if the passed in file name is null, pick a file     
277        if(strlen(fileStr)==0)
278                toLoadStr = DoOpenFileDialog("Pick the binary file to load")
279                if(strlen(toLoadStr)==0)
280                        return(0)
281                endif
282        endif
283        fileStr=toLoadStr
284
285        // make sure that N and T are correct, just read in the note
286        String noteStr=LoadNoteFunc("home",fileStr)
287        String abortStr
288       
289        // current values
290        NVAR FFT_T = root:FFT_T
291        NVAR FFT_N = root:FFT_N
292        NVAR FFT_SolventSLD = root:FFT_SolventSLD
293
294// possible set to load
295        Variable test_T,test_N,test_SolventSLD
296        test_T = NumberByKey("FFT_T", noteStr, "=" ,";")
297        test_N = NumberByKey("FFT_N", noteStr, "=" ,";")
298        test_SolventSLD = NumberByKey("FFT_SolventSLD", noteStr, "=" ,";")
299       
300// if I got bad values, warn and abort                 
301        if(FFT_T != test_T)
302                abortStr = "Current T = "+num2str(FFT_T)+" and Selected T = "+num2str(test_T)+", aborting load"
303                Abort abortStr
304        endif
305        if(FFT_N != test_N)
306                abortStr = "Current N = "+num2str(FFT_N)+" and Selected N = "+num2str(test_N)+", aborting load"
307                Abort abortStr
308        endif
309        if(FFT_SolventSLD != test_SolventSLD)
310                abortStr = "Current SolventSLD = "+num2str(FFT_SolventSLD)+" and Selected SolventSLD = "+num2str(test_SolventSLD)+", aborting load"
311                Abort abortStr
312        endif   
313               
314        // OK, then load and add the matrix
315        // the loaded matrix is "mat" as usual, so make a copy of the current first
316        Duplicate/O mat tmpMat
317        LoadWave/H/M/O/W/P=home         fileStr         //will ask for a file, Igor Binary format is assumed here
318        if(V_flag == 0)
319                return(0)               //user cancel
320        endif
321       
322        Wave mat=root:mat
323        FastOp mat = mat + tmpMat
324       
325        KillWaves/Z tmpMat
326       
327        Print "Loaded matrix parameters = ",noteStr
328        Execute "NumberOfPoints()"
329       
330        ColorizeGizmo()
331               
332        return(0)
333end
334
335// utility function that reads the wave note from a binary file
336// where did I get this from? Igor Exchange? I didn't write this myself...
337//
338Function/S LoadNoteFunc(PName,FName[,FileRef])
339        String PName, FName
340        Variable FileRef
341       
342        Variable noteStart, noteLength, version, dependLength
343        String noteStr, typeStr = ".ibw"
344        if (ParamIsDefault(FileRef))
345                Open/R/P=$PName/T=typeStr fileRef, as FName     //open the file
346        endif
347        FSetPos fileRef, 0
348        FBinRead/F=2 fileRef, version
349       
350       
351        if (version == 5)
352       
353                FSetPos fileRef, 4
354                Make/N=(3)/I/Free SizeWave
355                FBinRead FileRef,SizeWave
356                noteStart = SizeWave[0]
357                DependLength = SizeWave[1]
358                NoteLength = SizeWave[2]
359                noteStart += dependLength+64
360               
361        elseif (version == 2)
362       
363                FBinRead/F=3 fileRef, noteStart
364//              FBinRead/F=4 fileRef, dependLength
365                FBinRead/F=3 fileRef, noteLength
366                noteStart += 16
367               
368        else
369       
370                if (ParamIsDefault(FileRef))
371                        Close(FileRef)          //close the file
372                endif
373                return ""
374       
375        endif
376        if (!NoteLength)
377                if (ParamIsDefault(FileRef))
378                        Close(FileRef)          //close the file
379                endif
380                return("")
381        endif
382        FSetPos fileRef, noteStart
383        NoteStr = PadString("",NoteLength,0)
384        FBinRead FileRef,NoteStr
385       
386        if (ParamIsDefault(FileRef))
387                Close(FileRef)          //close the file
388        endif
389        return noteStr
390       
391End //LoadNoteFunc
392
393
394Function FFTHelpButton(ba) : ButtonControl
395        STRUCT WMButtonAction &ba
396       
397        String win = ba.win
398
399        switch (ba.eventCode)
400                case 2:
401                        // click code here
402                        DisplayHelpTopic/Z/K=1 "Real-Space Modeling of SANS Data"
403                        if(V_flag !=0)
404                                DoAlert 0,"The Real-Space Modeling Help file could not be found"
405                        endif
406                        break
407        endswitch
408
409        return 0
410End
411
412
413Function Interp2DSliceButton(ba) : ButtonControl
414        STRUCT WMButtonAction &ba
415
416        switch( ba.eventCode )
417                case 2: // mouse up
418                        // click code here
419                        String folderStr=""
420                        Prompt folderStr, "Pick a 2D data folder",popup,W_DataSetPopupList()            // Set prompt for x param
421                        DoPrompt "Enter data folder", folderStr
422                        if (V_Flag || strlen(folderStr) == 0)
423                                return -1                                                               // User canceled, null string entered
424                        endif                   
425       
426                        Interpolate2DSliceToData(folderStr)
427                       
428                        //display the 2D plane
429                        DoWindow FFT_Interp2D_log
430                        if(V_flag == 0)
431                                Execute "Display2DInterpSlice_log()"
432                        endif
433                       
434                        break
435        endswitch
436
437        return 0
438End
439
440
441
442Function FFT_Get2DSlice(ctrlName) : ButtonControl
443        String ctrlName
444               
445        WAVE/Z dmat = root:dmat
446        if(WaveExists(dmat)==1)
447                get2DSlice(dmat)
448        endif
449       
450        //display the 2D plane
451        DoWindow FFT_Slice2D
452        if(V_flag == 0)
453                Execute "Display2DSlice()"
454        endif
455       
456        //display the 2D plane
457        DoWindow FFT_Slice2D_log
458        if(V_flag == 0)
459                Execute "Display2DSlice_log()"
460        endif
461       
462        //display the 1D binning of the 2D plane
463        DoWindow FFT_Slice1D
464        if(V_flag == 0)
465                Execute "Slice2_1D()"
466        endif
467       
468        return(0)       
469End
470
471Proc Display2DSlice()
472        PauseUpdate; Silent 1           // building window...
473        Display /W=(1038,44,1404,403)
474        DoWindow/C FFT_Slice2D
475        AppendImage/T detPlane
476        ModifyImage detPlane ctab= {*,*,YellowHot,0}
477        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
478        ModifyGraph mirror=2
479        ModifyGraph nticks=4
480        ModifyGraph minor=1
481        ModifyGraph fSize=9
482        ModifyGraph standoff=0
483        ModifyGraph tkLblRot(left)=90
484        ModifyGraph btLen=3
485        ModifyGraph tlOffset=-2
486//      SetAxis/A/R left
487End
488
489Proc Display2DInterpSlice_log()
490        PauseUpdate; Silent 1           // building window...
491        Display /W=(1038,44,1404,403)
492        DoWindow/C FFT_Interp2D_log
493        AppendImage/T interp2DSlice_log
494        ModifyImage interp2DSlice_log ctab= {*,*,YellowHot,0}
495        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
496        ModifyGraph mirror=2
497        ModifyGraph nticks=4
498        ModifyGraph minor=1
499        ModifyGraph fSize=9
500        ModifyGraph standoff=0
501        ModifyGraph tkLblRot(left)=90
502        ModifyGraph btLen=3
503        ModifyGraph tlOffset=-2
504//      SetAxis/A/R left
505End
506
507Proc Display2DSlice_log()
508        PauseUpdate; Silent 1           // building window...
509        Display /W=(1038,44,1404,403)
510        DoWindow/C FFT_Slice2D_log
511        AppendImage/T logP
512        ModifyImage logP ctab= {*,*,YellowHot,0}
513        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
514        ModifyGraph mirror=2
515        ModifyGraph nticks=4
516        ModifyGraph minor=1
517        ModifyGraph fSize=9
518        ModifyGraph standoff=0
519        ModifyGraph tkLblRot(left)=90
520        ModifyGraph btLen=3
521        ModifyGraph tlOffset=-2
522//      SetAxis/A/R left
523End
524Proc Slice2_1D()
525        PauseUpdate; Silent 1           // building window...
526        Display /W=(1034,425,1406,763) iBin_2d vs qBin_2d
527        DoWindow/C FFT_Slice1D
528        ModifyGraph mode=4
529        ModifyGraph marker=19
530        ModifyGraph msize=2
531        ModifyGraph grid=1
532        ModifyGraph log=1
533        ModifyGraph mirror=2
534        Legend
535EndMacro
536
537
538Function FFT_TransposeMat(ctrlName) : ButtonControl
539        String ctrlName
540
541        Variable mode=1
542        Prompt mode,"Transform XYZ to:",popup,"XZY;ZXY;ZYX;YZX;YXZ;"
543        DoPrompt "Transform mode",mode
544        if (V_Flag)
545                return 0                                                                        // user canceled
546        endif
547       
548        fFFT_TransposeMat(mode)
549       
550        return (0)
551End
552
553// starting in XYZ
554// mode 1;2;3;4;5;
555//correspond to
556// "XZY;ZXY;ZYX;YZX;YXZ;"
557//
558Function fFFT_TransposeMat(mode)
559        Variable mode
560
561        WAVE/Z mat=mat
562        ImageTransform /G=(mode) transposeVol mat
563        WAVE M_VolumeTranspose=M_VolumeTranspose
564        Duplicate/O M_VolumeTranspose mat
565       
566        return(0)       
567End
568
569Function FFT_RotateMat(ctrlName) : ButtonControl
570        String ctrlName
571       
572        Variable angleX=45,angleY=0,angleZ=0
573        Prompt angleX, "Degrees of rotation around X-axis:"
574        Prompt angleY, "Degrees of rotation around Y-axis:"
575        Prompt angleZ, "Degrees of rotation around Z-axis:"
576        DoPrompt "Enter angles for rotation", angleX,angleY,angleZ
577       
578        if (V_Flag)
579                return 0                                                                        // user canceled
580        endif
581
582        XYZRotate(angleX,angleY,angleZ)
583
584
585//////////// old way, only one angle
586//      Variable degree=45,sense=1
587//      Prompt degree, "Degrees of rotation around Z-axis:"
588////    Prompt sense, "Direction of rotation:",popup,"CW;CCW;"
589//      DoPrompt "Enter parameters for rotation", degree
590//     
591//      if (V_Flag)
592//              return 0                                                                        // user canceled
593//      endif
594       
595
596        // old way using ImageRotate that interpolates, and is only around Z-axis
597//      fFFT_RotateMat(degree)
598        return(0)
599End
600
601// note that the rotation is not perfect. if the rotation produces an
602// odd number of pixels, then the object will "walk" one pixel. Also, some
603// small artifacts are introduced, simply due to the voxelization of the object
604// as it is rotated. rotating a cylinder 10 then 350 shows a few extra "bumps" on
605// the surface, but the calculated scattering is virtually identical.
606//
607// these issues, may be correctable, if needed
608//
609Function fFFT_RotateMat(degree)
610        Variable degree
611       
612        Variable fill=0
613        NVAR solventSLD = root:FFT_SolventSLD
614        fill = solventSLD
615       
616        WAVE mat = root:mat
617        Variable dx,dy,dz,nx,ny,nz
618        dx = DimSize(mat,0)
619        dy = DimSize(mat,1)
620        dz = DimSize(mat,2)
621               
622        //? the /W and /C flags seem to be special cases for 90 deg rotations
623        ImageRotate /A=(degree)/E=(fill) mat            //mat is not overwritten, unless /O
624       
625        nx = DimSize(M_RotatedImage,0)
626        ny = DimSize(M_RotatedImage,1)
627        nz = DimSize(M_RotatedImage,2)
628//      Print "rotated dims = ",dx,dy,dz,nx,ny,nz
629        Variable delx,dely,odd=0
630        delx = nx-dx
631        dely = ny-dy
632       
633        if(mod(delx,2) != 0)            //sometimes the new x,y dims are odd!
634                odd = 1
635        endif
636//      delx = (delx + odd)/2
637//      dely = (dely + odd)/2
638        delx = trunc(delx/2)
639        dely = trunc(dely/2)
640
641        //+odd removes an extra point if there is one
642        Duplicate/O/R=[delx+odd,nx-delx-1][dely+odd,ny-dely-1][0,nz-1] M_RotatedImage mat
643
644// - not sure why the duplicate operation changes the start and delta of mat, but I
645// need to reset the start and delta to be sure that the display is correct, and that the scaling
646// is read correctly later
647        SetScale/P x 0,1,"", mat
648        SetScale/P y 0,1,"", mat
649       
650        nx = DimSize(mat,0)
651        ny = DimSize(mat,1)
652        nz = DimSize(mat,2)
653//      Print "redim = ",dx,dy,dz,nx,ny,nz
654       
655        KillWaves/Z M_RotatedImage
656       
657End
658
659// also look for ImageTransform operations that will allow translation of the objects.
660// then simple objects could be oriented @0,0,0, and translated to the correct position (and added)
661//
662Function FFT_AddRotatedObject(ctrlName) : ButtonControl
663        String ctrlName
664
665        Print "Not yet implemented"
666End
667
668Function FFT_ChangeMatrixValuesButton(ctrlName)
669        String ctrlName
670       
671        Execute "ChangeMatrixValues()"
672
673end
674
675Function FFT_ReplaceSolventButton(ctrlname)
676        String ctrlName
677       
678        Execute "ReplaceSolvent()"
679end
680
681Function FFT_MakeMatrixButtonProc(ctrlName) : ButtonControl
682        String ctrlName
683
684        NVAR nn=root:FFT_N
685        if(mod(nn, 2 ) ==1)             //force an even number for FFT
686                nn +=1
687        endif
688        MakeMatrix("mat",nn,nn,nn)
689End
690
691Function FFTMakeGizmoButtonProc(ctrlName) : ButtonControl
692        String ctrlName
693       
694        DoWindow/F Gizmo_VoxelMat
695        if(V_flag==0)
696                Execute "Gizmo_VoxelMat()"
697        endif
698        ColorizeGizmo()
699End
700
701Function FFTDrawSphereButtonProc(ctrlName)  : ButtonControl
702        String ctrlName
703       
704        Execute "FFTDrawSphereProc()"
705End
706
707Proc FFTDrawSphereProc(matStr,rad,xc,yc,zc,fill,periodic)
708        String matStr="mat"
709        Variable rad=25,xc=50,yc=50,zc=50,fill=10,periodic=1
710        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
711        Prompt rad,"enter real radius (A)"
712        Prompt xc,"enter the X-center"
713        Prompt yc,"enter the Y-center"
714        Prompt zc,"enter the Z-center"
715        Prompt fill,"fill SLD value"
716        Prompt periodic,"enter 1 for periodic, 0 for non-periodic fill"
717       
718        Variable grid=root:FFT_T
719       
720        if(periodic)
721                FillSphereRadiusPeriodic($matStr,grid,rad,xc,yc,zc,fill)
722        else
723                FillSphereRadius($matStr,grid,rad,xc,yc,zc,fill)
724        endif
725End
726
727Function FFTDrawZCylinderButtonProc(ctrlName)  : ButtonControl
728        String ctrlName
729       
730        Execute "FFTDrawCylinder()"
731End
732
733Proc FFTDrawCylinder(direction,matStr,rad,len,xc,yc,zc,fill)
734        String direction
735        String matStr="mat"
736        Variable rad=25,len=300,xc=50,yc=50,zc=50,fill=10
737        Prompt direction, "Direction", popup "X;Y;Z;"
738        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
739        Prompt rad,"enter real radius (A)"
740        Prompt len,"enter length (A)"
741        Prompt xc,"enter the X-center"
742        Prompt yc,"enter the Y-center"
743        Prompt zc,"enter the Z-center"
744        Prompt fill,"fill SLD value"
745       
746       
747        Variable grid=root:FFT_T
748       
749        if(cmpstr(direction,"X")==0)
750                FillXCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
751        endif
752        if(cmpstr(direction,"Y")==0)
753                FillYCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
754        endif
755        if(cmpstr(direction,"Z")==0)
756                FillZCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
757        endif
758       
759End
760
761
762Function DoTheFFT_ButtonProc(ctrlName) : ButtonControl
763        String ctrlName
764
765        Calc_IQ_FFT()
766//      Execute "DoFFT()"
767End
768
769Function FFT_PlotResultsButtonProc(ctrlName) : ButtonControl
770        String ctrlName
771
772        Variable first=0
773        DoWindow/F FFT_IQ
774        if(V_flag==0)
775                first = 1
776                Display /W=(295,44,627,302)
777                DoWindow/C FFT_IQ
778        Endif
779       
780        // append the desired data, if it's not already there
781        // FFTButton_4 = FFT            = iBin
782        // FFTButton_7a = binned = _XOP
783        // FFTButton_8a = Debye = _full
784        // FFTButton_14a = SLD = _SLD
785        // 17 = iso USANS
786        // 18 = Anisotropic USANS
787        //
788        strswitch(ctrlName)     
789                case "FFTButton_4":
790                        if(!isTraceOnGraph("iBin","FFT_IQ") && exists("iBin")==1)               //only append if it's not already there
791                                AppendToGraph /W=FFT_IQ iBin vs qBin
792                                ModifyGraph mode=4,marker=19,msize=2,rgb(iBin)=(65535,0,0)
793                        endif
794                        break
795                case "FFTButton_7a":
796                        if(!isTraceOnGraph("ival_XOP","FFT_IQ") && exists("ival_XOP")==1)               //only append if it's not already there
797                                AppendToGraph /W=FFT_IQ ival_XOP vs qval_XOP
798                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_XOP)=(1,12815,52428)
799                        endif           
800                        break
801                case "FFTButton_8a":
802                        if(!isTraceOnGraph("ival_full","FFT_IQ") && exists("ival_full")==1)             //only append if it's not already there
803                                AppendToGraph /W=FFT_IQ ival_full vs qval_full
804                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_full)=(0,0,0)
805                        endif           
806                        break
807                case "FFTButton_14a":
808                        if(!isTraceOnGraph("ival_SLD","FFT_IQ") && exists("ival_SLD")==1)               //only append if it's not already there
809                                AppendToGraph /W=FFT_IQ ival_SLD vs qval_SLD
810                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_SLD)=(2,39321,1)
811                        endif           
812                        break
813                case "FFTButton_14a":
814                        if(!isTraceOnGraph("ival_SLD","FFT_IQ") && exists("ival_SLD")==1)               //only append if it's not already there
815                                AppendToGraph /W=FFT_IQ ival_SLD vs qval_SLD
816                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_SLD)=(2,39321,1)
817                        endif           
818                        break
819                case "FFTButton_17":
820                        if(!isTraceOnGraph("FFT_iUSANS_i","FFT_IQ") && exists("FFT_iUSANS_i")==1)               //only append if it's not already there
821                                AppendToGraph /W=FFT_IQ FFT_iUSANS_i vs FFT_iUSANS_q
822                                ModifyGraph mode=4,marker=19,msize=2,rgb(FFT_iUSANS_i)=(39321,1,31457)
823                        endif           
824                        break
825                case "FFTButton_18":
826                        if(!isTraceOnGraph("FFT_aUSANS_i","FFT_IQ") && exists("FFT_aUSANS_i")==1)               //only append if it's not already there
827                                AppendToGraph /W=FFT_IQ FFT_aUSANS_i vs FFT_aUSANS_q
828                                ModifyGraph mode=4,marker=19,msize=2,rgb(FFT_aUSANS_i)=(52428,34958,1)
829                        endif           
830                        break
831                       
832                       
833        endswitch
834       
835        if(first)
836                ModifyGraph mode=4
837                ModifyGraph marker=19
838                ModifyGraph msize=2
839                ModifyGraph gaps=0
840                ModifyGraph grid=1
841                ModifyGraph log=1
842                ModifyGraph mirror=2
843                Legend
844        endif
845       
846        return(0)
847End
848
849Function isTraceOnGraph(traceStr,winStr)
850        String traceStr,winStr
851       
852        Variable isOn=0
853        String str
854        str = TraceNameList(winStr,";",1)               //only normal traces
855        isOn = strsearch(str,traceStr,0)                //is the trace there?
856        isOn = isOn == -1 ? 0 : 1                       // return 0 if not there, 1 if there
857
858        return(isOn)
859End
860
861Function FFTEraseMatrixButtonProc(ctrlName) : ButtonControl
862        String ctrlName
863       
864        Wave mat=root:mat
865        FastOp mat=0
866        return(0)
867End
868
869Function FFTFillSolventMatrixProc(ctrlName) : ButtonControl
870        String ctrlName
871       
872        Wave mat=root:mat
873        NVAR val=root:FFT_SolventSLD
874        FastOp mat=(val)
875        return(0)
876End
877
878Function FFT_AltiSpheresButtonProc(ctrlName) : ButtonControl
879        String ctrlName
880
881        Execute "DoSpheresCalcFFTPanel()"
882End
883
884Function FFT_BinnedSpheresButtonProc(ctrlName) : ButtonControl
885        String ctrlName
886
887        Execute "DoBinnedSpheresCalcFFTPanel()"
888End
889
890Function FFT_BinnedSLDButtonProc(ctrlName) : ButtonControl
891        String ctrlName
892
893        Execute "DoBinnedSLDCalcFFTPanel()"
894End
895
896Function FFT_Iso2USANS(ctrlName) : ButtonControl
897        String ctrlName
898
899        Execute "Isotropic_FFT_to_USANS()"
900        FFT_PlotResultsButtonProc(ctrlName)
901End
902
903Function FFT_Aniso2USANS(ctrlName) : ButtonControl
904        String ctrlName
905
906        Execute "Anisotropic_FFT_to_USANS()"
907        FFT_PlotResultsButtonProc(ctrlName)
908End
909
910
911
912
913/////UTILITIES
914
915// inverts the values in the matrix 0<->1
916Function InvertMatrixFill(mat)
917        Wave mat
918       
919        mat = (mat==1) ? 0 : 1
920End
921
922// replaces specified values
923Proc ChangeMatrixValues(old,new)
924        Variable old,new
925       
926        mat = (mat==old) ? new : mat
927        // sequence of steps to get the gizmo to update the display correctly
928        RemoveFromGizmo/N=Gizmo_VoxelMat object=Voxelgram0
929        RemoveFromGizmo/N=Gizmo_VoxelMat displayItem=axes0
930        AppendToGizmo/N=Gizmo_VoxelMat voxelgram=root:mat,name=voxelgram0
931        ModifyGizmo/N=Gizmo_VoxelMat  setDisplayList=-1, object=voxelgram0
932        ModifyGizmo/N=Gizmo_VoxelMat  setDisplayList=-1, object=axes0           //so that the axes are drawn last
933        ModifyGizmo ModifyObject=voxelgram0 property={ pointSize,3}
934
935       
936        ColorizeGizmo()
937End
938
939
940
941// replaces the solvent value and updates the global
942Proc ReplaceSolvent(newSolv)
943        Variable newSolv
944       
945        Variable solv = root:FFT_SolventSLD
946       
947        mat = (mat==solv) ? newSolv : mat
948       
949        root:FFT_SolventSLD = newSolv
950
951// sequence of steps to get the gizmo to update the display correctly
952        RemoveFromGizmo/N=Gizmo_VoxelMat object=Voxelgram0
953        RemoveFromGizmo/N=Gizmo_VoxelMat displayItem=axes0
954        AppendToGizmo/N=Gizmo_VoxelMat voxelgram=root:mat,name=voxelgram0
955        ModifyGizmo/N=Gizmo_VoxelMat  setDisplayList=-1, object=voxelgram0
956        ModifyGizmo/N=Gizmo_VoxelMat  setDisplayList=-1, object=axes0           //so that the axes are drawn last
957        ModifyGizmo ModifyObject=voxelgram0 property={ pointSize,3}
958
959
960        ColorizeGizmo()
961End
962
963
964//overwrites any existing matrix
965// matrix is byte, to save space
966//forces the name to be "mat"
967//
968// switched to signed byte
969//
970Function MakeMatrix(nam,xd,yd,zd)
971        String nam
972        Variable xd,yd,zd
973       
974        nam="mat"
975        Print "Matrix has been created and named \"mat\""
976//      Make/O/U/B/N=(xd,yd,zd) $nam
977        Make/O/B/N=(xd,yd,zd) $nam
978End
979
980
981// calculate the average center-to-center distance between points
982// assumes evenly spaced on a cubic grid
983//
984Proc Center_to_Center(np)
985        Variable np = 100
986       
987        Variable Nedge = root:FFT_N
988        Variable Tscale = root:FFT_T
989       
990        Variable davg
991       
992        davg = (Nedge*Tscale)^3 / np
993        davg = davg^(1/3)
994        Print "Average separation (A) = ",davg
995       
996End
997
998// calculate the number of points required on a grid
999// to yield a given average center-to-center distance between points
1000// assumes evenly spaced on a cubic grid
1001//
1002// davg and Tscale are the same units, Angstroms
1003//
1004Proc Davg_to_Np(davg)
1005        Variable davg = 400
1006       
1007        Variable Nedge = root:FFT_N
1008        Variable Tscale = root:FFT_T
1009       
1010        Variable np
1011       
1012        np = (Nedge*Tscale)^3 / davg^3
1013        Print "Number of points required = ",np
1014       
1015End
1016
1017
1018
1019// The matrix is not necessarily 0|1, this reports the number of filled voxels
1020// - needed to estimate the time required for the AltiVec_Spheres calculation
1021Proc NumberOfPoints()
1022       
1023        Print "Number of points = ",NumberOfPoints_Occ(root:mat)
1024        Print "Fraction occupied = ",VolumeFraction_Occ(root:mat)
1025        Print "Overall Cube Edge [A] = ",root:FFT_T * root:FFT_N
1026        Print "Found values in matrix = ",ListOfValues(root:mat)
1027       
1028End
1029
1030Function NumberOfPoints_Occ(m)
1031        Wave m
1032       
1033        Variable num = NonZeroValues(m)
1034
1035        return(num)
1036End
1037
1038
1039Function VolumeFraction_Occ(m)
1040        Wave m
1041       
1042        Variable num = NonZeroValues(m)
1043        Variable dim = DimSize(m,0)
1044
1045        return(num/dim^3)
1046End
1047
1048//it's a byte wave, so I can't use the trick of setting the zero values to NaN
1049// since NaN can't be expressed as byte. So make it binary 0|1
1050//
1051// - the assumption here is that the defined solvent value is not part of the
1052// particle volume.
1053//
1054Function NonZeroValues(m)
1055        Wave m
1056       
1057        Variable num
1058        NVAR val = root:FFT_SolventSLD
1059       
1060//      Variable t1=ticks,dim
1061//      dim = DimSize(m, 0 )            //assume NxNxN
1062        Duplicate/O m,mz
1063       
1064        MultiThread mz = (m[p][q][r] != val) ? 1 : 0
1065       
1066        WaveStats/Q/M=1 mz              // NaN and Inf are not reported in V_npnts
1067       
1068        num = V_npnts*V_avg
1069       
1070        KillWaves/Z mz
1071//      Print "Number of points = ",num
1072//      Print "Fraction occupied = ",num/V_npnts
1073       
1074//      Print "time = ",(ticks-t1)/60.15
1075        return(num)
1076End
1077
1078//
1079// return a list of the different values of the voxels in the matrix
1080//
1081Function/S ListOfValues(m)
1082        Wave m
1083       
1084        String list=""
1085        Variable done
1086
1087        Duplicate/O m,mz
1088       
1089        done=0
1090        do
1091                WaveStats/Q/M=1 mz              // NaN and Inf are not reported in V_npnts
1092                if(V_max == V_min)
1093                        list += num2str(V_min) + ";"
1094                        done = 1
1095                else
1096                        list += num2str(V_max) + ";"
1097                        MultiThread mz = mz[p][q][r] == V_max ? V_min : mz[p][q][r]             // replace the max with min                     
1098                endif
1099        while(!done)   
1100       
1101//      Print "Found values in matrix = ",list
1102//      KillWaves/Z mz
1103
1104        return(list)
1105End
1106
1107
1108
1109// returns estimate in seconds
1110//
1111//              updated for quad-core iMac, 2010
1112//
1113// type 3 = binned distances and SLD
1114// type 2 = binned distances
1115// type 1 is rather meaningless
1116// type 0 = is the Debye, double sum (multithreaded)
1117//
1118//
1119// - types 2 and 3, the binned methods are inherently AAO, so there is no significant
1120// dependence on the number of q-values (unless it's  >> 1000). So values reported are
1121// for 100 q-values.
1122//
1123// There is a function Timing_Method(type) in FFT_ConcentratedSpheres.ipf to automate some of this
1124//
1125Function EstimatedTime(nx,nq,type)
1126        Variable nx,nq,type
1127       
1128        Variable est=0
1129       
1130        if(type==0)             // "full" XOP
1131                est = 0.4 + 1.622e-5*nx+4.86e-7*nx^2
1132                est /= 100
1133                est *= nq
1134        endif
1135       
1136        if(type==1)             //Igor function (much slower, 9x)
1137                est = (nx^2)*0.0000517          //empirical values (seconds) for igor code, 100 q-values
1138                est /= 100                                      // per q now...
1139                est *= nq
1140        endif
1141       
1142        if(type==2)             //binned distances, slow parts XOPed
1143                est = 0.0680 + 2.94e-6*nx+4.51e-9*nx^2                  //with threading
1144               
1145//              est /= 100                                      // per q now...
1146//              est *= nq
1147        endif
1148       
1149        if(type==3)             //binned distances AND sld, slow parts XOPed
1150                est = 0.576 + 4.22e-7*nx+1.76e-8*nx^2           //with threading
1151               
1152//              est /= 100                                      // per q now...
1153//              est *= nq
1154        endif
1155       
1156        return(est)
1157End
1158
1159////////////// my functions to rotate the matrix in a XYZ coordinate system
1160//
1161//
1162// definitely not generic, is expecting NxNxN volume
1163//
1164// not very friendly in that it "clips" anything that rotates out of the volume
1165//
1166// does translate the center of the box to 000, rotates, then translates back
1167//
1168// friendly in the sense that the rotated matrix is the same size as the original.
1169//  --this is important for my final application (FFT)
1170//
1171// does no interpolation of values, so be sure to keep a copy of the original
1172// -- multiple rotation steps are going to make a mess of things.
1173//
1174// The multi axis rotation is done as one step, and probably violates every conventional
1175//  coordinate system. The rotation is applied as RxRyRz, but this could easily be changed
1176//
1177// I just want it to be correct, so speed was not an issue.
1178// -- it's nested for loops.
1179// -- it's working with the full matrix, even when 99% is empty.
1180//
1181// 20 NOV 2013 SRK
1182//
1183//
1184
1185//
1186// mat is the input volume
1187// rotVol is the output rotated volume
1188Function XYZRotate(angleX,angleY,angleZ)
1189        Variable angleX,angleY,angleZ
1190       
1191        NVAR FFT_N=root:FFT_N
1192        WAVE mat=root:mat
1193        Variable dist=FFT_N/2
1194
1195
1196// convert the NxNxN into 3xN xyz locations + wave of "w" values named "values"
1197        fVolumeToXYZTriplet(mat,"trip")
1198        Wave trip=root:trip
1199
1200// translate to get the center of the xyz values to 0,0,0       
1201        fTranslateCoordinate(trip,dist)         //subtracts dist
1202
1203// do the rotation as a matrix multiplication   
1204// putting zero is no rotation around that axis
1205// the triplet wave "trip" is overwritten with the output
1206        DoRotation(trip,angleX,angleY,angleZ)
1207//      Wave rotated=root:rotated
1208
1209// translate back to a 0->N based coordinate
1210//      fTranslateCoordinate(rotated,-dist)
1211        fTranslateCoordinate(trip,-dist)
1212        Wave values=root:values
1213
1214// convert the triplet back to a volume
1215// this CLIPS anything that has rotated out of the NxNxN volume
1216//      fXYZTripletToVolume(rotated,values,"rotVol",FFT_N)
1217        fXYZTripletToVolume(trip,values,"rotVol",FFT_N)
1218
1219// clean up by killng the extra waves that were generated
1220//
1221        KillWaves/Z trip,rotated,values
1222       
1223        Wave rotVol=root:rotVol
1224        mat=rotVol
1225       
1226        return(0)
1227End
1228
1229
1230
1231Function fVolumeToXYZTriplet(matrixWave, outputName)
1232        Wave matrixWave
1233        String outputName       
1234       
1235        Variable dimx=DimSize(matrixWave,0)
1236        Variable dimy=DimSize(matrixWave,1)
1237        Variable dimz=DimSize(matrixWave,2)
1238        Variable rows=dimx*dimy*dimz
1239        Make/O/N=(3,rows) $outputName
1240        Make/O/N=(rows) values
1241        WAVE TripletWave= $outputName
1242        Wave values=values
1243       
1244
1245        Variable ii,jj,kk,count=0
1246        Variable xVal,yVal,zval
1247        for(kk=0;kk<dimz;kk+=1)                 // kk is z (layer)
1248                zval=kk
1249                for(jj=0;jj<dimy;jj+=1)         // jj is y (column)
1250                        yVal=jj
1251                        for(ii=0;ii<dimx;ii+=1) // ii is x (row)
1252                                xVal=ii
1253                                TripletWave[0][count]=xVal
1254                                TripletWave[1][count]=yVal
1255                                TripletWave[2][count]=zval
1256                                values[count]=matrixWave[ii][jj][kk] // value at [row][col][lay]
1257                                count+=1
1258                        endfor
1259                endfor
1260        endfor
1261       
1262        return(0)
1263End
1264
1265
1266Function fXYZTripletToVolume(triplet, values, outputName, outputDim)
1267        Wave triplet,values
1268        String outputName
1269        Variable outputDim
1270       
1271        Variable numPt=DimSize(triplet,1)
1272
1273        Variable num = outputDim
1274       
1275        Make/O/B/N=(num,num,num) $outputName
1276        WAVE newVol= $outputName
1277
1278        FastOp newVol = 0
1279       
1280        Variable ii,jj,kk,count=0
1281        Variable xVal,yVal,zval
1282        Variable xOK, yOK, zOK
1283
1284       
1285        for(ii=0;ii<numPt;ii+=1)
1286                xval = round(triplet[0][ii])
1287                yval = round(triplet[1][ii])
1288                zval = round(triplet[2][ii])
1289               
1290                // round and keep in bounds (returns truth)
1291                xOK = inRange(xval,0,num-1)
1292                yOK = inRange(yval,0,num-1)
1293                zOK = inRange(zval,0,num-1)
1294               
1295                if(xOK && yOK && zOK)
1296                        newVol[xval][yval][zval] = values[ii]
1297                endif
1298                       
1299        endfor
1300
1301
1302        return(0)
1303End
1304
1305// if val < lo or > hi, bad val
1306// fi both of these pass, pt is OK
1307ThreadSafe Static Function inRange(val, lo, hi)
1308        Variable val, lo, hi
1309 
1310        if(val < lo)
1311                return (0)
1312        endif
1313        if(val > hi)
1314                return (0)
1315        endif
1316       
1317        return(1)
1318 
1319End
1320
1321// Rotation is applied in the order Rx Ry Rz
1322Function DoRotation(triplet,angleX,angleY,angleZ)
1323        Wave triplet
1324        Variable angleX,angleY,angleZ
1325       
1326        Variable thetaX,thetaY,thetaZ
1327        thetaX = angleX*pi/180          // convert degrees to radians
1328        thetaY = angleY*pi/180          // convert degrees to radians
1329        thetaZ = angleZ*pi/180          // convert degrees to radians
1330       
1331        Make/O/D/N=(3,3) Rx,Ry,Rz
1332        Rx=0
1333        Ry=0
1334        Rz=0
1335       
1336        Rx[0][0] = 1
1337        Rx[1][1] = cos(thetaX)
1338        Rx[1][2] = -sin(thetaX)
1339        Rx[2][1] = sin(thetaX)
1340        Rx[2][2] = cos(thetaX)
1341       
1342        Ry[0][0] = cos(thetaY)
1343        Ry[0][2] = sin(thetaY)
1344        Ry[1][1] = 1
1345        Ry[2][0] = -sin(thetaY)
1346        Ry[2][2] = cos(thetaY)
1347       
1348        Rz[0][0] = cos(thetaZ)
1349        Rz[0][1] = -sin(thetaZ)
1350        Rz[1][0] = sin(thetaZ)
1351        Rz[1][1] = cos(thetaZ)
1352        Rz[2][2] = 1   
1353       
1354       
1355//      MatrixOp/O rotated = Rx x Ry x Rz x triplet
1356        MatrixOp/O triplet = Rx x Ry x Rz x triplet
1357       
1358
1359        return(0)
1360end
1361
1362
1363// the rotation matrix, as I copied from wikipedia, is a rotation around (0,0,0), not
1364// the center of the gizmo plot
1365Function fTranslateCoordinate(trip,dist)
1366        Wave trip
1367        Variable dist
1368       
1369//      MatrixOp/O trip = trip - dist
1370        trip = trip - dist
1371       
1372        return(0)
1373End
1374
Note: See TracBrowser for help on using the repository browser.