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

Last change on this file since 971 was 971, checked in by srkline, 7 years ago

Filling in missing "dummy" functions when XML xop is not present.

Added new V_ procedures for the beginning of the reduction steps - initialization, main panel, preferences, and a menu. These are lifted directly from the SANS routines and much still needs to be modified.
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="%d"
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] == V_max ? V_min : mz[p][q]           // 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.