Ignore:
Timestamp:
Jan 7, 2009 1:17:09 PM (14 years ago)
Author:
srkline
Message:

Final fixes, w/Charles Dewhurst, to finally calculate Qx and Qy correctly as vector components of the vector Q, as it is projected onto the XY planar detector. An additional function has been added to calculate Qz as well, but is not written out anywhere. This calculation only has an effect at the highest q-values, for QxQy? export. At a Q ~ 0.45, the correct Qx and Qy are about 2% different than the old calculation. |Q| has always been correct.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/RawWindowHook.ipf

    r459 r467  
    271271//returned magnitude of Q is in 1/Angstroms 
    272272// 
    273 // repaired the dumb error of incorrect qx and qy calculation 3 dec 08 SRK (Lionel...) 
    274 // 
    275273Function CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    276274        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
    277275         
    278         Variable dx,dy,qval,theta,dist 
     276        Variable dx,dy,qval,two_theta,dist 
    279277         
    280278        Variable pixSizeX=pixSize 
     
    286284        dist = sqrt(dx^2 + dy^2) 
    287285         
    288         theta = atan(dist/sdd) 
    289  
    290         qval = 4*Pi/lam*sin(theta/2) 
     286        two_theta = atan(dist/sdd) 
     287 
     288        qval = 4*Pi/lam*sin(two_theta/2) 
    291289         
    292290        return qval 
     
    301299//wavelength is in Angstroms 
    302300// 
    303 // repaired the dumb error of incorrect qx and qy calculation 3 dec 08 SRK (Lionel...) 
     301// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) 
     302// now properly accounts for qz 
    304303// 
    305304Function CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    306305        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
    307306 
    308         Variable qx,qval,phi,dx,dy 
     307        Variable qx,qval,phi,dx,dy,dist,two_theta 
    309308         
    310309        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     
    315314        phi = FindPhi(dx,dy) 
    316315         
    317         qx = qval*cos(phi) 
     316        //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
     317        dist = sqrt(dx^2 + dy^2) 
     318        two_theta = atan(dist/sdd) 
     319 
     320        qx = qval*cos(two_theta/2)*cos(phi) 
    318321         
    319322        return qx 
     
    327330//wavelength is in Angstroms 
    328331// 
    329 // repaired the dumb error of incorrect qx and qy calculation 3 dec 08 SRK (Lionel...) 
     332// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) 
     333// now properly accounts for qz 
    330334// 
    331335Function CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    332336        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
    333337         
    334         Variable dy,qval,dx,phi,qy 
     338        Variable dy,qval,dx,phi,qy,dist,two_theta 
    335339         
    336340        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     
    341345        phi = FindPhi(dx,dy) 
    342346         
    343         qy = qval*sin(phi) 
     347        //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
     348        dist = sqrt(dx^2 + dy^2) 
     349        two_theta = atan(dist/sdd) 
     350         
     351        qy = qval*cos(two_theta/2)*sin(phi) 
    344352         
    345353        return qy 
    346354End 
    347355 
     356//calculates just the z-component of the q-vector, not measured on the detector 
     357//input/output is the same as CalcQval() 
     358//ALL inputs are in detector coordinates 
     359//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
     360//sdd is in meters 
     361//wavelength is in Angstroms 
     362// 
     363// not actually used, but here for completeness if anyone asks 
     364// 
     365Function CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     366        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
     367         
     368        Variable dy,qval,dx,phi,qz,dist,two_theta 
     369         
     370        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     371         
     372        sdd *=100               //convert to cm 
     373         
     374        //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
     375        dx = (xaxval - xctr)*pixSize            //delta x in cm 
     376        dy = (yaxval - yctr)*pixSize            //delta y in cm 
     377        dist = sqrt(dx^2 + dy^2) 
     378        two_theta = atan(dist/sdd) 
     379         
     380        qz = qval*sin(two_theta/2) 
     381         
     382        return qz 
     383End 
    348384 
    349385//phi is defined from +x axis, proceeding CCW around [0,2Pi] 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/WriteQIS.ipf

    r459 r467  
    736736//      Endif 
    737737         
    738         Duplicate/O data,qx_val,qy_val,z_val 
     738        Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val 
    739739         
    740740//      Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val 
     
    743743        qx_val = CalcQx(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)          //+1 converts to detector coordinate system 
    744744        qy_val = CalcQy(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    745  
     745         
    746746        Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val 
     747         
     748        // un-comment these if you want to write out qz_val and qval too, then use the proper save command 
     749//      qval = CalcQval(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     750//      qz_val = CalcQz(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     751//      Redimension/N=(pixelsX*pixelsY) qz_val,qval 
    747752 
    748753        //not demo-compatible, but approx 8x faster!!    
    749754#if(cmpstr(stringbykey("IGORKIND",IgorInfo(0),":",";"),"pro") == 0)      
    750         Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val as fullpath       // /M=termStr specifies terminator 
     755        Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val as fullpath       // /M=termStr specifies terminator       
     756//      Save/G/M="\r\n" labelWave,qx_val,qy_val,qz_val,qval,z_val as fullpath   // for debugging, write out everything 
    751757#else 
    752758        Open refNum as fullpath 
     
    757763#endif 
    758764         
    759         Killwaves/Z spWave,labelWave,qx_val,qy_val,z_val 
     765        Killwaves/Z spWave,labelWave,qx_val,qy_val,z_val,qval,qz_val 
    760766         
    761767        Print "QxQy_Export File written: ", GetFileNameFromPathNoSemi(fullPath) 
Note: See TracChangeset for help on using the changeset viewer.