Ignore:
Timestamp:
Feb 1, 2019 2:25:12 PM (4 years ago)
Author:
srkline
Message:

added procedures to output QxQy_ASCII data. Each panel is output into its own file. Output format is the same as for 2D SANS data, including the 2D resolution function. However, reading the data back in with the SANS analysis macros currently does not redimension the data back to the matrix correctly as it assumes a square detector.

I will need to add the X and Y dimensions of each panel into the header, and then make use of these values when they are read in. Also, writing the QxQy? data is quick for the M and F panels ( 0.8 s) but is extremely slow for the back, High-res panel (120 s) since there are 1.1.E6 points there vs. 6144 pts. I'll need to find a way to speed this operation up.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Utils.ipf

    r1117 r1119  
    14271427        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type) 
    14281428        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type) 
    1429                                  
    1430  
    1431         Variable ret1,ret2,ret3 
    1432         Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses 
    1433  
    1434 // TODO: check the units of all of the inputs 
    1435  
    1436 // lambda = wavelength [A] 
    1437         lambda = V_getWavelength(folderStr) 
    1438          
    1439 // lambdaWidth = [dimensionless] 
    1440         lambdaWidth = V_getWavelength_spread(folderStr) 
    1441          
    1442 // DDet = detector pixel resolution [cm]        **assumes square pixel 
    1443         // V_getDet_pixel_fwhm_x(folderStr,detStr) 
    1444         // V_getDet_pixel_fwhm_y(folderStr,detStr) 
    1445 //      DDet = 0.8              // TODO -- this is hard-wired 
    1446  
    1447         if(strlen(type) == 1) 
    1448                 // it's "B" 
    1449                 DDet = V_getDet_pixel_fwhm_x(folderStr,type)            // value is already in cm 
    1450         else 
    1451                 DDet = V_getDet_pixel_fwhm_x(folderStr,type[0,1])               // value is already in cm 
    1452         endif 
    1453                  
    1454 // apOff = sample aperture to sample distance [cm] 
    1455         apOff = 10              // TODO -- this is hard-wired 
    1456          
    1457 // S1 = source aperture diameter [mm] 
    1458 // may be either circle or rectangle 
    1459         String s1_shape="",bs_shape="" 
    1460         Variable width,height,equiv_S1,equiv_bs 
    1461          
    1462          
    1463         s1_shape = V_getSourceAp_shape(folderStr) 
    1464         if(cmpstr(s1_shape,"CIRCLE") == 0) 
    1465                 S1 = str2num(V_getSourceAp_size(folderStr)) 
    1466         else 
    1467                 S1 = V_getSourceAp_height(folderStr)            // TODO: need the width or at least an equivalent diameter 
    1468         endif 
    1469          
    1470          
    1471 // S2 = sample aperture diameter [cm] 
    1472 // as of 3/2018, the "internal" sample aperture is not in use, only the external 
    1473 // TODO : verify the units on the Ap2 (external) 
    1474 // sample aperture 1(internal) is set to report "12.7 mm" as a STRING 
    1475 // sample aperture 2(external) reports the number typed in... 
    1476 // 
    1477 // so I'm trusting [cm] is in the file 
    1478         S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [mm] 
    1479          
    1480 // L1 = source to sample distance [m]  
    1481         L1 = V_getSourceAp_distance(folderStr)/100 
    1482  
    1483 // L2 = sample to detector distance [m] 
    1484 // take the first two characters of the "type" to get the correct distance. 
    1485 // if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution 
    1486 // is not an issue for the slightly different distances. 
    1487         if(strlen(type) == 1) 
    1488                 // it's "B" 
    1489                 L2 = V_getDet_ActualDistance(folderStr,type)/100                //convert cm to m 
    1490         else 
    1491                 L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m 
    1492         endif 
    1493          
    1494 // BS = beam stop diameter [mm] 
    1495 //TODO:? which BS is in? carr2, carr3, none? 
    1496 // -- need to check the detector, num_beamstops field, then description, then shape/size or shape/height and shape/width 
    1497 // 
    1498 // TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!! 
    1499         BS = V_DeduceBeamstopDiameter(folderStr,type)           //returns diameter in [mm] 
    1500 //      BS = V_getBeamStopC2_size(folderStr)            // Units are [mm]  
    1501 //      BS = 25.4                       //TODO hard-wired value 
    1502          
    1503 //      bs_shape = V_getBeamStopC2_shape(folderStr) 
    1504 //      if(cmpstr(s1_shape,"CIRCLE") == 0) 
    1505 //              bs = V_getBeamStopC2_size(folderStr) 
    1506 //      else 
    1507 //              bs = V_getBeamStopC2_height(folderStr)   
    1508 //      endif 
    1509  
    1510  
    1511          
    1512 // del_r = step size [mm] = binWidth*(mm/pixel)  
    1513         del_r = 1*DDet*10               // TODO: this is probably not the correct value 
    1514  
    1515 // usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam 
    1516         usingLenses = 0 
    1517  
    1518 if(cmpstr(detStr,"FL")==0) 
    1519         Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses" 
    1520         Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses 
    1521 endif 
    1522  
    1523  
    1524 // TODO: 
    1525 // this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc) 
     1429        Variable ret1,ret2,ret3                  
     1430 
     1431// all of the different collimation conditions are handled within the V_getResolution function 
     1432// which is responsible for switching based on the different collimation types (white beam, slit, Xtal, etc) 
    15261433// to calculate the correct resolution, or fill the waves with the correct "flags" 
    15271434// 
    … …  
    15421449// narrowSlit 
    15431450// narrowSlit_whiteBeam 
    1544  
    1545         if(cmpstr(collimationStr,"pinhole") == 0) 
    1546  
    1547                 ii=0 
    1548                 do 
    1549                         V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3) 
    1550                         sigmaq[ii] = ret1        
    1551                         qbar[ii] = ret2  
    1552                         fsubs[ii] = ret3         
    1553                         ii+=1 
    1554                 while(ii<nq) 
    1555          
    1556         endif 
    1557          
    1558  
    1559         if(cmpstr(collimationStr,"pinhole_whiteBeam") == 0) 
    1560  
    1561 //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution. 
    1562 // the white beam distribution will need to be flagged some other way 
    1563 // 
    1564                 lambdaWidth = 0 
    1565                  
    1566                 ii=0 
    1567                 do 
    1568                         V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3) 
    1569                         sigmaq[ii] = ret1        
    1570                         qbar[ii] = ret2  
    1571                         fsubs[ii] = ret3         
    1572                         ii+=1 
    1573                 while(ii<nq) 
    1574          
    1575         endif 
    1576  
    1577         if(cmpstr(collimationStr,"convergingPinholes") == 0) 
    1578  
    1579 //              set usingLenses == 1 so that the Gaussian resolution calculation will be for a focus condition 
    1580 // 
    1581                 usingLenses = 1 
    1582                  
    1583                 ii=0 
    1584                 do 
    1585                         V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3) 
    1586                         sigmaq[ii] = ret1        
    1587                         qbar[ii] = ret2  
    1588                         fsubs[ii] = ret3         
    1589                         ii+=1 
    1590                 while(ii<nq) 
    1591          
    1592         endif 
    1593  
    1594  
    1595 // should not end up here, except for odd testing cases 
    1596         if(cmpstr(collimationStr,"narrowSlit") == 0) 
    1597  
    1598                 Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???" 
    1599                 ii=0 
    1600                 do 
    1601                         V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3) 
    1602                         sigmaq[ii] = ret1        
    1603                         qbar[ii] = ret2  
    1604                         fsubs[ii] = ret3         
    1605                         ii+=1 
    1606                 while(ii<nq) 
    1607          
    1608         endif 
    1609          
    1610 // should not end up here, except for odd testing cases 
    1611         if(cmpstr(collimationStr,"narrowSlit_whiteBeam") == 0) 
    1612  
    1613 //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution. 
    1614 // the white beam distribution will need to be flagged some other way 
    1615 // 
    1616                 Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???" 
    1617  
    1618                 lambdaWidth = 0 
    1619                  
    1620                 ii=0 
    1621                 do 
    1622                         V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3) 
    1623                         sigmaq[ii] = ret1        
    1624                         qbar[ii] = ret2  
    1625                         fsubs[ii] = ret3         
    1626                         ii+=1 
    1627                 while(ii<nq) 
    1628          
    1629         endif 
    1630  
    1631  
    1632  
    1633                  
     1451// 
     1452        ii=0 
     1453        do 
     1454                V_getResolution(qBin_qxqy[ii],folderStr,type,collimationStr,ret1,ret2,ret3) 
     1455                sigmaq[ii] = ret1        
     1456                qbar[ii] = ret2  
     1457                fsubs[ii] = ret3         
     1458                ii+=1 
     1459        while(ii<nq) 
     1460 
     1461 
    16341462        SetDataFolder root: 
    16351463         

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_BeamCenter.ipf

    r1118 r1119  
    4747        PopupMenu popup_0,mode=1,popvalue="FL",value= #"\"FL;FR;FT;FB;ML;MR;MT;MB;B;\"" 
    4848        PopupMenu popup_1,pos={200,20},size={157,20},proc=V_DetModelPopMenuProc,title="Model Function" 
    49         PopupMenu popup_1,mode=1,popvalue="BroadPeak",value= #"\"BroadPeak;other;\"" 
     49        PopupMenu popup_1,mode=1,popvalue="BroadPeak",value= #"\"BroadPeak;BroadPeak_constrained;PowerLaw;\"" 
    5050        PopupMenu popup_2,pos={20,20},size={109,20},title="Data Source"//,proc=SetFldrPopMenuProc 
    5151        PopupMenu popup_2,mode=1,popvalue="RAW",value= #"\"RAW;SAM;VCALC;\"" 
    … …  
    5757        Button button_4,pos={615,400},size={110,20},proc=V_CtrTableButtonProc,title="Ctr table" 
    5858        Button button_5,pos={730,440},size={110,20},proc=V_WriteCtrTableButtonProc,title="Write table",disable=2 
     59 
     60 
     61        Button button_6,pos={615,470},size={110,20},proc=V_MaskBeforeFitButtonProc,title="Mask Panel" 
     62        Button button_7,pos={730,470},size={110,20},proc=V_ConvertFitPix2cmButtonProc,title="Convert Pix2Cm" 
     63        Button button_8,pos={615,500},size={110,20},proc=V_GizmoFitButtonProc,title="Gizmo" 
    5964 
    6065 
    … …  
    352357// 
    353358// TODO - make a better guess (how?) 
     359// TODO - make the guess appropriate for the fitted model 
    354360// 
    355361Function V_DetFitGuessButtonProc(ba) : ButtonControl 
    … …  
    376382End 
    377383 
     384 
     385Function V_GizmoFitButtonProc(ba) : ButtonControl 
     386        STRUCT WMButtonAction &ba 
     387 
     388        switch( ba.eventCode ) 
     389                case 2: // mouse up 
     390                        // click code here 
     391 
     392                        Execute "V_Gizmo_PeakFit()" 
     393                         
     394                        break 
     395                case -1: // control being killed 
     396                        break 
     397        endswitch 
     398 
     399        return 0 
     400End 
     401 
     402 
     403Function V_MaskBeforeFitButtonProc(ba) : ButtonControl 
     404        STRUCT WMButtonAction &ba 
     405 
     406        switch( ba.eventCode ) 
     407                case 2: // mouse up 
     408                        // click code here 
     409 
     410                        Execute "V_NaN_BeforeFit()" 
     411                         
     412                        break 
     413                case -1: // control being killed 
     414                        break 
     415        endswitch 
     416 
     417        return 0 
     418End 
     419 
     420 
     421// TODO - read from the popup for the panel string 
     422// TODO - read the appropriate coeffients for xy, depending on the model selected 
     423// 
     424Function V_ConvertFitPix2cmButtonProc(ba) : ButtonControl 
     425        STRUCT WMButtonAction &ba 
     426 
     427        switch( ba.eventCode ) 
     428                case 2: // mouse up 
     429                        // click code here 
     430 
     431                        String detStr 
     432                        Variable xPix,yPix 
     433//                      Prompt detStr, "enter the panel string" 
     434//                      Prompt xPix, "enter the x pixel center" 
     435//                      Prompt yPix, "enter the y pixel center" 
     436//                      DoPrompt "enter the values",detStr,xPix,yPix 
     437 
     438                        ControlInfo popup_0 
     439                        detStr = S_Value 
     440                        Wave coefW=root:coef_PeakPix2D 
     441         
     442                        xPix = coefW[9] 
     443                        yPix = coefW[10] 
     444                         
     445                        V_Convert_FittedPix_2_cm(detStr,xPix,yPix) 
     446                         
     447                        break 
     448                case -1: // control being killed 
     449                        break 
     450        endswitch 
     451 
     452        return 0 
     453End 
     454 
     455 
    378456// 
    379457// TODO -- currently hard-wired for coefficients from the only fit function 
    … …  
    430508                        Wave coefW=root:coef_PeakPix2D 
    431509                         
    432                         FuncFitMD/H="11000111100"/NTHR=0 V_BroadPeak_Pix2D coefW  dispW /D                       
     510                        FuncFitMD/H="11000101100"/NTHR=0/M=2 V_BroadPeak_Pix2D coefW  dispW /D                   
    433511                         
    434512                        Wave ws=W_sigma 
    … …  
    811889// in V_Marquee_Operation.ipf 
    812890// 
    813 // ** updated these values for the FRONT only with fitted arcs of AgBeh (March 2018 data, run 4494) 
     891// ** updated these values for the FRONT only with fitted arcs of AgBeh (Dec 2018 data, multiple runs) 
    814892// 
    815893Proc V_fDeriveBeamCenters_VelSel(x_FrontReference,y_FrontReference,x_MiddleReference,y_MiddleReference) 
    … …  
    821899        newYCtr_cm[1] = y_FrontReference 
    822900        // FL 
    823 //      newXCtr_cm[0] = x_FrontReference - (0.03 + 0.03)/2              //OLD, pre Nov 2018 
     901//      newXCtr_cm[0] = x_FrontReference - (0.03 + 0.03)/2              //OLD, pre Dec 2018 
    824902//      newYCtr_cm[0] = y_FrontReference + (0.34 + 0.32)/2 
    825         newXCtr_cm[0] = x_FrontReference + 0.26 
    826         newYCtr_cm[0] = y_FrontReference + 0.33 
     903        newXCtr_cm[0] = x_FrontReference + 0.13                         //NEW Dec 2018 
     904        newYCtr_cm[0] = y_FrontReference + 0.35 
    827905        // FB 
    828 //      newXCtr_cm[3] = x_FrontReference - (2.02 + 2.06)/2              // OLD, pre Nov 2018 
     906//      newXCtr_cm[3] = x_FrontReference - (2.02 + 2.06)/2              // OLD, pre Dec 2018 
    829907//      newYCtr_cm[3] = y_FrontReference - (0.12 + 0.19)/2              // (-) is correct here 
    830         newXCtr_cm[3] = x_FrontReference - 1.84 
    831         newYCtr_cm[3] = y_FrontReference + 0.41 
    832         // FT (not a duplicate of FB anymore) 
    833         newXCtr_cm[2] = x_FrontReference - 1.33 
    834         newYCtr_cm[2] = y_FrontReference - 0.30 
     908        newXCtr_cm[3] = x_FrontReference + 0.95                                 // NEW Dec 2018 
     909        newYCtr_cm[3] = y_FrontReference + 0.77 
     910        // FT  
     911//      newXCtr_cm[2] = newXCtr_cm[3]                           // OLD, pre Dec 2018 
     912//      newYCtr_cm[2] = newYCtr_cm[3] 
     913        newXCtr_cm[2] = x_FrontReference + 1.59                         // NEW Dec 2018 (not a duplicate of FB anymore) 
     914        newYCtr_cm[2] = y_FrontReference + 0.09 
    835915         
    836916        // MR 
    … …  
    838918        newYCtr_cm[5] = y_MiddleReference 
    839919        // ML 
    840         newXCtr_cm[4] = x_MiddleReference - (0.06 + 0.05)/2 
    841         newYCtr_cm[4] = y_MiddleReference + (0.14 + 0.01)/2 
     920//      newXCtr_cm[4] = x_MiddleReference - (0.06 + 0.05)/2 
     921//      newYCtr_cm[4] = y_MiddleReference + (0.14 + 0.01)/2 
     922        newXCtr_cm[4] = x_MiddleReference + 0.26 
     923        newYCtr_cm[4] = y_MiddleReference - 0.16 
    842924        // MB 
    843         newXCtr_cm[7] = x_MiddleReference - (0.51 + 0.62)/2 
    844         newYCtr_cm[7] = y_MiddleReference + (0.79 + 0.74)/2 
    845         // MT (duplicate MB) 
    846         newXCtr_cm[6] = newXCtr_cm[7] 
    847         newYCtr_cm[6] = newYCtr_cm[7]    
     925//      newXCtr_cm[7] = x_MiddleReference - (0.51 + 0.62)/2 
     926//      newYCtr_cm[7] = y_MiddleReference + (0.79 + 0.74)/2 
     927        newXCtr_cm[7] = x_MiddleReference - 0.89 
     928        newYCtr_cm[7] = y_MiddleReference + 0.96 
     929        // MT  
     930        newXCtr_cm[6] = x_MiddleReference - 0.28 
     931        newYCtr_cm[6] = y_MiddleReference + 0.60 
    848932         
    849933         
    … …  
    912996End 
    913997 
     998 
     999 
     1000 
     1001 
     1002Window V_Gizmo_PeakFit() : GizmoPlot 
     1003        PauseUpdate; Silent 1           // building window... 
     1004        // Building Gizmo 7 window... 
     1005        NewGizmo/W=(35,45,550,505) 
     1006        ModifyGizmo startRecMacro=700 
     1007        ModifyGizmo scalingOption=63 
     1008        AppendToGizmo Surface=root:Packages:NIST:VSANS:Globals:BeamCenter:curDispPanel,name=surface0 
     1009        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ srcMode,0} 
     1010        ModifyGizmo ModifyObject=surface0,objectType=surface,property={ surfaceCTab,Rainbow} 
     1011        AppendToGizmo Axes=boxAxes,name=axes0 
     1012        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1} 
     1013        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1} 
     1014        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,3} 
     1015        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,3} 
     1016        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,3} 
     1017        ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0} 
     1018        AppendToGizmo Surface=root:PeakPix2D_mat,name=surface1 
     1019        ModifyGizmo ModifyObject=surface1,objectType=surface,property={ srcMode,0} 
     1020        ModifyGizmo ModifyObject=surface1,objectType=surface,property={ surfaceCTab,Grays} 
     1021        ModifyGizmo setDisplayList=0, object=surface0 
     1022        ModifyGizmo setDisplayList=1, object=axes0 
     1023        ModifyGizmo setDisplayList=2, object=surface1 
     1024        ModifyGizmo autoscaling=1 
     1025        ModifyGizmo currentGroupObject="" 
     1026        ModifyGizmo showInfo 
     1027        ModifyGizmo infoWindow={551,23,1368,322} 
     1028        ModifyGizmo endRecMacro 
     1029        ModifyGizmo idleEventQuaternion={1.38005e-05,-1.48789e-05,-6.11841e-06,1} 
     1030EndMacro 
     1031 
     1032 
     1033Proc V_NaN_BeforeFit(x1,x2,y1,y2) 
     1034        Variable x1,x2,y1,y2 
     1035         
     1036        root:Packages:NIST:VSANS:Globals:BeamCenter:curDispPanel[x1,x2][y1,y2] = NaN 
     1037End 
     1038 
     1039 
     1040Function V_Convert_FittedPix_2_cm(panel,xPix,yPix) 
     1041        String panel 
     1042        Variable xPix,yPix 
     1043         
     1044        Make/O/D/N=128 tmpTube 
     1045         
     1046        String pathStr = "root:Packages:NIST:VSANS:RAW:entry:instrument:detector_" 
     1047        Variable x_cm,y_cm 
     1048         
     1049         
     1050        Wave xW = $(pathStr+panel+":data_realDistX") 
     1051        Wave yW = $(pathStr+panel+":data_realDistY") 
     1052 
     1053        strswitch(panel)        // string switch 
     1054                case "FL": 
     1055                case "ML": 
     1056                        // for Left/Right 
     1057                        tmpTube = yW[0][p]                       
     1058                        // for left 
     1059                        x_cm = (xW[47][0] + (xPix-47)*8.4)/10 
     1060                        y_cm = tmpTube[yPix]/10  
     1061         
     1062                        break            
     1063                case "FR":       
     1064                case "MR": 
     1065                        // for Left/Right 
     1066                        tmpTube = yW[0][p]                       
     1067                        // for right 
     1068                        x_cm = (xW[0][0] + xPix*8.4)/10 
     1069                        y_cm = tmpTube[yPix]/10 
     1070                         
     1071                        break 
     1072                case "FT":       
     1073                case "MT": 
     1074                        // for Top/Bottom 
     1075                        tmpTube = xW[p][0] 
     1076                         
     1077                        x_cm = tmpTube[xPix]/10 
     1078                        y_cm = (yW[0][0] + yPix*8.4)/10 
     1079                         
     1080                        break            
     1081                case "FB":       
     1082                case "MB": 
     1083                        // for Top/Bottom 
     1084                        tmpTube = xW[p][0] 
     1085                         
     1086                        x_cm = tmpTube[xPix]/10 
     1087                        y_cm = (yW[0][47] + (yPix-47)*8.4)/10 
     1088                                                 
     1089                        break 
     1090                default:                        // optional default expression executed 
     1091                        Print "No case matched in V_Convert_FittedPix_2_cm" 
     1092                        return(1) 
     1093        endswitch 
     1094 
     1095         
     1096        Print "Converted Center = ",x_cm,y_cm 
     1097        return(0) 
     1098end 
     1099 
     1100Function V_MakeCorrelationMatrix() 
     1101         
     1102        Wave M_Covar=M_Covar 
     1103        Duplicate M_Covar, CorMat        // You can use any name instead of CorMat 
     1104        CorMat = M_Covar[p][q]/sqrt(M_Covar[p][p]*M_Covar[q][q]) 
     1105        Edit/K=0 root:CorMat 
     1106 
     1107        return(0) 
     1108End 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_BroadPeak_Pix_2D.ipf

    r1079 r1119  
    234234        if(ratio > 1) 
    235235        //      qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2)                    // use if the pixels are square 
    236                 qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2/ratio)                      // use for LR panels where the y pixels are half the size of x   
     236                qval = sqrt((xw-xCtr)^2+((yw-yCtr)^2)/ratio)                    // use for LR panels where the y pixels are half the size of x   
    237237        else 
    238                 qval = sqrt((xw-xCtr)^2*ratio+(yw-yCtr)^2)                      // use for TB panels where the y pixels are twice the size of x  
     238                qval = sqrt(((xw-xCtr)^2)*ratio+(yw-yCtr)^2)                    // use for TB panels where the y pixels are twice the size of x  
    239239        endif 
    240240 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Instrument_Resolution.ipf

    r1081 r1119  
    3838 
    3939 
     40 
     41 
     42//********************** 
     43// Resolution calculation - used by the averaging routines 
     44// to calculate the resolution function at each q-value 
     45// - the return value is not used 
     46// 
     47// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR 
     48// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114 
     49// 
     50// - 21 MAR 07 uses projected BS diameter on the detector 
     51// - APR 07 still need to add resolution with lenses. currently there is no flag in the  
     52//          raw data header to indicate the presence of lenses. 
     53// 
     54// - Aug 07 - added input to switch calculation based on lenses (==1 if in) 
     55// 
     56// - SANS -- called by CircSectAvg.ipf and RectAnnulAvg.ipf 
     57// 
     58// - VSANS -- called in VC_fDoBinning_QxQy2D(folderStr, binningType) 
     59// 
     60// DDet is the detector pixel resolution 
     61// apOff is the offset between the sample aperture and the sample position 
     62// 
     63// 
     64// INPUT: 
     65// inQ = q-value [1/A] 
     66// folderStr = folder with the current reduction step 
     67// type = binning type (not the same as the detStr) 
     68// collimationStr = collimation type, to switch for lenses, etc. 
     69 
     70// READ/DERIVED within the function 
     71// lambda = wavelength [A] 
     72// lambdaWidth = [dimensionless] 
     73// DDet = detector pixel resolution [cm]        **assumes square pixel 
     74// apOff = sample aperture to sample distance [cm] 
     75// S1 = source aperture diameter [mm] 
     76// S2 = sample aperture diameter [mm] 
     77// L1 = source to sample distance [m]  
     78// L2 = sample to detector distance [m] 
     79// BS = beam stop diameter [mm] 
     80// del_r = step size [mm] = binWidth*(mm/pixel)  
     81// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam 
     82// 
     83// OUPUT: 
     84// SigmaQ 
     85// QBar 
     86// fSubS 
     87// 
     88// 
     89Function V_getResolution(inQ,folderStr,type,collimationStr,SigmaQ,QBar,fSubS) 
     90        Variable inQ 
     91        String folderStr,type,collimationStr 
     92        Variable &SigmaQ, &QBar, &fSubS         //these are the output quantities at the input Q value 
     93 
     94 
     95        Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses 
     96                 
     97        //lots of calculation variables 
     98        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g 
     99        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r 
     100 
     101        //Constants 
     102        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron 
     103        Variable g = 981.0                              //gravity acceleration [cm/s^2] 
     104 
     105///////// get all of the values from the header 
     106// TODO: check the units of all of the inputs 
     107// lambda = wavelength [A] 
     108        lambda = V_getWavelength(folderStr) 
     109         
     110// lambdaWidth = [dimensionless] 
     111        lambdaWidth = V_getWavelength_spread(folderStr) 
     112         
     113// DDet = detector pixel resolution [cm]        **assumes square pixel 
     114        // V_getDet_pixel_fwhm_x(folderStr,detStr) 
     115        // V_getDet_pixel_fwhm_y(folderStr,detStr) 
     116//      DDet = 0.8              // TODO -- this is hard-wired 
     117 
     118        if(strlen(type) == 1) 
     119                // it's "B" 
     120                DDet = V_getDet_pixel_fwhm_x(folderStr,type)            // value is already in cm 
     121        else 
     122                DDet = V_getDet_pixel_fwhm_x(folderStr,type[0,1])               // value is already in cm 
     123        endif 
     124                 
     125// apOff = sample aperture to sample distance [cm] 
     126        apOff = 10              // TODO -- this is hard-wired 
     127 
     128 
     129         
     130// S1 = source aperture diameter [mm] 
     131// may be either circle or rectangle 
     132        String s1_shape="",bs_shape="" 
     133        Variable width,height,equiv_S1,equiv_bs 
     134         
     135         
     136        s1_shape = V_getSourceAp_shape(folderStr) 
     137        if(cmpstr(s1_shape,"CIRCLE") == 0) 
     138                S1 = str2num(V_getSourceAp_size(folderStr)) 
     139        else 
     140                S1 = V_getSourceAp_height(folderStr)            // TODO: need the width or at least an equivalent diameter 
     141        endif 
     142         
     143         
     144// S2 = sample aperture diameter [cm] 
     145// as of 3/2018, the "internal" sample aperture is not in use, only the external 
     146// TODO : verify the units on the Ap2 (external) 
     147// sample aperture 1(internal) is set to report "12.7 mm" as a STRING 
     148// sample aperture 2(external) reports the number typed in... 
     149// 
     150// so I'm trusting [cm] is in the file 
     151        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [mm] 
     152         
     153// L1 = source to sample distance [m]  
     154        L1 = V_getSourceAp_distance(folderStr)/100 
     155 
     156// L2 = sample to detector distance [m] 
     157// take the first two characters of the "type" to get the correct distance. 
     158// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution 
     159// is not an issue for the slightly different distances. 
     160        if(strlen(type) == 1) 
     161                // it's "B" 
     162                L2 = V_getDet_ActualDistance(folderStr,type)/100                //convert cm to m 
     163        else 
     164                L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m 
     165        endif 
     166         
     167// BS = beam stop diameter [mm] 
     168//TODO:? which BS is in? carr2, carr3, none? 
     169// -- need to check the detector, num_beamstops field, then description, then shape/size or shape/height and shape/width 
     170// 
     171// TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!! 
     172        BS = V_DeduceBeamstopDiameter(folderStr,type)           //returns diameter in [mm] 
     173//      BS = V_getBeamStopC2_size(folderStr)            // Units are [mm]  
     174//      BS = 25.4                       //TODO hard-wired value 
     175         
     176//      bs_shape = V_getBeamStopC2_shape(folderStr) 
     177//      if(cmpstr(s1_shape,"CIRCLE") == 0) 
     178//              bs = V_getBeamStopC2_size(folderStr) 
     179//      else 
     180//              bs = V_getBeamStopC2_height(folderStr)   
     181//      endif 
     182 
     183 
     184         
     185// del_r = step size [mm] = binWidth*(mm/pixel)  
     186        del_r = 1*DDet*10               // TODO: this is probably not the correct value 
     187 
     188// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam 
     189        usingLenses = 0 
     190 
     191//if(cmpstr(type[0,1],"FL")==0) 
     192//      Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses" 
     193//      Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses 
     194//endif 
     195 
     196 
     197 
     198// TODO: 
     199// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc) 
     200// to calculate the correct resolution, or fill the waves with the correct "flags" 
     201// 
     202 
     203// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other  
     204//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution, 
     205//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral, 
     206//  with the inner(lambda) calculated first, then the outer(geometry). 
     207// 
     208 
     209// possible values are: 
     210// 
     211// pinhole 
     212// pinhole_whiteBeam 
     213// convergingPinholes 
     214// 
     215// *slit data should be reduced using the slit routine, not here, proceed but warn 
     216// narrowSlit 
     217// narrowSlit_whiteBeam 
     218 
     219 
     220        if(cmpstr(collimationStr,"pinhole") == 0) 
     221                //nothing to change      
     222        endif 
     223 
     224        if(cmpstr(collimationStr,"pinhole_whiteBeam") == 0) 
     225                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution. 
     226                // the white beam distribution will need to be flagged some other way 
     227                // 
     228                lambdaWidth = 0 
     229        endif 
     230 
     231        if(cmpstr(collimationStr,"convergingPinholes") == 0) 
     232 
     233                //              set usingLenses == 1 so that the Gaussian resolution calculation will be for a focus condition 
     234                usingLenses = 1 
     235        endif 
     236 
     237 
     238// should not end up here, except for odd testing cases 
     239        if(cmpstr(collimationStr,"narrowSlit") == 0) 
     240 
     241                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???" 
     242        endif 
     243         
     244// should not end up here, except for odd testing cases 
     245        if(cmpstr(collimationStr,"narrowSlit_whiteBeam") == 0) 
     246 
     247                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution. 
     248                // the white beam distribution will need to be flagged some other way 
     249                // 
     250                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???" 
     251 
     252                lambdaWidth = 0 
     253        endif 
     254 
     255 
     256///////////////////////////// 
     257///////////////////////////// 
     258// do the calculation 
     259        S1 *= 0.5*0.1                   //convert to radius and [cm] 
     260        S2 *= 0.5*0.1 
     261 
     262        L1 *= 100.0                     // [cm] 
     263        L1 -= apOff                             //correct the distance 
     264 
     265        L2 *= 100.0 
     266        L2 += apOff 
     267        del_r *= 0.1                            //width of annulus, convert mm to [cm] 
     268         
     269        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm] 
     270        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture 
     271        Variable LB 
     272        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical) 
     273        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax 
     274         
     275        //Start resolution calculation 
     276        a2 = S1*L2/L1 + S2*(L1+L2)/L1 
     277        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2) 
     278        lp = 1.0/( 1.0/L1 + 1.0/L2) 
     279 
     280        v_lambda = lambdaWidth^2/6.0 
     281         
     282//      if(usingLenses==1)                      //SRK 2007 
     283        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header 
     284                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term 
     285        else 
     286                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form 
     287        endif 
     288         
     289        v_d = (DDet/2.3548)^2 + del_r^2/12.0                    //the 2.3548 is a conversion from FWHM->Gauss, see http://mathworld.wolfram.com/GaussianFunction.html 
     290        vz = vz_1 / lambda 
     291        yg = 0.5*g*L2*(L1+L2)/vz^2 
     292        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007 
     293 
     294        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) )) 
     295        delta = 0.5*(BS - r0)^2/v_d 
     296 
     297        if (r0 < BS)  
     298                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta)) 
     299        else 
     300                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta)) 
     301        endif 
     302 
     303        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) ) 
     304        if (fSubS <= 0.0)  
     305                fSubS = 1.e-10 
     306        endif 
     307        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi)) 
     308        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d 
     309 
     310        rmd = fr*r0 
     311        v_r1 = v_b + fv*v_d +v_g 
     312 
     313        rm = rmd + 0.5*v_r1/rmd 
     314        v_r = v_r1 - 0.5*(v_r1/rmd)^2 
     315        if (v_r < 0.0)  
     316                v_r = 0.0 
     317        endif 
     318        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2)) 
     319        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda) 
     320 
     321 
     322// more readable method for calculating the variance in Q 
     323// EXCEPT - this is calculated for Qo, NOT qBar 
     324// (otherwise, they are nearly equivalent, except for close to the beam stop) 
     325//      Variable kap,a_val,a_val_l2,m_h 
     326//      g = 981.0                               //gravity acceleration [cm/s^2] 
     327//      m_h     = 252.8                 // m/h [=] s/cm^2 
     328//       
     329//      kap = 2*pi/lambda 
     330//      a_val = L2*(L1+L2)*g/2*(m_h)^2 
     331//      a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2 
     332// 
     333//      sigmaQ = 0 
     334//      sigmaQ = 3*(S1/L1)^2 
     335//       
     336//      if(usingLenses != 0) 
     337//              sigmaQ += 2*(S2/lp)^2*(lambdaWidth)^2   //2nd term w/ lenses 
     338//      else     
     339//              sigmaQ += 2*(S2/lp)^2                                           //2nd term w/ no lenses 
     340//      endif 
     341//       
     342//      sigmaQ += (del_r/L2)^2 
     343//      sigmaQ += 2*(r0/L2)^2*(lambdaWidth)^2 
     344//      sigmaQ += 4*(a_val_l2)^2*lambda^4*(lambdaWidth)^2 
     345//       
     346//      sigmaQ *= kap^2/12 
     347//      sigmaQ = sqrt(sigmaQ) 
     348 
     349 
     350        Return (0) 
     351End 
     352 
     353 
     354// 
     355//********************** 
     356// 2D resolution function calculation - ***NOT*** in terms of X and Y 
     357// but written in terms of Parallel and perpendicular to the Q vector at each point 
     358// 
     359// -- it is more naturally written this way since the 2D function is an ellipse with its major 
     360// axis pointing in the direction of Q_parallel. Hence there is no way to properly define the  
     361// elliptical gaussian in terms of sigmaX and sigmaY 
     362// 
     363// For a full description of the gravity effect on the resolution, see: 
     364// 
     365//      "The effect of gravity on the resolution of small-angle neutron diffraction peaks" 
     366//      D.F.R Mildner, J.G. Barker & S.R. Kline J. Appl. Cryst. (2011). 44, 1127-1129. 
     367//      [ doi:10.1107/S0021889811033322 ] 
     368// 
     369//              2/17/12 SRK 
     370//              NOTE: the first 2/3 of this code is the 1D code, copied here just to have the beam stop 
     371//                              calculation here, if I decide to implement it. The real calculation is all at the  
     372//                              bottom and is quite compact 
     373// 
     374// 
     375// 
     376// 
     377// - 21 MAR 07 uses projected BS diameter on the detector 
     378// - APR 07 still need to add resolution with lenses. currently there is no flag in the  
     379//          raw data header to indicate the presence of lenses. 
     380// 
     381// - Aug 07 - added input to switch calculation based on lenses (==1 if in) 
     382// 
     383// passed values are read from RealsRead 
     384// except DDet and apOff, which are set from globals before passing 
     385// 
     386// phi is the azimuthal angle, CCW from +x axis 
     387// r_dist is the real-space distance from ctr of detector to QxQy pixel location 
     388// 
     389// MAR 2011 - removed the del_r terms, they don't apply since no binning is done to the 2D data 
     390// 
     391Function V_get2DResolution(inQ,phi,r_dist,folderStr,type,collimationStr,SigmaQX,SigmaQY,fSubS) 
     392        Variable inQ,phi,r_dist 
     393        String folderStr,type,collimationStr 
     394        Variable &SigmaQX,&SigmaQY,&fSubS               //these are the output quantities at the input Q value 
     395//      Variable SigmaQX,SigmaQY,fSubS          //these are the output quantities at the input Q value 
     396 
     397         
     398        Variable lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses 
     399 
     400        //      phi = FindPhi( pixSize*((p+1)-xctr) , pixSize*((q+1)-yctr)+(2)*yg_d)            //(dx,dy+yg_d) 
     401        //      r_dist = sqrt(  (pixSize*((p+1)-xctr))^2 +  (pixSize*((q+1)-yctr)+(2)*yg_d)^2 )         //radial distance from ctr to pt 
     402 
     403///////// get all of the values from the header 
     404// TODO: check the units of all of the inputs 
     405// lambda = wavelength [A] 
     406        lambda = V_getWavelength(folderStr) 
     407         
     408// lambdaWidth = [dimensionless] 
     409        lambdaWidth = V_getWavelength_spread(folderStr) 
     410         
     411// DDet = detector pixel resolution [cm]        **assumes square pixel 
     412        // V_getDet_pixel_fwhm_x(folderStr,detStr) 
     413        // V_getDet_pixel_fwhm_y(folderStr,detStr) 
     414//      DDet = 0.8              // TODO -- this is hard-wired 
     415 
     416        if(strlen(type) == 1) 
     417                // it's "B" 
     418                DDet = V_getDet_pixel_fwhm_x(folderStr,type)            // value is already in cm 
     419        else 
     420                DDet = V_getDet_pixel_fwhm_x(folderStr,type[0,1])               // value is already in cm 
     421        endif 
     422                 
     423// apOff = sample aperture to sample distance [cm] 
     424        apOff = 10              // TODO -- this is hard-wired 
     425 
     426 
     427// S1 = source aperture diameter [mm] 
     428// may be either circle or rectangle 
     429        String s1_shape="",bs_shape="" 
     430        Variable width,height,equiv_S1,equiv_bs 
     431         
     432        s1_shape = V_getSourceAp_shape(folderStr) 
     433        if(cmpstr(s1_shape,"CIRCLE") == 0) 
     434                S1 = str2num(V_getSourceAp_size(folderStr)) 
     435        else 
     436                S1 = V_getSourceAp_height(folderStr)            // TODO: need the width or at least an equivalent diameter 
     437        endif 
     438         
     439         
     440// S2 = sample aperture diameter [cm] 
     441// as of 3/2018, the "internal" sample aperture is not in use, only the external 
     442// TODO : verify the units on the Ap2 (external) 
     443// sample aperture 1(internal) is set to report "12.7 mm" as a STRING 
     444// sample aperture 2(external) reports the number typed in... 
     445// 
     446// so I'm trusting [cm] is in the file 
     447        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [mm] 
     448         
     449// L1 = source to sample distance [m]  
     450        L1 = V_getSourceAp_distance(folderStr)/100 
     451 
     452// L2 = sample to detector distance [m] 
     453// take the first two characters of the "type" to get the correct distance. 
     454// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution 
     455// is not an issue for the slightly different distances. 
     456        if(strlen(type) == 1) 
     457                // it's "B" 
     458                L2 = V_getDet_ActualDistance(folderStr,type)/100                //convert cm to m 
     459        else 
     460                L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m 
     461        endif 
     462         
     463// BS = beam stop diameter [mm] 
     464//TODO:? which BS is in? carr2, carr3, none? 
     465// -- need to check the detector, num_beamstops field, then description, then shape/size or shape/height and shape/width 
     466// 
     467// TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!! 
     468        BS = V_DeduceBeamstopDiameter(folderStr,type)           //returns diameter in [mm] 
     469//      BS = V_getBeamStopC2_size(folderStr)            // Units are [mm]  
     470//      BS = 25.4                       //TODO hard-wired value 
     471         
     472//      bs_shape = V_getBeamStopC2_shape(folderStr) 
     473//      if(cmpstr(s1_shape,"CIRCLE") == 0) 
     474//              bs = V_getBeamStopC2_size(folderStr) 
     475//      else 
     476//              bs = V_getBeamStopC2_height(folderStr)   
     477//      endif 
     478 
     479 
     480         
     481// del_r = step size [mm] = binWidth*(mm/pixel)  
     482        del_r = 1*DDet*10               // TODO: this is probably not the correct value 
     483 
     484// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam 
     485        usingLenses = 0 
     486 
     487//if(cmpstr(type[0,1],"FL")==0) 
     488//      Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses" 
     489//      Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses 
     490//endif 
     491 
     492 
     493 
     494// TODO: 
     495// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc) 
     496// to calculate the correct resolution, or fill the waves with the correct "flags" 
     497// 
     498 
     499// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other  
     500//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution, 
     501//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral, 
     502//  with the inner(lambda) calculated first, then the outer(geometry). 
     503// 
     504 
     505// possible values are: 
     506// 
     507// pinhole 
     508// pinhole_whiteBeam 
     509// convergingPinholes 
     510// 
     511// *slit data should be reduced using the slit routine, not here, proceed but warn 
     512// narrowSlit 
     513// narrowSlit_whiteBeam 
     514 
     515 
     516        if(cmpstr(collimationStr,"pinhole") == 0) 
     517                //nothing to change      
     518        endif 
     519 
     520        if(cmpstr(collimationStr,"pinhole_whiteBeam") == 0) 
     521                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution. 
     522                // the white beam distribution will need to be flagged some other way 
     523                // 
     524                lambdaWidth = 0 
     525        endif 
     526 
     527        if(cmpstr(collimationStr,"convergingPinholes") == 0) 
     528 
     529                //              set usingLenses == 1 so that the Gaussian resolution calculation will be for a focus condition 
     530                usingLenses = 1 
     531        endif 
     532 
     533 
     534// should not end up here, except for odd testing cases 
     535        if(cmpstr(collimationStr,"narrowSlit") == 0) 
     536 
     537                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???" 
     538        endif 
     539         
     540// should not end up here, except for odd testing cases 
     541        if(cmpstr(collimationStr,"narrowSlit_whiteBeam") == 0) 
     542 
     543                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution. 
     544                // the white beam distribution will need to be flagged some other way 
     545                // 
     546                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???" 
     547 
     548                lambdaWidth = 0 
     549        endif 
     550 
     551 
     552         
     553        //lots of calculation variables 
     554        Variable a2, lp, v_lambda, v_b, v_d, vz, yg, v_g 
     555        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r 
     556 
     557        //Constants 
     558        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron 
     559        Variable g = 981.0                              //gravity acceleration [cm/s^2] 
     560        Variable m_h    = 252.8                 // m/h [=] s/cm^2 
     561 
     562 
     563        S1 *= 0.5*0.1                   //convert to radius and [cm] 
     564        S2 *= 0.5*0.1 
     565 
     566        L1 *= 100.0                     // [cm] 
     567        L1 -= apOff                             //correct the distance 
     568 
     569        L2 *= 100.0 
     570        L2 += apOff 
     571        del_r *= 0.1                            //width of annulus, convert mm to [cm] 
     572         
     573        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm] 
     574        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture 
     575        Variable LB 
     576        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical) 
     577        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax 
     578         
     579        //Start resolution calculation 
     580        a2 = S1*L2/L1 + S2*(L1+L2)/L1 
     581        lp = 1.0/( 1.0/L1 + 1.0/L2) 
     582 
     583        v_lambda = lambdaWidth^2/6.0 
     584         
     585//      if(usingLenses==1)                      //SRK 2007 
     586        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header 
     587                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term 
     588        else 
     589                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form 
     590        endif 
     591         
     592        v_d = (DDet/2.3548)^2 + del_r^2/12.0 
     593        vz = vz_1 / lambda 
     594        yg = 0.5*g*L2*(L1+L2)/vz^2 
     595        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007 
     596 
     597        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) )) 
     598        delta = 0.5*(BS - r0)^2/v_d 
     599 
     600        if (r0 < BS)  
     601                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta)) 
     602        else 
     603                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta)) 
     604        endif 
     605 
     606        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) ) 
     607        if (fSubS <= 0.0)  
     608                fSubS = 1.e-10 
     609        endif 
     610//      fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi)) 
     611//      fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d 
     612// 
     613//      rmd = fr*r0 
     614//      v_r1 = v_b + fv*v_d +v_g 
     615// 
     616//      rm = rmd + 0.5*v_r1/rmd 
     617//      v_r = v_r1 - 0.5*(v_r1/rmd)^2 
     618//      if (v_r < 0.0)  
     619//              v_r = 0.0 
     620//      endif 
     621 
     622        Variable kap,a_val,a_val_L2,proj_DDet 
     623         
     624        kap = 2*pi/lambda 
     625        a_val = L2*(L1+L2)*g/2*(m_h)^2 
     626        a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2 
     627 
     628 
     629        // the detector pixel is square, so correct for phi 
     630        proj_DDet = DDet*cos(phi) + DDet*sin(phi) 
     631 
     632 
     633///////// OLD - don't use --- 
     634//in terms of Q_parallel ("x") and Q_perp ("y") - this works, since parallel is in the direction of Q and I 
     635// can calculate that from the QxQy (I just need the projection) 
     636//// for test case with no gravity, set a_val = 0 
     637//// note that gravity has no wavelength dependence. the lambda^4 cancels out. 
     638//// 
     639////    a_val = 0 
     640////    a_val_l2 = 0 
     641// 
     642//       
     643//      // this is really sigma_Q_parallel 
     644//      SigmaQX = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (sin(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2) 
     645//      SigmaQX += inQ*inQ*v_lambda 
     646// 
     647//      //this is really sigma_Q_perpendicular 
     648//      proj_DDet = DDet*sin(phi) + DDet*cos(phi)               //not necessary, since DDet is the same in both X and Y directions 
     649// 
     650//      SigmaQY = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (cos(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2) 
     651//       
     652//      SigmaQX = sqrt(SigmaQX) 
     653//      SigmaQy = sqrt(SigmaQY) 
     654//       
     655 
     656///////////////////////////////////////////////// 
     657/////    
     658//      ////// this is all new, inclusion of gravity effect into the parallel component 
     659//       perpendicular component is purely geometric, no gravity component 
     660// 
     661// the shadow factor is calculated as above -so keep the above calculations, even though 
     662// most of them are redundant. 
     663// 
     664         
     665////    // 
     666        Variable yg_d,acc,sdd,ssd,lambda0,DL_L,sig_l 
     667        Variable var_qlx,var_qly,var_ql,qx,qy,sig_perp,sig_para, sig_para_new 
     668         
     669        G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2 
     670        acc = vz_1              //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC 
     671        SDD = L2                //1317 
     672        SSD = L1                //1627          //cm 
     673        lambda0 = lambda                //              15 
     674        DL_L = lambdaWidth              //0.236 
     675        SIG_L = DL_L/sqrt(6) 
     676        YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2 
     677/////   Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG 
     678//              Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd 
     679         
     680        sig_perp = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2) 
     681        sig_perp = sqrt(sig_perp) 
     682         
     683// TODO -- not needed???         
     684//      FindQxQy(inQ,phi,qx,qy) 
     685 
     686 
     687// missing a factor of 2 here, and the form is different than the paper, so re-write     
     688//      VAR_QLY = SIG_L^2 * (QY+4*PI*YG_d/(2*SDD*LAMBDA0))^2 
     689//      VAR_QLX = (SIG_L*QX)^2 
     690//      VAR_QL = VAR_QLY + VAR_QLX  //! WAVELENGTH CONTRIBUTION TO VARIANCE 
     691//      sig_para = (sig_perp^2 + VAR_QL)^0.5 
     692         
     693        // r_dist is passed in, [=]cm 
     694        // from the paper 
     695        a_val = 0.5*G*SDD*(SSD+SDD)*m_h^2 * 1e-16               //units now are cm /(A^2) 
     696         
     697        var_QL = 1/6*(kap/SDD)^2*(DL_L)^2*(r_dist^2 - 4*r_dist*a_val*lambda0^2*sin(phi) + 4*a_val^2*lambda0^4) 
     698        sig_para_new = (sig_perp^2 + VAR_QL)^0.5 
     699         
     700         
     701///// return values PBR  
     702        SigmaQX = sig_para_new 
     703        SigmaQy = sig_perp 
     704         
     705////     
     706 
     707        Return (0) 
     708End 
    40709 
    41710 
    … …  
    82751// 
    83752// 
    84 Function V_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,SigmaQ,QBar,fSubS) 
     753Function V_getResolution_old(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,SigmaQ,QBar,fSubS) 
    85754        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses 
    86755        Variable &SigmaQ, &QBar, &fSubS         //these are the output quantities at the input Q value 
    … …  
    190859 
    191860 
    192 // 
    193 //********************** 
    194 // 2D resolution function calculation - ***NOT*** in terms of X and Y 
    195 // but written in terms of Parallel and perpendicular to the Q vector at each point 
    196 // 
    197 // -- it is more naturally written this way since the 2D function is an ellipse with its major 
    198 // axis pointing in the direction of Q_parallel. Hence there is no way to properly define the  
    199 // elliptical gaussian in terms of sigmaX and sigmaY 
    200 // 
    201 // For a full description of the gravity effect on the resolution, see: 
    202 // 
    203 //      "The effect of gravity on the resolution of small-angle neutron diffraction peaks" 
    204 //      D.F.R Mildner, J.G. Barker & S.R. Kline J. Appl. Cryst. (2011). 44, 1127-1129. 
    205 //      [ doi:10.1107/S0021889811033322 ] 
    206 // 
    207 //              2/17/12 SRK 
    208 //              NOTE: the first 2/3 of this code is the 1D code, copied here just to have the beam stop 
    209 //                              calculation here, if I decide to implement it. The real calculation is all at the  
    210 //                              bottom and is quite compact 
    211 // 
    212 // 
    213 // 
    214 // 
    215 // - 21 MAR 07 uses projected BS diameter on the detector 
    216 // - APR 07 still need to add resolution with lenses. currently there is no flag in the  
    217 //          raw data header to indicate the presence of lenses. 
    218 // 
    219 // - Aug 07 - added input to switch calculation based on lenses (==1 if in) 
    220 // 
    221 // passed values are read from RealsRead 
    222 // except DDet and apOff, which are set from globals before passing 
    223 // 
    224 // phi is the azimuthal angle, CCW from +x axis 
    225 // r_dist is the real-space distance from ctr of detector to QxQy pixel location 
    226 // 
    227 // MAR 2011 - removed the del_r terms, they don't apply since no bining is done to the 2D data 
    228 // 
    229 Function V_get2DResolution(inQ,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,SigmaQX,SigmaQY,fSubS) 
    230         Variable inQ, phi,lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses,r_dist 
    231         Variable &SigmaQX,&SigmaQY,&fSubS               //these are the output quantities at the input Q value 
    232          
    233         //lots of calculation variables 
    234         Variable a2, lp, v_lambda, v_b, v_d, vz, yg, v_g 
    235         Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r 
    236  
    237         //Constants 
    238         Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron 
    239         Variable g = 981.0                              //gravity acceleration [cm/s^2] 
    240         Variable m_h    = 252.8                 // m/h [=] s/cm^2 
    241  
    242  
    243         S1 *= 0.5*0.1                   //convert to radius and [cm] 
    244         S2 *= 0.5*0.1 
    245  
    246         L1 *= 100.0                     // [cm] 
    247         L1 -= apOff                             //correct the distance 
    248  
    249         L2 *= 100.0 
    250         L2 += apOff 
    251         del_r *= 0.1                            //width of annulus, convert mm to [cm] 
    252          
    253         BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm] 
    254         // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture 
    255         Variable LB 
    256         LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical) 
    257         BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax 
    258          
    259         //Start resolution calculation 
    260         a2 = S1*L2/L1 + S2*(L1+L2)/L1 
    261         lp = 1.0/( 1.0/L1 + 1.0/L2) 
    262  
    263         v_lambda = lambdaWidth^2/6.0 
    264          
    265 //      if(usingLenses==1)                      //SRK 2007 
    266         if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header 
    267                 v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term 
    268         else 
    269                 v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form 
    270         endif 
    271          
    272         v_d = (DDet/2.3548)^2 + del_r^2/12.0 
    273         vz = vz_1 / lambda 
    274         yg = 0.5*g*L2*(L1+L2)/vz^2 
    275         v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007 
    276  
    277         r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) )) 
    278         delta = 0.5*(BS - r0)^2/v_d 
    279  
    280         if (r0 < BS)  
    281                 inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta)) 
    282         else 
    283                 inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta)) 
    284         endif 
    285  
    286         fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) ) 
    287         if (fSubS <= 0.0)  
    288                 fSubS = 1.e-10 
    289         endif 
    290 //      fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi)) 
    291 //      fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d 
    292 // 
    293 //      rmd = fr*r0 
    294 //      v_r1 = v_b + fv*v_d +v_g 
    295 // 
    296 //      rm = rmd + 0.5*v_r1/rmd 
    297 //      v_r = v_r1 - 0.5*(v_r1/rmd)^2 
    298 //      if (v_r < 0.0)  
    299 //              v_r = 0.0 
    300 //      endif 
    301  
    302         Variable kap,a_val,a_val_L2,proj_DDet 
    303          
    304         kap = 2*pi/lambda 
    305         a_val = L2*(L1+L2)*g/2*(m_h)^2 
    306         a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2 
    307  
    308  
    309         // the detector pixel is square, so correct for phi 
    310         proj_DDet = DDet*cos(phi) + DDet*sin(phi) 
    311  
    312  
    313 ///////// OLD - don't use --- 
    314 //in terms of Q_parallel ("x") and Q_perp ("y") - this works, since parallel is in the direction of Q and I 
    315 // can calculate that from the QxQy (I just need the projection) 
    316 //// for test case with no gravity, set a_val = 0 
    317 //// note that gravity has no wavelength dependence. the lambda^4 cancels out. 
    318 //// 
    319 ////    a_val = 0 
    320 ////    a_val_l2 = 0 
    321 // 
    322 //       
    323 //      // this is really sigma_Q_parallel 
    324 //      SigmaQX = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (sin(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2) 
    325 //      SigmaQX += inQ*inQ*v_lambda 
    326 // 
    327 //      //this is really sigma_Q_perpendicular 
    328 //      proj_DDet = DDet*sin(phi) + DDet*cos(phi)               //not necessary, since DDet is the same in both X and Y directions 
    329 // 
    330 //      SigmaQY = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (cos(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2) 
    331 //       
    332 //      SigmaQX = sqrt(SigmaQX) 
    333 //      SigmaQy = sqrt(SigmaQY) 
    334 //       
    335  
    336 ///////////////////////////////////////////////// 
    337 /////    
    338 //      ////// this is all new, inclusion of gravity effect into the parallel component 
    339 //       perpendicular component is purely geometric, no gravity component 
    340 // 
    341 // the shadow factor is calculated as above -so keep the above calculations, even though 
    342 // most of them are redundant. 
    343 // 
    344          
    345 ////    // 
    346         Variable yg_d,acc,sdd,ssd,lambda0,DL_L,sig_l 
    347         Variable var_qlx,var_qly,var_ql,qx,qy,sig_perp,sig_para, sig_para_new 
    348          
    349         G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2 
    350         acc = vz_1              //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC 
    351         SDD = L2                //1317 
    352         SSD = L1                //1627          //cm 
    353         lambda0 = lambda                //              15 
    354         DL_L = lambdaWidth              //0.236 
    355         SIG_L = DL_L/sqrt(6) 
    356         YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2 
    357 /////   Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG 
    358 //              Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd 
    359          
    360         sig_perp = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2) 
    361         sig_perp = sqrt(sig_perp) 
    362          
    363 // TODO -- not needed???         
    364 //      FindQxQy(inQ,phi,qx,qy) 
    365  
    366  
    367 // missing a factor of 2 here, and the form is different than the paper, so re-write     
    368 //      VAR_QLY = SIG_L^2 * (QY+4*PI*YG_d/(2*SDD*LAMBDA0))^2 
    369 //      VAR_QLX = (SIG_L*QX)^2 
    370 //      VAR_QL = VAR_QLY + VAR_QLX  //! WAVELENGTH CONTRIBUTION TO VARIANCE 
    371 //      sig_para = (sig_perp^2 + VAR_QL)^0.5 
    372          
    373         // r_dist is passed in, [=]cm 
    374         // from the paper 
    375         a_val = 0.5*G*SDD*(SSD+SDD)*m_h^2 * 1e-16               //units now are cm /(A^2) 
    376          
    377         var_QL = 1/6*(kap/SDD)^2*(DL_L)^2*(r_dist^2 - 4*r_dist*a_val*lambda0^2*sin(phi) + 4*a_val^2*lambda0^4) 
    378         sig_para_new = (sig_perp^2 + VAR_QL)^0.5 
    379          
    380          
    381 ///// return values PBR  
    382         SigmaQX = sig_para_new 
    383         SigmaQy = sig_perp 
    384          
    385 ////     
    386  
    387         Return (0) 
    388 End 
    389  
    390  
    391  

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Marquee_Operations.ipf

    r1118 r1119  
    367367// these corrections are exactly the opposite of what is done in V_fDeriveBeamCenters(xFR,yFR,xMR,yMR) 
    368368                if(cmpstr(detStr,"FL") == 0) 
    369                         Print "FRONT Reference X-center (cm) (Velocity Selector) = ",x_mm/10 - 0.26 
    370                         Print "FRONT Reference Y-center (cm) (Velocity Selector) = ",y_mm/10 - 0.33 
    371                          
    372                         Print "FRONT Reference X-center (cm) (Graphite) = ",x_mm/10 + 0.03 
    373                         Print "FRONT Reference Y-center (cm) (Graphite) = ",y_mm/10 - 0.28 
     369//                      Print "FRONT Reference X-center (cm) (Velocity Selector) = ",x_mm/10 + (0.03 + 0.03)/2   // OLD pre Dec 2018 
     370//                      Print "FRONT Reference Y-center (cm) (Velocity Selector) = ",y_mm/10 - (0.34 + 0.32)/2 
     371                        Print "FRONT Reference X-center (cm) (Velocity Selector) = ",x_mm/10 - 0.13     // NEW Dec 2018 values 
     372                        Print "FRONT Reference Y-center (cm) (Velocity Selector) = ",y_mm/10 - 0.35 
     373                                                 
     374        //              Print "FRONT Reference X-center (cm) (Graphite) = ",x_mm/10 + 0.03 
     375        //              Print "FRONT Reference Y-center (cm) (Graphite) = ",y_mm/10 - 0.28 
    374376                endif 
    375377                 
    376378                if(cmpstr(detStr,"ML") == 0) 
    377                         Print "MIDDLE Reference X-center (cm) (Velocity Selector) = ",x_mm/10 + (0.06 + 0.05)/2 
    378                         Print "MIDDLE Reference Y-center (cm) (Velocity Selector) = ",y_mm/10 - (0.14 + 0.01)/2 
    379                          
    380                         Print "MIDDLE Reference X-center (cm) (Graphite) = ",x_mm/10 + (0.06 + 0.05)/2 
    381                         Print "MIDDLE Reference Y-center (cm) (Graphite) = ",y_mm/10 - (0.14 + 0.01)/2 
     379//                      Print "MIDDLE Reference X-center (cm) (Velocity Selector) = ",x_mm/10 + (0.06 + 0.05)/2 
     380//                      Print "MIDDLE Reference Y-center (cm) (Velocity Selector) = ",y_mm/10 - (0.14 + 0.01)/2 
     381                        Print "MIDDLE Reference X-center (cm) (Velocity Selector) = ",x_mm/10 - 0.26 
     382                        Print "MIDDLE Reference Y-center (cm) (Velocity Selector) = ",y_mm/10 + 0.16 
     383                                         
     384        //              Print "MIDDLE Reference X-center (cm) (Graphite) = ",x_mm/10 + (0.06 + 0.05)/2 
     385        //              Print "MIDDLE Reference Y-center (cm) (Graphite) = ",y_mm/10 - (0.14 + 0.01)/2 
    382386                endif 
    383387        endif 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Menu.ipf

    r1118 r1119  
    2929                "-" 
    3030                "Derive Beam Centers - VelSel",V_DeriveBeamCenters_VelSel() 
    31                 "Derive Beam Centers - Graphite",V_DeriveBeamCenters_Graphite() 
     31//              "Derive Beam Centers - Graphite",V_DeriveBeamCenters_Graphite() 
    3232                "-" 
    3333                "Back Detector Saturation",Vm_NumberSaturated() 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Protocol_Reduction.ipf

    r1117 r1119  
    13401340 
    13411341//      Prompt av_typ, "Type of Average",popup,"Circular;Sector;Rectangular;Annular;2D_ASCII;QxQy_ASCII;PNG_Graphic;Sector_PlusMinus;" 
    1342         Prompt av_typ, "Type of Average",popup,"Circular;Narrow_Slit;Annular;" 
     1342        Prompt av_typ, "Type of Average",popup,"Circular;Narrow_Slit;Annular;QxQy_ASCII;" 
    13431343 
    13441344// comment out above line in DEMO_MODIFIED version, and uncomment the line below (to disable PNG save) 
    … …  
    24042404        endif 
    24052405         
     2406        gAbsStr = RemoveFromList("ask",gAbsStr)         //now that values are added, remove "ask" 
    24062407 
    24072408        SetDataFolder root: 
    … …  
    33593360                                break 
    33603361                        case "QxQy_ASCII": 
    3361 //                              QxQy_Export(activeType,fullPath,dialog) 
     3362                                fullPath = S_Path + newFileName+".DAT" 
     3363                                V_QxQy_Export(activeType,fullPath,newFileName,dialog) 
    33623364                                break 
    33633365                        case "PNG_Graphic": 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_WorkFolderUtils.ipf

    r1117 r1119  
    561561                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 
    562562                        detStr = StringFromList(ii, ksDetectorListAll, ";") 
    563                         Wave w = V_getDetectorDataW(fname,detStr) 
    564                         Wave w_err = V_getDetectorDataErrW(fname,detStr) 
    565                          
    566                         V_DIVCorrection(w,w_err,detStr,newType) 
     563                         
     564                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1) 
     565                                // do nothing 
     566                        else 
     567                                Wave w = V_getDetectorDataW(fname,detStr) 
     568                                Wave w_err = V_getDetectorDataErrW(fname,detStr) 
     569                                 
     570                                V_DIVCorrection(w,w_err,detStr,newType) 
     571                        endif 
    567572                endfor 
    568573        else 
    … …  
    632637                         
    633638                endfor 
    634                  
    635                 //"B" is separate 
    636                 detStr = "B" 
    637                 Wave w = V_getDetectorDataW(fname,detStr) 
    638                 Wave cal_x = V_getDet_cal_x(fname,detStr) 
    639                 Wave cal_y = V_getDet_cal_y(fname,detStr) 
    640                  
    641                 V_NonLinearCorrection_B(fname,w,cal_x,cal_y,detStr,destPath) 
    642  
    643 // "B" is always naturally defined in terms of a pixel center. This can be converted to mm,  
    644 // but the experiment will measure pixel x,y - just like ordela detectors. 
    645                  
    646 //              if(kBCTR_CM) 
    647 // 
    648 //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm") 
    649 //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm") 
    650 //                      WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm") 
    651 //                      WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm") 
    652 //                      x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm 
    653 //                      y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm 
    654 //                       
    655 //                      // now I need to convert the beam center in mm to pixels 
    656 //                      // and have some rational place to look for it... 
    657 //                      V_ConvertBeamCtr_to_pixB(fname,detStr,destPath) 
    658 //              else 
    659                         // beam center is in pixels, so use the old routine 
    660                         V_ConvertBeamCtrPix_to_mmB(fname,detStr,destPath) 
    661  
    662 //              endif 
    663                 V_Detector_CalcQVals(fname,detStr,destPath) 
     639 
     640                if(gIgnoreDetB==1)               
     641                        //"B" is separate 
     642                        detStr = "B" 
     643                        Wave w = V_getDetectorDataW(fname,detStr) 
     644                        Wave cal_x = V_getDet_cal_x(fname,detStr) 
     645                        Wave cal_y = V_getDet_cal_y(fname,detStr) 
     646                         
     647                        V_NonLinearCorrection_B(fname,w,cal_x,cal_y,detStr,destPath) 
     648         
     649        // "B" is always naturally defined in terms of a pixel center. This can be converted to mm,  
     650        // but the experiment will measure pixel x,y - just like ordela detectors. 
     651                         
     652        //              if(kBCTR_CM) 
     653        // 
     654        //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm") 
     655        //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm") 
     656        //                      WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm") 
     657        //                      WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm") 
     658        //                      x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm 
     659        //                      y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm 
     660        //                       
     661        //                      // now I need to convert the beam center in mm to pixels 
     662        //                      // and have some rational place to look for it... 
     663        //                      V_ConvertBeamCtr_to_pixB(fname,detStr,destPath) 
     664        //              else 
     665                                // beam center is in pixels, so use the old routine 
     666                                V_ConvertBeamCtrPix_to_mmB(fname,detStr,destPath) 
     667         
     668        //              endif 
     669                        V_Detector_CalcQVals(fname,detStr,destPath) 
     670                 
     671                endif 
    664672                 
    665673        else 
    … …  
    686694                        ctTime = V_getCount_time(fname) 
    687695 
    688                         if(cmpstr(detStr,"B") == 0) 
    689                                 Variable b_dt = V_getDetector_deadtime_B(fname,detStr) 
    690                                 // do the correction for the back panel 
    691                                 countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts 
    692  
    693                                 w = w/(1-countRate*b_dt) 
    694                                 w_err = w_err/(1-countRate*b_dt) 
    695                                                                  
     696                        if(cmpstr(detStr,"B") == 0 ) 
     697                                if(gIgnoreDetB == 0) 
     698                                        Variable b_dt = V_getDetector_deadtime_B(fname,detStr) 
     699                                        // do the correction for the back panel 
     700                                        countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts 
     701         
     702                                        w = w/(1-countRate*b_dt) 
     703                                        w_err = w_err/(1-countRate*b_dt) 
     704                                endif                    
    696705                        else 
    697706                                // do the corrections for 8 tube panels 
    … …  
    720729                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 
    721730                        detStr = StringFromList(ii, ksDetectorListAll, ";") 
    722                         Wave w = V_getDetectorDataW(fname,detStr) 
    723                         Wave w_err = V_getDetectorDataErrW(fname,detStr) 
    724                         // any other dimensions to pass in? 
    725                         V_SolidAngleCorrection(w,w_err,fname,detStr,destPath) 
     731                         
     732                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1) 
     733                                // do nothing 
     734                        else 
     735                                Wave w = V_getDetectorDataW(fname,detStr) 
     736                                Wave w_err = V_getDetectorDataErrW(fname,detStr) 
     737                                // any other dimensions to pass in? 
     738                                V_SolidAngleCorrection(w,w_err,fname,detStr,destPath) 
     739                        endif 
    726740                         
    727741                endfor 
    … …  
    756770                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 
    757771                        detStr = StringFromList(ii, ksDetectorListAll, ";") 
    758                         Wave w = V_getDetectorDataW(fname,detStr) 
    759                         Wave w_err = V_getDetectorDataErrW(fname,detStr) 
    760                          
    761                         V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath) 
    762                          
     772                         
     773                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1) 
     774                                // do nothing 
     775                        else 
     776                                Wave w = V_getDetectorDataW(fname,detStr) 
     777                                Wave w_err = V_getDetectorDataErrW(fname,detStr) 
     778                                 
     779                                V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath) 
     780                        endif 
    763781                endfor 
    764782        else 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Write_VSANS_QIS.ipf

    r1108 r1119  
    340340        V_SaveIQ_ButtonProc(ba) 
    341341end 
     342 
     343 
     344///////// QxQy Export  ////////// 
     345// 
     346// (see the similar-named SANS routine for additonal steps - like resolution, etc.) 
     347//ASCII export of data as 8-columns qx-qy-Intensity-err-qz-sigmaQ_parall-sigmaQ_perp-fShad 
     348// + limited header information 
     349// 
     350//      Jan 2019 -- first version, simply exports the basic matrix of data with no resolution information 
     351// 
     352// 
     353Function V_QxQy_Export(type,fullpath,newFileName,dialog) 
     354        String type,fullpath,newFileName 
     355        Variable dialog         //=1 will present dialog for name 
     356         
     357        String typeStr="" 
     358        Variable refnum 
     359        String detStr="",detSavePath 
     360 
     361        SVAR gProtoStr = root:Packages:NIST:VSANS:Globals:Protocols:gProtoStr 
     362         
     363        // declare, or make a fake protocol if needed (if the export type is RAW) 
     364        String rawTag="" 
     365        if(cmpstr(type,"RAW")==0) 
     366                Make/O/T/N=(kNumProtocolSteps) proto 
     367                RawTag = "RAW Data File: "       
     368        else 
     369                Wave/T proto=$("root:Packages:NIST:VSANS:Globals:Protocols:"+gProtoStr)  
     370        endif 
     371         
     372        SVAR samFiles = $("root:Packages:NIST:VSANS:"+type+":gFileList") 
     373         
     374        //check each wave - MUST exist, or will cause a crash 
     375//      If(!(WaveExists(data))) 
     376//              Abort "data DNExist QxQy_Export()" 
     377//      Endif 
     378 
     379        if(dialog) 
     380                PathInfo/S catPathName 
     381                fullPath = DoSaveFileDialog("Save data as") 
     382                If(cmpstr(fullPath,"")==0) 
     383                        //user cancel, don't write out a file 
     384                        Close/A 
     385                        Abort "no data file was written" 
     386                Endif 
     387                //Print "dialog fullpath = ",fullpath 
     388        Endif 
     389 
     390        // data values to populate the file header 
     391        String fileName,fileDate,fileLabel 
     392        Variable monCt,lambda,offset,dist,trans,thick 
     393        Variable bCentX,bCentY,a2,a1a2_dist,deltaLam,bstop 
     394        String a1Str 
     395        Variable pixX,pixY 
     396        Variable numTextLines=19,ii,jj,kk 
     397 
     398        Make/O/T/N=(numTextLines) labelWave 
     399 
     400        // 
     401         
     402        //loop over all of the detector panels 
     403        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB 
     404 
     405        String detList 
     406        if(gIgnoreDetB) 
     407                detList = ksDetectorListNoB 
     408        else 
     409                detList = ksDetectorListAll 
     410        endif 
     411         
     412        for(kk=0;kk<ItemsInList(detList);kk+=1) 
     413 
     414                detStr = StringFromList(kk, detList, ";") 
     415                detSavePath = fullPath + "_" + detStr 
     416                 
     417                pixX = V_getDet_pixel_num_x(type,detStr) 
     418                pixY = V_getDet_pixel_num_y(type,detStr) 
     419                 
     420                fileName = newFileName 
     421                fileDate = V_getDataStartTime(type)             // already a string 
     422                fileLabel = V_getSampleDescription(type) 
     423                 
     424                monCt = V_getBeamMonNormData(type) 
     425                lambda = V_getWavelength(type) 
     426         
     427        // TODO - switch based on panel type 
     428        //      V_getDet_LateralOffset(fname,detStr) 
     429        //      V_getDet_VerticalOffset(fname,detStr) 
     430                offset = V_getDet_LateralOffset(type,detStr) 
     431         
     432                dist = V_getDet_ActualDistance(type,detStr) 
     433                trans = V_getSampleTransmission(type) 
     434                thick = V_getSampleThickness(type) 
     435                 
     436                bCentX = V_getDet_beam_center_x(type,detStr) 
     437                bCentY = V_getDet_beam_center_y(type,detStr) 
     438                a1Str = V_getSourceAp_size(type)                //already a string 
     439                a2 = V_getSampleAp2_size(type) 
     440                a1a2_dist = V_getSourceAp_distance(type) 
     441                deltaLam = V_getWavelength_spread(type) 
     442                // TODO -- decipher which beamstop, if any is actually in place 
     443        // or -- V_getBeamStopC3_size(type) 
     444                bstop = V_getBeamStopC2_size(type) 
     445                 
     446        ///////// 
     447                labelWave[0] = "FILE: "+fileName+"   CREATED: "+fileDate 
     448                labelWave[1] = "LABEL: "+fileLabel 
     449                labelWave[2] = "MON CNT   LAMBDA (A)  DET_OFF(cm)   DET_DIST(cm)   TRANS   THICK(cm)" 
     450                labelWave[3] = num2str(monCt)+"  "+num2str(lambda)+"       "+num2str(offset)+"     "+num2str(dist) 
     451                labelWave[3] += "     "+num2str(trans)+"     "+num2str(thick) 
     452                labelWave[4] = "BCENT(X,Y)(cm)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)" 
     453                labelWave[5] = num2str(bCentX)+"  "+num2str(bCentY)+"  "+a1Str+"    "+num2str(a2)+"    " 
     454                labelWave[5] += num2str(a1a2_dist)+"    "+num2str(deltaLam)+"    "+num2str(bstop) 
     455                labelWave[6] =  "SAM: "+rawTag+samFiles 
     456                labelWave[7] =  "BGD: "+proto[0] 
     457                labelWave[8] =  "EMP: "+proto[1] 
     458                labelWave[9] =  "DIV: "+proto[2] 
     459                labelWave[10] =  "MASK: "+proto[3] 
     460                labelWave[11] =  "ABS Parameters (3-6): "+proto[4] 
     461                labelWave[12] = "Average Choices: "+proto[5] 
     462                labelWave[13] = "Collimation type: "+proto[9] 
     463                labelWave[14] = "" 
     464                labelWave[15] = "*** Data written from "+type+" folder and may not be a fully corrected data file ***" 
     465//              labelWave[16] = "Data columns are Qx - Qy - Qz - I(Qx,Qy) - Err I(Qx,Qy)" 
     466        //      labelWave[16] = "Data columns are Qx - Qy - I(Qx,Qy) - Qz - SigmaQ_parall - SigmaQ_perp - fSubS(beam stop shadow)" 
     467                labelWave[16] = "Data columns are Qx - Qy - I(Qx,Qy) - err(I) - Qz - SigmaQ_parall - SigmaQ_perp - fSubS(beam stop shadow)" 
     468                labelWave[17] = "The error wave may not be properly propagated (1/2019)" 
     469                labelWave[18] = "ASCII data created " +date()+" "+time() 
     470                //strings can be too long to print-- must trim to 255 chars 
     471                for(jj=0;jj<numTextLines;jj+=1) 
     472                        labelWave[jj] = (labelWave[jj])[0,240] 
     473                endfor 
     474         
     475         
     476        // get the data waves for output 
     477        // QxQyQz have already been calculated for VSANS data 
     478                 
     479                WAVE data = V_getDetectorDataW(type,detStr) 
     480                WAVE data_err = V_getDetectorDataErrW(type,detStr) 
     481                 
     482                // TOOD - replace hard wired paths with Read functions 
     483                // hard-wired 
     484                Wave qx_val = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qx_"+detStr) 
     485                Wave qy_val = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qy_"+detStr) 
     486                Wave qz_val = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qz_"+detStr) 
     487                Wave qTot = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qTot_"+detStr) 
     488                 
     489        ///// calculation of the resolution function (2D) 
     490 
     491        // 
     492                Variable acc,ssd,lambda0,yg_d,qstar,g,L1,L2,vz_1,sdd 
     493                // L1 = source to sample distance [cm]  
     494                L1 = V_getSourceAp_distance(type) 
     495         
     496        // L2 = sample to detector distance [cm] 
     497                L2 = V_getDet_ActualDistance(type,detStr)               //cm 
     498 
     499        //               
     500                G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2 
     501                vz_1 =  3.956E5 //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC 
     502                acc = vz_1 
     503                SDD = L2                //1317 
     504                SSD = L1                //1627          //cm 
     505                lambda0 = lambda                //              15 
     506                YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2 
     507                Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG_d 
     508        ////            Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd 
     509                qstar = -2*pi/lambda0*2*yg_d/sdd 
     510        //       
     511        // 
     512        //// the gravity center is not the resolution center 
     513        //// gravity center = beam center 
     514        //// resolution center = offset y = dy + (2)*yg_d 
     515        /////************ 
     516        //// do everything to write out the resolution too 
     517        //      // un-comment these if you want to write out qz_val and qval too, then use the proper save command 
     518        //      qval = CalcQval(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     519                Duplicate/O qTot,phi,r_dist 
     520                Variable pixSizeX,pixSizeY 
     521                pixSizeX = V_getDet_x_pixel_size(type,detStr) 
     522                pixSizeY = V_getDet_y_pixel_size(type,detStr) 
     523 
     524                phi = V_FindPhi( pixSizeX*((p+1)-bCentX) , pixSizeY*((q+1)-bCentY)+(2)*yg_d)            //(dx,dy+yg_d) 
     525                r_dist = sqrt(  (pixSizeX*((p+1)-bCentX))^2 +  (pixSizeY*((q+1)-bCentY)+(2)*yg_d)^2 )           //radial distance from ctr to pt 
     526         
     527                //make everything in 1D now 
     528                Duplicate/O qTot SigmaQX,SigmaQY,fsubS,qval      
     529                Redimension/N=(pixX*pixY) SigmaQX,SigmaQY,fsubS,qval,phi,r_dist 
     530 
     531                Variable ret1,ret2,ret3,nq 
     532                String collimationStr 
     533                 
     534                 
     535                collimationStr = proto[9] 
     536                 
     537                nq = pixX*pixY 
     538                ii=0 
     539 
     540s_tic()  
     541// TODO 
     542// this loop is the slow step. it takes Å 0.7 s for F or M panels, and Å 120 s for the Back panel (6144 pts vs. 1.12e6 pts) 
     543// find some way to speed this up! 
     544// MultiThreading will be difficult as it requires all the dependent functions (HDF5 reads, etc.) to be threadsafe as well 
     545// and there are a lot of them... 
     546                //type = work folder 
     547                 
     548//              (this doesn't work...and isn't any faster) 
     549//              Duplicate/O qval dum 
     550//              dum = V_get2DResolution(qval,phi,r_dist,type,detStr,collimationStr,SigmaQX,SigmaQY,fsubS) 
     551 
     552                do 
     553                        V_get2DResolution(qval[ii],phi[ii],r_dist[ii],type,detStr,collimationStr,ret1,ret2,ret3) 
     554                        SigmaQX[ii] = ret1       
     555                        SigmaQY[ii] = ret2       
     556                        fsubs[ii] = ret3         
     557                        ii+=1 
     558                while(ii<nq)     
     559s_toc() 
     560         
     561        ////*********************        
     562                Duplicate/O qx_val,qx_val_s 
     563                Duplicate/O qy_val,qy_val_s 
     564                Duplicate/O qz_val,qz_val_s 
     565                Duplicate/O data,z_val_s 
     566                Duplicate/O SigmaQx,sigmaQx_s 
     567                Duplicate/O SigmaQy,sigmaQy_s 
     568                Duplicate/O fSubS,fSubS_s 
     569                Duplicate/O data_err,sw_s 
     570                 
     571                //so that double precision data is not written out 
     572                Redimension/S qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s 
     573                Redimension/S SigmaQx_s,SigmaQy_s,fSubS_s 
     574         
     575                Redimension/N=(pixX*pixY) qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s 
     576                 
     577                //not demo-compatible, but approx 8x faster!!    
     578#if(strsearch(stringbykey("IGORKIND",IgorInfo(0),":",";"), "demo", 0 ) == -1) 
     579                 
     580//              Save/O/G/M="\r\n" labelWave,qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s as detSavePath      // without resolution 
     581                Save/O/G/M="\r\n" labelWave,qx_val_s,qy_val_s,z_val_s,sw_s,qz_val_s,SigmaQx_s,SigmaQy_s,fSubS_s as detSavePath  // write out the resolution information 
     582#else 
     583                Open refNum as detSavePath 
     584                wfprintf refNum,"%s\r\n",labelWave 
     585                fprintf refnum,"\r\n" 
     586//              wfprintf refNum,"%8g\t%8g\t%8g\t%8g\t%8g\r\n",qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s 
     587                wfprintf refNum,"%8g\t%8g\t%8g\t%8g\t%8g\t%8g\t%8g\t%8g\r\n",qx_val,qy_val,z_val,sw,qz_val,SigmaQx,SigmaQy,fSubS 
     588                Close refNum 
     589#endif 
     590                 
     591                KillWaves/Z qx_val_s,qy_val_s,z_val_s,qz_val_s,SigmaQx_s,SigmaQy_s,fSubS_s,sw,sw_s 
     592                 
     593                Killwaves/Z qx_val,qy_val,z_val,qval,qz_val,sigmaQx,SigmaQy,fSubS 
     594                 
     595                Print "QxQy_Export File written: ", V_GetFileNameFromPathNoSemi(detSavePath) 
     596         
     597        endfor 
     598         
     599        KillWaves/Z labelWave,dum 
     600        return(0) 
     601End 
     602 
     603// this assumes that: 
     604// --QxQy data was written out in the format specified by the Igor macros, that is the x varies most rapidly 
     605// 
     606// TODO -- this needs to be made generic for reading in different panels with different XY dimensions 
     607// -- add the XY dimensions to the QxQyASCII file header somewhere so that it can be read in and used here 
     608// 
     609// the SANS analysis 2D loader assumes that the matrix is square, mangling the VSANS data. 
     610// the column data (for fitting) is still fine, but the matrix representation is incorrect. 
     611// 
     612Function V_ConvertQxQy2Mat(Qx,Qy,inten,matStr) 
     613        Wave Qx,Qy,inten 
     614        String matStr 
     615         
     616        String folderStr=GetWavesDataFolder(Qx,1) 
     617         
     618        Variable numX,numY 
     619        numX=48 
     620        numY=128 
     621        Make/O/D/N=(numX,numY) $(folderStr + matStr) 
     622        Wave mat=$(folderStr + matStr) 
     623         
     624        WaveStats/Q Qx 
     625        SetScale/I x, V_min, V_max, "", mat 
     626        WaveStats/Q Qy 
     627        SetScale/I y, V_min, V_max, "", mat 
     628         
     629        Variable xrows=numX 
     630         
     631        mat = inten[q*xrows+p] 
     632         
     633        return(0) 
     634End 
     635 
     636 
Note: See TracChangeset for help on using the changeset viewer.