Ignore:
Timestamp:
Nov 19, 2013 9:50:29 AM (9 years ago)
Author:
srkline
Message:

Adding changes to the real space modeling that allow calculation of the USANS (slit-smeared) intensity from an anisotropic structure that has been generated using a real-space construction. The detector plane from the FFT result is converted to USANS, on absolute scale. Will be extended as it gets some use, but it is functional and is found from the "3D Examples" Macros submenu. Documentation is on the way too.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Cubes.ipf

    r798 r925  
    33 
    44// FFT and DebyeSphere calculations are now properly normalized to the conditions of: 
    5 // deltaRho = 1e-6 (set as a global) 
     5// deltaRho = 1e-7 (set as a global) 
    66// volume fraction = occupied fraction 
    77// 
     
    518518        ImageTransform/P=0 getPlane M_VolumeTranspose 
    519519        WAVE M_ImagePlane 
     520         
    520521        Duplicate/O M_ImagePlane detPlane 
    521         Duplicate/O M_ImagePlane logP 
    522          
     522         
     523        //get the scaling right 
    523524        NVAR FFT_N=root:FFT_N 
     525        NVAR FFT_T=root:FFT_T 
     526        NVAR delRho=root:FFT_delRho 
     527                 
     528        detPlane *= delRho*delRho 
     529        detPlane *= 1e8 
     530        detPlane *= (FFT_T/FFT_N)^3 
     531        Duplicate/O detPlane logP 
     532         
     533 
    524534        detPlane[FFT_N/2][FFT_N/2] = NaN                        //hopefully this trims out the singularity at the center 
    525535        logP = log(detPlane) 
     
    563573        return(0) 
    564574end 
     575 
     576 
     577////////// Routines to take the FFT result and convert it to USANS data 
     578// 
     579// - this is done by taking the 2D detector image and summing the columns (of corrrect Dq resolution) 
     580// 
     581// OCT/NOV 2013 
     582 
     583 
     584// 
     585// This is for when an anisotropic 2D result is desired 
     586// 
     587// these steps calculate the slit-smeared USANS from the FFT 
     588// 
     589// I think that the scaling of the USANS result is correct now 
     590// 
     591// -- for USANS length scales, need to set N ~ 256 and T ~ 500 to get the q-range into 
     592//    the proper range (deltaQ and qMax). 256*500 = 128,000 A = 12.8 microns per edge of the box 
     593// 
     594// 1D output is named "FFT_aUSANS_" (i and q) to signify "Anisotropic" USANS 
     595// 
     596// 
     597// SRK 18 OCT 2013 
     598// 
     599Function Anisotropic_FFT_to_USANS() 
     600 
     601//figure out the proper Qmin, Qmax for the N and T 
     602//want deltaQ = 2e-5 to match the USANS instrumental resolution 
     603 
     604        Variable maxQx,Qy,num,deltaQ 
     605        String str,qmStr 
     606         
     607        NVAR FFT_N = root:FFT_N 
     608        NVAR FFT_T = root:FFT_T 
     609        NVAR FFT_QMaxReal = root:FFT_QMaxReal 
     610        NVAR delRho = root:FFT_delRho 
     611         
     612        deltaQ = 2e-5 
     613        num = round(FFT_QMaxReal/deltaQ) 
     614        if(num>1000) 
     615                Print "num = ",num 
     616                Abort "N and T are not set correctly for USANS. Reset these (N>=256 and T>=500) for an appropriate range." 
     617        endif 
     618         
     619//      print "num = ",num 
     620         
     621        Qmstr = num2str(FFT_QMaxReal-2*deltaQ)          //using FFT_QMaxReal gives NaN at edges of interp2Dslice 
     622        str = "FakeQxQy(0,"+qmstr+",-"+qmstr+","+qmstr+","+num2str(num)+",\"fakeHalf\",1,0)"            //overwrites the data folder, no plot 
     623        Execute str 
     624 
     625// do the FFT 
     626        Execute "DoFFT()" 
     627// get the detector slice -- the slice is already on absolute scale 
     628        FFT_Get2DSlice("") 
     629// interpolate to the USANS-sized detector 
     630        Interpolate2DSliceToData("fakeHalf") 
     631        Wave interp2Dslice = $("root:interp2Dslice") 
     632 
     633// sum the rows  
     634        MatrixOp/O FFT_aUSANS_i=sumRows(interp2Dslice)  //automatically generates the destination 
     635         
     636        Make/O/D/N=(num) FFT_aUSANS_q 
     637        FFT_aUSANS_q = FFT_QMaxReal/num*x 
     638         
     639        //now get the scaling correct 
     640        // q-integration (rectangular), matrixOp simply summed, so I need to multiply by dy 
     641        FFT_aUSANS_i *= dimdelta(interp2dslice,1) 
     642 
     643        FFT_aUSANS_i *= 4               //why the factor of 4??? 
     644         
     645        // not needed -- occupancy and the normal FFT scaling are already done -- the slice is already on absolute scale 
     646        Variable occ = VolumeFraction_Occ(root:mat) 
     647        print "Fraction occupied = ",occ 
     648//      FFT_USANS_i *= occ 
     649 
     650        return(0) 
     651End 
     652 
     653// 
     654// 
     655//// this is for when an isotropic 2D result is desired 
     656// 
     657// FIRST -- do the FFT - this calculates the 1D averaged result 
     658// SECOND -- get the 2D slice. this sets up a 2D matrix appropriately scaled to the q-space. 
     659// then - // 
     660// --Take the 1D FFT result - which is orientationally averaged, and convert it to a 2D plot 
     661// use the detPlane slice to get the 2D scaling - otherwise there's no need for detPlane 
     662// and I could just ask for qMin qMax in XY directions. 
     663// 
     664// --(optional) after this is called -- call fDoBinning_Scaled2D(avg2d,0) to confirm that it returns 
     665// back to the correct, orientationally averaged 1D 
     666// 
     667// --then to plot -- DisplayFFT_to2DAverage() 
     668// 
     669// 
     670// finally, you can use the 2D result to either interpolate to a real 2D range, or more useful, convert 
     671// the 2D results to a USANS result. Then the case of an oddly shaped object can be converted to USANS. Previously 
     672// only the FFT slice (which is directional) could be converted to USANS. This limited the calculations to only 
     673// symmetric structures, which was not very useful. 
     674// 
     675// 1D output is named "FFT_iUSANS_" (i and q) to signify "isotropic" USANS 
     676// 
     677Function Isotropic_FFT_to_USANS() 
     678 
     679////////////// 
     680        Variable maxQx,Qy,num,deltaQ 
     681        String str,qmStr 
     682         
     683        NVAR FFT_N = root:FFT_N 
     684        NVAR FFT_T = root:FFT_T 
     685        NVAR FFT_QMaxReal = root:FFT_QMaxReal 
     686        NVAR delRho = root:FFT_delRho 
     687         
     688        deltaQ = 2e-5 
     689        num = round(FFT_QMaxReal/deltaQ) 
     690        if(num>1000) 
     691                Print "num = ",num 
     692                Abort "N and T are not set correctly for USANS. Reset these (N>=256 and T>=500) for an appropriate range." 
     693        endif 
     694         
     695//      print "num = ",num 
     696         
     697        Qmstr = num2str(FFT_QMaxReal-2*deltaQ)          //using FFT_QMaxReal gives NaN at edges of interp2Dslice 
     698        str = "FakeQxQy(0,"+qmstr+",-"+qmstr+","+qmstr+","+num2str(num)+",\"fakeHalf\",1,0)"            //overwrites the data folder, no plot 
     699        Execute str 
     700 
     701// do the FFT 
     702        Execute "DoFFT()" 
     703// get the detector slice -- the slice is already on absolute scale 
     704        FFT_Get2DSlice("") 
     705 
     706 
     707         
     708// now get the isotropically averaged I(Q) = iBin        
     709 
     710        WAVE iBin = iBin 
     711        WAVE qBin = qBin 
     712 
     713        Wave detPlane = root:detPlane           //just for the 2D scaling 
     714         
     715        Duplicate/O detPlane avg2d 
     716        avg2d = 0 
     717         
     718        Variable rowOff,colOff,rowDel,colDel 
     719        rowOff = DimOffset(avg2d,0) 
     720        colOff = DimOffset(avg2d,1) 
     721        rowDel = DimDelta(avg2d,0) 
     722        colDel = DimDelta(avg2d,1) 
     723 
     724// fill the 2D plane (qx,qy) with I(q) from the FFT 
     725         
     726        avg2d = interp( sqrt((rowOff + p*rowDel)^2 + (colOff + q*colDel)^2), qBin, iBin ) 
     727//      avg2d = interp( sqrt((rowOff + p*rowDel)^2 + (colOff + q*colDel)^2), qval_XOP, ival_XOP ) 
     728         
     729        Duplicate/O avg2d logAvg2d 
     730        logAvg2d = log(avg2d) 
     731 
     732         
     733 
     734         
     735// interpolate to the USANS-sized detector, this is equivalent to: 
     736//      Interpolate2DSliceToData("fakeHalf")            //except using avg2d rather than "calc = detPlane" (the anisotropic result) 
     737        String folderStr = "fakeHalf" 
     738        Wave data = $("root:"+folderStr+":"+folderStr+"_lin") 
     739         
     740        Duplicate/O data,interp2DSlice                  //keeps the scaling of the data 
     741         
     742        rowOff = DimOffset(data,0) 
     743        colOff = DimOffset(data,1) 
     744        rowDel = DimDelta(data,0) 
     745        colDel = DimDelta(data,1) 
     746 
     747        interp2DSlice = Interp2D (avg2d, rowOff + p*rowDel, colOff + q*colDel ) 
     748 
     749        Duplicate/O interp2DSlice interp2DSlice_log 
     750        interp2DSlice_log = log(interp2DSlice) 
     751 
     752 
     753 
     754// sum the rows  
     755        MatrixOp/O FFT_iUSANS_i=sumRows(interp2Dslice)  //automatically generates the destination 
     756         
     757        Make/O/D/N=(num) FFT_iUSANS_q 
     758        FFT_iUSANS_q = FFT_QMaxReal/num*x 
     759         
     760        //now get the scaling correct 
     761        // q-integration (rectangular), matrixOp simply summed, so I need to multiply by dy 
     762        FFT_iUSANS_i *= dimdelta(interp2dslice,1) 
     763 
     764        FFT_iUSANS_i *= 4               //why the factor of 4??? 
     765 
     766        Execute "DisplayFFT_to2DAverage()" 
     767 
     768// if you want to scale this to match a model calculation, you'll need the volume fraction 
     769        Variable occ = VolumeFraction_Occ(root:mat) 
     770        print "Fraction occupied = ",occ 
     771                 
     772        return(0) 
     773End 
     774 
     775Proc DisplayFFT_to2DAverage() 
     776 
     777        DoWindow FFT_Avg2D 
     778        if(V_flag == 0) 
     779                PauseUpdate; Silent 1           // building window... 
     780                Display /W=(1038,44,1404,403) 
     781                DoWindow/C FFT_Avg2D 
     782                AppendImage/T avg2d 
     783                ModifyImage avg2d ctab= {*,*,YellowHot,0} 
     784                ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 
     785                ModifyGraph mirror=2 
     786                ModifyGraph nticks=4 
     787                ModifyGraph minor=1 
     788                ModifyGraph fSize=9 
     789                ModifyGraph standoff=0 
     790                ModifyGraph tkLblRot(left)=90 
     791                ModifyGraph btLen=3 
     792                ModifyGraph tlOffset=-2 
     793                SetAxis/A/R left 
     794        endif 
     795         
     796        DoWindow FFT_logAvg2D 
     797        if(V_flag == 0) 
     798                Display /W=(1038,44,1404,403) 
     799                DoWindow/C FFT_logAvg2D 
     800                AppendImage/T logavg2d 
     801                ModifyImage logavg2d ctab= {*,*,YellowHot,0} 
     802                ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 
     803                ModifyGraph mirror=2 
     804                ModifyGraph nticks=4 
     805                ModifyGraph minor=1 
     806                ModifyGraph fSize=9 
     807                ModifyGraph standoff=0 
     808                ModifyGraph tkLblRot(left)=90 
     809                ModifyGraph btLen=3 
     810                ModifyGraph tlOffset=-2 
     811                SetAxis/A/R left 
     812        endif 
     813         
     814End 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Cubes_Includes.ipf

    r847 r925  
    4040                "Put X-axis Cylinders Hexagonal Grid",PutXAxisCylindersHexagonal() 
    4141                "Core-shell Cylinders Hex Grid",PutXAxisCoreShellCyl_HexGrid() 
     42                "-" 
     43                "Anisotropic_FFT_to_USANS" 
     44                "Isotropic_FFT_to_USANS" 
    4245        end 
    4346        Submenu "Matrix Viewing" 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Panel.ipf

    r923 r925  
    7979        SetVariable FFTSetVar_0,limits={50,512,2},value= FFT_N,live= 1 
    8080        SetVariable FFTSetVar_1,pos={7,33},size={150,15},title="Length per Cell (T)" 
    81         SetVariable FFTSetVar_1,limits={1,500,0.2},value= FFT_T,live= 1 
     81        SetVariable FFTSetVar_1,limits={1,5000,0.2},value= FFT_T,live= 1 
    8282        SetVariable FFTSetVar_2,pos={183,26},size={120,15},title="Real Qmax" 
    8383        SetVariable FFTSetVar_2,limits={0,0,0},value= FFT_QmaxReal,noedit= 1,live= 1 
     
    210210 
    211211// utility function that reads the wave note from a binary file 
     212// where did I get this from? Igor Exchange? I didn't write this myself... 
    212213// 
    213214Function/S LoadNoteFunc(PName,FName[,FileRef]) 
     
    293294                        // click code here 
    294295                        String folderStr="" 
    295                         Prompt folderStr, "Pick a 2D data folder"               // Set prompt for x param 
     296                        Prompt folderStr, "Pick a 2D data folder",popup,W_DataSetPopupList()            // Set prompt for x param 
    296297                        DoPrompt "Enter data folder", folderStr 
    297298                        if (V_Flag || strlen(folderStr) == 0) 
     
    448449        Prompt degree, "Degrees of rotation around Z-axis:" 
    449450//      Prompt sense, "Direction of rotation:",popup,"CW;CCW;" 
    450         DoPrompt "Enter paramters for rotation", degree 
     451        DoPrompt "Enter parameters for rotation", degree 
    451452         
    452453        if (V_Flag) 
     
    457458        return(0) 
    458459End 
     460 
    459461// note that the rotation is not perfect. if the rotation produces an 
    460462// odd number of pixels, then the object will "walk" one pixel. Also, some 
     
    469471         
    470472        Variable fill=0 
     473        NVAR solventSLD = root:FFT_SolventSLD 
     474        fill = solventSLD 
    471475         
    472476        WAVE mat = mat 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_VoxelDisplay.ipf

    r924 r925  
    1313        ModifyGizmo startRecMacro 
    1414        AppendToGizmo voxelgram=root:mat,name=voxelgram0 
    15         ModifyGizmo ModifyObject=voxelgram0 property={ valueRGBA,0,10,1.5259e-05,0.195544,0.8,0.3} 
     15        ModifyGizmo ModifyObject=voxelgram0 property={ valueRGBA,0,10,1.5259e-05,0.195544,0.8,0.15} 
    1616        ModifyGizmo ModifyObject=voxelgram0 property={ mode,0} 
    1717        ModifyGizmo ModifyObject=voxelgram0 property={ pointSize,3} 
     
    6767        Variable ii,ind 
    6868//// now loop through and colorize the different parts 
    69         Make/O/D/N=5 rr,gg,bb,aa 
    70         rr = {0,1,4.57771e-05,0.8,0} 
    71         gg = {0.244434,0,0.6,0.799954,0} 
    72         bb = {1,0,0,0,0} 
    73         aa = 0.15                               //alpha = transparency 
     69        Make/O/D/N=5 rr_vox,gg_vox,bb_vox,aa_vox 
     70        rr_vox = {0,1,4.57771e-05,0.8,0} 
     71        gg_vox = {0.244434,0,0.6,0.799954,0} 
     72        bb_vox = {1,0,0,0,0} 
     73        aa_vox = 0.15                           //alpha = transparency 
    7474         
    7575        for(ii=0;ii<numDiffSLD;ii+=1) 
     
    8080                str += num2istr(ind)+"," 
    8181                str += StringFromList(ind, listStr ,";")   +"," 
    82                 str += num2str(rr[ind])+"," 
    83                 str += num2str(gg[ind])+"," 
    84                 str += num2str(bb[ind])+"," 
    85                 str += num2str(aa[ind])+"}" 
     82                str += num2str(rr_vox[ind])+"," 
     83                str += num2str(gg_vox[ind])+"," 
     84                str += num2str(bb_vox[ind])+"," 
     85                str += num2str(aa_vox[ind])+"}" 
    8686                         
    8787//              ModifyGizmo/N=Gizmo_VoxelMat ModifyObject=voxelgram0 property={ valueRGBA,ii,SLD_id[ii],rr[ii],gg[ii],bb[ii],aa[ii]} 
Note: See TracChangeset for help on using the changeset viewer.