Ignore:
Timestamp:
Jun 16, 2010 2:10:47 PM (12 years ago)
Author:
srkline
Message:

SA_Includes_v410 : now include Smear_2D

PeakGauss_2D, Sphere_2D : included threaded resolution smearing calculations for testing

DataSetHandling? : Included a quick and dirty batch converter for XML->6-col. See the top
of the file for the command to run

GaussUtils? : re-define the ResSemear_2D_AAOStruct. Relocated q-value and phi calculations from
RawWindowHook? to this file so they would be available for reduction and analysis

Smear_2D : now has a generic (non-threaded) smearing routine. Threading must be done in
individual functions since FUNCREF can't be passed to threads (plus a few other issues)

PlotUtils_2D : updated loader for new QxQy? columns. Fixes to Wrapper_2D to enable smeared fits

RawWindowHook? : removed Q-value calculation functions and moved these to GaussUtils?

WriteQIS : now writes out 8-columns for QxQy? data, defining the resolution
function in terms of directions parallel and perpendicular to Q. TEMPORARILY in the data
file an error in intensity is generated that is SQRT(I), being careful to
replace any NaN, inf, or zero with an average error value

MultiScatter_MonteCarlo_2D : 4-processor aware

NCNR_Utils : 2D resolution calculation is now in terms of parallel and perpendicular
rather than x and y. Gravity is included in the y-component

Location:
sans/Dev/trunk/NCNR_User_Procedures/Common
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Common/DataSetHandling.ipf

    r681 r708  
    11#pragma rtGlobals=1             // Use modern global access method. 
     2 
     3 
     4///////// SRK - VERY SIMPLE batch converter has been added: see 
     5// 
     6//      Function batchXML26ColConvert() 
     7// 
     8 
     9 
    210 
    311// Functions and interfaces to manage datasets now that they are in data folders 
     
    15091517// 
    15101518//End 
     1519 
     1520 
     1521 
     1522///////// SRK - VERY SIMPLE batch converter 
     1523// no header information is preserved 
     1524// file names are partially preserved 
     1525// 
     1526 
     1527/// to use this: 
     1528// -open the Plot Manager and set the path 
     1529// -run this function 
     1530// 
     1531// it doesn't matter if the XML ouput flag is set - this overrides. 
     1532Function batchXML26ColConvert() 
     1533 
     1534        String list, item,path,fname 
     1535        Variable num,ii 
     1536         
     1537        PathInfo CatPathName 
     1538        path = S_Path 
     1539 
     1540        list = A_ReducedDataFileList("") 
     1541        num = itemsInList(list) 
     1542        Print num 
     1543        for(ii=0;ii<num;ii+=1) 
     1544                item = StringFromList(ii, list ,";") 
     1545                fname=path + item 
     1546                Execute "A_LoadOneDDataWithName(\""+fname+"\",0)"               //won't plot 
     1547//              easier to load all, then write out, since the name will be changed 
     1548        endfor 
     1549         
     1550         
     1551        list = DM_DataSetPopupList() 
     1552 
     1553        num = itemsInList(list) 
     1554        Print num 
     1555        for(ii=0;ii<num;ii+=1) 
     1556                item = StringFromList(ii, list ,";") 
     1557                fReWrite1DData_noPrompt(item,"tab","CR") 
     1558        endfor 
     1559         
     1560End 
     1561 
     1562// quick version (copied from fReWrite1DData() that NEVER asks for a new fileName 
     1563// - and right now, always expect 6-column data, either SANS or USANS (re-writes -dQv) 
     1564// - AJJ Nov 2009 : better make sure we always fake 6 columns on reading then.... 
     1565Function fReWrite1DData_noPrompt(folderStr,delim,term) 
     1566        String folderStr,delim,term 
     1567         
     1568        String formatStr="",fullpath="" 
     1569        Variable refnum,dialog=1 
     1570         
     1571        String dataSetFolderParent,basestr 
     1572         
     1573        //setup delimeter and terminator choices 
     1574        If(cmpstr(delim,"tab")==0) 
     1575                //tab-delimeted 
     1576                formatStr="%15.8g\t%15.8g\t%15.8g\t%15.8g\t%15.8g\t%15.8g" 
     1577        else 
     1578                //use 3 spaces 
     1579                formatStr="%15.8g   %15.8g   %15.8g   %15.8g   %15.8g   %15.8g" 
     1580        Endif 
     1581        If(cmpstr(term,"CR")==0) 
     1582                formatStr += "\r" 
     1583        Endif 
     1584        If(cmpstr(term,"LF")==0) 
     1585                formatStr += "\n" 
     1586        Endif 
     1587        If(cmpstr(term,"CRLF")==0) 
     1588                formatStr += "\r\n" 
     1589        Endif 
     1590         
     1591        //Abuse ParseFilePath to get path without folder name 
     1592        dataSetFolderParent = ParseFilePath(1,folderStr,":",1,0) 
     1593        //Abuse ParseFilePath to get basestr 
     1594        basestr = ParseFilePath(0,folderStr,":",1,0) 
     1595         
     1596        //make sure the waves exist 
     1597        SetDataFolder $(dataSetFolderParent+basestr) 
     1598        WAVE/Z qw = $(baseStr+"_q") 
     1599        WAVE/Z iw = $(baseStr+"_i") 
     1600        WAVE/Z sw = $(baseStr+"_s") 
     1601        WAVE/Z resw = $(baseStr+"_res") 
     1602         
     1603        if(WaveExists(qw) == 0) 
     1604                Abort "q is missing" 
     1605        endif 
     1606        if(WaveExists(iw) == 0) 
     1607                Abort "i is missing" 
     1608        endif 
     1609        if(WaveExists(sw) == 0) 
     1610                Abort "s is missing" 
     1611        endif 
     1612        if(WaveExists(resw) == 0) 
     1613                Abort "Resolution information is missing." 
     1614        endif 
     1615         
     1616        Duplicate/O qw qbar,sigQ,fs 
     1617        if(dimsize(resW,1) > 4) 
     1618                //it's USANS put -dQv back in the last 3 columns 
     1619                NVAR/Z dQv = USANS_dQv 
     1620                if(NVAR_Exists(dQv) == 0) 
     1621                        Abort "It's USANS data, and I don't know what the slit height is." 
     1622                endif 
     1623                sigQ = -dQv 
     1624                qbar = -dQv 
     1625                fs = -dQv 
     1626        else 
     1627                //it's SANS 
     1628                sigQ = resw[p][0] 
     1629                qbar = resw[p][1] 
     1630                fs = resw[p][2] 
     1631        endif 
     1632         
     1633        dialog=0 
     1634        if(dialog) 
     1635                PathInfo/S catPathName 
     1636//              fullPath = DoSaveFileDialog("Save data as",fname=baseStr+".txt") 
     1637                fullPath = DoSaveFileDialog("Save data as",fname=baseStr[0,strlen(BaseStr)-2]) 
     1638                Print fullPath 
     1639                If(cmpstr(fullPath,"")==0) 
     1640                        //user cancel, don't write out a file 
     1641                        Close/A 
     1642                        Abort "no data file was written" 
     1643                Endif 
     1644                //Print "dialog fullpath = ",fullpath 
     1645        Endif 
     1646        PathInfo catPathName 
     1647        fullPath = S_Path + baseStr[0,strlen(BaseStr)-2] 
     1648 
     1649        Open refnum as fullpath 
     1650         
     1651        fprintf refnum,"Modified data written from folder %s on %s\r\n",baseStr,(date()+" "+time()) 
     1652        wfprintf refnum,formatStr,qw,iw,sw,sigQ,qbar,fs 
     1653        Close refnum 
     1654         
     1655        KillWaves/Z sigQ,qbar,fs 
     1656         
     1657        SetDataFolder root: 
     1658        return(0) 
     1659End 
     1660 
     1661///////////end SRK 
  • sans/Dev/trunk/NCNR_User_Procedures/Common/GaussUtils_v40.ipf

    r570 r708  
    6262 
    6363 
    64 // tentative pass at 2D resolution smearing 
     64//// tentative pass at 2D resolution smearing 
     65//// 
     66//Structure ResSmear_2D_AAOStruct 
     67//      Wave coefW 
     68//      Wave zw                 //answer 
     69//      Wave qy                 // q-value 
     70//      Wave qx 
     71//      Wave qz 
     72//      Wave sQpl               //resolution parallel to Q 
     73//      Wave sQpp               //resolution perpendicular to Q 
     74//      Wave fs 
     75//      String info 
     76//EndStructure 
     77 
     78// reformat the structure ?? WM Fit compatible 
     79// -- 2D resolution smearing 
    6580// 
    6681Structure ResSmear_2D_AAOStruct 
    6782        Wave coefW 
    6883        Wave zw                 //answer 
    69         Wave qy                 // q-value 
    70         Wave qx 
     84        Wave xw[2]              // qx-value is [0], qy is xw[1] 
    7185        Wave qz 
    72         Wave sigQx              //resolution 
    73         Wave sigQy 
     86        Wave sQpl               //resolution parallel to Q 
     87        Wave sQpp               //resolution perpendicular to Q 
    7488        Wave fs 
    7589        String info 
     
    10771091 
    10781092// prototype function for 2D smearing routine 
    1079 Function SANS_2D_ModelAAO_proto(w,zw,xw,yw) 
     1093ThreadSafe Function SANS_2D_ModelAAO_proto(w,zw,xw,yw) 
    10801094        Wave w,zw,xw,yw 
    10811095         
     
    13461360        return(0) 
    13471361end 
     1362 
     1363 
     1364//// 
     1365//// moved from RawWindowHook to here - where the Q calculations are available to 
     1366//   reduction and analysis 
     1367// 
     1368 
     1369//phi is defined from +x axis, proceeding CCW around [0,2Pi] 
     1370Threadsafe Function FindPhi(vx,vy) 
     1371        variable vx,vy 
     1372         
     1373        variable phi 
     1374         
     1375        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2 
     1376         
     1377        // special cases 
     1378        if(vx==0 && vy > 0) 
     1379                return(pi/2) 
     1380        endif 
     1381        if(vx==0 && vy < 0) 
     1382                return(3*pi/2) 
     1383        endif 
     1384        if(vx >= 0 && vy == 0) 
     1385                return(0) 
     1386        endif 
     1387        if(vx < 0 && vy == 0) 
     1388                return(pi) 
     1389        endif 
     1390         
     1391         
     1392        if(vx > 0 && vy > 0) 
     1393                return(phi) 
     1394        endif 
     1395        if(vx < 0 && vy > 0) 
     1396                return(phi + pi) 
     1397        endif 
     1398        if(vx < 0 && vy < 0) 
     1399                return(phi + pi) 
     1400        endif 
     1401        if( vx > 0 && vy < 0) 
     1402                return(phi + 2*pi) 
     1403        endif 
     1404         
     1405        return(phi) 
     1406end 
     1407 
     1408         
     1409//function to calculate the overall q-value, given all of the necesary trig inputs 
     1410//NOTE: detector locations passed in are pixels = 0.5cm real space on the detector 
     1411//and are in detector coordinates (1,128) rather than axis values 
     1412//the pixel locations need not be integers, reals are ok inputs 
     1413//sdd is in meters 
     1414//wavelength is in Angstroms 
     1415// 
     1416//returned magnitude of Q is in 1/Angstroms 
     1417// 
     1418Function CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     1419        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
     1420         
     1421        Variable dx,dy,qval,two_theta,dist 
     1422         
     1423        Variable pixSizeX=pixSize 
     1424        Variable pixSizeY=pixSize 
     1425         
     1426        sdd *=100               //convert to cm 
     1427        dx = (xaxval - xctr)*pixSizeX           //delta x in cm 
     1428        dy = (yaxval - yctr)*pixSizeY           //delta y in cm 
     1429        dist = sqrt(dx^2 + dy^2) 
     1430         
     1431        two_theta = atan(dist/sdd) 
     1432 
     1433        qval = 4*Pi/lam*sin(two_theta/2) 
     1434         
     1435        return qval 
     1436End 
     1437 
     1438//calculates just the q-value in the x-direction on the detector 
     1439//input/output is the same as CalcQval() 
     1440//ALL inputs are in detector coordinates 
     1441// 
     1442//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
     1443//sdd is in meters 
     1444//wavelength is in Angstroms 
     1445// 
     1446// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) 
     1447// now properly accounts for qz 
     1448// 
     1449Function CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     1450        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
     1451 
     1452        Variable qx,qval,phi,dx,dy,dist,two_theta 
     1453         
     1454        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     1455         
     1456        sdd *=100               //convert to cm 
     1457        dx = (xaxval - xctr)*pixSize            //delta x in cm 
     1458        dy = (yaxval - yctr)*pixSize            //delta y in cm 
     1459        phi = FindPhi(dx,dy) 
     1460         
     1461        //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
     1462        dist = sqrt(dx^2 + dy^2) 
     1463        two_theta = atan(dist/sdd) 
     1464 
     1465        qx = qval*cos(two_theta/2)*cos(phi) 
     1466         
     1467        return qx 
     1468End 
     1469 
     1470//calculates just the q-value in the y-direction on the detector 
     1471//input/output is the same as CalcQval() 
     1472//ALL inputs are in detector coordinates 
     1473//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
     1474//sdd is in meters 
     1475//wavelength is in Angstroms 
     1476// 
     1477// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) 
     1478// now properly accounts for qz 
     1479// 
     1480Function CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     1481        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
     1482         
     1483        Variable dy,qval,dx,phi,qy,dist,two_theta 
     1484         
     1485        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     1486         
     1487        sdd *=100               //convert to cm 
     1488        dx = (xaxval - xctr)*pixSize            //delta x in cm 
     1489        dy = (yaxval - yctr)*pixSize            //delta y in cm 
     1490        phi = FindPhi(dx,dy) 
     1491         
     1492        //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
     1493        dist = sqrt(dx^2 + dy^2) 
     1494        two_theta = atan(dist/sdd) 
     1495         
     1496        qy = qval*cos(two_theta/2)*sin(phi) 
     1497         
     1498        return qy 
     1499End 
     1500 
     1501//calculates just the z-component of the q-vector, not measured on the detector 
     1502//input/output is the same as CalcQval() 
     1503//ALL inputs are in detector coordinates 
     1504//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
     1505//sdd is in meters 
     1506//wavelength is in Angstroms 
     1507// 
     1508// not actually used, but here for completeness if anyone asks 
     1509// 
     1510Function CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     1511        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
     1512         
     1513        Variable dy,qval,dx,phi,qz,dist,two_theta 
     1514         
     1515        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
     1516         
     1517        sdd *=100               //convert to cm 
     1518         
     1519        //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
     1520        dx = (xaxval - xctr)*pixSize            //delta x in cm 
     1521        dy = (yaxval - yctr)*pixSize            //delta y in cm 
     1522        dist = sqrt(dx^2 + dy^2) 
     1523        two_theta = atan(dist/sdd) 
     1524         
     1525        qz = qval*sin(two_theta/2) 
     1526         
     1527        return qz 
     1528End 
     1529 
     1530//for command-line testing, replace the function declaration 
     1531//Function FindQxQy(qq,phi) 
     1532//      Variable qq,phi 
     1533//      Variable qx,qy 
     1534// 
     1535// 
     1536ThreadSafe Function FindQxQy(qq,phi,qx,qy) 
     1537        Variable qq,phi,&qx,&qy 
     1538 
     1539        qx = sqrt(qq^2/(1+tan(phi)*tan(phi))) 
     1540        qy = qx*tan(phi) 
     1541         
     1542        if(phi >= 0 && phi <= pi/2) 
     1543                qx = abs(qx) 
     1544                qy = abs(qy) 
     1545        endif 
     1546         
     1547        if(phi > pi/2 && phi <= pi) 
     1548                qx = -abs(qx) 
     1549                qy = abs(qy) 
     1550        endif 
     1551         
     1552        if(phi > pi && phi <= pi*3/2) 
     1553                qx = -abs(qx) 
     1554                qy = -abs(qy) 
     1555        endif 
     1556         
     1557        if(phi > pi*3/2 && phi < 2*pi) 
     1558                qx = abs(qx) 
     1559                qy = -abs(qy) 
     1560        endif    
     1561         
     1562         
     1563//      Print "recalculated qx,qy,q = ",qx,qy,sqrt(qx*qx+qy*qy) 
     1564         
     1565        return(0) 
     1566end 
     1567         
  • sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtils2D_v40.ipf

    r665 r708  
    44 
    55// 
    6 // WARNING: there are a lot of assumptions that went into this code at the first pass 
    7 // that will need to be made generic before a release to the wild 
    8 // 
    9 // 
    10 // - some manipulations assume a square matrix 
    11 // - naming schemes are assumed (to be documented) 
    12 // 
    13 // 
    14  
    15  
    16 // 
    17 // changed May 2009 to read in the resolution information (now 7 columns) 
     6// TO DO: 
     7// 
     8// - more intelligent beam stop masking 
     9// - constraints 
     10// - interactive masking 
     11// 
     12// 
     13// 
     14 
     15 
     16// 
     17// changed June 2010 to read in the resolution information (now 8! columns) 
    1818// -- subject to change -- 
    1919// 
     
    2525        Variable numCols = V_flag 
    2626 
    27         String w0,w1,w2,n0,n1,n2 
    28         String w3,w4,w5,w6,n3,n4,n5,n6 
     27        String w0,w1,w2,w3,w4,w5,w6,w7 
     28        String n0,n1,n2,n3,n4,n5,n6,n7 
    2929                 
    3030        // put the names of the three loaded waves into local names 
     
    3636        n5 = StringFromList(5, S_waveNames ,";" ) 
    3737        n6 = StringFromList(6, S_waveNames ,";" ) 
     38        n7 = StringFromList(7, S_waveNames ,";" ) 
    3839         
    3940        //remove the semicolon AND period from files from the VAX 
     
    4142        w1 = CleanupName((S_fileName + "_qy"),0) 
    4243        w2 = CleanupName((S_fileName + "_i"),0) 
    43         w3 = CleanupName((S_fileName + "_qz"),0) 
    44         w4 = CleanupName((S_fileName + "_sigQx"),0) 
    45         w5 = CleanupName((S_fileName + "_sigQy"),0) 
    46         w6 = CleanupName((S_fileName + "_fs"),0) 
     44        w3 = CleanupName((S_fileName + "_iErr"),0) 
     45        w4 = CleanupName((S_fileName + "_qz"),0) 
     46        w5 = CleanupName((S_fileName + "_sQpl"),0) 
     47        w6 = CleanupName((S_fileName + "_sQpp"),0) 
     48        w7 = CleanupName((S_fileName + "_fs"),0) 
    4749 
    4850        String baseStr=w1[0,strlen(w1)-4] 
     
    5153                        if(V_flag==2)   //user selected No, don't load the data 
    5254                                SetDataFolder root: 
    53                                 KillWaves $n0,$n1,$n2,$n3,$n4,$n5,$n6           // kill the default waveX that were loaded 
     55                                KillWaves $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7               // kill the default waveX that were loaded 
    5456                                return          //quits the macro 
    5557                        endif 
     
    8486        Duplicate/O $("root:"+n5), $w5 
    8587        Duplicate/O $("root:"+n6), $w6 
     88        Duplicate/O $("root:"+n7), $w7 
    8689         
    8790        Variable/G gIsLogScale = 0 
     
    9497         
    9598        PlotQxQy(baseStr)               //this sets the data folder back to root:!! 
    96          
    97         //don't create the triplet - not really of much use 
    98 //      Make/D/O/N=(num,3) $(baseStr+"_tri") 
    99 //      $(baseStr+"_tri")[][0] = $w0[p]         // Qx 
    100 //      $(baseStr+"_tri")[][1] = $w1[p]         // Qy 
    101 //      $(baseStr+"_tri")[][2] = $w2[p]         // Intensity 
    10299 
    103100        //clean up               
    104101        SetDataFolder root: 
    105         KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6 
     102        KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7 
    106103EndMacro 
    107104 
     
    547544        String DF="root:"+folderStr+":"  
    548545         
    549         Struct ResSmearAAOStruct fs 
    550 //      WAVE resW = $(DF+folderStr+"_res")               
    551 //      WAVE fs.resW =  resW 
    552546        WAVE inten=$(DF+folderStr+"_i") 
    553         WAVE Qx=$(DF+folderStr+"_qx") 
    554         WAVE Qy=$(DF+folderStr+"_qy") 
    555 //      Wave fs.coefW = cw 
    556 //      Wave fs.yW = yw 
    557 //      Wave fs.xW = xw 
    558          
    559         // generate my own error wave for I(qx,qy) 
    560         Duplicate/O inten sw 
    561         sw = sqrt(sw)           //assumes Poisson statistics for each cell (counter) 
    562 //      sw = 0.05*sw            // uniform 5% error? tends to favor the low intensity too strongly 
     547        WAVE sw=$(DF+folderStr+"_iErr") 
     548        WAVE qx=$(DF+folderStr+"_qx") 
     549        WAVE qy=$(DF+folderStr+"_qy") 
     550        WAVE qz=$(DF+folderStr+"_qz") 
     551        WAVE sQpl=$(DF+folderStr+"_sQpl") 
     552        WAVE sQpp=$(DF+folderStr+"_sQpp") 
     553        WAVE shad=$(DF+folderStr+"_fs") 
     554 
     555//just a dummy - I shouldn't need this 
     556        Duplicate/O qx resultW 
     557        resultW=0 
     558         
     559        STRUCT ResSmear_2D_AAOStruct s 
     560        WAVE s.coefW = cw        
     561        WAVE s.zw = resultW      
     562        WAVE s.xw[0] = qx 
     563        WAVE s.xw[1] = qy 
     564        WAVE s.qz = qz 
     565        WAVE s.sQpl = sQpl 
     566        WAVE s.sQpp = sQpp 
     567        WAVE s.fs = shad 
     568         
    563569 
    564570        // generate my own mask wave - as a matrix first, then redimension to N*N vector 
     
    571577        Endif 
    572578         
    573  
    574          
    575 //      Duplicate/O yw $(DF+"FitYw") 
    576 //      WAVE fitYw = $(DF+"FitYw") 
    577 //      fitYw = NaN 
     579         
     580        Duplicate/O inten inten_masked 
     581        inten_masked = (mask[p][q] == 0) ? NaN : inten[p][q] 
     582         
    578583 
    579584        //for now, use res will always be 0 for 2D functions     
     
    582587                useRes=1 
    583588        endif 
    584          
    585         // do not construct constraints for any of the coefficients that are being held 
    586         // -- this will generate an "unknown error" from the curve fitting 
    587         Make/O/T/N=0 constr 
     589 
     590        // can't use constraints in this way for multivariate fits. See the curve fitting help file 
     591        // and "Contraint Matrix and Vector" 
    588592        if(useConstr) 
    589                 String constraintExpression 
    590                 Variable i, nPnts=DimSize(lolim, 0),nextRow=0 
    591                 for (i=0; i < nPnts; i += 1) 
    592                         if (strlen(lolim[i]) > 0 && hold[i] == 0) 
    593                                 InsertPoints nextRow, 1, constr 
    594                                 sprintf constraintExpression, "K%d > %s", i, lolim[i] 
    595                                 constr[nextRow] = constraintExpression 
    596                                 nextRow += 1 
    597                         endif 
    598                         if (strlen(hilim[i]) > 0 && hold[i] == 0) 
    599                                 InsertPoints nextRow, 1, constr 
    600                                 sprintf constraintExpression, "K%d < %s", i, hilim[i] 
    601                                 constr[nextRow] = constraintExpression 
    602                                 nextRow += 1 
    603                         endif 
    604                 endfor 
    605         endif 
     593                Print "Constraints not yet implemented" 
     594                useConstr = 0 
     595        endif    
     596         
     597//      // do not construct constraints for any of the coefficients that are being held 
     598//      // -- this will generate an "unknown error" from the curve fitting 
     599//      Make/O/T/N=0 constr 
     600//      if(useConstr) 
     601//              String constraintExpression 
     602//              Variable i, nPnts=DimSize(lolim, 0),nextRow=0 
     603//              for (i=0; i < nPnts; i += 1) 
     604//                      if (strlen(lolim[i]) > 0 && hold[i] == 0) 
     605//                              InsertPoints nextRow, 1, constr 
     606//                              sprintf constraintExpression, "K%d > %s", i, lolim[i] 
     607//                              constr[nextRow] = constraintExpression 
     608//                              nextRow += 1 
     609//                      endif 
     610//                      if (strlen(hilim[i]) > 0 && hold[i] == 0) 
     611//                              InsertPoints nextRow, 1, constr 
     612//                              sprintf constraintExpression, "K%d < %s", i, hilim[i] 
     613//                              constr[nextRow] = constraintExpression 
     614//                              nextRow += 1 
     615//                      endif 
     616//              endfor 
     617//      endif 
    606618 
    607619///// NO CURSORS for 2D waves 
     
    631643// dispatch the fit 
    632644        //      FuncFit/H="11110111111"/NTHR=0 Cylinder2D_D :cyl2d_c_txt:coef_Cyl2D_D  :cyl2d_c_txt:cyl2d_c_txt_i /X={:cyl2d_c_txt:cyl2d_c_txt_qy,:cyl2d_c_txt:cyl2d_c_txt_qx} /W=:cyl2d_c_txt:sw /I=1 /M=:cyl2d_c_txt:mask /D  
     645        Variable t1=StopMSTimer(-2) 
     646 
     647 
     648// /NTHR=1 means just one thread for the fit (since the function itself is threaded) 
    633649 
    634650        do 
    635651                if(useRes && useEps && useCursors && useConstr)         //do it all 
    636                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps /C=constr 
     652                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s 
    637653                        break 
    638654                endif 
    639655                 
    640656                if(useRes && useEps && useCursors)              //no constr 
    641                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps 
     657                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /STRC=s 
    642658                        break 
    643659                endif 
    644660                 
    645661                if(useRes && useEps && useConstr)               //no crsr 
    646                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps /C=constr 
     662                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s 
    647663                        break 
    648664                endif 
    649665                 
    650666                if(useRes && useCursors && useConstr)           //no eps 
    651                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D /C=constr 
     667                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /C=constr /STRC=s 
    652668                        break 
    653669                endif 
    654670                 
    655671                if(useRes && useCursors)                //no eps, no constr 
    656                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D 
     672                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /STRC=s 
    657673                        break 
    658674                endif 
    659675                 
    660676                if(useRes && useEps)            //no crsr, no constr 
    661                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps 
     677                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /STRC=s 
    662678                        break 
    663679                endif 
    664680         
    665681                if(useRes && useConstr)         //no crsr, no eps 
    666                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D /C=constr 
     682                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /C=constr /STRC=s 
    667683                        break 
    668684                endif 
    669685                 
    670686                if(useRes)              //just res 
    671                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D 
     687                        Print "useRes only" 
     688                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /STRC=s 
     689//                      FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /W=sw /I=1 /STRC=s 
     690//                      FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 
     691//                      FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten_masked /X={Qx,Qy} /W=sw /I=1 /STRC=s 
    672692                        break 
    673693                endif 
     
    675695/////   same as above, but all without useRes (no /STRC flag) 
    676696                if(useEps && useCursors && useConstr)           //do it all 
    677                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps /C=constr 
     697                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr 
    678698                        break 
    679699                endif 
    680700                 
    681701                if(useEps && useCursors)                //no constr 
    682                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps 
     702                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps 
    683703                        break 
    684704                endif 
     
    686706                 
    687707                if(useEps && useConstr)         //no crsr 
    688                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps /C=constr 
     708                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr 
    689709                        break 
    690710                endif 
    691711                 
    692712                if(useCursors && useConstr)             //no eps 
    693                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D /C=constr 
     713                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /C=constr 
    694714                        break 
    695715                endif 
    696716                 
    697717                if(useCursors)          //no eps, no constr 
    698                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 /D 
     718                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten[pcsr(A),pcsr(B)] /X={Qx,Qy} /M=mask /W=sw /I=1 
    699719                        break 
    700720                endif 
    701721                 
    702722                if(useEps)              //no crsr, no constr 
    703                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D /E=eps 
     723                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps 
    704724                        break 
    705725                endif 
    706726         
    707727                if(useConstr)           //no crsr, no eps 
    708                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D /C=constr 
     728                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /C=constr 
    709729                        break 
    710730                endif 
    711731                 
    712732                //just a plain vanilla fit 
    713                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D 
     733                print "unsmeared vanilla" 
     734                        FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 
     735//                      FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten /X={Qx,Qy} /W=sw /I=1 
     736//                      FuncFit/H=getHStr(hold) /NTHR=1 $funcStr cw, inten_masked /X={Qx,Qy} /W=sw /I=1 
    714737         
    715738        while(0) 
    716739         
     740        Print "elapsed 2D fit time  = ",(StopMSTimer(-2) - t1)/1e6," s = ",(StopMSTimer(-2) - t1)/1e6/60," min" 
     741 
    717742        // append the fit 
    718743        // need to manage duplicate copies 
     
    754779         
    755780        if(yesReport) 
    756                 String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //this is *hopefully* one wave 
     781                String parStr = getFunctionParams(funcStr) 
     782//              String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //old way, but doesn't work in 2D folders 
    757783                String topGraph= TopGizmoWindow()       //this is the topmost Gizmo (XOP) window 
    758784                 
  • sans/Dev/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf

    r570 r708  
    22#pragma IgorVersion=6.1 
    33 
    4 Function TestSmear_2D() 
    5  
    6         Variable DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA 
    7         DX = 0.5 
    8         num = 128 
    9 //      x0 = 64 
    10         x0 = 114 
    11         y0 = 64 
    12         L1 = 300                //units of cm ?? 
    13         L2 = 130 
    14         s1 = 5.0/2 
    15         s2 = 1.27/2 
    16         sig_det = 0.5                   //not sure about this 
    17         dlamb = 0.15 
    18         lambda = 6 
    19          
    20         Duplicate/O root:no_gravity_dat:no_gravity_dat_mat root:no_gravity_dat:John_mat 
    21         Wave data=root:no_gravity_dat:John_mat 
    22          
    23         SUB_SMEAR_2D(DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA,DATA) 
    24          
    25         Duplicate/O root:no_gravity_dat:John_mat root:no_gravity_dat:John_mat_log 
    26         Wave log_data=root:no_gravity_dat:John_mat_log 
    27          
    28         log_data = log(data) 
    29          
    30 end 
    31  
    32 // I have changed the array indexing to [0,], so subtract 1 from the x0,Y0 center  
    33 // to shift from detector coordinates to Igor array index 
    34 // 
    35 // 
    36 // !! the wi values do not match what is written in John's notebook. Are these the  
    37 // correct values for hermite integration?? 
    38 // 
    39 Function SUB_SMEAR_2D(DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA,DATA)               //,Q_MODEL_NAME) 
    40         Variable DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA 
    41         Wave data 
    42          
    43         Variable I,J,KI,KJ              //integers 
    44         Variable SUMm,THET0,Q0,R_PL,R_PD,Q0_PL,Q0_PD,LP,V_R,V_L 
    45         Variable PHI,R0,SIGQ_R,SIGQ_A,Q_PL,Q_PD,DIF_PD_I 
    46         Variable RES_I,RES_J,RES,DIF_PL_J,DIF_PD_J,DIF_PL_I 
    47 //      DIMENSION DATA(128,128),XI(3),WI(3) 
    48 //      EXTERNAL Q_MODEL_NAME 
    49 //      PARAMETER PI = 3.14159265 
    50         Variable N_QUAD = 3 
    51         Make/O/D xi_h = {.6167065887,1.889175877,3.324257432} 
    52         Make/O/D wi_h = {.72462959522,.15706732032,.45300099055E-2} 
    53          
    54 //C     DATA XI/.4360774119,1.3358490740,2.3506049736/ 
    55 //      DATA XI/.6167065887,1.889175877,3.324257432/ 
    56 //      DATA WI/.72462959522,.15706732032,.45300099055E-2/ 
    57 //C     DX :    PIXEL SIZE, CM 
    58 //C     NUM:    NUMBER OF PIXELS ACROSS DETECTOR. (128) 
    59 //C     X0,Y0:  BEAM CENTER, IN UNITS OF PIXELS. 
    60 //C     L1:     SOURCE TO SAMPLE DISTANCE. 
    61 //C     L2:     SAMPLE TO DETECTOR DISTANCE. 
    62 //C     S1:     SOURCE APERTURE RADIUS. 
    63 //C     S2:     SAMPLE APERTURE RADIUS. 
    64 //C     SIG_DET:STANDARD DEVIATION OF DETECTOR SPATIAL RESOLUTION. 
    65 //C     DLAMB:  FWHM WAVLENGTH RESOLUTION. 
    66 //C     LAMBDA: MEAN WAVELENGTH. 
    67 //C     DATA:   OUTPUT SMEARED ARRAY (NUM,NUM) 
    68  
    69         Make/O/D/N=(128,128) sigQR, sigQA 
    70  
    71  
    72         LP = 1 / ( 1/L1 + 1/L2 ) 
    73         V_R = 0.25*(S1/L1)^2 + 0.25*(S2/LP)^2 + (SIG_DET/L2)^2 
    74         V_L = DLAMB^2/6. 
    75         for(i=0;i<num;i+=1) 
    76           R_PL = DX*(I-X0) 
    77           for(j=0;j<num;j+=1) 
    78             R_PD = DX*(J-Y0) 
    79             PHI = ATAN(R_PD/R_PL)               //do I need atan2 here? 
    80             R0 = SQRT(R_PL^2+R_PD^2) 
    81             THET0 = ATAN(R0/L2) 
    82             Q0 = 4*PI*SIN(0.5*THET0)/LAMBDA 
    83 //C     DETERMINE Q VECTOR, CARTESIAN REPRESENTATION. 
    84             Q0_PL = Q0*COS(PHI) 
    85             Q0_PD = Q0*SIN(PHI) 
    86 //C     DETERMINE SIGMA'S FOR RESOLUTION FUNCTION, RADIALLY, AZIMUTHAL 
    87             SIGQ_R = Q0*SQRT(V_R+V_L) 
    88             SIGQ_A = Q0*SQRT(V_R) 
    89              
    90                  sigQR[i][j] = sigq_R 
    91                  sigQA[i][j] = sigq_A 
    92  
    93             SUMm = 0.0 
    94             for(KI=0;ki<N_quad;ki+=1) 
    95               DIF_PL_I = SIGQ_R*COS(PHI)*xi_h[ki] 
    96               DIF_PD_I = SIGQ_R*SIN(PHI)*xi_h[ki] 
    97               for( KJ=0;kj<N_QUAD;kj+=1) 
    98                                 DIF_PL_J = SIGQ_A*SIN(PHI)*xi_h[kj] 
    99                                 DIF_PD_J = SIGQ_A*COS(PHI)*xi_h[kj] 
    100                 //C             -,- 
    101                                 Q_PL = Q0_PL - DIF_PL_I - DIF_PL_J 
    102                                 Q_PD = Q0_PD - DIF_PD_I - DIF_PD_J 
    103                                 SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD)              //,Q_MODEL_NAME)                 
    104                 //C             -,+ 
    105                                 Q_PL = Q0_PL - DIF_PL_I + DIF_PL_J 
    106                                 Q_PD = Q0_PD - DIF_PD_I + DIF_PD_J 
    107                                 SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD)              //,Q_MODEL_NAME)                 
    108                 //C             +,- 
    109                                 Q_PL = Q0_PL + DIF_PL_I - DIF_PL_J 
    110                                 Q_PD = Q0_PD + DIF_PD_I - DIF_PD_J 
    111                                 SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD)              //,Q_MODEL_NAME)                 
    112                 //C             +,+ 
    113                                 Q_PL = Q0_PL + DIF_PL_I + DIF_PL_J 
    114                                 Q_PD = Q0_PD + DIF_PD_I + DIF_PD_J 
    115                                 SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD)              //,Q_MODEL_NAME)                 
    116                 endfor  //   KJ 
    117             endfor  // KI 
    118             DATA[i][j] = SUMm / PI 
    119           endfor  //   J 
    120         endfor  // I 
    121          
    122         RETURN(0) 
    123  
    124 END 
    125  
    126 /// --- either way, same to machine precision 
    127 Function I_MACRO(Q_PL,Q_PD)             //,Q_MODEL_NAME) 
    128         Variable Q_PL,Q_PD 
    129          
    130         Variable I_MACRO,Q,PHI,PHI_MODEL,NU 
    131          
    132         //Q_MODEL_NAME 
    133         //eccentricity factor for ellipse in John's code... 
    134         NU = 1 
    135          
    136 //C     PHI = ATAN(Q_PD/Q_PL) 
    137  
    138         Q = SQRT((NU*Q_PD)^2+Q_PL^2) 
    139          
    140         WAVE cw = $"root:coef_Peak_Gauss" 
    141          
    142         I_MACRO = Peak_Gauss_modelX(cw,Q) 
    143 //      I_MACRO = Q_MODEL_NAME(Q) 
    144          
    145         RETURN(I_MACRO) 
    146 END 
    147  
    148 //Function I_MACRO(Q_PL,Q_PD)           //,Q_MODEL_NAME) 
    149 //      Variable Q_PL,Q_PD 
    150 //       
    151 //      Variable I_MACRO 
    152 //       
    153 //      Make/O/D/N=1 fcnRet,xptW,yPtw 
    154 //      xptw[0] = q_pl 
    155 //      yptw[0] = q_pd 
    156 // 
    157 //      WAVE cw = $"root:coef_sf" 
    158 // 
    159 //      I_MACRO = Sphere_2DX(cw,xptw,yptw) 
    160 //       
    161 //      RETURN(I_MACRO) 
    162 //END 
    163  
    164 ////Structure ResSmear_2D_AAOStruct 
    165 ////    Wave coefW 
    166 ////    Wave zw                 //answer 
    167 ////    Wave qy                 // q-value 
    168 ////    Wave qx 
    169 ////    Wave qz 
    170 ////    Wave sigQx              //resolution 
    171 ////    Wave sigQy 
    172 ////    Wave fs 
    173 ////    String info 
    174 ////EndStructure 
    175 // 
    176 Function Smear_2DModel_5(fcn,s) 
     4 
     5 
     6// there are utlity functions here (commented out) for plotting the 2D resolution function 
     7// at any point on the detector 
     8 
     9 
     10// I can only calculate (nord) points in an AAO fashion, so threading the function calculation 
     11// doesn't help. Even if I can rewrite this to calculate nord*nord AAO, that will typically be 10*10=100 
     12// and that is not enough points to thread with any benefit 
     13 
     14// but the calculation of each q value is independent, so I can split the total number of q-values between processors 
     15// 
     16// -- but I must be careful to pass this a function that is not already threaded! 
     17// --- BUT - I can't pass either a function reference, OR a structure to a thread! 
     18// ---- So sadly, each 2D resolution calculation must be threaded by hand. 
     19// 
     20// See the Sphere_2D function for an example of how to thread the smearing 
     21// 
     22// 
     23// SRK May 2010 
     24 
     25// 
     26// NOTE: there is a definition of FindTheta = -1* FindPhi() that is a duplicate of what is in RawWindowHook.ipf 
     27// and should eventually be in a common location for analysis and reduction packages 
     28// 
     29 
     30 
     31//// this is the completely generic 2D smearing, not threaded, but takes FUNCREF and STRUCT parameters 
     32// 
     33// uses resolution ellipse defined perpendicular "Y" and parallel "X", 
     34// rotating the ellipse into its proper orientaiton based on qx,qy 
     35// 
     36// 5 gauss points is not enough - it gives artifacts as a function of phi 
     37// 10 gauss points is minimally sufficient 
     38// 20 gauss points are needed if lots of oscillations (just like in 1D) 
     39// even more may be necessary for highly peaked functions 
     40// 
     41// 
     42Function Smear_2DModel_PP(fcn,s,nord) 
    17743        FUNCREF SANS_2D_ModelAAO_proto fcn 
    17844        Struct ResSmear_2D_AAOStruct &s 
    179          
    180         String weightStr="gauss5wt",zStr="gauss5z" 
    181         Variable nord=5 
    182  
    183 //      if wt,z waves don't exist, create them (only check for weight, should really check for both) 
    184         if (WaveExists($weightStr) == 0) // wave reference is not valid,  
    185                 Make/D/N=(nord) $weightStr,$zStr 
    186                 Wave wt = $weightStr 
    187                 Wave xi = $zStr         // wave references to pass 
    188                 Make5GaussPoints(wt,xi)  
    189         else 
    190                 if(exists(weightStr) > 1)  
    191                          Abort "wave name is already in use"            //executed only if name is in use elsewhere 
    192                 endif 
    193                 Wave wt = $weightStr 
    194                 Wave xi = $zStr         // create the wave references 
    195         endif 
    196          
    197         Variable ii,jj,kk,ax,bx,ay,by,num 
    198         Variable qx,qy,qz,qval,sqx,sqy,fs 
    199         Variable qy_pt,qx_pt,res_x,res_y,res_tot,answer,sumIn,sumOut 
    200         num=numpnts(s.qx) 
    201          
    202         // end points of integration 
    203         // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero 
    204         // +/- 3 sigq catches 99.73% of distrubution 
    205         // change limits (and spacing of zi) at each evaluation based on R() 
    206         //integration from va to vb 
    207         Make/O/D/N=1 fcnRet,xptW,yPtw 
     45        Variable nord 
     46         
     47        String weightStr,zStr 
     48 
     49        Variable ii,jj,kk,num 
     50        Variable qx,qy,qz,qval,fs 
     51        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut 
     52         
     53        Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3 
     54         
     55        switch(nord)     
     56                case 5:          
     57                        weightStr="gauss5wt" 
     58                        zStr="gauss5z" 
     59                        if (WaveExists($weightStr) == 0) 
     60                                Make/O/D/N=(nord) $weightStr,$zStr 
     61                                Make5GaussPoints($weightStr,$zStr)       
     62                        endif 
     63                        break                            
     64                case 10:                 
     65                        weightStr="gauss10wt" 
     66                        zStr="gauss10z" 
     67                        if (WaveExists($weightStr) == 0) 
     68                                Make/O/D/N=(nord) $weightStr,$zStr 
     69                                Make10GaussPoints($weightStr,$zStr)      
     70                        endif 
     71                        break                            
     72                case 20:                 
     73                        weightStr="gauss20wt" 
     74                        zStr="gauss20z" 
     75                        if (WaveExists($weightStr) == 0) 
     76                                Make/O/D/N=(nord) $weightStr,$zStr 
     77                                Make20GaussPoints($weightStr,$zStr)      
     78                        endif 
     79                        break 
     80                default:                                                         
     81                        Abort "Smear_2DModel_PP called with invalid nord value"                                  
     82        endswitch 
     83         
     84        Wave/Z wt = $weightStr 
     85        Wave/Z xi = $zStr 
     86         
     87/// keep these waves local 
     88//      Make/O/D/N=1 yPtW 
     89        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW 
     90         
     91        // now just loop over the points as specified 
     92        num=numpnts(s.xw[0]) 
    20893         
    20994        answer=0 
     95         
     96        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     97        Variable t1=StopMSTimer(-2) 
     98 
    21099        //loop over q-values 
    211100        for(ii=0;ii<num;ii+=1) 
    212                 qx = s.qx[ii] 
    213                 qy = s.qy[ii] 
     101         
     102//              if(mod(ii, 1000 ) == 0) 
     103//                      Print "ii= ",ii 
     104//              endif 
     105                 
     106                qx = s.xw[0][ii] 
     107                qy = s.xw[1][ii] 
    214108                qz = s.qz[ii] 
    215109                qval = sqrt(qx^2+qy^2+qz^2) 
    216                 sqx = s.sigQx[ii] 
    217                 sqy = s.sigQy[ii] 
     110                spl = s.sQpl[ii] 
     111                spp = s.sQpp[ii] 
    218112                fs = s.fs[ii] 
    219113                 
    220                 ax = -3*sqx + qx                //qx integration limits 
    221                 bx = 3*sqx + qx 
    222                 ay = -3*sqy + qy                //qy integration limits 
    223                 by = 3*sqy + qy 
    224                  
    225                 // 5-pt quadrature loops 
     114                normFactor = 2*pi*spl*spp 
     115                 
     116                phi = -1*FindTheta(qx,qy)               //Findtheta is an exact duplicate of FindPhi() * -1 
     117                 
     118                apl = -numStdDev*spl + qval             //parallel = q integration limits 
     119                bpl = numStdDev*spl + qval 
     120                app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
     121                bpp = numStdDev*spp + phi 
     122                 
     123                //make sure the limits are reasonable. 
     124                if(apl < 0) 
     125                        apl = 0 
     126                endif 
     127                // do I need to specially handle limits when phi ~ 0? 
     128         
     129                 
    226130                sumOut = 0 
    227                 for(jj=0;jj<nord;jj+=1)         // call qy the "outer' 
    228                         qy_pt = (xi[jj]*(by-ay)+ay+by)/2 
    229                         res_y = exp((-1*(qy - qy_pt)^2)/(2*sqy*sqy)) 
     131                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
     132                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
    230133 
    231134                        sumIn=0 
    232                         for(kk=0;kk<nord;kk+=1) 
    233  
    234                                 qx_pt = (xi[kk]*(bx-ax)+ax+bx)/2 
    235                                 res_x = exp((-1*(qx - qx_pt)^2)/(2*sqx*sqx)) 
     135                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
     136 
     137                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
    236138                                 
    237                                 res_tot = res_x*res_y/(2*pi*sqx*sqy) 
    238                                 xptw[0] = qx_pt 
    239                                 yptw[0] = qy_pt 
    240                                 fcn(s.coefW,fcnRet,xptw,yptw)                   //fcn passed in is an AAO 
    241                                 sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
     139                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     140                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
     141                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
     142                                 
     143                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     144                                res_tot[kk] /= normFactor 
     145//                              res_tot[kk] *= fs 
     146 
    242147                        endfor 
    243                         answer += (bx-ax)/2.0*sumIn             //this is NOT the right normalization 
     148                         
     149                        fcn(s.coefW,fcnRet,xptw,yptw)                   //calculate nord pts at a time 
     150                         
     151                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
     152                        fcnRet *= wt[jj]*wt*res_tot 
     153                        // 
     154                        answer += (bpl-apl)/2.0*sum(fcnRet)             // 
    244155                endfor 
    245156 
    246                 answer *= (by-ay)/2.0 
     157                answer *= (bpp-app)/2.0 
    247158                s.zw[ii] = answer 
    248 //              s.zw[ii] = sumIn 
    249159        endfor 
    250160         
     161        Variable elap = (StopMSTimer(-2) - t1)/1e6 
     162        Print "elapsed time = ",elap 
    251163         
    252164        return(0) 
    253165end 
    254166 
    255 Function Smear_2DModel_20(fcn,s) 
    256         FUNCREF SANS_2D_ModelAAO_proto fcn 
    257         Struct ResSmear_2D_AAOStruct &s 
    258          
    259         String weightStr="gauss20wt",zStr="gauss20z" 
    260         Variable nord=20 
    261  
    262 //      if wt,z waves don't exist, create them (only check for weight, should really check for both) 
    263         if (WaveExists($weightStr) == 0) // wave reference is not valid,  
    264                 Make/D/N=(nord) $weightStr,$zStr 
    265                 Wave wt = $weightStr 
    266                 Wave xi = $zStr         // wave references to pass 
    267                 Make20GaussPoints(wt,xi)         
    268         else 
    269                 if(exists(weightStr) > 1)  
    270                          Abort "wave name is already in use"            //executed only if name is in use elsewhere 
    271                 endif 
    272                 Wave wt = $weightStr 
    273                 Wave xi = $zStr         // create the wave references 
    274         endif 
    275          
    276         Variable ii,jj,kk,ax,bx,ay,by,num 
    277         Variable qx,qy,qz,qval,sqx,sqy,fs 
    278         Variable qy_pt,qx_pt,res_x,res_y,res_tot,answer,sumIn,sumOut 
    279         num=numpnts(s.qx) 
    280          
    281         // end points of integration 
    282         // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero 
    283         // +/- 3 sigq catches 99.73% of distrubution 
    284         // change limits (and spacing of zi) at each evaluation based on R() 
    285         //integration from va to vb 
    286         Make/O/D/N=1 fcnRet,xptW,yPtw 
    287          
    288         answer=0 
    289         //loop over q-values 
    290         for(ii=0;ii<num;ii+=1) 
    291                 qx = s.qx[ii] 
    292                 qy = s.qy[ii] 
    293                 qz = s.qz[ii] 
    294                 qval = sqrt(qx^2+qy^2+qz^2) 
    295                 sqx = s.sigQx[ii] 
    296                 sqy = s.sigQy[ii] 
    297                 fs = s.fs[ii] 
    298                  
    299                 ax = -3*sqx + qx                //qx integration limits 
    300                 bx = 3*sqx + qx 
    301                 ay = -3*sqy + qy                //qy integration limits 
    302                 by = 3*sqy + qy 
    303                  
    304                 // 20-pt quadrature loops 
    305                 sumOut = 0 
    306                 for(jj=0;jj<nord;jj+=1)         // call qy the "outer' 
    307                         qy_pt = (xi[jj]*(by-ay)+ay+by)/2 
    308                         res_y = exp((-1*(qy - qy_pt)^2)/(2*sqy*sqy)) 
    309  
    310                         sumIn=0 
    311                         for(kk=0;kk<nord;kk+=1) 
    312  
    313                                 qx_pt = (xi[kk]*(bx-ax)+ax+bx)/2 
    314                                 res_x = exp((-1*(qx - qx_pt)^2)/(2*sqx*sqx)) 
    315                                  
    316                                 res_tot = res_x*res_y/(2*pi*sqx*sqy) 
    317                                 xptw[0] = qx_pt 
    318                                 yptw[0] = qy_pt 
    319                                 fcn(s.coefW,fcnRet,xptw,yptw)                   //fcn passed in is an AAO 
    320                                 sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
    321                         endfor 
    322                         answer += (bx-ax)/2.0*sumIn             //this is NOT the right normalization 
    323                 endfor 
    324  
    325                 answer *= (by-ay)/2.0 
    326                 s.zw[ii] = answer 
    327 //              s.zw[ii] = sumIn 
    328         endfor 
    329          
    330          
    331         return(0) 
     167 
     168 
     169 
     170 
     171//phi is defined from +x axis, proceeding CCW around [0,2Pi] 
     172//rotate the resolution function by theta,  = -phi 
     173// 
     174// this is only different by (-1) from FindPhi 
     175// I'd just call FindPhi, but it's awkward to include 
     176// 
     177Threadsafe Function FindTheta(vx,vy) 
     178        variable vx,vy 
     179 
     180         
     181        return(-1 * FindPhi(vx,vy)) 
    332182end 
     183 
     184 
     185 
     186         
     187////////////////////////////////////////////////////////// 
     188////////////////////////////////////////////////////////// 
     189/// utility functions to test the calculation of the resolution function and  
     190//       plotting of it for visualization 
     191// 
     192// -- these need to have the reduction package included for it to compile 
     193// 
     194// this works with the raw 2D data, calculating resolution in place 
     195// type is "RAW" 
     196// xx,yy are in detector coordinates 
     197// 
     198// call as:  
     199// PlotResolution_atPixel("RAW",pcsr(A),qcsr(A)) 
     200// 
     201// then: 
     202// Display;AppendImage res;AppendMatrixContour res;ModifyContour res labels=0,autoLevels={*,*,3} 
     203//  
     204 
     205//Function PlotResolution_atPixel(type,xx,yy) 
     206//      String type 
     207//      Variable xx,yy 
     208//       
     209/////from QxQyExport 
     210//      String destStr="",typeStr="" 
     211//      Variable step=1,refnum 
     212//      destStr = "root:Packages:NIST:"+type 
     213//       
     214//      //must select the linear_data to export 
     215//      NVAR isLog = $(destStr+":gIsLogScale") 
     216//      if(isLog==1) 
     217//              typeStr = ":linear_data" 
     218//      else 
     219//              typeStr = ":data" 
     220//      endif 
     221//       
     222//      NVAR pixelsX = root:myGlobals:gNPixelsX 
     223//      NVAR pixelsY = root:myGlobals:gNPixelsY 
     224//       
     225//      Wave data=$(destStr+typeStr) 
     226//      WAVE intw=$(destStr + ":integersRead") 
     227//      WAVE rw=$(destStr + ":realsRead") 
     228//      WAVE/T textw=$(destStr + ":textRead") 
     229//       
     230////    Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist 
     231//      Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist 
     232// 
     233//      Variable xctr,yctr,sdd,lambda,pixSize 
     234//      xctr = rw[16] 
     235//      yctr = rw[17] 
     236//      sdd = rw[18] 
     237//      lambda = rw[26] 
     238//      pixSize = rw[13]/10             //convert mm to cm (x and y are the same size pixels) 
     239//       
     240////    qx_val = CalcQx(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)          //+1 converts to detector coordinate system 
     241////    qy_val = CalcQy(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     242//       
     243//      qx_val = CalcQx(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)                //+1 converts to detector coordinate system 
     244//      qy_val = CalcQy(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     245//       
     246////    Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val 
     247// 
     248/////************ 
     249//// do everything to write out the resolution too 
     250//      // un-comment these if you want to write out qz_val and qval too, then use the proper save command 
     251//      qval = CalcQval(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     252//      qz_val = CalcQz(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     253//      phi = FindPhi( pixSize*((xx+1)-xctr) , pixSize*((yy+1)-yctr))           //(dx,dy) 
     254//      r_dist = sqrt(  (pixSize*((xx+1)-xctr))^2 +  (pixSize*((yy+1)-yctr))^2 )                //radial distance from ctr to pt 
     255////    Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist 
     256//      //everything in 1D now 
     257////    Duplicate/O qval SigmaQX,SigmaQY,fsubS 
     258//      Variable SigmaQX,SigmaQY,fsubS 
     259// 
     260//      Variable L2 = rw[18] 
     261//      Variable BS = rw[21] 
     262//      Variable S1 = rw[23] 
     263//      Variable S2 = rw[24] 
     264//      Variable L1 = rw[25] 
     265//      Variable lambdaWidth = rw[27]    
     266//      Variable usingLenses = rw[28]           //new 2007 
     267// 
     268//      //Two parameters DDET and APOFF are instrument dependent.  Determine 
     269//      //these from the instrument name in the header. 
     270//      //From conversation with JB on 01.06.99 these are the current good values 
     271//      Variable DDet 
     272//      NVAR apOff = root:myGlobals:apOff               //in cm 
     273//      DDet = rw[10]/10                        // header value (X) is in mm, want cm here 
     274// 
     275//      Variable ret1,ret2,ret3,del_r 
     276//      del_r = rw[10] 
     277//       
     278//      get2DResolution(qval,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,ret1,ret2,ret3) 
     279//      SigmaQX = ret1   
     280//      SigmaQY = ret2   
     281//      fsubs = ret3     
     282/////// 
     283// 
     284////    Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0 
     285//      Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor 
     286//      Variable qx_ret,qy_ret 
     287//       
     288////    theta = FindPhi(qx_val,qy_val) 
     289//// need to rotate properly - theta is defined a starting from +y axis, moving CW 
     290//// we define phi starting from +x and moving CCW 
     291//      theta = -phi                    //seems to give the right behavior... 
     292//       
     293//       
     294//      Print qx_val,qy_val,qval 
     295//      Print "phi, theta",phi,theta 
     296//       
     297//      FindQxQy(qval,phi,qx_ret,qy_ret) 
     298//       
     299//      sx = SigmaQx 
     300//      sy = sigmaQy 
     301//      x0 = qx_val 
     302//      y0 = qy_val 
     303//       
     304//      a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) 
     305//      b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) 
     306//      c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) 
     307//       
     308//      normFactor = pi/sqrt(a*c-b*b) 
     309// 
     310//      Make/O/D/N=(num,num) res 
     311//      // so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle 
     312//      maxSig = max(sx,sy) 
     313//      Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res 
     314//      Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res 
     315////    Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res 
     316////    Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res 
     317//       
     318//      Variable xPt,yPt,delx,dely,offx,offy 
     319//      delx = DimDelta(res,0) 
     320//      dely = DimDelta(res,1) 
     321//      offx = DimOffset(res,0) 
     322//      offy = DimOffset(res,1) 
     323// 
     324//      Print "sx,sy = ",sx,sy 
     325//      for(ii=0;ii<num;ii+=1) 
     326//              xPt = offx + ii*delx 
     327//              for(jj=0;jj<num;jj+=1) 
     328//                      yPt = offy + jj*dely 
     329//                      res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2)) 
     330//              endfor 
     331//      endfor   
     332//      res /= normFactor 
     333//       
     334//      //Print sum(res,-inf,inf)*delx*dely 
     335//      if(WaveExists($"coef")==0) 
     336//              Make/O/D/N=6 coef 
     337//      endif 
     338//      Wave coef=coef 
     339//      coef[0] = 1 
     340//      coef[1] = qx_val 
     341//      coef[2] = qy_val 
     342//      coef[3] = sx 
     343//      coef[4] = sy 
     344//      coef[5] = theta 
     345// 
     346////    Variable t1=StopMSTimer(-2) 
     347// 
     348//// 
     349////    do2dIntegrationGauss(-nStdDev*maxSig+x0,nStdDev*maxSig+x0,-nStdDev*maxSig+y0,nStdDev*maxSig+y0) 
     350//// 
     351// 
     352////    Variable elap = (StopMSTimer(-2) - t1)/1e6 
     353////    Print "elapsed time = ",elap 
     354////    Print "time for 16384 = (minutes)",16384*elap/60 
     355//      return(0) 
     356//End 
     357 
     358 
     359//// this is called each time to integrate the gaussian 
     360//Function do2dIntegrationGauss(xMin,xMax,yMin,yMax) 
     361//      Variable xMin,xMax,yMin,yMax 
     362//       
     363//      Variable/G globalXmin=xMin 
     364//      Variable/G globalXmax=xMax 
     365//      Variable/G globalY 
     366//                       
     367//      Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5)         
     368//      KillVariables/z globalXmax,globalXmin,globalY 
     369//      print "integration of 2D = ",result 
     370//End 
     371// 
     372//Function Gauss2DFuncOuter(inY) 
     373//      Variable inY 
     374//       
     375//      NVAR globalXmin,globalXmax,globalY 
     376//      globalY=inY 
     377//       
     378//      return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5)           
     379//End 
     380// 
     381//Function Gauss2DFuncInner(inX) 
     382//      Variable inX 
     383//       
     384//      NVAR globalY 
     385//      Wave coef=coef 
     386//       
     387//      return Gauss2D_theta(coef,inX,GlobalY) 
     388//End 
     389// 
     390//Function Gauss2D_theta(w,x,y) 
     391//      Wave w 
     392//      Variable x,y 
     393//       
     394//      Variable val,a,b,c 
     395//      Variable scale,x0,y0,sx,sy,theta,normFactor 
     396//       
     397//      scale = w[0] 
     398//      x0 = w[1] 
     399//      y0 = w[2] 
     400//      sx = w[3] 
     401//      sy = w[4] 
     402//      theta = w[5] 
     403//       
     404//      a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) 
     405//      b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) 
     406//      c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) 
     407//       
     408//      val = exp(-1*(a*(x-x0)^2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)^2)) 
     409//       
     410//      normFactor = pi/sqrt(a*c-b*b) 
     411//       
     412//      return(scale*val/normFactor) 
     413//end 
     414// 
     415//////////////////////////////////////////////////////////// 
     416//////////////////////////////////////////////////////////// 
Note: See TracChangeset for help on using the changeset viewer.