Changeset 708


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
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/PeakGauss2D_v40.ipf

    r570 r708  
    157157 
    158158//non-threaded version of the function 
    159 Function PeakGauss2D_noThread(cw,zw,xw,yw) 
     159ThreadSafe Function PeakGauss2D_noThread(cw,zw,xw,yw) 
    160160        WAVE cw,zw, xw,yw 
    161161         
    162162#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed 
    163          
    164163        zw= PeakGauss2DX(cw,xw,yw) 
    165164#else 
     
    175174        Struct ResSmear_2D_AAOStruct &s 
    176175         
    177 //      Smear_2DModel_20(PeakGauss2D_noThread,s) 
    178         Smear_2DModel_5(PeakGauss2D_noThread,s) 
     176        //no threads 
     177//      Smear_2DModel_PP(PeakGauss2D_noThread,s,10) 
     178 
     179        //or threaded... 
     180        SmearPeakGauss2D_THR(s,10) 
     181         
    179182        return(0) 
    180183end 
     
    190193        WAVE qy = $(DF+str+"_qy") 
    191194        WAVE qz = $(DF+str+"_qz") 
    192         WAVE sigQx = $(DF+str+"_sigQx") 
    193         WAVE sigQy = $(DF+str+"_sigQy") 
     195        WAVE sQpl = $(DF+str+"_sQpl") 
     196        WAVE sQpp = $(DF+str+"_sQpp") 
    194197        WAVE shad = $(DF+str+"_fs") 
    195198         
     
    197200        WAVE s.coefW = coefW     
    198201        WAVE s.zw = resultW      
    199         WAVE s.qx = qx 
    200         WAVE s.qy = qy 
     202        WAVE s.xw[0] = qx 
     203        WAVE s.xw[1] = qy 
    201204        WAVE s.qz = qz 
    202         WAVE s.sigQx = sigQx 
    203         WAVE s.sigQy = sigQy 
     205        WAVE s.sQpl = sQpl 
     206        WAVE s.sQpp = sQpp 
    204207        WAVE s.fs = shad 
    205208         
     
    220223//      retval = fPeak_Gauss_model(w,qval) 
    221224         
    222         RETURN(retVal) 
    223 END 
     225        return(retVal) 
     226End 
     227 
     228/////////////// 
     229// 
     230// this is the threaded version, that dispatches the calculation out to threads 
     231// 
     232// must be written specific to each 2D function 
     233//  
     234Function SmearPeakGauss2D_THR(s,nord) 
     235        Struct ResSmear_2D_AAOStruct &s 
     236        Variable nord 
     237         
     238        String weightStr,zStr 
     239         
     240// create all of the necessary quadrature waves here - rather than inside a threadsafe function 
     241        switch(nord)     
     242                case 5:          
     243                        weightStr="gauss5wt" 
     244                        zStr="gauss5z" 
     245                        if (WaveExists($weightStr) == 0) 
     246                                Make/O/D/N=(nord) $weightStr,$zStr 
     247                                Make5GaussPoints($weightStr,$zStr)       
     248                        endif 
     249                        break                            
     250                case 10:                 
     251                        weightStr="gauss10wt" 
     252                        zStr="gauss10z" 
     253                        if (WaveExists($weightStr) == 0) 
     254                                Make/O/D/N=(nord) $weightStr,$zStr 
     255                                Make10GaussPoints($weightStr,$zStr)      
     256                        endif 
     257                        break                            
     258                case 20:                 
     259                        weightStr="gauss20wt" 
     260                        zStr="gauss20z" 
     261                        if (WaveExists($weightStr) == 0) 
     262                                Make/O/D/N=(nord) $weightStr,$zStr 
     263                                Make20GaussPoints($weightStr,$zStr)      
     264                        endif 
     265                        break 
     266                default:                                                         
     267                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                         
     268        endswitch 
     269         
     270        Wave/Z wt = $weightStr 
     271        Wave/Z xi = $zStr               // wave references to pass 
     272 
     273        Variable npt=numpnts(s.xw[0]) 
     274        Variable i,nthreads= ThreadProcessorCount 
     275        variable mt= ThreadGroupCreate(nthreads) 
     276 
     277        Variable t1=StopMSTimer(-2) 
     278         
     279        for(i=0;i<nthreads;i+=1) 
     280        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     281                ThreadStart mt,i,SmearPeakGauss2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     282        endfor 
     283 
     284        do 
     285                variable tgs= ThreadGroupWait(mt,100) 
     286        while( tgs != 0 ) 
     287 
     288        variable dummy= ThreadGroupRelease(mt) 
     289         
     290// comment out the threading + uncomment this for testing to make sure that the single thread works 
     291//      nThreads=1 
     292//      SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     293 
     294        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
     295         
     296        return(0) 
     297end 
     298 
     299// 
     300// - worker function for threads of PeakGauss2D 
     301// 
     302ThreadSafe Function SmearPeakGauss2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord) 
     303        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi 
     304        Variable pt1,pt2,nord 
     305         
     306// now passed in.... 
     307//      Wave wt = $weightStr 
     308//      Wave xi = $zStr          
     309 
     310        Variable ii,jj,kk,num 
     311        Variable qx,qy,qz,qval,sx,sy,fs 
     312        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut 
     313         
     314        Variable normFactor,phi,theta,maxSig,numStdDev=3 
     315         
     316/// keep these waves local 
     317        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW 
     318         
     319        // now just loop over the points as specified 
     320         
     321        answer=0 
     322         
     323        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     324 
     325        //loop over q-values 
     326        for(ii=pt1;ii<(pt2+1);ii+=1) 
     327                 
     328                qx = qxw[ii] 
     329                qy = qyw[ii] 
     330                qz = qzw[ii] 
     331                qval = sqrt(qx^2+qy^2+qz^2) 
     332                spl = sxw[ii] 
     333                spp = syw[ii] 
     334                fs = fsw[ii] 
     335                 
     336                normFactor = 2*pi*spl*spp 
     337                 
     338                phi = FindPhi(qx,qy) 
     339                 
     340                apl = -numStdDev*spl + qval             //parallel = q integration limits 
     341                bpl = numStdDev*spl + qval 
     342                app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
     343                bpp = numStdDev*spp + phi 
     344                 
     345                //make sure the limits are reasonable. 
     346                if(apl < 0) 
     347                        apl = 0 
     348                endif 
     349                // do I need to specially handle limits when phi ~ 0? 
     350         
     351                 
     352                sumOut = 0 
     353                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
     354                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     355 
     356                        sumIn=0 
     357                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
     358 
     359                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
     360                                 
     361                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     362                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
     363                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
     364                                 
     365                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     366                                res_tot[kk] /= normFactor 
     367//                              res_tot[kk] *= fs 
     368 
     369                        endfor 
     370                         
     371                        PeakGauss2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO 
     372                         
     373                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
     374                        fcnRet *= wt[jj]*wt*res_tot 
     375                        // 
     376                        answer += (bpl-apl)/2.0*sum(fcnRet)             // 
     377                endfor 
     378 
     379                answer *= (bpp-app)/2.0 
     380                zw[ii] = answer 
     381        endfor 
     382         
     383        return(0) 
     384end 
     385 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Sphere_2D_v40.ipf

    r570 r708  
    108108End 
    109109 
    110  
    111 //threaded version of the function 
    112 ThreadSafe Function Sphere2D_T(cw,zw,xw,yw,p1,p2) 
    113         WAVE cw,zw,xw,yw 
    114         Variable p1,p2 
    115          
    116 #if exists("Sphere_2DX")                        //to hide the function if XOP not installed 
    117          
    118         zw[p1,p2]= Sphere_2DX(cw,xw,yw) 
    119  
    120 #endif 
    121  
    122         return 0 
    123 End 
    124  
    125110// 
    126111//  Fit function that is actually a wrapper to dispatch the calculation to N threads 
     
    140125 
    141126#endif 
     127 
     128        return(0) 
     129End 
     130 
     131 
     132/// 
     133//// keep this section as an example 
     134// 
    142135 
    143136//      Variable npt=numpnts(yw) 
     
    160153////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
    161154//       
    162         return(0) 
    163 End 
    164  
    165  
    166 //non-threaded version of the function 
    167 Function Sphere2D_noThread(cw,zw,xw,yw) 
     155// 
     156////// end example of threading 
     157 
     158 
     159 
     160//threaded version of the function 
     161//ThreadSafe Function Sphere2D_T(cw,zw,xw,yw,p1,p2) 
     162//      WAVE cw,zw,xw,yw 
     163//      Variable p1,p2 
     164//       
     165//#if exists("Sphere_2DX")                      //to hide the function if XOP not installed 
     166//       
     167//      zw[p1,p2]= Sphere_2DX(cw,xw,yw) 
     168// 
     169//#endif 
     170// 
     171//      return 0 
     172//End 
     173 
     174 
     175//non-threaded version of the function, necessary for the smearing calculation 
     176// -- the smearing calculation can only calculate (nord) points at a time. 
     177// 
     178ThreadSafe Function Sphere2D_noThread(cw,zw,xw,yw) 
    168179        WAVE cw,zw, xw,yw 
    169180         
    170181#if exists("Sphere_2DX")                        //to hide the function if XOP not installed 
    171          
    172182        zw= Sphere_2DX(cw,xw,yw) 
    173  
    174183#endif 
    175184 
     
    177186End 
    178187 
    179 // I think I did this because when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading 
    180 // would just be a slowdown 
     188 
     189//  when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading 
     190// is just a slowdown (by more than a factor of 2 !!) 
     191// 
     192// - the threading gives a clean speedup of 2 for N=2, even for this simple calculation 
     193// 
     194// nord = 5,10,20 allowed 
     195// 
    181196Function SmearedSphere2D(s) 
    182197        Struct ResSmear_2D_AAOStruct &s 
    183198         
    184         Smear_2DModel_5(Sphere2D_noThread,s) 
    185 //      Smear_2DModel_20(Sphere2D_noThread,s) 
     199//// non-threaded, but generic calculation 
     200//// the last param is nord      
     201//      Smear_2DModel_PP(Sphere2D_noThread,s,10) 
     202 
     203//// the threaded version must be specifically written, since 
     204//// FUNCREF can't be passed into  a threaded calc (or structures) 
     205//// the last param is nord 
     206        SmearSphere2D_THR(s,10)          
     207 
    186208        return(0) 
    187209end 
     
    197219        WAVE qy = $(DF+str+"_qy") 
    198220        WAVE qz = $(DF+str+"_qz") 
    199         WAVE sigQx = $(DF+str+"_sigQx") 
    200         WAVE sigQy = $(DF+str+"_sigQy") 
     221        WAVE sQpl = $(DF+str+"_sQpl") 
     222        WAVE sQpp = $(DF+str+"_sQpp") 
    201223        WAVE shad = $(DF+str+"_fs") 
    202224         
     
    204226        WAVE s.coefW = coefW     
    205227        WAVE s.zw = resultW      
    206         WAVE s.qx = qx 
    207         WAVE s.qy = qy 
     228        WAVE s.xw[0] = qx 
     229        WAVE s.xw[1] = qy 
    208230        WAVE s.qz = qz 
    209         WAVE s.sigQx = sigQx 
    210         WAVE s.sigQy = sigQy 
     231        WAVE s.sQpl = sQpl 
     232        WAVE s.sQpp = sQpp 
    211233        WAVE s.fs = shad 
    212234         
     
    216238        return (0) 
    217239End 
     240 
     241 
     242 
     243 
     244// 
     245// this is the threaded version, that dispatches the calculation out to threads 
     246// 
     247// must be written specific to each 2D function 
     248//  
     249Function SmearSphere2D_THR(s,nord) 
     250        Struct ResSmear_2D_AAOStruct &s 
     251        Variable nord 
     252         
     253        String weightStr,zStr 
     254         
     255// create all of the necessary quadrature waves here - rather than inside a threadsafe function 
     256        switch(nord)     
     257                case 5:          
     258                        weightStr="gauss5wt" 
     259                        zStr="gauss5z" 
     260                        if (WaveExists($weightStr) == 0) 
     261                                Make/O/D/N=(nord) $weightStr,$zStr 
     262                                Make5GaussPoints($weightStr,$zStr)       
     263                        endif 
     264                        break                            
     265                case 10:                 
     266                        weightStr="gauss10wt" 
     267                        zStr="gauss10z" 
     268                        if (WaveExists($weightStr) == 0) 
     269                                Make/O/D/N=(nord) $weightStr,$zStr 
     270                                Make10GaussPoints($weightStr,$zStr)      
     271                        endif 
     272                        break                            
     273                case 20:                 
     274                        weightStr="gauss20wt" 
     275                        zStr="gauss20z" 
     276                        if (WaveExists($weightStr) == 0) 
     277                                Make/O/D/N=(nord) $weightStr,$zStr 
     278                                Make20GaussPoints($weightStr,$zStr)      
     279                        endif 
     280                        break 
     281                default:                                                         
     282                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                         
     283        endswitch 
     284         
     285        Wave/Z wt = $weightStr 
     286        Wave/Z xi = $zStr               // wave references to pass 
     287 
     288        Variable npt=numpnts(s.xw[0]) 
     289        Variable i,nthreads= ThreadProcessorCount 
     290        variable mt= ThreadGroupCreate(nthreads) 
     291 
     292        Variable t1=StopMSTimer(-2) 
     293         
     294        for(i=0;i<nthreads;i+=1) 
     295//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     296                ThreadStart mt,i,SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     297        endfor 
     298 
     299        do 
     300                variable tgs= ThreadGroupWait(mt,100) 
     301        while( tgs != 0 ) 
     302 
     303        variable dummy= ThreadGroupRelease(mt) 
     304         
     305// comment out the threading + uncomment this for testing to make sure that the single thread works 
     306//      nThreads=1 
     307//      SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     308 
     309        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
     310         
     311        return(0) 
     312end 
     313 
     314// 
     315// - worker function for threads of Sphere2D 
     316// 
     317ThreadSafe Function SmearSphere2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord) 
     318        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi 
     319        Variable pt1,pt2,nord 
     320         
     321// now passed in.... 
     322//      Wave wt = $weightStr 
     323//      Wave xi = $zStr          
     324 
     325        Variable ii,jj,kk,num 
     326        Variable qx,qy,qz,qval,sx,sy,fs 
     327        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut 
     328         
     329        Variable normFactor,phi,theta,maxSig,numStdDev=3 
     330         
     331/// keep these waves local 
     332        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW 
     333         
     334        // now just loop over the points as specified 
     335         
     336        answer=0 
     337         
     338        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     339 
     340        //loop over q-values 
     341        for(ii=pt1;ii<(pt2+1);ii+=1) 
     342                 
     343                qx = qxw[ii] 
     344                qy = qyw[ii] 
     345                qz = qzw[ii] 
     346                qval = sqrt(qx^2+qy^2+qz^2) 
     347                spl = sxw[ii] 
     348                spp = syw[ii] 
     349                fs = fsw[ii] 
     350                 
     351                normFactor = 2*pi*spl*spp 
     352                 
     353                phi = FindPhi(qx,qy) 
     354                 
     355                apl = -numStdDev*spl + qval             //parallel = q integration limits 
     356                bpl = numStdDev*spl + qval 
     357                app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
     358                bpp = numStdDev*spp + phi 
     359                 
     360                //make sure the limits are reasonable. 
     361                if(apl < 0) 
     362                        apl = 0 
     363                endif 
     364                // do I need to specially handle limits when phi ~ 0? 
     365         
     366                 
     367                sumOut = 0 
     368                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
     369                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     370 
     371                        sumIn=0 
     372                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
     373 
     374                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
     375                                 
     376                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     377                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
     378                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
     379                                 
     380                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     381                                res_tot[kk] /= normFactor 
     382//                              res_tot[kk] *= fs 
     383 
     384                        endfor 
     385                         
     386                        Sphere2D_noThread(coef,fcnRet,xptw,yptw)                        //fcn passed in is an AAO 
     387                         
     388                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
     389                        fcnRet *= wt[jj]*wt*res_tot 
     390                        // 
     391                        answer += (bpl-apl)/2.0*sum(fcnRet)             // 
     392                endfor 
     393 
     394                answer *= (bpp-app)/2.0 
     395                zw[ii] = answer 
     396        endfor 
     397         
     398        return(0) 
     399end 
     400 
     401 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Packages/Wrapper_v40.ipf

    r704 r708  
    10441044         
    10451045        if(yesReport) 
    1046                 String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //this is *hopefully* one wave 
     1046                String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              // this is *hopefully* one wave 
    10471047                String topGraph= WinName(0,1)   //this is the topmost graph 
    10481048         
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/SA_includes_v410.ipf

    r693 r708  
    2828 
    2929#include "DataSetHandling"                                      //added Nov 2009 AJJ 
     30#include "Smear_2D"                                             //for 2D resolution smearing, May 2010 
    3031 
    3132Menu "SANS Models" 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/WriteModelData_v40.ipf

    r607 r708  
    193193                PathInfo/S catPathName 
    194194                fullPath = DoSaveFileDialog("Save data as",fname=baseStr+".txt") 
     195                Print fullPath 
    195196                If(cmpstr(fullPath,"")==0) 
    196197                        //user cancel, don't write out a file 
     
    200201                //Print "dialog fullpath = ",fullpath 
    201202        Endif 
    202          
     203 
    203204        Open refnum as fullpath 
    204205         
  • 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//////////////////////////////////////////////////////////// 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r685 r708  
    7676//   should always default to... 
    7777// 
     78 
     79// - to add --- 
     80// -- wavelength distribution = another RNG to select the wavelength 
     81// -- quartz windows (an empirical model?? or measure some real data - power Law + background) 
     82// -- blocked beam (measure this too, and have some empirical model for this too - Broad Peak) 
    7883// 
    7984 
     
    143148                        retWave0[0] = -1*trunc(datetime-gInitTime)              //to initialize ran3 
    144149                        ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0) 
    145                         Print "started thread 0" 
     150                        Print "started thread 1" 
    146151                endif 
    147152                if(i==1) 
     
    150155                        retWave1[0] = -1*trunc(datetime-gInitTime-2)            //to initialize ran1 
    151156                        ThreadStart mt,i,Monte_SANS_W2(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1) 
    152                         Print "started thread 1" 
     157                        Print "started thread 2" 
    153158                endif 
    154159                if(i==2) 
     
    156161                        retWave2[0] = -1*trunc(datetime-gInitTime-3)            //to initialize ran3a 
    157162                        ThreadStart mt,i,Monte_SANS_W3(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2) 
     163                        Print "started thread 3" 
    158164                endif 
    159165                if(i==3) 
     
    161167                        retWave3[0] = -1*trunc(datetime-gInitTime-4)            //to initialize ran1a 
    162168                        ThreadStart mt,i,Monte_SANS_W4(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3) 
     169                        Print "started thread 4" 
    163170                endif 
    164171        endfor 
     
    12921299        t0 = (stopMSTimer(-2) - t0)*1e-6 
    12931300        t0 *= imon/1000/ThreadProcessorCount                    //projected time, in seconds (using threads for the calculation) 
    1294  
     1301        t0 *= 2         //empirical correction 
    12951302        Print "Estimated Simulation time (s) = ",t0 
    12961303         
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NCNR_Utils.ipf

    r683 r708  
    153153End 
    154154 
    155 // 
    156 // 
    157 //      this should end up in FACILITY_Utils once it's fully functional 
    158 // 
    159 // lots of unnecessary junk... 
     155 
    160156// 
    161157//********************** 
    162 // 2D resolution function calculation - in terms of X and Y 
     158// 2D resolution function calculation - ***NOT*** in terms of X and Y 
     159// but written in terms of Parallel and perpendicular to the Q vector at each point 
     160// 
     161// -- it must be written this way since the 2D function is an ellipse with its major 
     162// axis pointing in the direction of Q_parallel. Hence there is no way to properly define the  
     163// elliptical gaussian in terms of sigmaX and sigmaY 
    163164// 
    164165// based on notes from David Mildner, 2008 
     166// 
    165167// 
    166168// - 21 MAR 07 uses projected BS diameter on the detector 
     
    187189        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron 
    188190        Variable g = 981.0                              //gravity acceleration [cm/s^2] 
    189         Variable h_m = 3995                             // h/m [=] A*m/s 
     191        Variable m_h    = 252.8                 // m/h [=] s/cm^2 
    190192 
    191193        String results 
     
    251253//      endif 
    252254 
    253         Variable kap,a_val 
     255        Variable kap,a_val,a_val_L2 
    254256         
    255257        kap = 2*pi/lambda 
    256         a_val = (L1+L2)*g/2/(h_m)^2 
    257  
    258 //      lambdaWidth = 0.5        
    259         SigmaQX = 3*(S1/L1)^2 + 3*(S2/LP)^2 + 2*(DDet/L2)^2 + 2*(r_dist/L2)^2*(lambdaWidth)^2*(cos(phi))^2 
    260  
    261         SigmaQY = 3*(S1/L1)^2 + 3*(S2/LP)^2 + 2*(DDet/L2)^2 + 2*(r_dist/L2)^2*(lambdaWidth)^2*(sin(phi))^2 + 8*(a_val/L2)^2*lambda^4*lambdaWidth^2 
    262 //      SigmaQY = 3*(S1/L1)^2 + 3*(S2/LP)^2 + 2*(DDet/L2)^2 + 2*(r_dist/L2)^2*(lambdaWidth)^2*(sin(phi))^2 
    263  
    264         SigmaQX = sqrt(kap*kap/12*SigmaQX) 
    265         SigmaQy = sqrt(kap*kap/12*SigmaQY) 
    266  
     258        a_val = L2*(L1+L2)*g/2*(m_h)^2 
     259        a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2 
     260 
     261 
     262///////// in terms of Q_parallel ("x") and Q_perp ("y") - this works, since parallel is in the direction of Q and I 
     263// can calculate that from the QxQy (I just need the projection) 
     264 
     265 
     266// for test case with no gravity, set a_val = 0 
     267// note that gravity has no wavelength dependence. the lambda^4 cancels out. 
     268// 
     269//      a_val = 0 
     270 
     271        // the detector pixel is square, so correct for phi 
     272        DDet = DDet*cos(phi) + DDet*sin(phi) 
     273         
     274        // this is really sigma_Q_parallel 
     275        SigmaQX = kap*kap/12* (3*(S1/L1)^2 + 3*(S2/LP)^2 + (DDet/L2)^2 + (del_r/L2)^2 + (sin(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2) 
     276        SigmaQX += inQ*inQ*v_lambda 
     277 
     278        //this is really sigma_Q_perpendicular 
     279        SigmaQY = 3*(S1/L1)^2 + 3*(S2/LP)^2 + (DDet/L2)^2 + (del_r/L2)^2 + (cos(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2 
     280        SigmaQY *= kap*kap/12 
     281         
     282        SigmaQX = sqrt(SigmaQX) 
     283        SigmaQy = sqrt(SigmaQY) 
     284         
     285         
    267286        results = "success" 
    268287        Return results 
     
    13701389        Return filename 
    13711390End 
     1391 
     1392 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/RawWindowHook.ipf

    r665 r708  
    261261        return 0 
    262262end 
    263          
    264 //function to calculate the overall q-value, given all of the necesary trig inputs 
    265 //NOTE: detector locations passed in are pixels = 0.5cm real space on the detector 
    266 //and are in detector coordinates (1,128) rather than axis values 
    267 //the pixel locations need not be integers, reals are ok inputs 
    268 //sdd is in meters 
    269 //wavelength is in Angstroms 
    270 // 
    271 //returned magnitude of Q is in 1/Angstroms 
    272 // 
    273 Function CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    274         Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
    275          
    276         Variable dx,dy,qval,two_theta,dist 
    277          
    278         Variable pixSizeX=pixSize 
    279         Variable pixSizeY=pixSize 
    280          
    281         sdd *=100               //convert to cm 
    282         dx = (xaxval - xctr)*pixSizeX           //delta x in cm 
    283         dy = (yaxval - yctr)*pixSizeY           //delta y in cm 
    284         dist = sqrt(dx^2 + dy^2) 
    285          
    286         two_theta = atan(dist/sdd) 
    287  
    288         qval = 4*Pi/lam*sin(two_theta/2) 
    289          
    290         return qval 
    291 End 
    292  
    293 //calculates just the q-value in the x-direction on the detector 
    294 //input/output is the same as CalcQval() 
    295 //ALL inputs are in detector coordinates 
    296 // 
    297 //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
    298 //sdd is in meters 
    299 //wavelength is in Angstroms 
    300 // 
    301 // repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) 
    302 // now properly accounts for qz 
    303 // 
    304 Function CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    305         Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
    306  
    307         Variable qx,qval,phi,dx,dy,dist,two_theta 
    308          
    309         qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    310          
    311         sdd *=100               //convert to cm 
    312         dx = (xaxval - xctr)*pixSize            //delta x in cm 
    313         dy = (yaxval - yctr)*pixSize            //delta y in cm 
    314         phi = FindPhi(dx,dy) 
    315          
    316         //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
    317         dist = sqrt(dx^2 + dy^2) 
    318         two_theta = atan(dist/sdd) 
    319  
    320         qx = qval*cos(two_theta/2)*cos(phi) 
    321          
    322         return qx 
    323 End 
    324  
    325 //calculates just the q-value in the y-direction on the detector 
    326 //input/output is the same as CalcQval() 
    327 //ALL inputs are in detector coordinates 
    328 //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
    329 //sdd is in meters 
    330 //wavelength is in Angstroms 
    331 // 
    332 // repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) 
    333 // now properly accounts for qz 
    334 // 
    335 Function CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    336         Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
    337          
    338         Variable dy,qval,dx,phi,qy,dist,two_theta 
    339          
    340         qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    341          
    342         sdd *=100               //convert to cm 
    343         dx = (xaxval - xctr)*pixSize            //delta x in cm 
    344         dy = (yaxval - yctr)*pixSize            //delta y in cm 
    345         phi = FindPhi(dx,dy) 
    346          
    347         //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
    348         dist = sqrt(dx^2 + dy^2) 
    349         two_theta = atan(dist/sdd) 
    350          
    351         qy = qval*cos(two_theta/2)*sin(phi) 
    352          
    353         return qy 
    354 End 
    355  
    356 //calculates just the z-component of the q-vector, not measured on the detector 
    357 //input/output is the same as CalcQval() 
    358 //ALL inputs are in detector coordinates 
    359 //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
    360 //sdd is in meters 
    361 //wavelength is in Angstroms 
    362 // 
    363 // not actually used, but here for completeness if anyone asks 
    364 // 
    365 Function CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    366         Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize 
    367          
    368         Variable dy,qval,dx,phi,qz,dist,two_theta 
    369          
    370         qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) 
    371          
    372         sdd *=100               //convert to cm 
    373          
    374         //get scattering angle to project onto flat detector => Qr = qval*cos(theta) 
    375         dx = (xaxval - xctr)*pixSize            //delta x in cm 
    376         dy = (yaxval - yctr)*pixSize            //delta y in cm 
    377         dist = sqrt(dx^2 + dy^2) 
    378         two_theta = atan(dist/sdd) 
    379          
    380         qz = qval*sin(two_theta/2) 
    381          
    382         return qz 
    383 End 
    384  
    385 //phi is defined from +x axis, proceeding CCW around [0,2Pi] 
    386 Threadsafe Function FindPhi(vx,vy) 
    387         variable vx,vy 
    388          
    389         variable phi 
    390          
    391         phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2 
    392          
    393         // special cases 
    394         if(vx==0 && vy > 0) 
    395                 return(pi/2) 
    396         endif 
    397         if(vx==0 && vy < 0) 
    398                 return(3*pi/2) 
    399         endif 
    400         if(vx >= 0 && vy == 0) 
    401                 return(0) 
    402         endif 
    403         if(vx < 0 && vy == 0) 
    404                 return(pi) 
    405         endif 
    406          
    407          
    408         if(vx > 0 && vy > 0) 
    409                 return(phi) 
    410         endif 
    411         if(vx < 0 && vy > 0) 
    412                 return(phi + pi) 
    413         endif 
    414         if(vx < 0 && vy < 0) 
    415                 return(phi + pi) 
    416         endif 
    417         if( vx > 0 && vy < 0) 
    418                 return(phi + 2*pi) 
    419         endif 
    420          
    421         return(phi) 
    422 end 
    423  
    424263 
    425264//function to set the q-axis scaling after the data has been read in 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/WriteQIS.ipf

    r703 r708  
    641641 
    642642// NEW additions - May 2009 
    643 //ASCII export of data as 7-columns qx-qy-Intensity-qz-sigmaQx-sigmaQy-fShad 
     643//ASCII export of data as 7-columns qx-qy-Intensity-qz-sigmaQ_parall-sigmaQ_perp-fShad 
    644644//limited header information? 
    645645// 
     
    649649// -- when the Qz and resolution are written, be sure to change the tw[15] in the header back to the  
    650650//              proper labels 
    651 // 
     651// - May 2010: 
     652// now the smearing is correct, and is now defined in terms of Q_parallel and Q_perpendicular 
     653// 
     654// - June 2010: 
     655// TEMPORARY - I've added "fake" error that is sqrt(value). It really needs to be propogated 
     656//  in a more correct way, but this is at least a placeholder for the error column 
    652657// 
    653658// - creates the qx and qy data here, based on the data and header information 
     
    736741        labelWave[13] = "" 
    737742        labelWave[14] = "*** Data written from "+type+" folder and may not be a fully corrected data file ***" 
    738         labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy)" 
    739 //      labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy) - Qz - SigmaQx - SigmaQy - fSubS(beam stop shadow)" 
    740         labelWave[16] = "" 
     743//      labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy)" 
     744//      labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy) - Qz - SigmaQ_parall - SigmaQ_perp - fSubS(beam stop shadow)" 
     745        labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy) - err(I) - Qz - SigmaQ_parall - SigmaQ_perp - fSubS(beam stop shadow)" 
     746        labelWave[16] = "ERROR WAVE IS ONLY AN ESTIMATE  - 6/2010" 
    741747        labelWave[17] = "ASCII data created " +date()+" "+time() 
    742748        //strings can be too long to print-- must trim to 255 chars 
     
    812818//********************* 
    813819 
     820        // generate my own error wave for I(qx,qy) 
     821        Duplicate/O z_val sw 
     822        sw = sqrt(z_val)                //assumes Poisson statistics for each cell (counter) 
     823        //      sw = 0.05*sw            // uniform 5% error? tends to favor the low intensity too strongly 
     824        // get rid of the "bad" errorsby replacing the NaN, Inf, and zero with V_avg 
     825        // THIS IS EXTREMEMLY IMPORTANT - if this is not done, there are some "bad" values in the  
     826        // error wave (things that are not numbers) - and this wrecks the smeared model fitting. 
     827        // It appears to have no effect on the unsmeared model. 
     828        WaveStats/Q sw 
     829        sw = numtype(sw[p]) == 0 ? sw[p] : V_avg 
     830        sw = sw[p] != 0 ? sw[p] : V_avg 
     831         
     832 
    814833        //not demo-compatible, but approx 8x faster!!    
    815 #if(cmpstr(stringbykey("IGORKIND",IgorInfo(0),":",";"),"pro") == 0)      
    816         Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val as fullpath       // without resolution 
    817 //      Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val,qz_val,SigmaQx,SigmaQy,fSubS as fullpath  // write out the resolution information 
     834#if(cmpstr(stringbykey("IGORKIND",IgorInfo(0),":",";"),"pro") == 0) 
     835        Duplicate/O qx_val,qx_val_s 
     836        Duplicate/O qy_val,qy_val_s 
     837        Duplicate/O qz_val,qz_val_s 
     838        Duplicate/O z_val,z_val_s 
     839        Duplicate/O SigmaQx,sigmaQx_s 
     840        Duplicate/O SigmaQy,sigmaQy_s 
     841        Duplicate/O fSubS,fSubS_s 
     842        Duplicate/O sw,sw_s 
     843         
     844        //so that double precision is not written ou 
     845        Redimension/S qx_val_s,qy_val_s,qz_val_s,z_val_s,sigmaQx_s,sigmaQy_s,fSubS_s,sw_s 
     846         
     847//      Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val as fullpath       // without resolution 
     848        Save/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 fullpath       // write out the resolution information 
    818849#else 
    819850        Open refNum as fullpath 
    820851        wfprintf refNum,"%s\r\n",labelWave 
    821852        fprintf refnum,"\r\n" 
    822         wfprintf refNum,"%8g\t%8g\t%8g\r\n",qx_val,qy_val,z_val 
     853//      wfprintf refNum,"%8g\t%8g\t%8g\r\n",qx_val,qy_val,z_val 
     854        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 
    823855        Close refNum 
    824856#endif 
     857         
     858        KillWaves/Z qx_val_s,qy_val_s,z_val_s,qz_val_s,SigmaQx_s,SigmaQy_s,fSubS_s 
    825859         
    826860        Killwaves/Z spWave,labelWave,qx_val,qy_val,z_val,qval,qz_val,sigmaQx,SigmaQy,fSubS,phi,r_dist 
Note: See TracChangeset for help on using the changeset viewer.