Ignore:
Timestamp:
Nov 24, 2010 6:19:12 PM (12 years ago)
Author:
ajj
Message:

Merging Dev changes r704:765 into Release trunk

Location:
sans/Release/trunk/NCNR_User_Procedures
Files:
60 edited
3 copied

Legend:

Unmodified
Added
Removed
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/CoreShellCylinder_v40.ipf

    r633 r766  
    5252        // Setup parameter table for model function 
    5353        make/o/d smear_coef_cscyl =  {1.,20.,10.,400,1.0e-6,4.0e-6,1.0e-6,0.01} 
    54         make/o/t smear_parameters_cscyl = {"scale","core radius (A)","shell radius (A)","length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"} 
     54        make/o/t smear_parameters_cscyl = {"scale","core radius (A)","shell THICKNESS (A)","CORE length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"} 
    5555        Edit smear_parameters_cscyl,smear_coef_cscyl 
    5656         
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/CoreShellCyl2D_v40.ipf

    r570 r766  
    124124// 
    125125 
    126 //threaded version of the function 
    127 ThreadSafe Function CoreShellCylinder2D_T(cw,zw,xw,yw,p1,p2) 
    128         WAVE cw,zw,xw,yw 
    129         Variable p1,p2 
    130          
    131 #if exists("CoreShellCylinder_2DX")                     //to hide the function if XOP not installed 
    132  
    133         Make/O/D/N=15 CSCyl2D_tmp 
    134         CSCyl2D_tmp = cw 
    135         CSCyl2D_tmp[14] = 25 
    136         CSCyl2D_tmp[7] = 0              //send a zero background to the calculation, add it in later 
    137  
    138         zw[p1,p2]= CoreShellCylinder_2DX(CSCyl2D_tmp,xw,yw) + cw[7] 
    139          
    140 #endif 
    141  
    142         return 0 
    143 End 
     126////threaded version of the function 
     127//ThreadSafe Function CoreShellCylinder2D_T(cw,zw,xw,yw,p1,p2) 
     128//      WAVE cw,zw,xw,yw 
     129//      Variable p1,p2 
     130//       
     131//#if exists("CoreShellCylinder_2DX")                   //to hide the function if XOP not installed 
     132// 
     133//      Make/O/D/N=15 CSCyl2D_tmp 
     134//      CSCyl2D_tmp = cw 
     135//      CSCyl2D_tmp[14] = 25 
     136//      CSCyl2D_tmp[7] = 0              //send a zero background to the calculation, add it in later 
     137// 
     138//      zw[p1,p2]= CoreShellCylinder_2DX(CSCyl2D_tmp,xw,yw) + cw[7] 
     139//       
     140//#endif 
     141// 
     142//      return 0 
     143//End 
     144// 
     145//// 
     146////  Fit function that is actually a wrapper to dispatch the calculation to N threads 
     147//// 
     148//// nthreads is 1 or an even number, typically 2 
     149//// it doesn't matter if npt is odd. In this case, fractional point numbers are passed 
     150//// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points 
     151//// and the points "2.5" and "3.5" evaluate correctly as 2 and 3 
     152//// 
     153//Function CoreShellCylinder2D(cw,zw,xw,yw) : FitFunc 
     154//      Wave cw,zw,xw,yw 
     155//       
     156//      Variable npt=numpnts(yw) 
     157//      Variable i,nthreads= ThreadProcessorCount 
     158//      variable mt= ThreadGroupCreate(nthreads) 
     159// 
     160//      for(i=0;i<nthreads;i+=1) 
     161//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     162//              ThreadStart mt,i,CoreShellCylinder2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1)) 
     163//      endfor 
     164// 
     165//      do 
     166//              variable tgs= ThreadGroupWait(mt,100) 
     167//      while( tgs != 0 ) 
     168// 
     169//      variable dummy= ThreadGroupRelease(mt) 
     170//       
     171//      return(0) 
     172//End 
    144173 
    145174// 
     
    154183        Wave cw,zw,xw,yw 
    155184         
    156         Variable npt=numpnts(yw) 
    157         Variable i,nthreads= ThreadProcessorCount 
    158         variable mt= ThreadGroupCreate(nthreads) 
     185#if exists("CoreShellCylinder_2DX")                     //to hide the function if XOP not installed 
    159186 
    160         for(i=0;i<nthreads;i+=1) 
    161         //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
    162                 ThreadStart mt,i,CoreShellCylinder2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1)) 
    163         endfor 
     187        Make/O/D/N=15 CSCyl2D_tmp 
     188        CSCyl2D_tmp = cw 
     189        CSCyl2D_tmp[14] = 25 
     190        CSCyl2D_tmp[7] = 0              //send a zero background to the calculation, add it in later 
    164191 
    165         do 
    166                 variable tgs= ThreadGroupWait(mt,100) 
    167         while( tgs != 0 ) 
     192        MultiThread zw= CoreShellCylinder_2DX(CSCyl2D_tmp,xw,yw) + cw[7] 
     193         
     194#endif 
    168195 
    169         variable dummy= ThreadGroupRelease(mt) 
    170          
    171196        return(0) 
    172197End 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Cylinder_2D_v40.ipf

    r570 r766  
    122122// 
    123123 
    124 //threaded version of the function 
    125 ThreadSafe Function Cylinder2D_T(cw,zw,xw,yw,p1,p2) 
    126         WAVE cw,zw,xw,yw 
    127         Variable p1,p2 
     124////threaded version of the function 
     125//ThreadSafe Function Cylinder2D_T(cw,zw,xw,yw,p1,p2) 
     126//      WAVE cw,zw,xw,yw 
     127//      Variable p1,p2 
     128//       
     129//#if exists("Cylinder_2DX")                    //to hide the function if XOP not installed 
     130// 
     131//      Make/O/D/N=12 Cyl2D_tmp                         // there seems to be no speed penalty for doing this... 
     132//      Cyl2D_tmp = cw 
     133//      Cyl2D_tmp[11] = 25                                      // hard-wire the number of integration points 
     134//      Cyl2D_tmp[5] = 0                                                // send a background of zero 
     135//       
     136//      zw[p1,p2]= Cylinder_2DX(Cyl2D_tmp,xw,yw) + cw[5]                //add in the proper background here 
     137// 
     138//#endif 
     139// 
     140//      return 0 
     141//End 
     142// 
     143//// 
     144////  Fit function that is actually a wrapper to dispatch the calculation to N threads 
     145//// 
     146//// nthreads is 1 or an even number, typically 2 
     147//// it doesn't matter if npt is odd. In this case, fractional point numbers are passed 
     148//// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points 
     149//// and the points "2.5" and "3.5" evaluate correctly as 2 and 3 
     150//// 
     151//Function Cylinder2D(cw,zw,xw,yw) : FitFunc 
     152//      Wave cw,zw,xw,yw 
     153//       
     154//      Variable npt=numpnts(yw) 
     155//      Variable ii,nthreads= ThreadProcessorCount 
     156//      variable mt= ThreadGroupCreate(nthreads) 
     157// 
     158////    Variable t1=StopMSTimer(-2) 
     159//       
     160//      for(ii=0;ii<nthreads;ii+=1) 
     161//      //      Print (ii*npt/nthreads),((ii+1)*npt/nthreads-1) 
     162//              ThreadStart mt,ii,Cylinder2D_T(cw,zw,xw,yw,(ii*npt/nthreads),((ii+1)*npt/nthreads-1)) 
     163//      endfor 
     164// 
     165//      do 
     166//              variable tgs= ThreadGroupWait(mt,100) 
     167//      while( tgs != 0 ) 
     168// 
     169//      variable dummy= ThreadGroupRelease(mt) 
     170//       
     171////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
     172//       
     173//      return(0) 
     174//End 
     175 
     176 
     177/// now using the MultiThread keyword. as of Igor 6.20, the manual threading 
     178// as above gives a wave read error (index out of range). Same code works fine in Igor 6.12 
     179Function Cylinder2D(cw,zw,xw,yw) : FitFunc 
     180        Wave cw,zw,xw,yw 
    128181         
     182        //      Variable t1=StopMSTimer(-2) 
     183 
    129184#if exists("Cylinder_2DX")                      //to hide the function if XOP not installed 
    130185 
     
    134189        Cyl2D_tmp[5] = 0                                                // send a background of zero 
    135190         
    136         zw[p1,p2]= Cylinder_2DX(Cyl2D_tmp,xw,yw) + cw[5]                //add in the proper background here 
     191        MultiThread zw = Cylinder_2DX(Cyl2D_tmp,xw,yw) + cw[5]          //add in the proper background here 
    137192 
    138193#endif 
    139194 
    140         return 0 
    141 End 
    142  
    143 // 
    144 //  Fit function that is actually a wrapper to dispatch the calculation to N threads 
    145 // 
    146 // nthreads is 1 or an even number, typically 2 
    147 // it doesn't matter if npt is odd. In this case, fractional point numbers are passed 
    148 // and the wave indexing works just fine - I tested this with test waves of 7 and 8 points 
    149 // and the points "2.5" and "3.5" evaluate correctly as 2 and 3 
    150 // 
    151 Function Cylinder2D(cw,zw,xw,yw) : FitFunc 
    152         Wave cw,zw,xw,yw 
    153          
    154         Variable npt=numpnts(yw) 
    155         Variable i,nthreads= ThreadProcessorCount 
    156         variable mt= ThreadGroupCreate(nthreads) 
    157  
    158 //      Variable t1=StopMSTimer(-2) 
    159          
    160         for(i=0;i<nthreads;i+=1) 
    161         //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
    162                 ThreadStart mt,i,Cylinder2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1)) 
    163         endfor 
    164  
    165         do 
    166                 variable tgs= ThreadGroupWait(mt,100) 
    167         while( tgs != 0 ) 
    168  
    169         variable dummy= ThreadGroupRelease(mt) 
    170          
    171195//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
    172196         
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Ellipsoid2D_v40.ipf

    r570 r766  
    123123//      return(0) 
    124124//End 
     125//// 
    125126// 
    126  
    127 //threaded version of the function 
    128 ThreadSafe Function Ellipsoid2D_T(cw,zw,xw,yw,p1,p2) 
    129         WAVE cw,zw,xw,yw 
    130         Variable p1,p2 
    131          
    132 #if exists("Ellipsoid_2DX")                     //to hide the function if XOP not installed 
    133  
    134         Make/O/D/N=13 Ellip2D_tmp 
    135         Ellip2D_tmp = cw 
    136         Ellip2D_tmp[12] = 25 
    137         Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later 
    138          
    139         zw[p1,p2]= Ellipsoid_2DX(Ellip2D_tmp,xw,yw) + cw[5] 
    140          
    141 #endif 
    142  
    143         return 0 
    144 End 
     127////threaded version of the function 
     128//ThreadSafe Function Ellipsoid2D_T(cw,zw,xw,yw,p1,p2) 
     129//      WAVE cw,zw,xw,yw 
     130//      Variable p1,p2 
     131//       
     132//#if exists("Ellipsoid_2DX")                   //to hide the function if XOP not installed 
     133// 
     134//      Make/O/D/N=13 Ellip2D_tmp 
     135//      Ellip2D_tmp = cw 
     136//      Ellip2D_tmp[12] = 25 
     137//      Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later 
     138//       
     139//      zw[p1,p2]= Ellipsoid_2DX(Ellip2D_tmp,xw,yw) + cw[5] 
     140//       
     141//#endif 
     142// 
     143//      return 0 
     144//End 
     145// 
     146//// 
     147////  Fit function that is actually a wrapper to dispatch the calculation to N threads 
     148//// 
     149//// nthreads is 1 or an even number, typically 2 
     150//// it doesn't matter if npt is odd. In this case, fractional point numbers are passed 
     151//// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points 
     152//// and the points "2.5" and "3.5" evaluate correctly as 2 and 3 
     153//// 
     154//Function Ellipsoid2D(cw,zw,xw,yw) : FitFunc 
     155//      Wave cw,zw,xw,yw 
     156//       
     157//      Variable npt=numpnts(yw) 
     158//      Variable i,nthreads= ThreadProcessorCount 
     159//      variable mt= ThreadGroupCreate(nthreads) 
     160// 
     161//      for(i=0;i<nthreads;i+=1) 
     162//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     163//              ThreadStart mt,i,Ellipsoid2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1)) 
     164//      endfor 
     165// 
     166//      do 
     167//              variable tgs= ThreadGroupWait(mt,100) 
     168//      while( tgs != 0 ) 
     169// 
     170//      variable dummy= ThreadGroupRelease(mt) 
     171//       
     172//      return(0) 
     173//End 
    145174 
    146175// 
     
    155184        Wave cw,zw,xw,yw 
    156185         
    157         Variable npt=numpnts(yw) 
    158         Variable i,nthreads= ThreadProcessorCount 
    159         variable mt= ThreadGroupCreate(nthreads) 
     186#if exists("Ellipsoid_2DX")                     //to hide the function if XOP not installed 
    160187 
    161         for(i=0;i<nthreads;i+=1) 
    162         //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
    163                 ThreadStart mt,i,Ellipsoid2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1)) 
    164         endfor 
    165  
    166         do 
    167                 variable tgs= ThreadGroupWait(mt,100) 
    168         while( tgs != 0 ) 
    169  
    170         variable dummy= ThreadGroupRelease(mt) 
     188        Make/O/D/N=13 Ellip2D_tmp 
     189        Ellip2D_tmp = cw 
     190        Ellip2D_tmp[12] = 25 
     191        Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later 
     192         
     193        MultiThread zw= Ellipsoid_2DX(Ellip2D_tmp,xw,yw) + cw[5] 
     194         
     195#endif 
    171196         
    172197        return(0) 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/EllipticalCylinder2D_v40.ipf

    r570 r766  
    126126// 
    127127 
    128 //threaded version of the function 
    129 ThreadSafe Function EllipticalCylinder2D_T(cw,zw,xw,yw,p1,p2) 
    130         WAVE cw,zw,xw,yw 
    131         Variable p1,p2 
    132          
    133 #if exists("EllipticalCylinder_2DX")                    //to hide the function if XOP not installed 
    134  
    135         Make/O/D/N=15 EllCyl2D_tmp 
    136         EllCyl2D_tmp = cw 
    137         EllCyl2D_tmp[14] = 25 
    138         EllCyl2D_tmp[6] = 0             //pass in a zero background and add it in later 
    139          
    140         zw[p1,p2]= EllipticalCylinder_2DX(EllCyl2D_tmp,xw,yw) + cw[6] 
    141          
    142 #endif 
    143  
    144         return 0 
    145 End 
     128////threaded version of the function 
     129//ThreadSafe Function EllipticalCylinder2D_T(cw,zw,xw,yw,p1,p2) 
     130//      WAVE cw,zw,xw,yw 
     131//      Variable p1,p2 
     132//       
     133//#if exists("EllipticalCylinder_2DX")                  //to hide the function if XOP not installed 
     134// 
     135//      Make/O/D/N=15 EllCyl2D_tmp 
     136//      EllCyl2D_tmp = cw 
     137//      EllCyl2D_tmp[14] = 25 
     138//      EllCyl2D_tmp[6] = 0             //pass in a zero background and add it in later 
     139//       
     140//      zw[p1,p2]= EllipticalCylinder_2DX(EllCyl2D_tmp,xw,yw) + cw[6] 
     141//       
     142//#endif 
     143// 
     144//      return 0 
     145//End 
     146// 
     147//// 
     148////  Fit function that is actually a wrapper to dispatch the calculation to N threads 
     149//// 
     150//// nthreads is 1 or an even number, typically 2 
     151//// it doesn't matter if npt is odd. In this case, fractional point numbers are passed 
     152//// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points 
     153//// and the points "2.5" and "3.5" evaluate correctly as 2 and 3 
     154//// 
     155//Function EllipticalCylinder2D(cw,zw,xw,yw) : FitFunc 
     156//      Wave cw,zw,xw,yw 
     157//       
     158//      Variable npt=numpnts(yw) 
     159//      Variable i,nthreads= ThreadProcessorCount 
     160//      variable mt= ThreadGroupCreate(nthreads) 
     161// 
     162//      for(i=0;i<nthreads;i+=1) 
     163//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     164//              ThreadStart mt,i,EllipticalCylinder2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1)) 
     165//      endfor 
     166// 
     167//      do 
     168//              variable tgs= ThreadGroupWait(mt,100) 
     169//      while( tgs != 0 ) 
     170// 
     171//      variable dummy= ThreadGroupRelease(mt) 
     172//       
     173//      return(0) 
     174//End 
    146175 
    147176// 
     
    156185        Wave cw,zw,xw,yw 
    157186         
    158         Variable npt=numpnts(yw) 
    159         Variable i,nthreads= ThreadProcessorCount 
    160         variable mt= ThreadGroupCreate(nthreads) 
     187#if exists("EllipticalCylinder_2DX")                    //to hide the function if XOP not installed 
    161188 
    162         for(i=0;i<nthreads;i+=1) 
    163         //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
    164                 ThreadStart mt,i,EllipticalCylinder2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1)) 
    165         endfor 
    166  
    167         do 
    168                 variable tgs= ThreadGroupWait(mt,100) 
    169         while( tgs != 0 ) 
    170  
    171         variable dummy= ThreadGroupRelease(mt) 
     189        Make/O/D/N=15 EllCyl2D_tmp 
     190        EllCyl2D_tmp = cw 
     191        EllCyl2D_tmp[14] = 25 
     192        EllCyl2D_tmp[6] = 0             //pass in a zero background and add it in later 
     193         
     194        MultiThread zw= EllipticalCylinder_2DX(EllCyl2D_tmp,xw,yw) + cw[6] 
     195         
     196#endif 
    172197         
    173198        return(0) 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/PeakGauss2D_v40.ipf

    r570 r766  
    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/Release/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Sphere_2D_v40.ipf

    r570 r766  
    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 
     190//// the threaded version must be specifically written, since 
     191//// FUNCREF can't be passed into  a threaded calc (structures can't be passed either) 
     192// so in this implementation, the smearing is dispatched as threads to a function that 
     193// can calculate the function for a range of points in the input qxqyqz. It is important 
     194// that the worker calls the un-threaded model function (so write one) and that in the (nord x nord) 
     195// loop, vectors of length (nord) are calculated rather than pointwise, since the model 
     196// function is AAO. 
     197// -- makes things rather messy to code individual functions, but I really see no other way 
     198// given the restrictions of what can be passed to threaded functions. 
     199// 
     200// 
     201// The smearing is handled this way since 1D smearing is 20 x 200 pts = 4000 evaluations 
     202// and the 2D is (10 x 10) x 16000 pts = 1,600,000 evaluations (if it's done like the 1D, it's 4000x slower) 
     203// 
     204// 
     205// - the threading gives a clean speedup of 2 for N=2, even for this simple calculation 
     206//  -- 4.8X speedup for N=8 (4 real cores + 4 virtual cores) 
     207// 
     208// nord = 5,10,20 allowed 
     209// 
    181210Function SmearedSphere2D(s) 
    182211        Struct ResSmear_2D_AAOStruct &s 
    183212         
    184         Smear_2DModel_5(Sphere2D_noThread,s) 
    185 //      Smear_2DModel_20(Sphere2D_noThread,s) 
     213//// non-threaded, but generic calculation 
     214//// the last param is nord      
     215//      Smear_2DModel_PP(Sphere2D_noThread,s,10) 
     216 
     217 
     218//// the last param is nord 
     219        SmearSphere2D_THR(s,10)          
     220 
    186221        return(0) 
    187222end 
     
    197232        WAVE qy = $(DF+str+"_qy") 
    198233        WAVE qz = $(DF+str+"_qz") 
    199         WAVE sigQx = $(DF+str+"_sigQx") 
    200         WAVE sigQy = $(DF+str+"_sigQy") 
     234        WAVE sQpl = $(DF+str+"_sQpl") 
     235        WAVE sQpp = $(DF+str+"_sQpp") 
    201236        WAVE shad = $(DF+str+"_fs") 
    202237         
     
    204239        WAVE s.coefW = coefW     
    205240        WAVE s.zw = resultW      
    206         WAVE s.qx = qx 
    207         WAVE s.qy = qy 
     241        WAVE s.xw[0] = qx 
     242        WAVE s.xw[1] = qy 
    208243        WAVE s.qz = qz 
    209         WAVE s.sigQx = sigQx 
    210         WAVE s.sigQy = sigQy 
     244        WAVE s.sQpl = sQpl 
     245        WAVE s.sQpp = sQpp 
    211246        WAVE s.fs = shad 
    212247         
     
    216251        return (0) 
    217252End 
     253 
     254 
     255 
     256 
     257// 
     258// this is the threaded version, that dispatches the calculation out to threads 
     259// 
     260// must be written specific to each 2D function 
     261//  
     262Function SmearSphere2D_THR(s,nord) 
     263        Struct ResSmear_2D_AAOStruct &s 
     264        Variable nord 
     265         
     266        String weightStr,zStr 
     267         
     268// create all of the necessary quadrature waves here - rather than inside a threadsafe function 
     269        switch(nord)     
     270                case 5:          
     271                        weightStr="gauss5wt" 
     272                        zStr="gauss5z" 
     273                        if (WaveExists($weightStr) == 0) 
     274                                Make/O/D/N=(nord) $weightStr,$zStr 
     275                                Make5GaussPoints($weightStr,$zStr)       
     276                        endif 
     277                        break                            
     278                case 10:                 
     279                        weightStr="gauss10wt" 
     280                        zStr="gauss10z" 
     281                        if (WaveExists($weightStr) == 0) 
     282                                Make/O/D/N=(nord) $weightStr,$zStr 
     283                                Make10GaussPoints($weightStr,$zStr)      
     284                        endif 
     285                        break                            
     286                case 20:                 
     287                        weightStr="gauss20wt" 
     288                        zStr="gauss20z" 
     289                        if (WaveExists($weightStr) == 0) 
     290                                Make/O/D/N=(nord) $weightStr,$zStr 
     291                                Make20GaussPoints($weightStr,$zStr)      
     292                        endif 
     293                        break 
     294                default:                                                         
     295                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                         
     296        endswitch 
     297         
     298        Wave/Z wt = $weightStr 
     299        Wave/Z xi = $zStr               // wave references to pass 
     300 
     301        Variable npt=numpnts(s.xw[0]) 
     302        Variable i,nthreads= ThreadProcessorCount 
     303        variable mt= ThreadGroupCreate(nthreads) 
     304 
     305        Variable t1=StopMSTimer(-2) 
     306         
     307        for(i=0;i<nthreads;i+=1) 
     308//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     309                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) 
     310        endfor 
     311 
     312        do 
     313                variable tgs= ThreadGroupWait(mt,100) 
     314        while( tgs != 0 ) 
     315 
     316        variable dummy= ThreadGroupRelease(mt) 
     317         
     318// comment out the threading + uncomment this for testing to make sure that the single thread works 
     319//      nThreads=1 
     320//      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) 
     321 
     322        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
     323         
     324        return(0) 
     325end 
     326 
     327// 
     328// - worker function for threads of Sphere2D 
     329// 
     330ThreadSafe Function SmearSphere2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord) 
     331        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi 
     332        Variable pt1,pt2,nord 
     333         
     334// now passed in.... 
     335//      Wave wt = $weightStr 
     336//      Wave xi = $zStr          
     337 
     338        Variable ii,jj,kk,num 
     339        Variable qx,qy,qz,qval,sx,sy,fs 
     340        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut 
     341         
     342        Variable normFactor,phi,theta,maxSig,numStdDev=3 
     343         
     344/// keep these waves local 
     345        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW 
     346         
     347        // now just loop over the points as specified 
     348         
     349        answer=0 
     350         
     351        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     352 
     353        //loop over q-values 
     354        for(ii=pt1;ii<(pt2+1);ii+=1) 
     355                 
     356                qx = qxw[ii] 
     357                qy = qyw[ii] 
     358                qz = qzw[ii] 
     359                qval = sqrt(qx^2+qy^2+qz^2) 
     360                spl = sxw[ii] 
     361                spp = syw[ii] 
     362                fs = fsw[ii] 
     363                 
     364                normFactor = 2*pi*spl*spp 
     365                 
     366                phi = FindPhi(qx,qy) 
     367                 
     368                apl = -numStdDev*spl + qval             //parallel = q integration limits 
     369                bpl = numStdDev*spl + qval 
     370                app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
     371                bpp = numStdDev*spp + phi 
     372                 
     373                //make sure the limits are reasonable. 
     374                if(apl < 0) 
     375                        apl = 0 
     376                endif 
     377                // do I need to specially handle limits when phi ~ 0? 
     378         
     379                 
     380                sumOut = 0 
     381                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
     382                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     383 
     384                        sumIn=0 
     385                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
     386 
     387                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
     388                                 
     389                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     390                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
     391                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
     392                                 
     393                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     394                                res_tot[kk] /= normFactor 
     395//                              res_tot[kk] *= fs 
     396 
     397                        endfor 
     398                         
     399                        Sphere2D_noThread(coef,fcnRet,xptw,yptw)                        //fcn passed in is an AAO 
     400                         
     401                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
     402                        fcnRet *= wt[jj]*wt*res_tot 
     403                        // 
     404                        answer += (bpl-apl)/2.0*sum(fcnRet)             // 
     405                endfor 
     406 
     407                answer *= (bpp-app)/2.0 
     408                zw[ii] = answer 
     409        endfor 
     410         
     411        return(0) 
     412end 
     413 
     414 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/LamellarPS_HG_v40.ipf

    r570 r766  
    206206        ii=0 
    207207        Sq = 0 
    208         for(ii=1;ii<(NN-1);ii+=1) 
     208        for(ii=1;ii<=(NN-1);ii+=1) 
    209209                temp = 0 
    210210                alpha = Cp/4/pi/pi*(ln(pi*ii) + Euler) 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/LamellarPS_v40.ipf

    r570 r766  
    198198        ii=0 
    199199        Sq = 0 
    200         for(ii=1;ii<(NN-1);ii+=1) 
     200        for(ii=1;ii<=(NN-1);ii+=1) 
    201201                temp = 0 
    202202                alpha = Cp/4/pi/pi*(ln(pi*ii) + Euler) 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/Fractal_PolyCore_v40.ipf

    r570 r766  
    160160        //calculate S(q) 
    161161        sq = Df*exp(gammln(Df-1))*sin((Df-1)*atan(x*corr)) 
    162         sq /= (x*r0)^Df * (1 + 1/(x*corr)^2)^((Df-1)/2) 
     162        sq /= (x*(r0+thick))^Df * (1 + 1/(x*corr)^2)^((Df-1)/2) 
    163163        sq += 1 
    164164        //combine and return 
    165165        ans = pq*sq + bkg 
    166  
     166        //ans = pq +bkg 
     167         
    167168        return (ans) 
    168169End 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/FuzzySpheres_v40.ipf

    r690 r766  
    100100        Variable scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,lor_sf,lor_len 
    101101        Variable I0, L 
    102         I0 = w[0] 
    103         L = w[1] 
     102        //I0 = w[0] 
     103        //L = w[1] 
    104104 
    105105        //set up the coefficient values 
     
    246246        return(retval) 
    247247End 
     248 
     249Proc PlotFuzzySphereSLD() 
     250 
     251        Variable scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,lor_sf,lor_len 
     252        Variable I0, L 
     253        scale=coef_fuzz[0] 
     254        rad=coef_fuzz[1] 
     255        pd=coef_fuzz[2] 
     256        sig=pd*rad 
     257        sig_surf = coef_fuzz[3] 
     258        rho=coef_fuzz[4] 
     259        rhos=coef_fuzz[5] 
     260        delrho=rho-rhos 
     261        I0 = coef_fuzz[6]               //for the Lorentzian 
     262        L = coef_fuzz[7] 
     263        bkg=coef_fuzz[8] 
     264 
     265        Make/O/N=1000 fuzzy_rwave 
     266        Make/O/N=1000 fuzzy_sldwave 
     267         
     268        fuzzy_rwave = p*(rad+4*sig_surf)/numpnts(fuzzy_rwave) 
     269         
     270         
     271        fuzzy_sldwave = erf(fuzzy_rwave) 
     272         
     273        Display fuzzy_sldwave vs fuzzy_rwave 
     274 
     275End 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Packages/GeneticOptimization/NCNR_GenFitUtils.ipf

    r570 r766  
    2929// X add the XOP to the distribution, or instructions (better to NOT bundle...) 
    3030// X odd bug where first "fit" fails on AppendToGraph as GenCurveFit is started. May need to disable /D flag 
    31 // 
    32 // -- need to pass back the chi-squared and number of points. "V_" globals don't appear 
     31// X- need to pass back the chi-squared and number of points. "V_" globals don't appear 
     32// -- add the flags to calculate the residuals (/R), and then manually append the residual wave to the graph. 
     33// 
    3334// 
    3435// for the speed test. try writing my own wrapper for an unsmeared calculation, and see if it's still dog-slow. 
     
    154155// need to pass back the chi-squared and number of points. "V_" globals don't appear 
    155156// 
     157// function will only compile if the GenCurveFit XOP is installed - else it will return zero 
     158// 
    156159Function DoGenCurveFit(useRes,useCursors,sw,fitYw,fs,funcStr,holdStr,val,lolim,hilim,pt1,pt2) 
    157160        Variable useRes,useCursors 
     
    162165        WAVE/T lolim,hilim 
    163166        Variable pt1,pt2                //already sorted if cursors are needed 
     167 
     168#if exists("GenCurveFit") 
    164169 
    165170        // initialise the structure you will use 
     
    254259        // other useful flags:  /N /DUMP /TOL /D=fitYw 
    255260        //  /OPT=1 seems to make no difference whether it's used or not 
    256          
    257 #if exists("GenCurveFit") 
     261        // /N does not speed anything up at all, and you get no feedback 
     262        // /R does work, generating "res_" wave, but the append must be "manual" 
     263         
    258264        // append the fit 
    259265        //do this only because GenCurveFit tries to append too quickly? 
     
    297303                 
    298304        while(0)         
    299 #endif 
    300305         
    301306         
     
    306311 
    307312        return(chi) 
     313#else 
     314        return(0) 
     315#endif   
    308316end 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Packages/ModelPicker/SANSModelPicker_v40.ipf

    r663 r766  
    5454        NewDataFolder/O root:Packages 
    5555        NewDataFolder/O root:Packages:NIST 
     56         
     57        if(cmpstr("Macintosh",IgorInfo(2)) == 0) 
     58                String/G root:Packages:NIST:gAngstStr = num2char(-127) 
     59//              Variable/G root:myGlobals:gIsMac = 1 
     60        else 
     61                //either Windows or Windows NT 
     62                String/G root:Packages:NIST:gAngstStr = num2char(-59) 
     63//              Variable/G root:myGlobals:gIsMac = 0 
     64                //SetIgorOption to keep some PC's (graphics cards?) from smoothing the 2D image 
     65                Execute "SetIgorOption WinDraw,forceCOLORONCOLOR=1" 
     66        endif 
    5667         
    5768        if(!DataFolderExists("root:Packages:NIST:FileList")) 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/Packages/Wrapper_v40.ipf

    r693 r766  
    7979        PopupMenu popup_2,pos={30,93},size={123,20},title="Coefficients" 
    8080        PopupMenu popup_2,mode=1,value= #"W_CoefPopupList()",proc=Coef_PopMenuProc 
    81         CheckBox check_0,pos={440,19},size={79,14},title="Use Cursors?",value= 0 
     81        CheckBox check_0,pos={430,19},size={79,14},title="Use Cursors?",value= 0 
    8282        CheckBox check_0,proc=UseCursorsWrapperProc 
    83         CheckBox check_1,pos={440,42},size={74,14},title="Use Epsilon?",value= 0 
    84         CheckBox check_2,pos={440,65},size={95,14},title="Use Constraints?",value= 0 
     83        CheckBox check_1,pos={430,42},size={74,14},title="Use Epsilon?",value= 0 
     84        CheckBox check_2,pos={430,65},size={95,14},title="Use Constraints?",value= 0 
    8585        CheckBox check_3,pos={530,18},size={72,14},title="2D Functions?",value= 0 
    8686        CheckBox check_3 proc=Toggle2DControlsCheckProc 
    87         CheckBox check_4,pos={440,85},size={72,14},title="Report?",value= 0 
    88         CheckBox check_5,pos={454,103},size={72,14},title="Save it?",value= 0 
     87        CheckBox check_4,pos={430,85},size={72,14},title="Report?",value= 0 
     88        CheckBox check_5,pos={444,103},size={72,14},title="Save it?",value= 0 
     89        CheckBox check_6,pos={530,42},size={72,14},title="Use Residuals?",value= 0 
     90        CheckBox check_6,proc=UseResidualsCheckProc 
     91        CheckBox check_7,pos={530,65},size={72,14},title="Info Box?",value= 0 
     92        CheckBox check_7,proc=UseInfoTextBoxCheckProc 
    8993        //change draw order to put button over text of checkbox 
    9094        Button button_0,pos={520,93},size={100,20},proc=DoTheFitButton,title="Do 1D Fit" 
     
    679683 
    680684        String folderStr,funcStr,coefStr 
    681         Variable useCursors,useEps,useConstr 
     685        Variable useCursors,useEps,useConstr,useResiduals,useTextBox 
    682686         
    683687        switch( ba.eventCode ) 
     
    699703                        useConstr=V_Value 
    700704                         
     705                        ControlInfo/W=WrapperPanel check_6 
     706                        useResiduals=V_Value 
     707                         
     708                        ControlInfo/W=WrapperPanel check_7 
     709                        useTextBox = V_Value 
     710                         
    701711                        if(!CheckFunctionAndCoef(funcStr,coefStr)) 
    702712                                DoAlert 0,"The coefficients and function type do not match. Please correct the selections in the popup menus." 
     
    704714                        endif 
    705715                         
    706                         FitWrapper(folderStr,funcStr,coefStr,useCursors,useEps,useConstr) 
     716                        FitWrapper(folderStr,funcStr,coefStr,useCursors,useEps,useConstr,useResiduals,useTextBox) 
    707717                         
    708718                        //      DoUpdate (does not work!) 
     
    743753// 
    744754// 
    745 Function FitWrapper(folderStr,funcStr,coefStr,useCursors,useEps,useConstr) 
     755Function FitWrapper(folderStr,funcStr,coefStr,useCursors,useEps,useConstr,useResiduals,useTextBox) 
    746756        String folderStr,funcStr,coefStr 
    747         Variable useCursors,useEps,useConstr 
     757        Variable useCursors,useEps,useConstr,useResiduals,useTextBox 
    748758 
    749759        String suffix=getModelSuffix(funcStr) 
     
    763773        Wave/T lolim=$("LoLim_"+suffix) 
    764774        Wave/T hilim=$("HiLim_"+suffix) 
    765         Wave eps=$("epsilon_"+suffix) 
    766          
     775         
     776        if(useEps) 
     777                Wave eps=$("epsilon_"+suffix) 
     778        endif 
    767779// fill a struct instance whether I need one or not 
    768780        String DF="root:"+folderStr+":"  
     
    782794        fitYw = NaN 
    783795         
    784         Variable useRes=0,isUSANS=0,val 
     796        Variable useResol=0,isUSANS=0,val 
    785797        if(stringmatch(funcStr, "Smear*"))              // if it's a smeared function, need a struct 
    786                 useRes=1 
     798                useResol=1 
    787799        endif 
    788800        if(dimsize(resW,1) > 4) 
    789801                isUSANS=1 
    790802        endif 
     803         
    791804        // do not construct constraints for any of the coefficients that are being held 
    792805        // -- this will generate an "unknown error" from the curve fitting 
    793         Make/O/T/N=0 constr 
     806        // -- if constraints are not used, the constr wave is killed. This apparently 
     807        // confuses the /NWOK flag, since there is not even a null reference present. So generate one. 
    794808        if(useConstr) 
     809                Make/O/T/N=0 constr 
    795810                String constraintExpression 
    796811                Variable i, nPnts=DimSize(lolim, 0),nextRow=0 
     
    809824                        endif 
    810825                endfor 
     826        else 
     827                Wave/T/Z constr = constr 
     828                KillWaves/Z constr 
     829                Wave/T/Z constr = constr                //this is intentionally a null reference 
    811830        endif 
    812831 
     
    850869                        pt2 = pcsr(B) 
    851870                endif 
     871        else 
     872                //if cursors are not being used, find the first and last points of the data set, and pass them 
     873                pt1 = 0 
     874                pt2 = numpnts(yw)-1 
    852875        endif 
    853876                 
     
    866889// The textwave would have to be parsed into a constraint matrix first, then passed as /C={cMat,cVec}. 
    867890// -- just something to watch out for. 
     891 
     892// now two more flags... ,useResiduals,useTextBox 
     893        Variable tb = 1+2+4+8+16+256+512                //See CurveFit docs for bit settings for /TBOX flag 
    868894 
    869895        do 
     
    881907                        // useEps and useConstr are not needed 
    882908                        // pass the structure to get the current waves, including the trimmed USANS matrix 
     909                        // 
     910                        // I don't know that GetCurveFit can do residuals, so I'm not passing that flag, or the text box flag 
     911                        // 
    883912                        Variable chi,pt 
    884913 
    885                         chi = DoGenCurveFit(useRes,useCursors,sw,fitYw,fs,funcStr,getHStr(hold),val,lolim,hilim,pt1,pt2) 
     914                        chi = DoGenCurveFit(useResol,useCursors,sw,fitYw,fs,funcStr,getHStr(hold),val,lolim,hilim,pt1,pt2) 
    886915                        pt = val 
    887916 
     
    890919                endif 
    891920                 
    892                  
    893                 if(useRes && useEps && useCursors && useConstr)         //do it all 
    894                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs 
    895                         break 
    896                 endif 
    897                  
    898                 if(useRes && useEps && useCursors)              //no constr 
    899                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /E=eps /D=fitYw /STRC=fs 
    900                         break 
    901                 endif 
    902                  
    903                 if(useRes && useEps && useConstr)               //no crsr 
    904                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs 
    905                         break 
    906                 endif 
    907                  
    908                 if(useRes && useCursors && useConstr)           //no eps 
    909                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /D=fitYw /C=constr /STRC=fs 
    910                         break 
    911                 endif 
    912                  
    913                 if(useRes && useCursors)                //no eps, no constr 
    914                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /D=fitYw /STRC=fs 
    915                         break 
    916                 endif 
    917                  
    918                 if(useRes && useEps)            //no crsr, no constr 
    919 //                      Print "timing test for Cylinder_PolyRadius --- the supposedly threaded version ---" 
     921                // now useCursors, useEps, and useConstr are all handled w/ /NWOK 
     922                // so there are only three conditions to test == 1 + 3 + 3 + 1 = 8 conditions 
     923                 
     924                if(useResol && useResiduals && useTextBox)              //do it all 
     925                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs /R /NWOK 
     926                        break 
     927                endif 
     928                 
     929                if(useResol && useResiduals)            //res + resid 
     930                        FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs /R /NWOK 
     931                        break 
     932                endif 
     933 
     934                 
     935                if(useResol && useTextBox)              //res + text 
     936                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs /NWOK 
     937                        break 
     938                endif 
     939                 
     940                if(useResol)            //res only 
     941//                      Print "timing test for Cylinder_PolyRadius---" 
    920942//                      Variable t0 = stopMStimer(-2) 
    921                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /E=eps /D=fitYw /STRC=fs 
     943 
     944                        FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs /NWOK 
     945                         
    922946//                      t0 = (stopMSTimer(-2) - t0)*1e-6 
    923947//                      Printf  "CylPolyRad fit time using res and eps and /NTHR=0 time = %g seconds\r\r",t0 
    924                          
    925948//                      cw[0] = .01 
    926949//                      cw[1] = 20 
     
    931954//                       
    932955//                      t0 = stopMSTimer(-2) 
    933 //                      FuncFit/H=getHStr(hold) $funcStr cw, yw /X=xw /W=sw /I=1 /E=eps /D=fitYw /STRC=fs 
     956//                      FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs /NWOK 
    934957//                      t0 = (stopMSTimer(-2) - t0)*1e-6 
    935958//                      Printf  "CylPolyRad fit time using res and eps and NO THREADING time = %g seconds\r\r",t0 
    936959                        break 
    937960                endif 
    938          
    939                 if(useRes && useConstr)         //no crsr, no eps 
    940                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /D=fitYw /C=constr /STRC=fs 
    941                         break 
    942                 endif 
    943                  
    944                 if(useRes)              //just res 
    945                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /D=fitYw /STRC=fs 
    946                         break 
    947                 endif 
    948                  
    949 /////   same as above, but all without useRes (no /STRC flag) 
    950                 if(useEps && useCursors && useConstr)           //do it all 
    951                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr 
    952                         break 
    953                 endif 
    954                  
    955                 if(useEps && useCursors)                //no constr 
    956                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /E=eps /D=fitYw 
    957                         break 
    958                 endif 
    959                  
    960                 if(useEps && useConstr)         //no crsr 
    961                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr 
    962                         break 
    963                 endif 
    964                  
    965                 if(useCursors && useConstr)             //no eps 
    966                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /D=fitYw /C=constr 
    967                         break 
    968                 endif 
    969                  
    970                 if(useCursors)          //no eps, no constr 
    971                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pcsr(A),pcsr(B)] /X=xw /W=sw /I=1 /D=fitYw 
    972                         break 
    973                 endif 
    974                  
    975                 if(useEps)              //no crsr, no constr 
    976                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /E=eps /D=fitYw 
    977                         break 
    978                 endif 
    979          
    980                 if(useConstr)           //no crsr, no eps 
    981                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /D=fitYw /C=constr 
     961                         
     962                 
     963                 
     964/////   same as above, but all without useResol (no /STRC flag) 
     965                if(useResiduals && useTextBox)          //resid+ text 
     966                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /R /NWOK 
     967                        break 
     968                endif 
     969                 
     970                if(useResiduals)                //resid 
     971                        FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /R /NWOK 
     972                        break 
     973                endif 
     974 
     975                 
     976                if(useTextBox)          //text 
     977                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /NWOK 
    982978                        break 
    983979                endif 
    984980                 
    985981                //just a plain vanilla fit 
    986                 FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw /X=xw /W=sw /I=1 /D=fitYw 
    987          
     982 
     983                FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /NWOK 
     984                 
    988985        while(0) 
    989986         
    990987//      t0 = (stopMSTimer(-2) - t0)*1e-6 
    991988//      Printf  "fit time = %g seconds\r\r",t0 
     989         
    992990         
    993991        // append the fit 
     
    10251023        print w_sigma 
    10261024        String resultStr="" 
    1027                  
     1025         
     1026        if(waveexists(W_sigma)) 
     1027                //append it to the table, if it's not already there 
     1028                CheckDisplayed/W=WrapperPanel#T0 W_sigma 
     1029                if(V_flag==0) 
     1030                        //not there, append it 
     1031                        AppendtoTable/W=wrapperPanel#T0 W_sigma 
     1032                else 
     1033                        //remove it, and put it back on to make sure it's the right one (do I need to do this?) 
     1034                        // -- not really, since any switch of the function menu takes W_Sigma off 
     1035                endif 
     1036        endif 
     1037         
    10281038        //now re-write the results 
    10291039        sprintf resultStr,"Chi^2 = %g  Sqrt(X^2/N) = %g",V_chisq,sqrt(V_chisq/V_Npnts) 
     
    10441054         
    10451055        if(yesReport) 
    1046                 String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //this is *hopefully* one wave 
     1056                String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              // this is *hopefully* one wave 
    10471057                String topGraph= WinName(0,1)   //this is the topmost graph 
    10481058         
    10491059                DoUpdate                //force an update of the graph before making a copy of it for the report 
    10501060 
    1051                 W_GenerateReport(funcStr,folderStr,$parStr,cw,yesSave,V_chisq,W_sigma,V_npnts,V_FitError,V_FitQuitReason,V_startRow,V_endRow,topGraph) 
     1061                if(useGenCurveFit)       
     1062                        W_GenerateReport(funcStr,folderStr,$parStr,cw,yesSave,V_chisq,W_sigma,V_npnts,V_FitError,V_FitQuitReason,pt1,pt2,topGraph) 
     1063                else 
     1064                        W_GenerateReport(funcStr,folderStr,$parStr,cw,yesSave,V_chisq,W_sigma,V_npnts,V_FitError,V_FitQuitReason,V_startRow,V_endRow,topGraph) 
     1065                endif 
    10521066        endif 
    10531067         
     
    11451159        // insert graphs 
    11461160        if(WaveExists(dataXw)) 
    1147                 Notebook $nb picture={$topGraph(0, 0, 400, 300), -5, 1}, text="\r" 
     1161//              Notebook $nb picture={$topGraph(0, 0, 400, 300), -5, 1}, text="\r" 
     1162                Notebook $nb scaling={50, 50}, picture={$topGraph(0, 0, 800, 600), -5, 1}, text="\r" 
    11481163        // 
    11491164        else            //must be 2D Gizmo 
    11501165                Execute "ExportGizmo Clip"                      //this ALWAYS is a PICT or BMP. Gizmo windows are different... 
    11511166                LoadPict/Q/O "Clipboard",tmp_Gizmo 
    1152                 Notebook $nb picture={tmp_Gizmo(0, 0, 400, 300), 0, 1}, text="\r" 
    1153         endif 
    1154         //Notebook Report picture={Table1, 0, 0}, text="\r" 
     1167                Notebook $nb picture={tmp_Gizmo(0, 0, 800, 600), 0, 1}, text="\r" 
     1168        endif 
    11551169         
    11561170        // show the top of the report 
     
    11601174        if(yesSave) 
    11611175                String nameStr=CleanupName(func,0) 
    1162                 nameStr = nameStr[0,8]  //shorten the name 
    11631176                nameStr += "_"+dataname 
     1177                nameStr = ReplaceString("Smeared",nameStr,"Sm_")                //if Smeared function, shorten the name 
    11641178                //make sure the name is no more than 31 characters 
    11651179                namestr = namestr[0,30]         //if shorter than 31, this will NOT pad to 31 characters 
    11661180                Print "file saved as ",nameStr 
    11671181                SaveNotebook /O/P=home/S=2 $nb as nameStr 
    1168                 //save the graph separately as a PICT file, 2x screen 
     1182                //save the graph separately as a PNG file, 2x screen 
    11691183                pictStr += nameStr 
    11701184                pictStr = pictStr[0,28]         //need a shorter name - why? 
     
    11741188///             SavePICT /E=-5/O/P=home /I/W=(0,0,3,3) as pictStr 
    11751189                if(WaveExists(dataXw)) 
    1176                         SavePICT /E=-5/O/P=home/WIN=$topGraph /W=(0,0,400,300) as pictStr 
     1190                        SavePICT /E=-5/O/P=home/WIN=$topGraph /W=(0,0,800,600) as pictStr 
    11771191                else 
    1178                         Execute "ExportGizmo /P=home as \""+pictStr+"\"" 
     1192                        Execute "ExportGizmo /P=home as \""+pictStr+"\""                //this won't be of very high quality 
    11791193                        //SavePICT /E=-5/O/P=home/WIN=$topGraph /W=(0,0,400,300) as pictStr 
    11801194                endif 
     
    12251239        endswitch 
    12261240end 
     1241 
     1242Function UseInfoTextBoxCheckProc(cba) : CheckBoxControl 
     1243        STRUCT WMCheckboxAction &cba 
     1244 
     1245        switch( cba.eventCode ) 
     1246                case 2: // mouse up 
     1247                        Variable checked = cba.checked 
     1248                        if(checked) 
     1249                                //print "checked, use textBox in the next fit" 
     1250                        else 
     1251                                //print "unchecked, ask to remove TextBox from the graph" 
     1252                                ControlInfo/W=WrapperPanel popup_0 
     1253                                RemoveTextBox(S_value) 
     1254                        endif 
     1255                        break 
     1256        endswitch 
     1257 
     1258        return 0 
     1259End 
     1260 
     1261//does not report an error if the text box is not there 
     1262// -- so I'll just be lazy and not check to see if it's there 
     1263// 
     1264Function RemoveTextBox(folderStr) 
     1265        String folderStr 
     1266         
     1267        DoAlert 1,"Remove the TextBox from the graph?" 
     1268        if(V_flag == 1) 
     1269                String str = "CF_"+folderStr+"_i" 
     1270                TextBox/K/N=$str 
     1271        endif 
     1272        return(0) 
     1273End 
     1274 
     1275Function UseResidualsCheckProc(cba) : CheckBoxControl 
     1276        STRUCT WMCheckboxAction &cba 
     1277 
     1278        switch( cba.eventCode ) 
     1279                case 2: // mouse up 
     1280                        Variable checked = cba.checked 
     1281                        if(checked) 
     1282                                //print "checked, use them in the next fit" 
     1283                        else 
     1284                                //print "unchecked, ask to remove residuals from the graph" 
     1285                                ControlInfo/W=WrapperPanel popup_0 
     1286                                RemoveResiduals(S_value) 
     1287                        endif 
     1288                        break 
     1289        endswitch 
     1290 
     1291        return 0 
     1292End 
     1293 
     1294// the default name from the /R flag is "Res_" + yWaveName 
     1295// 
     1296// better to find the wave that starts with "Res_" and remove that one in case the 
     1297// wave names get too long 
     1298// 
     1299// the difficulty now is that the residual wave ends up in root: and not with the data.... 
     1300// -- not really a problem, but adds to clutter 
     1301Function RemoveResiduals(folderStr) 
     1302        String folderStr 
     1303         
     1304        String list="",topWin="" 
     1305        Variable num,ii 
     1306        String str 
     1307 
     1308        DoAlert 1,"Remove the residuals from the graph?" 
     1309        if(V_flag == 1) 
     1310//              String topGraph= WinName(0,1)   //this is the topmost graph 
     1311                list=TraceNameList("", ";", 1 )         //"" as first parameter == look on the target graph 
     1312                num=ItemsInList(list) 
     1313                 
     1314                for(ii=0;ii<num;ii+=1) 
     1315                        str = StringFromList(ii, list ,";") 
     1316                        if(strsearch(str, "Res_", 0) != -1) 
     1317                                RemoveFromGraph $str 
     1318                        endif 
     1319                endfor 
     1320         
     1321                SetDataFolder root: 
     1322        endif 
     1323         
     1324        return(0) 
     1325End 
    12271326 
    12281327Function Toggle2DControlsCheckProc(cba) : CheckBoxControl 
     
    12441343                                Button button_2D_0,disable=0            //visible again, and enabled 
    12451344                                Button button_2D_1,disable=0 
     1345                                 
     1346                                CheckBox check_6,disable=1                      //info box and residual check, remove these from view 
     1347                                CheckBox check_7,disable=1 
    12461348                        else 
    12471349                                //print "unchecked, change them back to 1D" 
     
    12531355                                Button button_2D_0,disable=3    //hide the extra 2D buttons, and disable 
    12541356                                Button button_2D_1,disable=3 
     1357                                 
     1358                                CheckBox check_6,disable=0                      //info box and residual check, bring them back 
     1359                                CheckBox check_7,disable=0 
    12551360                        endif 
    12561361                        break 
     
    12741379                                CheckBox check_0,value=0 
    12751380                                return(0) 
     1381                        endif 
     1382                         
     1383                        //if in 2D mode, just exit 
     1384                        ControlInfo/W=WrapperPanel check_3 
     1385                        if(V_Value == 1) 
     1386                                return (0) 
    12761387                        endif 
    12771388                                 
     
    13071418                                endif 
    13081419 
    1309                                 HideInfo 
     1420                                HideInfo/W=$topGraph 
    13101421                                Cursor/K A 
    13111422                                Cursor/K B 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/SA_includes_v410.ipf

    r693 r766  
    2828 
    2929#include "DataSetHandling"                                      //added Nov 2009 AJJ 
     30#include "Smear_2D"                                             //for 2D resolution smearing, May 2010 
    3031 
    3132Menu "SANS Models" 
    3233        "Fit Manager", Init_WrapperPanel() 
    3334        "Load Model Functions",Execute/P "INSERTINCLUDE \"SANSModelPicker_v40\"";Execute/P "COMPILEPROCEDURES ";Execute/P "ModelPicker_Panel()" 
    34         "Load and Plot Manager", Show_Plot_Manager() 
    35         "Freeze Model" 
    36         "Write Model Data" 
    37         "ReWrite Experimental Data",ReWrite1DData() 
    38         "1D Arithmetic Panel",MakeDAPanel() 
     35        "-" 
     36        Submenu "1D Utilities" 
     37                "Load and Plot Manager", Show_Plot_Manager() 
     38                "Freeze Model" 
     39                "Write Model Data" 
     40                "ReWrite Experimental Data",MakeDMPanel()               //,ReWrite1DData()      // SRK SEP10 
     41                "1D Arithmetic Panel",MakeDAPanel() 
     42                "ReBin 1D Data",OpenRebin() 
     43        end 
    3944        "-" 
    4045        Submenu "Packages" 
     
    4348                "Simple Global Fitting",Execute/P "INSERTINCLUDE \"GlobalFit2_NCNR_v40\"";Execute/P "INSERTINCLUDE \"SimpleGlobalFit_NCNR_v40\"";Execute/P "COMPILEPROCEDURES ";Execute/P "Init_SimpleGlobalFit()" 
    4449                "Determine Invariant",Execute/P "INSERTINCLUDE \"Invariant_v40\"";Execute/P "COMPILEPROCEDURES ";Execute/P "Make_Invariant_Panel()" 
    45                 "Do Linear Fits",Execute/P "INSERTINCLUDE \"LinearizedFits_v40\"";Execute/P "COMPILEPROCEDURES ";Execute/P "A_OpenFitPanel()" 
     50                "Do Linear Fits",Execute/P "INSERTINCLUDE \"LinearizedFits_v40\"";Execute/P "COMPILEPROCEDURES ";Execute/P "OpenFitPanel()" 
    4651                GenOpFlagEnable()+"Genetic Optimization Enabled", Init_GenOp() 
    4752                GenOpFlagDisable()+"Genetic Optimization Disabled", UnSet_GenOp() 
  • sans/Release/trunk/NCNR_User_Procedures/Analysis/WriteModelData_v40.ipf

    r607 r766  
    178178                NVAR/Z dQv = USANS_dQv 
    179179                if(NVAR_Exists(dQv) == 0) 
     180                        SetDataFolder root: 
    180181                        Abort "It's USANS data, and I don't know what the slit height is." 
    181182                endif 
     
    193194                PathInfo/S catPathName 
    194195                fullPath = DoSaveFileDialog("Save data as",fname=baseStr+".txt") 
     196                Print fullPath 
    195197                If(cmpstr(fullPath,"")==0) 
    196198                        //user cancel, don't write out a file 
     
    200202                //Print "dialog fullpath = ",fullpath 
    201203        Endif 
    202          
     204 
    203205        Open refnum as fullpath 
    204206         
  • sans/Release/trunk/NCNR_User_Procedures/Common/DataSetHandling.ipf

    r681 r766  
    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// Function batchGrasp26ColConvert() 
     9// 
     10// these both need a real interface, and a way to better define the name of the 
     11// converted output file. And some retention of the header would be nice too... 
     12// 
     13 
    214 
    315// Functions and interfaces to manage datasets now that they are in data folders 
     
    5062        Button Rename_button,title="Rename",pos={75,170},size={150,20} 
    5163        Button Rename_button,proc=DM_RenameProc 
    52         Button  Duplicate_button,title="Duplicate",pos={275,170},size={150,20} 
     64        Button Duplicate_button,title="Duplicate",pos={275,170},size={150,20} 
    5365        Button Duplicate_button,proc=DM_DuplicateProc 
    5466 
     
    5971        Button SaveAsXML_button,proc=DMSaveAsXMLproc     
    6072 
    61         Button SaveAs6col_button,title="Save as NIST 6 column",pos={275,230},size={150,20} 
     73        Button SaveAs6col_button,title="Save as NIST 6 column",pos={275,230},size={160,20} 
    6274        Button SaveAs6col_button,proc=DMSaveAs6colproc   
    6375         
    64         Button BatchConvertData_button,title="Convert Format of 1D Data Files",pos={75,290},size={350,20} 
     76        Button BatchConvertData_button,title="Batch Convert Format of 1D Data Files",pos={75,290},size={350,20} 
    6577        Button BatchConvertData_button,proc=DMBatchConvertProc 
    6678                         
     
    7688        else 
    7789                //fake call to popup 
    78                 STRUCT WMPopupAction pa 
    79                 pa.win = "DataManagementPanel" 
    80                 pa.ctrlName = "DS_popup" 
    81                 pa.eventCode = 2 
    82                 DM_PopupProc(pa) 
    83         endif 
    84  
    85 End 
     90                //// -- can't use STRUCT in a Proc, only a function 
     91//              STRUCT WMPopupAction pa 
     92//              pa.win = "DataManagementPanel" 
     93//              pa.ctrlName = "DS_popup" 
     94//              pa.eventCode = 2 
     95//              DM_PopupProc(pa) 
     96        endif 
     97 
     98End 
     99 
     100Function DMSaveAs6colproc(ba) : ButtonControl 
     101        STRUCT WMButtonAction &ba 
     102         
     103        switch( ba.eventCode) 
     104                case 2: 
     105                        ControlInfo/W=DataManagementPanel DS_popup 
     106                        String folderName = S_Value 
     107                        fReWrite1DData(folderName,"tab","CRLF") 
     108                        break 
     109        endswitch 
     110         
     111        return 0 
     112End 
     113 
     114Function DMSaveAsXMLproc(ba) : ButtonControl 
     115        STRUCT WMButtonAction &ba 
     116         
     117        switch( ba.eventCode) 
     118                case 2: 
     119                        ControlInfo/W=DataManagementPanel DS_popup 
     120                        String folderName = S_Value 
     121                        ReWrite1DXMLData(folderName) 
     122                        break 
     123        endswitch 
     124         
     125        return 0 
     126End 
     127 
     128 
     129//Function SaveDataSetToFile(folderName) 
     130//      String folderName 
     131// 
     132//      String protoStr = "" 
     133//      //Check for history string in folder 
     134//      //If it doesn't exist then  
     135// 
     136//      //Do saving of data file. 
     137//       
     138//      NVAR gXML_Write = root:Packages:NIST:gXML_Write  
     139// 
     140//      if (gXML_Write == 1) 
     141//              ReWrite1DXMLData(folderName) 
     142//      else 
     143//              fReWrite1DData(folderName,"tab","CRLF") 
     144//      endif 
     145// 
     146//      //Include history string to record what was done? 
     147// 
     148//End 
     149 
     150 
     151 
     152 
     153 
    86154 
    87155Function DMBatchConvertProc(ba) : ButtonControl 
     
    90158        switch( ba.eventCode) 
    91159                case 2: 
    92                         Execute "MakeBatchConvertPanel()" 
     160                        Execute "MakeNCNRBatchConvertPanel()" 
    93161                        break 
    94162        endswitch 
     
    351419        Make/O/T/N=1 filewave="" 
    352420        Make/O/N=1 selWave=0 
    353         Variable/G ind=0 
     421        Variable/G ind=0,gRadioVal=1 
    354422        SetDataFolder root: 
    355423End 
     
    358426        PauseUpdate; Silent 1           // building window... 
    359427        NewPanel /W=(658,347,1018,737)/N=NCNRBatchConvertPanel/K=2 as "Batch Convert 1D Data" 
     428//      NewPanel /W=(658,347,1018,737)/N=NCNRBatchConvertPanel as "Batch Convert 1D Data" 
    360429        ModifyPanel cbRGB=(40000,50000,32896) 
    361430        ModifyPanel fixedSize=1 
     
    367436        Button button7,pos={238,20},size={100,20},proc=NCNRBatchConvertNewFolder,title="New Folder" 
    368437        Button button7,help={"Select a new data folder"} 
    369         TitleBox msg0,pos={238,140},size={100,30},title="\JCShift-click to\rselect multiple files" 
     438        TitleBox msg0,pos={238,160},size={100,30},title="\JCShift-click to\rselect multiple files" 
    370439        TitleBox msg0,frame=0,fixedSize=1 
    371440         
    372441        GroupBox filterGroup,pos={13,200},size={206,60},title="Filter list by input file type" 
    373          
    374         Button button8,pos={238,76},size={100,20},proc=NCNRBatchConvertHelpProc,title="Help" 
    375         Button button9,pos={238,48},size={100,20},proc=NCNRBatchConvertRefresh,title="Refresh List" 
    376         Button button0,pos={238,106},size={100,20},proc=NCNRBatchConvertDone,title="Done" 
     442        CheckBox filterCheck_1,pos={24,220},size={36,14},title="XML",value= 1,mode=1, proc=BC_filterCheckProc 
     443        CheckBox filterCheck_2,pos={24,239},size={69,14},title="ABS or AVE",value= 0,mode=1, proc=BC_filterCheckProc 
     444        CheckBox filterCheck_3,pos={100,220},size={69,14},title="none",value= 0,mode=1, proc=BC_filterCheckProc 
     445         
     446        Button button8,pos={238,75},size={100,20},proc=NCNRBatchConvertHelpProc,title="Help" 
     447        Button button9,pos={238,47},size={100,20},proc=NCNRBatchConvertRefresh,title="Refresh List" 
     448        Button button10,pos={238,102},size={100,20},proc=NCNRBatchConvertSelectAll,title="Select All" 
     449        Button button0,pos={238,130},size={100,20},proc=NCNRBatchConvertDone,title="Done" 
    377450 
    378451        GroupBox outputGroup,pos={13,270},size={206,60},title="Output File Type" 
    379  
     452        CheckBox outputCheck_1,pos={24,289},size={36,14},title="XML",value= 0,mode=1, proc=BC_outputCheckProc 
     453        CheckBox outputCheck_2,pos={24,309},size={69,14},title="ABS or AVE",value= 1,mode=1, proc=BC_outputCheckProc 
    380454         
    381455        Button button6,pos={13,350},size={206,20},proc=NCNRBatchConvertFiles,title="Convert File(s)" 
     
    386460End 
    387461 
     462Function BC_filterCheckProc(ctrlName,checked) : CheckBoxControl 
     463        String ctrlName 
     464        Variable checked 
     465 
     466        NVAR gRadioVal= root:Packages:NIST:BatchConvert:gRadioVal 
     467         
     468        strswitch (ctrlName) 
     469                case "filterCheck_1": 
     470                        gRadioVal= 1 
     471                        break 
     472                case "filterCheck_2": 
     473                        gRadioVal= 2 
     474                        break 
     475                case "filterCheck_3": 
     476                        gRadioVal= 3 
     477                        break 
     478        endswitch 
     479        CheckBox filterCheck_1,value= gRadioVal==1 
     480        CheckBox filterCheck_2,value= gRadioVal==2 
     481        CheckBox filterCheck_3,value= gRadioVal==3 
     482         
     483        NCNRBatchConvertGetList() 
     484          
     485        return(0) 
     486End 
     487 
     488Function BC_outputCheckProc(ctrlName,checked) : CheckBoxControl 
     489        String ctrlName 
     490        Variable checked 
     491 
     492        if(cmpstr("outputCheck_1",ctrlName)==0) 
     493                CheckBox outputCheck_1,value=checked 
     494                CheckBox outputCheck_2,value=!checked 
     495        else 
     496                CheckBox outputCheck_1,value=!checked 
     497                CheckBox outputCheck_2,value=checked 
     498        endif 
     499         
     500        return(0) 
     501End 
    388502 
    389503Function NCNRBatchConvertFiles(ba) : ButtonControl 
    390504                STRUCT WMButtonAction &ba 
    391505                 
     506                 
    392507                switch (ba.eventCode) 
    393508                        case 2: 
     509                         
     510                                //check the input/output as best I can (none may be the input filter) 
     511                                Variable inputType,outputType=1 
     512                                NVAR gRadioVal= root:Packages:NIST:BatchConvert:gRadioVal 
     513                                inputType = gRadioVal 
     514                                ControlInfo outputCheck_1 
     515                                if(V_value==1) 
     516                                        outputType = 1          //xml 
     517                                else 
     518                                        outputType = 2          //6-col 
     519                                endif 
     520                                 
     521                                if(inputType==outputType) 
     522                                        DoAlert 0,"Input and output types are the same. Nothing will be converted" 
     523                                        return(0) 
     524                                endif 
     525                         
     526                         
     527                                        // input and output are different, proceed 
     528 
    394529                                Wave/T fileWave=$"root:Packages:NIST:BatchConvert:fileWave" 
    395                                 Wave sel=$"root:BatchConvert:NIST:BatchConvert:selWave" 
     530                                Wave sel=$"root:Packages:NIST:BatchConvert:selWave" 
     531                                 
     532                                String fname="",pathStr="",newFileName="" 
     533                                Variable ii,num 
     534                                PathInfo catPathName                    //this is where the files are 
     535                                pathStr=S_path 
     536                                                         
     537                                // process the selected items 
     538                                num=numpnts(sel) 
     539                                ii=0 
     540                                do 
     541                                        if(sel[ii] == 1) 
     542                                                fname=pathStr + fileWave[ii] 
     543                                                 
     544                                                if(outputType == 1) 
     545                                                        convertNISTtoNISTXML(fname) 
     546                                                endif 
     547                                                 
     548                                                if(outputType == 2) 
     549                                                        convertNISTXMLtoNIST6Col(fname) 
     550                                                endif 
     551                                        endif 
     552                                        ii+=1 
     553                                while(ii<num) 
    396554                                 
    397555                                break 
     
    401559End 
    402560 
     561Function NCNRBatchConvertSelectAll(ba) : ButtonControl 
     562        STRUCT WMButtonAction &ba 
     563         
     564        switch (ba.eventcode) 
     565                case 2: 
     566                        Wave sel=$"root:Packages:NIST:BatchConvert:selWave" 
     567                        sel = 1 
     568                        break 
     569        endswitch 
     570 
     571End 
    403572 
    404573Function NCNRBatchConvertNewFolder(ba) : ButtonControl 
     
    414583End 
    415584 
     585 
     586// filter is a bit harsh - will need to soften this by presenting an option to enter the suffix 
     587// 
    416588Function NCNRBatchConvertGetList() 
    417589 
     
    422594        Endif 
    423595         
    424         String newList = A_ReducedDataFileList("") 
     596        String newList = A_ReducedDataFileList(""),tmpList="" 
    425597        Variable num 
     598         
     599        NVAR gRadioVal= root:Packages:NIST:BatchConvert:gRadioVal 
     600        ControlInfo filterCheck_1 
     601        if(gRadioVal == 1) 
     602                //keep XML data 
     603                tmpList = ListMatch(newList, "*.ABSx" ,";")             //note that ListMatch is not case-sensitive 
     604                tmpList += ListMatch(newList, "*.AVEx" ,";") 
     605                tmpList += ListMatch(newList, "*.CORx" ,";") 
     606                tmpList += ListMatch(newList, "*.xml" ,";") 
     607        else 
     608                if(gRadioVal ==2) 
     609                        //keep ave, abs data 
     610                        tmpList = ListMatch(newList, "*.ABS" ,";") 
     611                        tmpList += ListMatch(newList, "*.AVE" ,";") 
     612                else 
     613                        //return everything 
     614                        tmpList = newList 
     615                endif 
     616        endif 
     617        newList = tmpList 
    426618         
    427619        num=ItemsInList(newlist,";") 
     
    432624        fileWave = StringFromList(p,newlist,";") 
    433625        Sort filewave,filewave 
     626         
     627        return 0 
    434628End 
    435629 
     
    583777        Checkbox DAPlot_log_cb, title="Log I(q)", pos={20,410},value=0 
    584778        Checkbox DAPlot_log_cb, proc=DALogLinIProc 
     779        Checkbox DAPlot_lin_cb, title="High Q Linear", pos={100,410},value=0 
     780        Checkbox DAPlot_lin_cb, proc=DAHighQLinProc 
    585781         
    586782End 
     
    750946                        SetDrawEnv fname="Times", fsize=22, fstyle= 2,textrgb= (0,0,0)                                   
    751947                        DrawText 105,30,")" 
    752                         break            
     948                        break    
    753949        endswitch 
    754950 
     
    10741270                         
    10751271                        ModifyGraph/W=DAPlotPanel#DAPlot log(left)=cba.checked 
    1076                  
    1077         endswitch 
    1078  
    1079  
    1080 End 
     1272                        ModifyGraph/W=DAPlotPanel#DAPlot log(bottom)=cba.checked 
     1273                        ModifyGraph/W=DAPlotPanel#DAPlot zero(left)=0 
     1274                        SetAxis/A/W=DAPlotPanel#DAPlot 
     1275                         
     1276                        if(cba.checked) 
     1277                                Checkbox DAPlot_lin_cb,value=0          //uncheck lin 
     1278                        endif 
     1279        endswitch 
     1280 
     1281 
     1282End 
     1283 
     1284Function DAHighQLinProc(cba) : CheckBoxControl 
     1285        STRUCT WMCheckBoxAction &cba 
     1286 
     1287        switch(cba.eventcode) 
     1288                case 2: 
     1289                        if(cba.checked) 
     1290                                ModifyGraph/W=DAPlotPanel#DAPlot log=0,zero(left)=1 
     1291                                SetAxis/W=DAPlotPanel#DAPlot left -0.1,0.1 
     1292                                SetAxis/W=DAPlotPanel#DAPlot bottom 0.1,* 
     1293                                SetAxis/W=DAPlotPanel#DAPlot left -0.02,0.02 
     1294                                 
     1295                                Checkbox DAPlot_log_cb,value=0          //uncheck the log 
     1296                        endif 
     1297        endswitch 
     1298 
     1299 
     1300End 
     1301 
    10811302 
    10821303Function DACursorButtonProc(ba) : ButtonControl 
     
    12901511        //If we are here, the folder (now) exists and the user has agreed to overwrite 
    12911512        //either in the function call or from the alert. 
     1513         
     1514        // here, GetIndexedObjectName copies all of the waves 
     1515        index = 0 
    12921516        do 
    12931517                objName = GetIndexedObjName(basestr,1,index) 
     
    12991523                index+=1 
    13001524        while(1) 
     1525 
     1526// -- for USANS data, we need the slit height. copy all of the "USANS_*" variables 
     1527// may need to augment this for other situations 
     1528        index = 0 
     1529        do 
     1530                objName = GetIndexedObjName(basestr,2,index) 
     1531                if (strlen(objName) == 0) 
     1532                        break 
     1533                endif 
     1534                if(stringmatch(objName,"USANS*") == 1) 
     1535                        objname = ":"+basestr+":"+objname 
     1536                        NVAR tmp = $objName 
     1537                        Variable/G $(ReplaceString(basestr,objName,newName))= tmp 
     1538                endif 
     1539                index+=1 
     1540        while(1) 
     1541         
    13011542 
    13021543        SetDataFolder root: 
     
    14891730        endif 
    14901731 
    1491 //      NISTSave1DData(folderName,protoStr) 
    1492  
    14931732        //Include history string to record what was done? 
    14941733 
     
    14971736 
    14981737 
    1499 //////////////////// Write data functions //////////////////////////////////// 
    1500 // AJJ Nov 2009 - should move towards unified writer 
     1738// still need to get the header information, and possibly the SASprocessnote from the XML load into the 6-column header 
    15011739// 
    1502 ////////////////////////////////////////////////////////////////////////// 
    1503  
    1504 //Function NISTSave1DData(folderPath,protoStr) 
    1505 //      String folderPath, protoStr 
    1506 //       
    1507 //      //protoStr is string that will be written to either  
    1508 //       
     1740// start by looking in:  
     1741//      String xmlReaderFolder = "root:Packages:CS_XMLreader:" 
     1742// for Title and Title_folder strings -> then the metadata (but the processnote is still not there 
    15091743// 
    1510 //End 
     1744// may need to get it directly using the filename 
     1745Function  convertNISTXMLtoNIST6Col(fname) 
     1746        String fname 
     1747 
     1748        String list, item,path 
     1749        Variable num,ii 
     1750         
     1751        //load the XML 
     1752         
     1753        LoadNISTXMLData(fname,"",0,0)           //no plot, no force overwrite 
     1754//      Execute "A_LoadOneDDataWithName(\""+fname+"\",0)"               //won't plot 
     1755 
     1756        // then rewrite what is in the data folder that was just loaded 
     1757        String basestr = ParseFilePath(0, fname, ":", 1, 0) 
     1758        baseStr = CleanupName(baseStr,0) 
     1759        print fname 
     1760        print basestr 
     1761 
     1762        fReWrite1DData_noPrompt(baseStr,"tab","CR") 
     1763 
     1764        return(0) 
     1765End 
     1766 
     1767 
     1768///////// SRK - VERY SIMPLE batch converter 
     1769// no header information is preserved 
     1770// file names are partially preserved 
     1771// 
     1772 
     1773/// to use this: 
     1774// -open the Plot Manager and set the path 
     1775// -run this function 
     1776// 
     1777// it doesn't matter if the XML ouput flag is set - this overrides. 
     1778Function batchXML26ColConvert() 
     1779 
     1780        String list, item,path,fname 
     1781        Variable num,ii 
     1782         
     1783        PathInfo CatPathName 
     1784        path = S_Path 
     1785 
     1786        list = A_ReducedDataFileList("") 
     1787        num = itemsInList(list) 
     1788        Print num 
     1789        for(ii=0;ii<num;ii+=1) 
     1790                item = StringFromList(ii, list ,";") 
     1791                fname=path + item 
     1792                Execute "A_LoadOneDDataWithName(\""+fname+"\",0)"               //won't plot 
     1793//              easier to load all, then write out, since the name will be changed 
     1794        endfor 
     1795         
     1796         
     1797        list = DM_DataSetPopupList() 
     1798 
     1799        num = itemsInList(list) 
     1800        Print num 
     1801        for(ii=0;ii<num;ii+=1) 
     1802                item = StringFromList(ii, list ,";") 
     1803                fReWrite1DData_noPrompt(item,"tab","CR") 
     1804        endfor 
     1805         
     1806End 
     1807 
     1808///////// SRK - VERY SIMPLE batch converter 
     1809// NO header information is preserved 
     1810// file names are partially preserved 
     1811// 
     1812 
     1813/// to use this: 
     1814// -open the Plot Manager and set the path 
     1815// -run this function 
     1816// 
     1817// it doesn't matter if the XML ouput flag is set - this overrides. 
     1818// 
     1819// The GRASP output data is 5-column Q-I-errI-sigQ-nCells 
     1820// which gets read in as q-i-s-ism-fit_ism (as if it was some wierd USANS data format) 
     1821// -- so fake the output... 
     1822// 
     1823Function batchGrasp26ColConvert() 
     1824 
     1825        String list, item,path,fname 
     1826        Variable num,ii,npt 
     1827         
     1828        PathInfo CatPathName 
     1829        path = S_Path 
     1830 
     1831        list = A_ReducedDataFileList("") 
     1832        num = itemsInList(list) 
     1833        Print num 
     1834        for(ii=0;ii<num;ii+=1) 
     1835                item = StringFromList(ii, list ,";") 
     1836                fname=path + item 
     1837                Execute "A_LoadOneDDataWithName(\""+fname+"\",0)"               //won't plot 
     1838//              easier to load all, then write out, since the name will be changed 
     1839        endfor 
     1840         
     1841         
     1842        list = DM_DataSetPopupList() 
     1843 
     1844        num = itemsInList(list) 
     1845         
     1846        Print num 
     1847        for(ii=0;ii<num;ii+=1) 
     1848                item = StringFromList(ii, list ,";") 
     1849                 
     1850                // fake the 6-column NIST data structure 
     1851                WAVE qw = $("root:"+item+":"+item+"_q") 
     1852                npt = numpnts(qw) 
     1853                Make/O/D/N=(npt,4) $("root:"+item+":"+item+"_res") 
     1854                WAVE res = $("root:"+item+":"+item+"_res") 
     1855                WAVE sigQ = $("root:"+item+":"+item+"_ism") 
     1856                res[][0] = sigQ[p]      // sigQ 
     1857                res[][1] = qw[p]                // qBar ~ q 
     1858                res[][2] = 1            //shadow 
     1859                res[][3] = qw[p]                // q 
     1860                 
     1861                 
     1862                fReWrite1DData_noPrompt(item,"tab","CR") 
     1863        endfor 
     1864         
     1865End 
     1866 
     1867// quick version (copied from fReWrite1DData() that NEVER asks for a new fileName 
     1868// - and right now, always expect 6-column data, either SANS or USANS (re-writes -dQv) 
     1869// - AJJ Nov 2009 : better make sure we always fake 6 columns on reading then.... 
     1870Function fReWrite1DData_noPrompt(folderStr,delim,term) 
     1871        String folderStr,delim,term 
     1872         
     1873        String formatStr="",fullpath="" 
     1874        Variable refnum,dialog=1 
     1875         
     1876        String dataSetFolderParent,basestr 
     1877         
     1878        //setup delimeter and terminator choices 
     1879        If(cmpstr(delim,"tab")==0) 
     1880                //tab-delimeted 
     1881                formatStr="%15.8g\t%15.8g\t%15.8g\t%15.8g\t%15.8g\t%15.8g" 
     1882        else 
     1883                //use 3 spaces 
     1884                formatStr="%15.8g   %15.8g   %15.8g   %15.8g   %15.8g   %15.8g" 
     1885        Endif 
     1886        If(cmpstr(term,"CR")==0) 
     1887                formatStr += "\r" 
     1888        Endif 
     1889        If(cmpstr(term,"LF")==0) 
     1890                formatStr += "\n" 
     1891        Endif 
     1892        If(cmpstr(term,"CRLF")==0) 
     1893                formatStr += "\r\n" 
     1894        Endif 
     1895         
     1896        //Abuse ParseFilePath to get path without folder name 
     1897        dataSetFolderParent = ParseFilePath(1,folderStr,":",1,0) 
     1898        //Abuse ParseFilePath to get basestr 
     1899        basestr = ParseFilePath(0,folderStr,":",1,0) 
     1900         
     1901        //make sure the waves exist 
     1902        SetDataFolder $(dataSetFolderParent+basestr) 
     1903        WAVE/Z qw = $(baseStr+"_q") 
     1904        WAVE/Z iw = $(baseStr+"_i") 
     1905        WAVE/Z sw = $(baseStr+"_s") 
     1906        WAVE/Z resw = $(baseStr+"_res") 
     1907         
     1908        if(WaveExists(qw) == 0) 
     1909                Abort "q is missing" 
     1910        endif 
     1911        if(WaveExists(iw) == 0) 
     1912                Abort "i is missing" 
     1913        endif 
     1914        if(WaveExists(sw) == 0) 
     1915                Abort "s is missing" 
     1916        endif 
     1917        if(WaveExists(resw) == 0) 
     1918                Abort "Resolution information is missing." 
     1919        endif 
     1920         
     1921        Duplicate/O qw qbar,sigQ,fs 
     1922        if(dimsize(resW,1) > 4) 
     1923                //it's USANS put -dQv back in the last 3 columns 
     1924                NVAR/Z dQv = USANS_dQv 
     1925                if(NVAR_Exists(dQv) == 0) 
     1926                        Abort "It's USANS data, and I don't know what the slit height is." 
     1927                endif 
     1928                sigQ = -dQv 
     1929                qbar = -dQv 
     1930                fs = -dQv 
     1931        else 
     1932                //it's SANS 
     1933                sigQ = resw[p][0] 
     1934                qbar = resw[p][1] 
     1935                fs = resw[p][2] 
     1936        endif 
     1937         
     1938        dialog=0 
     1939        if(dialog) 
     1940                PathInfo/S catPathName 
     1941//              fullPath = DoSaveFileDialog("Save data as",fname=baseStr+".txt") 
     1942                fullPath = DoSaveFileDialog("Save data as",fname=baseStr[0,strlen(BaseStr)-1]) 
     1943                Print fullPath 
     1944                If(cmpstr(fullPath,"")==0) 
     1945                        //user cancel, don't write out a file 
     1946                        Close/A 
     1947                        Abort "no data file was written" 
     1948                Endif 
     1949                //Print "dialog fullpath = ",fullpath 
     1950        Endif 
     1951        PathInfo catPathName 
     1952        fullPath = S_Path + baseStr[0,strlen(BaseStr)-1] 
     1953 
     1954        Open refnum as fullpath 
     1955         
     1956        fprintf refnum,"Modified data written from folder %s on %s\r\n",baseStr,(date()+" "+time()) 
     1957        wfprintf refnum,formatStr,qw,iw,sw,sigQ,qbar,fs 
     1958        Close refnum 
     1959         
     1960        KillWaves/Z sigQ,qbar,fs 
     1961         
     1962        SetDataFolder root: 
     1963        return(0) 
     1964End 
     1965 
     1966///////////end SRK 
     1967 
     1968 
     1969 
     1970///// SRK  
     1971///// Rebinning 
     1972 
     1973 
     1974// main entry procedure for subtraction panel 
     1975// re-initializes necessary folders and waves 
     1976Proc OpenRebin() 
     1977        DoWindow/F Rebin_Panel 
     1978        if(V_Flag==0) 
     1979                InitRebin() 
     1980                Rebin()         //the panel 
     1981//              Plot_Sub1D()            //the graph 
     1982        endif 
     1983         
     1984End 
     1985 
     1986Function InitRebin() 
     1987 
     1988        NewDataFolder/O/S root:Packages:NIST:Rebin 
     1989//      Variable/G gPtsBeg1=0 
     1990//      Variable/G gPtsEnd1=0 
     1991        Variable/G binning=2 
     1992 
     1993        SetDataFolder root: 
     1994End 
     1995 
     1996Proc Rebin()  
     1997 
     1998        PauseUpdate; Silent 1           // building window... 
     1999        NewPanel /W=(658,347,920,562)/K=1 
     2000        DoWindow/C Rebin_Panel 
     2001        ModifyPanel cbRGB=(65535,48662,45086),fixedSize=1 
     2002        Button button0,pos={161,178},size={50,20},proc=Rebin_Done,title="Done" 
     2003        PopupMenu popup0,pos={11,50},size={211,20},title="Data in Memory" 
     2004        PopupMenu popup0,mode=1,value= #"A_OneDDataInMemory()" 
     2005        Button button1,pos={60,14},size={150,20},proc=Rebin_Load,title="Load Data From File" 
     2006        Button button2,pos={118,84},size={100,20},proc=Rebin_Append,title="Append Data" 
     2007        Button button3,pos={11,84},size={80,20},proc=Rebin_newGraph,title="New Graph" 
     2008        Button button6,pos={11,138},size={70,20},proc=Rebin_by,title="Rebin by" 
     2009        Button button_8,pos={160,138},size={90,20},proc=SaveResultBin,title="Save Result" 
     2010         
     2011        SetVariable end_3,pos={97,140},size={40,14},title=" ",Font="Arial",fsize=10 
     2012        SetVariable end_3,limits={1,10,1},value= root:Packages:NIST:Rebin:binning 
     2013 
     2014EndMacro 
     2015 
     2016 
     2017Function Rebin_Done(ctrlName) : ButtonControl 
     2018        String ctrlName 
     2019        DoWindow/K Rebin_Panel 
     2020End 
     2021 
     2022Proc Rebin_Load(ctrlName) : ButtonControl 
     2023        String ctrlName 
     2024         
     2025        A_LoadOneDData() 
     2026        ControlUpdate/W=Rebin_Panel popup0 
     2027End 
     2028 
     2029Function Rebin_Append(ctrlName) : ButtonControl 
     2030        String ctrlName 
     2031        //appends data set from memory   /// joindre 
     2032        String iStr,qStr,eStr 
     2033        Variable rr,gg,bb 
     2034         
     2035        //get the current selection 
     2036        ControlInfo popup0 
     2037        if(strlen(S_Value)==0) 
     2038                Abort "You must load data from a file into memory before appending the data" 
     2039        Endif 
     2040         
     2041        A_PM_doAppend(S_Value) 
     2042         
     2043        DoWindow/F Rebin_Panel 
     2044End 
     2045 
     2046 
     2047Function Rebin_newGraph(ctrlName) : ButtonControl 
     2048        String ctrlName 
     2049 
     2050        //get the current selection 
     2051        ControlInfo popup0 
     2052        if(strlen(S_Value)==0 || cmpstr(S_Value,"No data loaded")==0) 
     2053                Abort "You must load data from a file into memory before plotting the data" 
     2054        Endif 
     2055         
     2056        A_PM_doNewGraph(S_Value) 
     2057         
     2058        DoWindow/F Rebin_Panel 
     2059End 
     2060 
     2061 
     2062 
     2063Function Rebin_by (ctrlName) : ButtonControl 
     2064        String ctrlName 
     2065         
     2066         
     2067        Variable len 
     2068        Variable ii,jj, helplp,kk 
     2069        String iStr,qStr,eStr,sqStr,qbStr,fsStr,wtStr,folderStr 
     2070        Variable rr,gg,bb 
     2071        NVAR binpts = root:Packages:NIST:Rebin:binning 
     2072         
     2073         
     2074         
     2075        ControlInfo popup0 
     2076        if(strlen(S_Value)==0) 
     2077                Abort "You must load data from a file into memory before plotting the data" 
     2078        Endif 
     2079        folderStr = S_Value      
     2080         
     2081        SetDataFolder $("root:"+folderStr) 
     2082         
     2083        Wave w0 = $(folderStr+"_q") 
     2084        Wave w1 = $(folderStr+"_i") 
     2085        Wave w2 = $(folderStr+"_s") 
     2086        Wave resW = $(folderStr+"_res") 
     2087         
     2088        len = numpnts(w0) 
     2089        Make/O/D/N=(len) w3,w4,w5  
     2090        w3 = resW[p][0]         //std dev of resolution fn 
     2091        w4 = resW[p][1]         //mean q-value 
     2092        w5 = resW[p][2]         //beamstop shadow factor 
     2093         
     2094        Make/O/D/N=(round((len-1)/binpts)) qbin,Ibin,sigbin,sqbin,qbbin,fsbin,helpTemp 
     2095 
     2096        qbin = 0 
     2097        Ibin = 0 
     2098        sigbin = 0 
     2099        sqbin = 0 
     2100        qbbin = 0 
     2101        fsbin = 0 
     2102        helptemp = 0 
     2103         
     2104        ii=0 
     2105        do 
     2106                qBin = qBin + w0[binpts*p+ii]/binpts 
     2107                iBin = iBin + w1[binpts*p+ii]/binpts 
     2108                helptemp = helptemp + w2[binpts*p+ii]^2 
     2109                sqBin = sqBin + w3[binpts*p+ii]/binpts 
     2110                qbBin = qbBin + w4[binpts*p+ii]/binpts 
     2111                fsBin = fsBin + w5[binpts*p+ii]/binpts 
     2112                 
     2113 
     2114      ii+=1 
     2115        while (ii<binpts) 
     2116       
     2117        sigBin  = sqrt(helptemp)/binpts 
     2118 
     2119        CheckDisplayed Ibin             //check top graph only 
     2120        if(V_flag==0) 
     2121                AppendToGraph Ibin vs qbin 
     2122                ModifyGraph log=1,mode(Ibin)=4,marker(Ibin)=29,msize(Ibin)=3,rgb(Ibin)=(65535,0,0) 
     2123                ErrorBars/T=0 ibin Y,wave=(sigbin,sigbin) 
     2124        endif 
     2125 
     2126        KillWaves/Z helpTemp,w3,w4,w5 
     2127         
     2128        SetDataFolder root: 
     2129End 
     2130 
     2131 
     2132// duplicates the binned data to its own data folder, fixing up the names 
     2133// as needed, then write out the file from the folder using the standard procedures 
     2134// then kill the data folder to clean up. 
     2135Function SaveResultBin(ctrlName) : ButtonControl 
     2136        String ctrlName 
     2137         
     2138        String finalname,folderStr 
     2139         
     2140        ControlInfo popup0 
     2141        folderStr=S_Value        
     2142         
     2143        SetDataFolder $("root:"+folderStr) 
     2144        WAVE qbin = qBin 
     2145        WAVE Ibin = iBin 
     2146        WAVE sigbin = sigBin 
     2147        WAVE sqbin = sqBin 
     2148        WAVE qbbin = qbBin 
     2149        WAVE fsbin = fsBin 
     2150         
     2151        NVAR/Z dQv = USANS_dQv 
     2152        Variable isSANS = (NVAR_Exists(dQv) == 0) 
     2153//      Print "isSANS = ",isSANS 
     2154        if(isSANS) 
     2155                // SANS data - make a smaller resWave 
     2156                Variable np=numpnts(qBin) 
     2157                Make/D/O/N=(np,4) $(folderStr+"_resBin") 
     2158                Wave res = $(folderStr+"_resBin") 
     2159                 
     2160                res[][0] = sqBin[p]             //sigQ 
     2161                res[][1] = qbBin[p]             //qBar 
     2162                res[][2] = fsBin[p]             //fShad 
     2163                res[][3] = qBin[p]              //Qvalues 
     2164         
     2165        endif 
     2166        finalname = folderStr+"bin" 
     2167         
     2168        SetDataFolder root: 
     2169        DuplicateDataSet(folderStr,finalname,1) 
     2170         
     2171        // rename the waves so that the correct ones are written out 
     2172        SetDataFolder $("root:"+finalname) 
     2173        Killwaves/Z $(finalname+"_q") 
     2174        Killwaves/Z $(finalname+"_i") 
     2175        Killwaves/Z $(finalname+"_s") 
     2176        if(isSANS)              //if USANS, keep the resolution waves, so that the writer will recognize the set as USANS 
     2177                Killwaves/Z $(finalname+"_res") 
     2178        endif 
     2179                 
     2180        WAVE qbin = qBin                //reset the wave references to the new folder (bin) 
     2181        WAVE Ibin = iBin 
     2182        WAVE sigbin = sigBin 
     2183        if(isSANS) 
     2184                Wave res = $(finalName+"_resBin") 
     2185        endif 
     2186         
     2187        Rename qbin $(finalname+"_q") 
     2188        Rename ibin $(finalname+"_i") 
     2189        Rename sigbin $(finalname+"_s") 
     2190        if(isSANS) 
     2191                Rename res $(finalname+"_res") 
     2192        endif 
     2193        SetDataFolder root: 
     2194         
     2195        // save out the data in the folder 
     2196        SaveDataSetToFile(finalName) 
     2197 
     2198        KillDataFolder/Z $("root:"+finalname) 
     2199 
     2200        SetDataFolder root: 
     2201        return(0) 
     2202End 
     2203 
     2204 
     2205 
     2206 
     2207/// end rebinning procedures 
  • sans/Release/trunk/NCNR_User_Procedures/Common/GaussUtils_v40.ipf

    r570 r766  
    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/Release/trunk/NCNR_User_Procedures/Common/NIST_XML_v40.ipf

    r690 r766  
    131131                                                        endif 
    132132                                                elseif(exists(xmlDataSetFolder+"dQl")) 
    133                                                         //USAS Data 
     133                                                        //USANS Data 
    134134                                                        Wave dQl = $(xmlDataSetFolder+"dQl") 
    135                                                         dQv = dQl[0] 
     135                                                        dQv = abs(dQl[0]) 
    136136                                                 
    137137                                                        USANS_CalcWeights(baseStr,dQv) 
     
    262262                                                endif 
    263263                                        elseif(exists(xmlDataFolder+"dQl")) 
    264                                                 //USAS Data 
     264                                                //USANS Data 
    265265                                                Wave dQl = $(xmlDataFolder+"dQl") 
    266                                                 dQv = - dQl[0]          //make it positive again 
     266                                                dQv = abs(dQl[0])               //make it positive again 
    267267                                         
    268268                                                USANS_CalcWeights(baseStr,dQv) 
     
    565565end 
    566566 
     567// 
     568// !!! nf.Sample_ID is not set correctly here, since it's not read in from the NIST 6-col data file 
     569// and SASprocessnote does not get set either! 
     570// 
    567571Function convertNISTtoNISTXML(fileStr) 
    568572        String fileStr 
     
    639643                         
    640644                endif 
    641                  
    642                 //Tidy up 
    643                 Variable i = 0 
    644                 do 
    645                         WAVE/Z wv= $(StringFromList(i,S_waveNames,";")) 
    646                         if( WaveExists(wv) == 0 ) 
    647                                 break 
    648                         endif 
    649                         KillWaves wv 
    650                         i += 1 
    651                 while (1)       // exit is via break statement 
    652  
    653645         
    654646        endif   //6-col data 
     
    666658 
    667659        writeNISTXML(outfileName, nf) 
     660 
     661        //Tidy up AFTER we're all done, since STRUCT points to wave0,wave1, etc. 
     662        Variable i = 0 
     663        do 
     664                WAVE/Z wv= $(StringFromList(i,S_waveNames,";")) 
     665                if( WaveExists(wv) == 0 ) 
     666                        break 
     667                endif 
     668                KillWaves wv 
     669                i += 1 
     670        while (1)       // exit is via break statement 
    668671 
    669672end 
     
    695698                if (stringmatch(buffer,"*FIRST File LABEL:*") == 1) 
    696699                        NISTfile.title = TrimWS(StringFromList(1,buffer, ":")) 
    697                 elseif(stringmatch(buffer,"*LABEL:*") == 1) 
     700                endif 
     701                if(stringmatch(buffer,"*LABEL:*") == 1) 
    698702                        NISTfile.title = TrimWS(StringFromList(1,buffer, ":")) 
     703                endif 
     704                if(stringmatch(buffer,"NSORT*") == 1) 
     705                        NISTfile.title = buffer 
    699706                endif 
    700707                 
     
    9941001 
    9951002 
    996  
    997  
    998 End 
    999  
    10001003/// See WriteModelData_v40.ipf for 6 column equivalent 
     1004// 
     1005// will abort if resolution wave is missing 
     1006// switches for USANS data if the proper global is found, otheriwse treats as SANS data 
     1007// 
    10011008Function ReWrite1DXMLData(folderStr) 
    10021009        String folderStr 
     
    10321039        endif 
    10331040         
    1034         Duplicate/O qw qbar,sigQ,fs 
    1035                  
    1036  
    1037                  
    1038         //Data 
    1039         Wave nf.Q = qw 
    1040         nf.unitsQ = "1/A" 
    1041         Wave nf.I = iw 
    1042         nf.unitsI = "1/cm" 
    1043         Wave nf.Idev = sw 
    1044         nf.unitsIdev = "1/cm" 
    1045         Wave nf.Qdev = sigQ 
    1046         nf.unitsQdev = "1/A" 
    1047         Wave nf.Qmean = qbar 
    1048         nf.unitsQmean = "1/A" 
    1049         Wave nf.Shadowfactor = fs 
    1050         nf.unitsShadowfactor = "none" 
    1051          
    1052          
    1053         //write out the standard header information 
    1054         //fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n",textw[0],textw[1] 
    1055          
    1056         //AJJ to fix with sensible values 
    1057         nf.run = "" 
    1058         nf.nameSASinstrument = "NIST IGOR Procedures" 
    1059         nf.SASnote = "" 
    1060         // 
    1061         nf.sample_ID = baseStr 
    1062         nf.title = baseStr 
    1063         nf.radiation = "neutron" 
    1064         //Do something with beamstop (rw[21]) 
    1065         nf.detector_name = "Re-written data" 
    1066  
    1067         nf.SASprocessnote =  "Modified data written from folder "+baseStr+" on "+(date()+" "+time()) 
    1068          
    1069         nf.nameSASProcess = "NIST IGOR" 
    1070  
    1071         //Close refnum 
     1041         
     1042        // if (USANS) 
     1043        // else (SANS is assumed) 
     1044        // endif 
     1045        NVAR/Z dQv = USANS_dQv          // in current DF 
     1046        if (NVAR_Exists(dQv)) 
     1047                //USANS data, proceed 
     1048                //Use the evil extra column for the resolution "information". Should probably switch to using slit_length in collimation. 
     1049                Duplicate/O qw,dumWave 
     1050                dumWave = dQv                   //written out as a positive value, since the column is identified by its label, dQl 
     1051                 
     1052                //Data 
     1053                Wave nf.Q = qw 
     1054                nf.unitsQ = "1/A" 
     1055                Wave nf.I = iw 
     1056                nf.unitsI = "1/cm" 
     1057                Wave nf.Idev = sw 
     1058                nf.unitsIdev = "1/cm" 
     1059                // for slit-smeared USANS, set only a 4th column to  -dQv 
     1060                Wave nf.dQl = dumWave 
     1061                nf.unitsdQl= "1/A" 
     1062         
     1063                //AJJ to fix with sensible values 
     1064                nf.run = "" 
     1065                nf.nameSASinstrument = "NIST IGOR Procedures" 
     1066                nf.SASnote = "" 
     1067                // 
     1068                nf.sample_ID = baseStr 
     1069                nf.title = baseStr 
     1070                nf.radiation = "neutron" 
     1071                //Do something with beamstop (rw[21]) 
     1072                nf.detector_name = "Re-written USANS data" 
     1073         
     1074                nf.SASprocessnote =  "Modified data written from folder "+baseStr+" on "+(date()+" "+time()) 
     1075                 
     1076                nf.nameSASProcess = "NIST IGOR" 
     1077                 
     1078        else 
     1079                //assume SANS data 
     1080                Duplicate/O qw qbar,sigQ,fs 
     1081                sigq = resw[p][0] 
     1082                qbar = resw[p][1] 
     1083                fs = resw[p][2] 
     1084         
     1085                         
     1086                //Data 
     1087                Wave nf.Q = qw 
     1088                nf.unitsQ = "1/A" 
     1089                Wave nf.I = iw 
     1090                nf.unitsI = "1/cm" 
     1091                Wave nf.Idev = sw 
     1092                nf.unitsIdev = "1/cm" 
     1093                Wave nf.Qdev = sigQ 
     1094                nf.unitsQdev = "1/A" 
     1095                Wave nf.Qmean = qbar 
     1096                nf.unitsQmean = "1/A" 
     1097                Wave nf.Shadowfactor = fs 
     1098                nf.unitsShadowfactor = "none" 
     1099                 
     1100                 
     1101                //write out the standard header information 
     1102                //fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n",textw[0],textw[1] 
     1103                 
     1104                //AJJ to fix with sensible values 
     1105                nf.run = "" 
     1106                nf.nameSASinstrument = "NIST IGOR Procedures" 
     1107                nf.SASnote = "" 
     1108                // 
     1109                nf.sample_ID = baseStr 
     1110                nf.title = baseStr 
     1111                nf.radiation = "neutron" 
     1112                //Do something with beamstop (rw[21]) 
     1113                nf.detector_name = "Re-written data" 
     1114         
     1115                nf.SASprocessnote =  "Modified data written from folder "+baseStr+" on "+(date()+" "+time()) 
     1116                 
     1117                nf.nameSASProcess = "NIST IGOR" 
     1118 
     1119        endif 
     1120 
    10721121         
    10731122        if(dialog) 
     
    10871136        Print "XML File written: ", GetFileNameFromPathNoSemi(fullPath) 
    10881137        KillWaves/Z tempShortProto 
     1138         
     1139        SetDataFolder root: 
     1140 
    10891141        Return(0) 
    10901142End 
     
    11421194        Variable fileref 
    11431195         
    1144         Open/R fileref as filestr 
     1196        Open/R/Z fileref as filestr 
    11451197        FReadLine fileref,  line 
    11461198        Close fileref 
  • sans/Release/trunk/NCNR_User_Procedures/Common/Packages/Invariant/Invariant_v40.ipf

    r639 r766  
    374374        endif 
    375375         
    376         Variable yesGuinier=0,nume,inv 
     376        Variable yesGuinier=0,nume,inv,scale 
    377377        // do the extrapolation of the correct type 
    378378        ControlInfo/W=Invariant_Panel check_0           //the Guinier box 
     
    397397         
    398398        if(yesGuinier) 
    399                 Make/O/D G_coef={1000,-1000}            //input 
     399//              Make/O/D G_coef={1000,-1000}            //input 
     400                WaveStats/M=1/Q/R=[0,(nbeg-1)] iw 
     401                scale = V_avg 
     402                Make/O/D G_coef={(scale),-1000}         //input -- with a better guess for the overall scale factor,  
    400403                FuncFit Guinier_Fit G_coef iw[0,(nbeg-1)] /I=1 /X=qw /W=sw /D  
    401404                extr_lqi= Guinier_Fit(G_coef,extr_lqq) 
  • sans/Release/trunk/NCNR_User_Procedures/Common/Packages/LinearizedFits/LinearizedFits_v40.ipf

    r616 r766  
    2020// - DID NOT prepend "A_" to the loader/unloader which is unique to this package 
    2121// 
     22// SEP 2010 - absorbed the duplicate procedure file FIT_Ops. 
     23// - removed "A_" prefix 
     24// - need to keep A_OpenFitPanel() in case old experiments are looking for this 
     25// - added the (completely unused) FITRPA procedures to this file 
    2226// 
    2327/////////////////////////////// 
     
    3337// of the FIT/RPA panel) 
    3438// 
    35 Proc A_OpenFitPanel() 
    36         If(WinType("A_FitPanel") == 0) 
     39Proc OpenFitPanel() 
     40        If(WinType("FitPanel") == 0) 
    3741                //create the necessary data folder 
    3842                NewDataFolder/O root:Packages 
     
    4751                Variable/G root:Packages:NIST:FIT:gBack = 0 
    4852                String/G root:Packages:NIST:FIT:gDataPopList = "none" 
    49                 A_FitPanel() 
     53                FitPanel() 
    5054        else 
    5155                //window already exists, just bring to front for update 
    52                 DoWindow/F A_FitPanel 
     56                DoWindow/F FitPanel 
    5357                CheckBox check0,value=0         //deselect the checkbox to use cursors 
    5458        endif 
    5559        //pop the file menu 
    56         A_FIT_FilePopMenuProc("",1,"") 
     60        FIT_FilePopMenuProc("",1,"") 
     61End 
     62 
     63Proc A_OpenFitPanel() 
     64        OpenFitPanel() 
    5765End 
    5866 
    5967//the actual window recreation macro to draw the fit panel. Globals and data folder must  
    6068// already be initialized 
    61 Window A_FitPanel() 
    62         //String angst = root:Packages:NIST:gAngstStr 
    63         String angst = "A" 
     69Window FitPanel() 
     70        String angst = StrVarOrDefault("root:Packages:NIST:gAngstStr", "A" ) 
     71//      String angst = "A" 
    6472        PauseUpdate; Silent 1           // building window... 
    6573        NewPanel /W=(461,46,735,455)/K=1 
     
    7482        PopupMenu ywave,pos={13,60},size={154,19},title="Data File" 
    7583        PopupMenu ywave,help={"Select the experimental intensity values"} 
    76         PopupMenu ywave,mode=1,value=root:Packages:NIST:FIT:gDataPopList,proc=A_FIT_FilePopMenuProc 
    77         Button loadButton,pos={13,92},size={130,20},proc=A_FIT_Load_Proc,title="Load and Plot File" 
     84        PopupMenu ywave,mode=1,value=root:Packages:NIST:FIT:gDataPopList,proc=FIT_FilePopMenuProc 
     85        Button loadButton,pos={13,92},size={130,20},proc=FIT_Load_Proc,title="Load and Plot File" 
    7886        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."} 
    79         Button helpButton,pos={237,28},size={25,20},proc=A_showFITHelp,title="?" 
     87        Button helpButton,pos={237,28},size={25,20},proc=showFITHelp,title="?" 
    8088        Button helpButton,help={"Show help file for linearized fitting"} 
    8189        PopupMenu ymodel,pos={20,281},size={76,19},title="y-axis" 
    8290        PopupMenu ymodel,help={"This popup selects how the y-axis will be linearized based on the chosen data"} 
    8391        PopupMenu ymodel,mode=1,value= #"\"I;log(I);ln(I);1/I;I^a;Iq^a;I^a q^b;1/sqrt(I);ln(Iq);ln(Iq^2)\"" 
    84         Button GoFit,pos={60,367},size={70,20},proc=A_DispatchModel,title="Do the Fit" 
     92        Button GoFit,pos={60,367},size={70,20},proc=DispatchModel,title="Do the Fit" 
    8593        Button GoFit,help={"This button will do the specified fit using the selections in this panel"} 
    86         Button DoneButton,pos={180,367},size={50,20},proc=A_FITDoneButton,title="Done" 
     94        Button DoneButton,pos={180,367},size={50,20},proc=FITDoneButton,title="Done" 
    8795        Button DoneButton,help={"This button will close the panel and the associated graph"} 
    8896        SetVariable lolim,pos={64,147},size={134,17},title="Lower Limit" 
     
    109117        SetVariable expc,help={"This sets the exponent \"c\" for some x-axis formats. The value is ignored if the model does not use \"c\" as an adjustable exponent"} 
    110118        SetVariable expc,limits={-10,10,0},value= root:Packages:NIST:FIT:gExpC 
    111         Button sh_all,pos={65,193},size={130,20},proc=A_ShowAllButtonProc,title="Show Full q-range" 
     119        Button sh_all,pos={65,193},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range" 
    112120        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."} 
    113121         
    114         Button FIT_PathButton,pos={10,28},size={80,20},proc=A_FIT_PickPathButtonProc,title="Pick Path" 
     122        Button FIT_PathButton,pos={10,28},size={80,20},proc=FIT_PickPathButtonProc,title="Pick Path" 
    115123 
    116124EndMacro 
    117125 
    118126 
    119 Proc A_FITDoneButton(ctrlName): ButtonControl 
     127Proc FITDoneButton(ctrlName): ButtonControl 
    120128        String ctrlName 
    121         DoWindow/K A_FitWindow 
    122         DoWindow/K A_FitPanel 
     129        DoWindow/K FitWindow 
     130        DoWindow/K FitPanel 
    123131end 
    124132 
    125 Proc A_showFITHelp(ctrlName): ButtonControl 
     133Proc showFITHelp(ctrlName): ButtonControl 
    126134        String ctrlName 
    127135        DisplayHelpTopic/Z/K=1 "Linearized Fits" 
     
    133141//Loads the selected file for fitting 
    134142//graphs the data as needed 
    135 Proc A_FIT_Load_Proc(ctrlName): ButtonControl 
     143Proc FIT_Load_Proc(ctrlName): ButtonControl 
    136144        String ctrlName 
    137145         
     
    149157        Endif 
    150158        //get a valid file based on this partialName and catPathName 
    151         tempName = A_FindValidFilename(partialName) 
     159        tempName = FindValidFilename(partialName) 
    152160 
    153161        //prepend path to tempName for read routine  
     
    157165         
    158166        //load in the data (into the root directory) 
    159         A_LoadOneDDataWithName(tempName,0)              //let A_Rescale_Data() do the plotting 
     167        A_LoadOneDDataWithName(tempName,0)              //let Rescale_Data() do the plotting 
    160168        //Print S_fileName 
    161169        //Print tempName 
     
    171179 
    172180        //Plot, and adjust the scaling to match the axis scaling set by the popups 
    173         A_Rescale_Data(dataStr) 
     181        Rescale_Data(dataStr) 
    174182         
    175183End 
     
    177185//gets a valid file list (simply not the files with ".SAn" in the name) 
    178186// 
    179 Function A_FIT_FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl 
     187Function FIT_FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl 
    180188        String ctrlName 
    181189        Variable popNum 
    182190        String popStr 
    183191         
    184         String tempStr=A_filterButtonProc(ctrlName) 
     192        String tempStr=filterButtonProc(ctrlName) 
    185193        if(strlen(tempStr)==0) 
    186194                tempStr = "Pick the data path" 
     
    200208//the fitted range, which is the default display after a fit is performed 
    201209// 
    202 Function A_ShowAllButtonProc(ctrlName) : ButtonControl 
     210Function ShowAllButtonProc(ctrlName) : ButtonControl 
    203211        String ctrlName 
    204212 
    205213        //bring the FitWindow to the front and Autoscale the axes 
    206         DoWindow/F A_FitWindow 
     214        DoWindow/F FitWindow 
    207215        SetAxis/A 
    208216End 
     
    213221// 
    214222// dataStr is the root:folder: of the data that was loaded 
    215 Function A_Rescale_Data(dataStr) 
     223Function Rescale_Data(dataStr) 
    216224        String dataStr 
    217225         
    218226        //Scaling exponents and background value 
    219227        Variable pow_a,pow_b,pow_c,bkg 
    220         ControlInfo/W=A_FitPanel expa 
     228        ControlInfo/W=FitPanel expa 
    221229        pow_a = V_value 
    222         ControlInfo/W=A_FitPanel expb 
     230        ControlInfo/W=FitPanel expb 
    223231        pow_b = V_value 
    224         ControlInfo/W=A_FitPanel expc 
     232        ControlInfo/W=FitPanel expc 
    225233        pow_c = V_value 
    226         ControlInfo/W=A_FitPanel back 
     234        ControlInfo/W=FitPanel back 
    227235        bkg = V_value 
    228236         
     
    241249        endif 
    242250        //if q^c is the x-scaling, c must be be within limits and also non-zero 
    243         ControlInfo/W=A_FitPanel xModel 
     251        ControlInfo/W=FitPanel xModel 
    244252        If (cmpstr("q^c",S_Value) == 0) 
    245253                if(pow_c == 0)  
     
    265273        String xlabel,ylabel,xstr,ystr 
    266274        //check for proper y-scaling selection, make the necessary waves 
    267         ControlInfo/W=A_FitPanel yModel 
     275        ControlInfo/W=FitPanel yModel 
    268276        ystr = S_Value 
    269277//      print "ystr = ",ystr 
     
    365373        Variable low,high 
    366374         
    367         ControlInfo/W=A_FitPanel lolim 
     375        ControlInfo/W=FitPanel lolim 
    368376        low = V_value 
    369         ControlInfo/W=A_FitPanel uplim 
     377        ControlInfo/W=FitPanel uplim 
    370378        high = V_value 
    371379        if ((high<low) || (high==low)) 
     
    374382        endif 
    375383         
    376         ControlInfo/W=A_FitPanel xModel 
     384        ControlInfo/W=FitPanel xModel 
    377385        xstr = S_Value 
    378386        do 
     
    422430         
    423431//      String cleanLastFileName = "root:"+CleanupName(gLastFileName,0) 
    424         If(WinType("A_FitWindow") == 0) 
     432        If(WinType("FitWindow") == 0) 
    425433                Display /W=(5,42,480,400)/K=1 yAxisWave vs xAxisWave 
    426434                ModifyGraph mode=3,standoff=0,marker=8,opaque=1 
    427435                ErrorBars/T=0 yAxisWave Y,wave=(yErrWave,yErrWave) 
    428                 DoWindow/C A_FitWindow 
     436                DoWindow/C FitWindow 
    429437        else 
    430438                //window already exists, just bring to front for update 
    431                 DoWindow/F A_FitWindow 
     439                DoWindow/F FitWindow 
    432440                // remove old text boxes 
    433441                TextBox/K/N=text_1 
     
    439447        TextBox/C/N=textLabel/A=RB "File = "+baseStr 
    440448        //clear the old fit from the window, if it exists 
    441         RemoveFromGraph/W=A_FitWindow/Z fit_yAxisWave 
     449        RemoveFromGraph/W=FitWindow/Z fit_yAxisWave 
    442450         
    443451        // add the cursors if desired...         
    444452        //see if the user wants to use the data specified by the cursors - else use numerical values 
    445453         
    446         ControlInfo/W=A_FitPanel check0         //V_value = 1 if it is checked, meaning yes, use cursors 
     454        ControlInfo/W=FitPanel check0           //V_value = 1 if it is checked, meaning yes, use cursors 
    447455        yes_cursors = V_value 
    448456 
    449         DoWindow/F A_FitWindow 
     457        DoWindow/F FitWindow 
    450458        ShowInfo 
    451459        if(yes_cursors) 
     
    461469                if(V_flag == 1)                 //level NOT found 
    462470                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value" 
    463                         //DoWindow/K A_FitWindow 
     471                        //DoWindow/K FitWindow 
    464472                        //Abort 
    465473                endif 
     
    469477                if(V_flag == 1) 
    470478                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value" 
    471                         //DoWindow/K A_FitWindow 
     479                        //DoWindow/K FitWindow 
    472480                        //Abort 
    473481                endif 
     
    489497//and the results are plotted 
    490498// function works in root level data folder (where the loaded 1-d data will be) 
    491 Function A_DispatchModel(GoFit) : ButtonControl 
     499Function DispatchModel(GoFit) : ButtonControl 
    492500        String GoFit 
    493501 
    494502        //check for the FitWindow - to make sure that there is data to fit 
    495         If(WinType("A_FitWindow") == 0)         //if the window doesn't exist 
     503        If(WinType("FitWindow") == 0)           //if the window doesn't exist 
    496504                Abort "You must Load and Plot a File before fitting the data" 
    497505        endif 
    498506         
    499507        // rescale the data, to make sure it's as selected on the panel 
    500         ControlInfo/W=A_FitPanel $"ywave" 
     508        ControlInfo/W=FitPanel $"ywave" 
    501509        String partialName = CleanupName(S_value,0) 
    502         A_Rescale_Data("root:"+partialName+":") 
     510        Rescale_Data("root:"+partialName+":") 
    503511         
    504512        // now go do the fit 
     
    507515        Variable low,high 
    508516         
    509         ControlInfo/W=A_FitPanel lolim 
     517        ControlInfo/W=FitPanel lolim 
    510518        low = V_value 
    511         ControlInfo/W=A_FitPanel uplim 
     519        ControlInfo/W=FitPanel uplim 
    512520        high = V_value 
    513521        if ((high<low) || (high==low)) 
     
    518526        //try including residuals on the graph /R=residWave, explicitly place on new axis 
    519527        //if only /R used, residuals are automatically placed on graph 
    520          
    521         CurveFit line yAxisWave(xcsr(A),xcsr(B)) /I=1 /X=xAxisWave /W=yWtWave /D   
     528        // -- NOTE that Rescale_Data() calculates the weighting wave as 1/err (like the old days) so the flag is correctly 
     529        // /I=0, not /I=1 
     530         
     531        CurveFit line yAxisWave(xcsr(A),xcsr(B)) /I=0 /X=xAxisWave /W=yWtWave /D   
    522532        //CurveFit line yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave  /R /D   
    523533        ModifyGraph rgb(fit_yAxisWave)=(0,0,0) 
     
    528538        //ystr and xstr are the axis strings - filter with a do-loop 
    529539        String ystr="",xstr="" 
    530         //ControlInfo/W=A_FitPanel ywave 
     540        //ControlInfo/W=FitPanel ywave 
    531541        Wave xw = $("root:"+partialName+":"+partialName + "_q") 
    532         ControlInfo/W=A_FitPanel yModel 
     542        ControlInfo/W=FitPanel yModel 
    533543        ystr = S_Value 
    534         ControlInfo/W=A_FitPanel xModel 
     544        ControlInfo/W=FitPanel xModel 
    535545        xstr = S_Value 
    536546         
     
    553563        do 
    554564                If (cmpstr("I",ystr) == 0) 
    555                         textstr_3 = "I(q=0) =  "  + num2str(W_coef[0]) 
     565                        textstr_3 = "I(q=0) =  "  + num2str(W_coef[0]) +" ± "+num2str(W_sigma[0]) 
    556566                        break    
    557567                endif 
    558568                If (cmpstr("ln(I)",ystr) == 0) 
    559                         textstr_3 = "I(q=0) =  "  + num2str(exp(W_coef[0])) 
     569                        izerr = abs(exp(W_coef[0]) - exp(W_coef[0]+W_sigma[0])) 
     570                        textstr_3 = "I(q=0) =  "  + num2str(exp(W_coef[0]))+" ± " + num2str(izerr) 
    560571                        if(cmpstr("q^2",xstr) == 0)     //then a Guinier plot for a sphere (3-d) 
    561572                                rg = sqrt(-3*W_coef[1]) 
     
    591602                endif 
    592603                If (cmpstr("1/sqrt(I)",ystr) == 0) 
    593                         textstr_3 = "I(q=0) =  "  + num2str((W_coef[0])^2) 
     604                        izerr = abs( (W_coef[0])^-2 - (W_coef[0]+W_sigma[0])^-2 ) 
     605                        textstr_3 = "I(q=0) =  "  + num2str((W_coef[0])^-2)+" ± " + num2str(izerr) 
    594606                        break 
    595607                endif 
     
    621633        while(0) 
    622634        //kill the old textboxes, if they exist 
    623         TextBox/W=A_FitWindow/K/N=text_1 
    624         TextBox/W=A_FitWindow/K/N=text_2 
    625         TextBox/W=A_FitWindow/K/N=text_3 
     635        TextBox/W=FitWindow/K/N=text_1 
     636        TextBox/W=FitWindow/K/N=text_2 
     637        TextBox/W=FitWindow/K/N=text_3 
    626638        // write the new text boxes 
    627         TextBox/W=A_FitWindow/N=text_1/A=LT textstr_1 
    628         TextBox/W=A_FitWindow/N=text_2/A=LC textstr_2 
     639        TextBox/W=FitWindow/N=text_1/A=LT textstr_1 
     640        TextBox/W=FitWindow/N=text_2/A=LC textstr_2 
    629641        If (cmpstr("",textstr_3) != 0)          //only display textstr_3 if it isn't null 
    630                 TextBox/W=A_FitWindow/N=text_3/A=LB textstr_3 
     642                TextBox/W=FitWindow/N=text_3/A=LB textstr_3 
    631643        endif 
    632644         
    633645        //adjust the plot range to reflect the actual fitted range 
    634646        //cursors are already on the graph, done by Rescale_Data() 
    635         A_AdjustAxisToCursors() 
     647        AdjustAxisToCursors() 
    636648         
    637649End 
     
    641653// 
    642654// will expand the scale to show an extra 5 points in each direction (if available) 
    643 Function A_AdjustAxisToCursors() 
    644  
    645         DoWindow/F A_FitWindow 
     655Function AdjustAxisToCursors() 
     656 
     657        DoWindow/F FitWindow 
    646658        WAVE xAxisWave = root:xAxisWave 
    647659        WAVE yAxisWave = root:yAxisWave 
     
    649661        Variable extraPts = 5, num=numpnts(xAxisWave) 
    650662         
    651         String csrA = CsrInfo(A ,"A_FitWindow") 
    652         String csrB = CsrInfo(B ,"A_FitWindow") 
     663        String csrA = CsrInfo(A ,"FitWindow") 
     664        String csrB = CsrInfo(B ,"FitWindow") 
    653665         
    654666        //x-levels, these are monotonic 
     
    718730// Windows, where the data file must be .txt (and possibly OSX) 
    719731// 
    720 Function/S A_filterButtonProc(ctrlName) 
     732Function/S filterButtonProc(ctrlName) 
    721733        String ctrlName 
    722734 
     
    744756        endfor 
    745757        //remove VAX version numbers 
    746         newList = A_RemoveVersNumsFromList(newList) 
     758        newList = RemoveVersNumsFromList(newList) 
    747759        //sort 
    748760        newList = SortList(newList,";",0) 
     
    753765 
    754766 
     767//////////////////////////////////////// FIT RPA /////////////////////////////////////////// 
     768//**************************************** 
     769//procedures for creating and initializing the FITRPA panel 
     770//global variables (numerical only) are kept in root:myGlobals:FITRPA folder, 
     771//created as needed 
     772// 
     773// very similar in function to the FIT panel 
     774// 
     775Proc OpenFitRPAPanel() 
     776        If(WinType("FitRPAPanel") == 0) 
     777                //create the necessary data folder 
     778                NewDataFolder/O root:myGlobals:FITRPA 
     779                //initialize the values 
     780                Variable/G root:myGlobals:FITRPA:gLolim = 0.02 
     781                Variable/G root:myGlobals:FITRPA:gUplim = 0.04 
     782                Variable/G root:myGlobals:FITRPA:gBack = 0 
     783                Variable/G root:myGlobals:FITRPA:gLambda = 6.0 
     784                PathInfo/S catPathName 
     785                String localpath = S_path 
     786                if (V_flag == 0) 
     787                //path does not exist - no folder selected 
     788                        String/G root:myGlobals:FITRPA:gPathStr = "no folder selected" 
     789                else 
     790                        String/G root:myGlobals:FITRPA:gPathStr = localpath 
     791                endif 
     792                String/G    root:myGlobals:FITRPA:gDataPopList = "none" 
     793                FitRPAPanel() 
     794        else 
     795                //window already exists, just bring to front for update 
     796                DoWindow/F FitRPAPanel 
     797        endif 
     798        //pop the menu 
     799        FilePopMenuProc("",1,"") 
     800End 
     801 
     802//used on the fit/rpa panel to select the path for the data  
     803// - automatically pops the file list after the new  path selection 
     804Function FITRPAPickPathButton(ctrlName) : ButtonControl 
     805        String ctrlName 
     806 
     807        Variable err = A_PickPath()             //sets global path value 
     808        SVAR pathStr = root:myGlobals:gCatPathStr 
     809         
     810        //set the global string for NSORT to the selected pathname 
     811        String/G root:myGlobals:FITRPA:gPathStr = pathStr 
     812         
     813        //call the popup menu proc's to re-set the menu choices 
     814        FilePopMenuProc("filePopup",1,"") 
     815         
     816End 
     817 
     818//gets a valid file list (simply not the files with ".SAn" in the name) 
     819// 
     820Function FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl 
     821        String ctrlName 
     822        Variable popNum 
     823        String popStr 
     824 
     825        String tempStr=ReducedDataFileList(ctrlName) 
     826        if(strlen(tempStr)==0) 
     827                tempStr = "Pick the data path" 
     828        Endif 
     829        String/G root:myGlobals:FITRPA:gDataPopList =   tempStr //function is in NSORT.ipf 
     830        ControlUpdate filePopup 
     831         
     832End 
     833 
     834 
     835 
     836// window recreation macro to draw the fit/rpa panel 
     837//globals and data folders must be present before drawing panel 
     838// 
     839Window FitRPAPanel()  
     840 
     841        String angst = root:myGlobals:gAngstStr 
     842        PauseUpdate; Silent 1           // building window... 
     843        NewPanel /W=(250,266,591,579)/K=1 
     844        ModifyPanel cbRGB=(32768,54528,65280) 
     845        ModifyPanel fixedSize=1 
     846        SetDrawLayer UserBack 
     847        SetDrawEnv fstyle= 1 
     848        DrawText 81,19,"Select Experimental Data" 
     849        SetDrawEnv fstyle= 1 
     850        DrawText 97,102,"q-range to fit ("+angst+"^-1)" 
     851        SetDrawEnv fstyle= 1 
     852        DrawText 87,239,"Select the fit parameters" 
     853        SetDrawEnv fillpat= 0 
     854        DrawRect 1,103,338,224 
     855        SetDrawEnv fillpat= 0 
     856        DrawRect 1,20,337,83 
     857        SetDrawEnv fillpat= 0 
     858        DrawRect 2,241,337,275 
     859//      Button PathButton,pos={6,26},size={80,20},proc=FitRPAPickPathButton,title="Pick Path" 
     860//      Button PathButton,help={"Select the local path to the folder containing your SANS data"} 
     861//      SetVariable setPath,pos={95,29},size={240,17},title="Path:" 
     862//      SetVariable setPath,help={"The current path to the local folder with SANS data"} 
     863//      SetVariable setPath,fSize=10 
     864//      SetVariable setPath,limits={0,0,0},value= root:myGlobals:FITRPA:gPathStr 
     865        PopupMenu filePopup,pos={8,30},size={96,21},proc=FilePopMenuProc,title="Files" 
     866        PopupMenu filePopup,help={"Select the data file to load."} 
     867        PopupMenu filePopup,mode=5,popvalue="none",value= #"root:myGlobals:FITRPA:gDataPopList" 
     868        SetVariable lambda,pos={111,250},size={120,18},title="Lambda ("+angst+")" 
     869        SetVariable lambda,help={"This sets the wavelength for the multiple scattering corrections."} 
     870        SetVariable lambda,limits={0,10,0},value= root:myGlobals:FITRPA:gLambda 
     871        Button GoFit,pos={60,286},size={80,20},proc=DoFITRPA,title="Do the Fit" 
     872        Button GoFit,help={"This button will do the specified fit using the selections in this panel"} 
     873        SetVariable lolim,pos={82,113},size={134,28},title="Lower Limit" 
     874        SetVariable lolim,help={"Enter the lower q-limit to perform the fit ("+angst+"^-1)"} 
     875        SetVariable lolim,limits={0,5,0},value= root:myGlobals:FITRPA:gLolim 
     876        SetVariable uplim,pos={80,140},size={134,28},title="Upper Limit" 
     877        SetVariable uplim,help={"Enter the upper q-limit to perform the fit ("+angst+"^-1)"} 
     878        SetVariable uplim,limits={0,5,0},value= root:myGlobals:FITRPA:gUplim 
     879        CheckBox RPA_check0,pos={64,198},size={190,20},title="Use cursor range from FitWindow" 
     880        CheckBox RPA_check0,help={"Checking this will perform a fit between the cursors on the graph in FitWindow and ignore the numerical limits typed above"},value=0 
     881        PopupMenu model,pos={3,249},size={101,21},title="Standard" 
     882        PopupMenu model,help={"This popup selects which standard should be used to fit this data"} 
     883        PopupMenu model,mode=1,popvalue="B",value= #"\"B;C;AS\"" 
     884        Button sh_all,pos={82,168},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range" 
     885        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."} 
     886        Button loadButton,pos={20,55},size={70,20},proc=FITRPA_Load_Proc,title="Load File" 
     887        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."} 
     888        Button helpButton,pos={270,55},size={25,20},proc=showFITHelp,title="?" 
     889        Button helpButton,help={"Show help file for RPA fitting"} 
     890        Button DoneButton,pos={200,286},size={50,20},proc=FITRPADoneButton,title="Done" 
     891        Button DoneButton,help={"This button will close the panel and the associated graph"} 
     892EndMacro 
     893 
     894 
     895Proc FITRPADoneButton(ctrlName) : ButtonControl 
     896        String ctrlName 
     897        DoWindow/K FitWindow 
     898        DoWindow/K FitRPAPanel 
     899end 
     900 
     901//dispatches the fit to the appropriate model 
     902//and the appropriate range, based on selections in the panel 
     903// 
     904Proc DoFITRPA(ctrlName) : ButtonControl 
     905        String ctrlName 
     906        // 
     907        String cleanLastFileName = CleanupName(root:Packages:NIST:gLastFileName,0) 
     908        String tmpStr = "root:"+cleanLastFileName+":"+cleanLastFileName 
     909 
     910        Duplicate/O $(tmpStr+"_q") xAxisWave 
     911        Duplicate/O $(tmpStr+"_i") yAxisWave 
     912        Duplicate/O $(tmpStr+"_s") yErrWave,yWtWave,residWave 
     913 
     914        yWtWave = 1/yErrWave 
     915         
     916        String xlabel = "q (A^-1)" 
     917        String ylabel = "Intensity" 
     918        //Check to see if the FitWindow exists 
     919        //Plot the data in a FitWindow 
     920        If(WinType("FitWindow") == 0) 
     921                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave 
     922                ModifyGraph mode=3,standoff=0,marker=8 
     923                ErrorBars/T=0 yAxisWave Y,wave=(yErrWave,yErrWave) 
     924                DoWindow/C FitWindow 
     925                ShowInfo 
     926        else 
     927                //window already exists, just bring to front for update 
     928                DoWindow/F FitWindow 
     929                // remove old text boxes 
     930                TextBox/K/N=text_1 
     931                TextBox/K/N=text_2 
     932                TextBox/K/N=text_3 
     933        endif 
     934         
     935        //see if the user wants to use the data specified by the cursors - else use numerical values 
     936        Variable xlow,xhigh,ylow,yhigh,yes_cursors 
     937        ControlInfo/W=FitRPAPanel RPA_check0            //V_value = 1 if it is checked, meaning yes, use cursors 
     938        yes_cursors = V_value 
     939 
     940        ControlInfo/W=FitRPAPanel lolim 
     941        xlow = V_value 
     942        ControlInfo/W=FitRPAPanel uplim 
     943        xhigh = V_value 
     944        if(yes_cursors) 
     945                xlow = xAxisWave[xcsr(A)] 
     946                xhigh = xAxisWave[xcsr(B)] 
     947                if(xlow > xhigh) 
     948                        xhigh = xlow 
     949                        xlow = xAxisWave[xcsr(B)] 
     950                endif 
     951//              Print "xlow,xhigh = ",xlow,xhigh 
     952                ylow = yAxisWave[xcsr(A)] 
     953                yhigh = yAxisWave[xcsr(B)] 
     954                if(ylow > yhigh) 
     955                        ylow=yhigh 
     956                        yhigh = yAxisWave[xcsr(A)] 
     957                endif 
     958        else 
     959                FindLevel/P/Q xAxisWave, xlow 
     960                if(V_flag == 1)                 //level NOT found 
     961                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value" 
     962                        DoWindow/K FitWindow 
     963                        Abort 
     964                endif 
     965                Cursor/P A, yAxisWave,trunc(V_LevelX)+1 
     966                ylow = yAxisWave[V_LevelX] 
     967                FindLevel/P/Q xAxisWave, xhigh 
     968                if(V_flag == 1) 
     969                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value" 
     970                        DoWindow/K FitWindow 
     971                        Abort 
     972                endif 
     973                Cursor/P B, yAxisWave,trunc(V_LevelX) 
     974                yhigh = yAxisWave[V_LevelX] 
     975                if(ylow > yhigh) 
     976                        yhigh=ylow 
     977                        ylow = yAxisWave[V_levelX] 
     978                endif 
     979        endif   //if(V_value) 
     980        SetAxis bottom,xlow,xhigh 
     981         
     982//      print "ylow,yhigh",ylow,yhigh 
     983         
     984        //Get the rest of the data from the panel 
     985        //such as which standard, the wavelength 
     986        ControlInfo/W=FitRPAPanel model 
     987 
     988        //find the model name 
     989        String modelName = S_value 
     990         
     991        Variable first_guess, seglength,iabs,iarb,thick 
     992        Make/D/O/N=2 fitParams 
     993         
     994        seglength = 6.8 
     995         
     996        first_guess = 1.0 
     997        fitParams[0] = first_guess 
     998        fitParams[1] = seglength 
     999 
     1000        If (cmpstr(modelName,"B")==0) 
     1001                iabs = BStandardFunction(fitParams,xlow) 
     1002                fitParams[0] = yhigh/iabs 
     1003//              Print fitParams[0],fitParams[1]  
     1004                FuncFit BStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D 
     1005                iarb = BStandardFunction(fitParams, 0.0) 
     1006                iabs = iarb/fitParams[0] 
     1007                thick = 0.153 
     1008        endif 
     1009        If (cmpstr(modelName,"C")==0) 
     1010                iabs = CStandardFunction(fitParams,xlow) 
     1011                fitParams[0] = yhigh/iabs 
     1012                FuncFit CStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D 
     1013                iarb = CStandardFunction(fitParams, 0.0) 
     1014                iabs = iarb/fitParams[0] 
     1015                thick= 0.153 
     1016        endif 
     1017        If (cmpstr(modelName,"AS")==0) 
     1018                iabs = ASStandardFunction(fitParams,xlow) 
     1019                fitParams[0] = yhigh/iabs 
     1020                FuncFit ASStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D 
     1021                iarb = ASStandardFunction(fitParams, 0.0) 
     1022                iabs = iarb/fitParams[0] 
     1023                thick = 0.1 
     1024        endif 
     1025        ModifyGraph rgb(fit_yAxisWave)=(0,0,0) 
     1026        Label left ylabel 
     1027        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead         
     1028 
     1029        ControlInfo/W=FitRPAPanel lambda 
     1030         
     1031        Variable cor_mult = 1.0 + 2.2e-4*V_Value^2 
     1032         
     1033        //WAVE W_coef=W_coef 
     1034        //WAVE W_sigma=W_sigma 
     1035        String textstr_1,textstr_2,textstr_3 = "" 
     1036        textstr_1 = "Scaling Parameter: "+num2str(fitParams[0])+" ± "+num2str(W_sigma[0]) 
     1037        textstr_1 += "\rSegment Length: "+num2str(fitParams[1])+" ± "+num2str(W_sigma[1]) 
     1038        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3)) 
     1039         
     1040        textstr_2 = "Cross section at q=0:  Iabs(0) = "+num2str(iabs)+"cm\S-1\M" 
     1041        textstr_2 += "\rData extrapolated to q=0: Im(0) = "+num2str(iarb)+" Counts/(10\S8\M  Mon cts)" 
     1042        textstr_2 += "\rData corrected for multiple scattering: I(0) = "+num2str(iarb/cor_mult)+" Counts/(10\S8\M  Mon cnts)" 
     1043         
     1044        textstr_3 = "In the ABS protocol, " 
     1045        textstr_3 += "\rStandard Thickness, d = "+num2str(thick)+"cm" 
     1046        textstr_3 += "\rI(0), Iarb(0) = "+num2str(iarb/cor_mult)+"Counts/(10\S8\M Mon cts)" 
     1047        textstr_3 += "\rStandard Cross Section, Iabs(0) = "+num2str(iabs)+"cm\S-1\M" 
     1048        TextBox/K/N=text_1 
     1049        TextBox/K/N=text_2 
     1050        TextBox/K/N=text_3 
     1051        TextBox/N=text_2/A=RT textstr_2 
     1052        TextBox/N=text_3/A=RC textstr_3 
     1053        TextBox/N=text_1/A=RB textstr_1 
     1054         
     1055End 
     1056 
     1057//loads the file selected in the popup for fitting with POL 
     1058//standard functions. Reads the wavelength from the header, using 
     1059//6 A as the default 
     1060//plots the data in FitWindow after reading the file 
     1061//updates lambda and full q-range  on the Panel 
     1062// 
     1063Proc FITRPA_Load_Proc(ctrlName): ButtonControl 
     1064        String ctrlName 
     1065        //Load the data 
     1066        String tempName="",partialName="" 
     1067        Variable err 
     1068        ControlInfo $"filePopup" 
     1069        //find the file from the partial filename 
     1070        If( (cmpstr(S_value,"")==0) || (cmpstr(S_value,"none")==0) ) 
     1071                //null selection, or "none" from any popup 
     1072                Abort "no file selected in popup menu" 
     1073        else 
     1074                //selection not null 
     1075                partialName = S_value 
     1076                //Print partialName 
     1077        Endif 
     1078        //get a valid file based on this partialName and catPathName 
     1079        tempName = FindValidFilename(partialName) 
     1080 
     1081        Variable lambdaFromFile=GetLambdaFromReducedData(tempName) 
     1082        Variable/G root:myGlobals:FITRPA:gLambda = lambdaFromFile 
     1083        Print "Lambda in file read as:", lambdaFromFile 
     1084         
     1085        //prepend path to tempName for read routine  
     1086        PathInfo catPathName 
     1087        tempName = S_path + tempName 
     1088         
     1089        //load in the data (into the root directory) 
     1090        LoadOneDDataWithName(tempName,0) 
     1091        //Print S_fileName 
     1092        //Print tempName 
     1093         
     1094        String cleanLastFileName = CleanupName(root:Packages:NIST:gLastFileName,0) 
     1095        String tmpStr = "root:"+cleanLastFileName+":"+cleanLastFileName 
     1096 
     1097        Duplicate/o $(tmpStr+"_q") xAxisWave 
     1098        Duplicate/o $(tmpStr+"_i") yAxisWave 
     1099        Duplicate/o $(tmpStr+"_s") yErrWave 
     1100         
     1101        Variable xmin, xmax 
     1102        WaveStats/Q xAxisWave 
     1103        root:myGlobals:FITRPA:gLolim=V_min 
     1104        root:myGlobals:FITRPA:gUplim=V_max 
     1105        ControlUpdate/W=FITRPAPanel/A 
     1106         
     1107        //Check to see if the FitWindow exists 
     1108        //Plot the data in a FitWindow 
     1109        If(WinType("FitWindow") == 0) 
     1110                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave 
     1111                ModifyGraph mode=3,standoff=0,marker=8 
     1112                ErrorBars/T=0 yAxisWave Y,wave=(yErrWave,yErrWave) 
     1113                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName 
     1114                DoWindow/C FitWindow 
     1115                ShowInfo 
     1116        else 
     1117                //window already exists, just bring to front for update 
     1118                DoWindow/F FitWindow 
     1119                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName 
     1120        endif 
     1121        // remove old text boxes 
     1122        TextBox/K/N=text_1 
     1123        TextBox/K/N=text_2 
     1124        TextBox/K/N=text_3 
     1125        RemoveFromGraph/W=fitWindow /Z fit_yAxisWave 
     1126        SetAxis/A 
     1127         
     1128        //put cursors on the graph at first and last points 
     1129        Cursor/P A  yAxisWave  0 
     1130        Cursor/P B yAxisWave (numpnts(yAxisWave) - 1) 
     1131End 
     1132 
     1133//Fitting function for the POL-B standards 
     1134// 
     1135Function BStandardFunction(parameterWave, x) 
     1136        Wave parameterWave; Variable x 
     1137         
     1138        //Model parameters 
     1139        Variable KN=4.114E-3,CH=9.613E4, CD=7.558E4, NH=1872, ND=1556 
     1140        Variable INC=0.32, CHIV=2.2E-6 
     1141        //Correction based on absolute flux measured 5/93 
     1142        Variable CORR = 1.1445 
     1143         
     1144        //Local variables 
     1145        Variable AP2,QRGH,QRGD,IABS_RPA,IABS 
     1146         
     1147        //Calculate the function here 
     1148        ap2=parameterWave[1]^2 
     1149        qrgh = x*sqrt(nh*ap2/6) 
     1150        qrgd = x*sqrt(nd*ap2/6) 
     1151        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv) 
     1152        iabs = corr*iabs_rpa + inc 
     1153         
     1154        //return the result 
     1155        return parameterWave[0]*iabs 
     1156         
     1157End 
     1158 
     1159//Fitting function for the POL-C standards 
     1160// 
     1161Function CStandardFunction(parameterWave, x) 
     1162        Wave parameterWave; Variable x 
     1163         
     1164        //Model parameters 
     1165        Variable KN=4.114E-3,CH=2.564E5, CD=1.912E5, NH=4993, ND=3937 
     1166        Variable INC=0.32, CHIV=2.2E-6 
     1167        //Correction based on absolute flux measured 5/93 
     1168        Variable CORR = 1.0944 
     1169         
     1170        //Local variables 
     1171        Variable AP2,QRGH,QRGD,IABS_RPA,IABS 
     1172         
     1173        //Calculate the function here 
     1174        ap2=parameterWave[1]^2 
     1175        qrgh = x*sqrt(nh*ap2/6) 
     1176        qrgd = x*sqrt(nd*ap2/6) 
     1177        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv) 
     1178        iabs = corr*iabs_rpa + inc 
     1179         
     1180        //return the result 
     1181        return parameterWave[0]*iabs 
     1182         
     1183End 
     1184 
     1185//fitting function for the POL-AS standards 
     1186// 
     1187Function ASStandardFunction(parameterWave, x) 
     1188        Wave parameterWave; Variable x 
     1189         
     1190        //Model parameters 
     1191        Variable KN=64.5,CH=1.0, CD=1.0, NH=766, ND=766 
     1192        Variable INC=0.32, CHIV=0.0 
     1193         
     1194        //Local variables 
     1195        Variable AP2,QRGH,QRGD,IABS_RPA,IABS 
     1196         
     1197        //Calculate the function here 
     1198        ap2=parameterWave[1]^2 
     1199        qrgh = x*sqrt(nh*ap2/6) 
     1200 
     1201//The following lines were commented out in the fortran function 
     1202        //qrgd = x*sqrt(nd*ap2/6) 
     1203        //iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv) 
     1204 
     1205        iabs_rpa = kn*FIT_dbf(qrgh) 
     1206        iabs = iabs_rpa + inc 
     1207         
     1208        //return the result 
     1209        return parameterWave[0]*iabs 
     1210         
     1211End 
     1212 
     1213//Debye Function used for polymer standards 
     1214Function FIT_dbf(rgq) 
     1215        Variable rgq 
     1216        Variable x 
     1217         
     1218        x=rgq*rgq 
     1219        if (x < 5.0E-3) 
     1220                return 1.0 - x/3 + x^2/12 
     1221        else 
     1222                return 2*(exp(-x) + x - 1)/x^2 
     1223        endif 
     1224End 
     1225 
     1226 
     1227//////////////////////////////// EMD FIT RPA ////////////////////////////////////// 
     1228 
     1229 
     1230 
     1231 
     1232 
    7551233 
    7561234 
    7571235 
    7581236//procedures to clean up after itself 
     1237// 
     1238// -- note that it's never going to unload if the SANS Reduction is open too - since 
     1239// the old FIT_Ops was always loaded, and now points to this file. 
     1240// 
    7591241Function UnloadLinFit() 
    7601242 
    761         if (WinType("A_FitPanel") == 7) 
    762                 DoWindow/K A_FitPanel 
    763         endif 
    764         if (WinType("A_FitWindow") != 0) 
    765                 DoWindow/K $"A_FitWindow" 
     1243        if (WinType("FitPanel") == 7) 
     1244                DoWindow/K FitPanel 
     1245        endif 
     1246        if (WinType("FitWindow") != 0) 
     1247                DoWindow/K $"FitWindow" 
    7661248        endif 
    7671249        SetDataFolder root: 
    7681250        Killwaves/Z xAxisWave,yAxisWave,yErrWave,residWave,yWtWave,fit_yAxisWave 
    7691251         
    770         SVAR fileVerExt=root:Packages:NIST:SANS_ANA_EXTENSION 
     1252        SVAR/Z fileVerExt=root:Packages:NIST:SANS_ANA_EXTENSION 
     1253        if(SVAR_Exists(fileVerExt) == 0) 
     1254                return 0 
     1255        endif 
     1256         
    7711257        String fname="LinearizedFits" 
    7721258        Execute/P "DELETEINCLUDE \""+fname+fileVerExt+"\"" 
     
    7741260end 
    7751261 
     1262#if(exists("root:Packages:NIST:SANS_ANA_EXTENSION") != 0) 
     1263//this keeps the SANS Models menu from appearing in the SANS Reduction package 
    7761264Menu "SANS Models" 
    7771265        Submenu "Packages" 
     
    7801268end 
    7811269 
    782 Function A_FIT_PickPathButtonProc(ctrlName) : ButtonControl 
     1270#endif 
     1271 
     1272Function FIT_PickPathButtonProc(ctrlName) : ButtonControl 
    7831273        String ctrlName 
    7841274 
    7851275        A_PickPath() 
    7861276        //pop the file menu 
    787         A_FIT_FilePopMenuProc("",1,"") 
    788 End 
     1277        FIT_FilePopMenuProc("",1,"") 
     1278End 
  • sans/Release/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotManager_v40.ipf

    r681 r766  
    4444        Button button7,pos={238,20},size={100,20},proc=A_OneDLoader_NewFolderButton,title="New Folder" 
    4545        Button button7,help={"Select a new data folder"} 
     46        Checkbox check0,pos={240,190},title="Plot data on loading?",noproc,value=1 
    4647        GroupBox group0,pos={222,127},size={50,4},title="Shift-click to load" 
    4748        GroupBox group0_1,pos={222,143},size={50,4},title="multiple files" 
     
    297298                if(sel[ii] == 1) 
    298299                        fname=pathStr + fileWave[ii] 
    299                         Execute "A_LoadOneDDataWithName(\""+fname+"\",1)" 
     300                        Execute "A_LoadOneDDataWithName(\""+fname+"\","+num2str(doGraph)+")" 
    300301                        cnt += 1        //a file was loaded 
    301302                endif 
     
    341342        list = RemoveFromList(ListMatch(list,"*.GSP",";"), list, ";", 0) 
    342343        list = RemoveFromList(ListMatch(list,"*.MASK",";"), list, ";", 0) 
     344        #if(exists("QUOKKA") == 6) 
     345        list = RemoveFromList(ListMatch(list,"*.nx.hdf",";"), list, ";", 0)      
     346        list = RemoveFromList(ListMatch(list,"*.bin",";"), list, ";", 0)         
     347        #endif 
     348         
    343349 
    344350        //sort 
  • sans/Release/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtils2D_v40.ipf

    r665 r766  
    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// - OCT 2010   - added fDoBinning_QxQy2D(folderStr) that takes the QxQyQz data and bins it to I(Q) 
     14//                                      1D result must still be plotted manually, can be automated later. 
     15//                                      - for Scaled Image data (QxQy), the function fDoBinning_Scaled2D (in FFT_Cubes) could be 
     16//                                              modified and added to this file. It was meant for an FFT slice, but applies to  
     17//                                              2D model calculations (model_lin) 2D data files as well. 
     18 
     19 
     20// 
     21// changed June 2010 to read in the resolution information (now 8! columns) 
    1822// -- subject to change -- 
    1923// 
     24// look for either the old-style 3-column (no resolution information) or the newer 8-column format 
    2025Proc LoadQxQy() 
    2126 
     
    2530        Variable numCols = V_flag 
    2631 
    27         String w0,w1,w2,n0,n1,n2 
    28         String w3,w4,w5,w6,n3,n4,n5,n6 
    29                  
    30         // put the names of the three loaded waves into local names 
    31         n0 = StringFromList(0, S_waveNames ,";" ) 
    32         n1 = StringFromList(1, S_waveNames ,";" ) 
    33         n2 = StringFromList(2, S_waveNames ,";" ) 
    34         n3 = StringFromList(3, S_waveNames ,";" ) 
    35         n4 = StringFromList(4, S_waveNames ,";" ) 
    36         n5 = StringFromList(5, S_waveNames ,";" ) 
    37         n6 = StringFromList(6, S_waveNames ,";" ) 
    38          
    39         //remove the semicolon AND period from files from the VAX 
    40         w0 = CleanupName((S_fileName + "_qx"),0) 
    41         w1 = CleanupName((S_fileName + "_qy"),0) 
    42         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) 
    47  
    48         String baseStr=w1[0,strlen(w1)-4] 
    49         if(DataFolderExists("root:"+baseStr)) 
    50                         DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?" 
    51                         if(V_flag==2)   //user selected No, don't load the data 
    52                                 SetDataFolder root: 
    53                                 KillWaves $n0,$n1,$n2,$n3,$n4,$n5,$n6           // kill the default waveX that were loaded 
    54                                 return          //quits the macro 
    55                         endif 
    56                         SetDataFolder $("root:"+baseStr) 
    57         else 
    58                 NewDataFolder/S $("root:"+baseStr) 
    59         endif 
    60          
    61         //read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank) 
    62         Make/O/T/N=18 header 
    63         Variable refnum,ii 
    64         string tmpStr="" 
    65         Open/R refNum  as (path+filename) 
    66         ii=0 
    67         do 
    68                 tmpStr = "" 
    69                 FReadLine refNum, tmpStr 
    70                 header[ii] = tmpStr 
    71                 ii+=1 
    72         while(ii < 18)           
    73         Close refnum 
    74          
    75         // ? parse to get nice variable names? put these in a structure just for fun? 
    76          
    77          
    78         ////overwrite the existing data, if it exists 
    79         Duplicate/O $("root:"+n0), $w0 
    80         Duplicate/O $("root:"+n1), $w1 
    81         Duplicate/O $("root:"+n2), $w2 
    82         Duplicate/O $("root:"+n3), $w3 
    83         Duplicate/O $("root:"+n4), $w4 
    84         Duplicate/O $("root:"+n5), $w5 
    85         Duplicate/O $("root:"+n6), $w6 
     32        String w0,w1,w2,w3,w4,w5,w6,w7 
     33        String n0,n1,n2,n3,n4,n5,n6,n7 
     34                 
     35        if(numCols == 8) 
     36                // put the names of the 8 loaded waves into local names 
     37                n0 = StringFromList(0, S_waveNames ,";" ) 
     38                n1 = StringFromList(1, S_waveNames ,";" ) 
     39                n2 = StringFromList(2, S_waveNames ,";" ) 
     40                n3 = StringFromList(3, S_waveNames ,";" ) 
     41                n4 = StringFromList(4, S_waveNames ,";" ) 
     42                n5 = StringFromList(5, S_waveNames ,";" ) 
     43                n6 = StringFromList(6, S_waveNames ,";" ) 
     44                n7 = StringFromList(7, S_waveNames ,";" ) 
     45                 
     46                //remove the semicolon AND period from file names 
     47                w0 = CleanupName((S_fileName + "_qx"),0) 
     48                w1 = CleanupName((S_fileName + "_qy"),0) 
     49                w2 = CleanupName((S_fileName + "_i"),0) 
     50                w3 = CleanupName((S_fileName + "_iErr"),0) 
     51                w4 = CleanupName((S_fileName + "_qz"),0) 
     52                w5 = CleanupName((S_fileName + "_sQpl"),0) 
     53                w6 = CleanupName((S_fileName + "_sQpp"),0) 
     54                w7 = CleanupName((S_fileName + "_fs"),0) 
     55         
     56                String baseStr=w1[0,strlen(w1)-4] 
     57                if(DataFolderExists("root:"+baseStr)) 
     58                                DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?" 
     59                                if(V_flag==2)   //user selected No, don't load the data 
     60                                        SetDataFolder root: 
     61                                        KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7             // kill the default waveX that were loaded 
     62                                        return          //quits the macro 
     63                                endif 
     64                                SetDataFolder $("root:"+baseStr) 
     65                else 
     66                        NewDataFolder/S $("root:"+baseStr) 
     67                endif 
     68                 
     69                //read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank) 
     70                Make/O/T/N=18 header 
     71                Variable refnum,ii 
     72                string tmpStr="" 
     73                Open/R refNum  as (path+filename) 
     74                ii=0 
     75                do 
     76                        tmpStr = "" 
     77                        FReadLine refNum, tmpStr 
     78                        header[ii] = tmpStr 
     79                        ii+=1 
     80                while(ii < 18)           
     81                Close refnum             
     82                 
     83                ////overwrite the existing data, if it exists 
     84                Duplicate/O $("root:"+n0), $w0 
     85                Duplicate/O $("root:"+n1), $w1 
     86                Duplicate/O $("root:"+n2), $w2 
     87                Duplicate/O $("root:"+n3), $w3 
     88                Duplicate/O $("root:"+n4), $w4 
     89                Duplicate/O $("root:"+n5), $w5 
     90                Duplicate/O $("root:"+n6), $w6 
     91                Duplicate/O $("root:"+n7), $w7 
     92         
     93        endif           //8-columns 
     94         
     95        if(numCols == 3) 
     96                // put the names of the 3 loaded waves into local names 
     97                n0 = StringFromList(0, S_waveNames ,";" ) 
     98                n1 = StringFromList(1, S_waveNames ,";" ) 
     99                n2 = StringFromList(2, S_waveNames ,";" ) 
     100 
     101                //remove the semicolon AND period from file names 
     102                w0 = CleanupName((S_fileName + "_qx"),0) 
     103                w1 = CleanupName((S_fileName + "_qy"),0) 
     104                w2 = CleanupName((S_fileName + "_i"),0) 
     105                w3 = CleanupName((S_fileName + "_iErr"),0)              //make a name for the error wave, to be generated here 
     106 
     107                String baseStr=w1[0,strlen(w1)-4] 
     108                if(DataFolderExists("root:"+baseStr)) 
     109                                DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?" 
     110                                if(V_flag==2)   //user selected No, don't load the data 
     111                                        SetDataFolder root: 
     112                                        KillWaves/Z $n0,$n1,$n2         // kill the default waveX that were loaded 
     113                                        return          //quits the macro 
     114                                endif 
     115                                SetDataFolder $("root:"+baseStr) 
     116                else 
     117                        NewDataFolder/S $("root:"+baseStr) 
     118                endif 
     119                 
     120                //read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank) 
     121                Make/O/T/N=18 header 
     122                Variable refnum,ii 
     123                string tmpStr="" 
     124                Open/R refNum  as (path+filename) 
     125                ii=0 
     126                do 
     127                        tmpStr = "" 
     128                        FReadLine refNum, tmpStr 
     129                        header[ii] = tmpStr 
     130                        ii+=1 
     131                while(ii < 18)           
     132                Close refnum             
     133                 
     134                ////overwrite the existing data, if it exists 
     135                Duplicate/O $("root:"+n0), $w0 
     136                Duplicate/O $("root:"+n1), $w1 
     137                Duplicate/O $("root:"+n2), $w2 
     138         
     139 
     140 
     141                // generate my own error wave for I(qx,qy). This is exactly the same procedure that is used in the QxQy_Export function 
     142                Duplicate/O $("root:"+n0), $w3 
     143                $w3 = sqrt($w2)         //assumes Poisson statistics for each cell (counter) 
     144                //      sw = 0.05*sw            // uniform 5% error? tends to favor the low intensity too strongly 
     145                // get rid of the "bad" errorsby replacing the NaN, Inf, and zero with V_avg 
     146                // THIS IS EXTREMEMLY IMPORTANT - if this is not done, there are some "bad" values in the  
     147                // error wave (things that are not numbers) - and this wrecks the smeared model fitting. 
     148                // It appears to have no effect on the unsmeared model. 
     149                WaveStats/Q $w3 
     150                $w3 = numtype($w3[p]) == 0 ? $w3[p] : V_avg 
     151                $w3 = $w3[p] != 0 ? $w3[p] : V_avg 
     152         
     153        endif           //3-columns 
     154         
     155         
     156        /// do this for all 2D data, whether or not resolution information was read in 
    86157         
    87158        Variable/G gIsLogScale = 0 
     
    94165         
    95166        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 
    102167 
    103168        //clean up               
    104169        SetDataFolder root: 
    105         KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6 
     170        KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7 
     171         
    106172EndMacro 
    107173 
     
    518584// 
    519585// Currently the limitations are: 
    520 // - I have no error waves for the intensity 
    521 // - There is no smeared model 
     586// - I have no error waves for the intensity (fixed 10/2010) 
     587// - There is no smeared model (coming soon after 10/2010) 
    522588// - Cursors can't be used 
    523 // - the report *probably* won't work 
     589// - the report works OK, but I have little control over the graphics 
     590// - the mask is generated here with a default radius of 8 pixels around the beam center 
     591// 
    524592Function FitWrapper2D(folderStr,funcStr,coefStr,useCursors,useEps,useConstr) 
    525593        String folderStr,funcStr,coefStr 
    526594        Variable useCursors,useEps,useConstr 
    527595 
     596        //These only make sense for the 1D fits, but put them here so keep the look of the dispatching the same 
     597        Variable useResiduals, useTextBox 
     598        useResiduals = 0 
     599        useTextBox = 0 
     600         
    528601        String suffix=getModelSuffix(funcStr) 
    529602         
     
    545618         
    546619// fill a struct instance whether I need one or not 
     620// note that the resolution waves may or may not exist, and may or may not be used in the fitting 
    547621        String DF="root:"+folderStr+":"  
    548622         
    549         Struct ResSmearAAOStruct fs 
    550 //      WAVE resW = $(DF+folderStr+"_res")               
    551 //      WAVE fs.resW =  resW 
    552623        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 
     624        WAVE sw=$(DF+folderStr+"_iErr") 
     625        WAVE qx=$(DF+folderStr+"_qx") 
     626        WAVE qy=$(DF+folderStr+"_qy") 
     627        WAVE/Z qz=$(DF+folderStr+"_qz") 
     628        WAVE/Z sQpl=$(DF+folderStr+"_sQpl") 
     629        WAVE/Z sQpp=$(DF+folderStr+"_sQpp") 
     630        WAVE/Z shad=$(DF+folderStr+"_fs") 
     631 
     632//just a dummy - I shouldn't need this 
     633        Duplicate/O qx resultW 
     634        resultW=0 
     635         
     636        STRUCT ResSmear_2D_AAOStruct s 
     637        WAVE s.coefW = cw        
     638        WAVE s.zw = resultW      
     639        WAVE s.xw[0] = qx 
     640        WAVE s.xw[1] = qy 
     641        WAVE/Z s.qz = qz 
     642        WAVE/Z s.sQpl = sQpl 
     643        WAVE/Z s.sQpp = sQpp 
     644        WAVE/Z s.fs = shad 
     645         
    563646 
    564647        // generate my own mask wave - as a matrix first, then redimension to N*N vector 
     
    571654        Endif 
    572655         
    573  
    574          
    575 //      Duplicate/O yw $(DF+"FitYw") 
    576 //      WAVE fitYw = $(DF+"FitYw") 
    577 //      fitYw = NaN 
     656         
     657        Duplicate/O inten inten_masked 
     658        inten_masked = (mask[p][q] == 0) ? NaN : inten[p][q] 
     659         
    578660 
    579661        //for now, use res will always be 0 for 2D functions     
    580         Variable useRes=0 
     662        Variable useResol=0 
    581663        if(stringmatch(funcStr, "Smear*"))              // if it's a smeared function, need a struct 
    582                 useRes=1 
     664                useResol=1 
    583665        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 
     666 
     667        // can't use constraints in this way for multivariate fits. See the curve fitting help file 
     668        // and "Contraint Matrix and Vector" 
    588669        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 
    606  
     670                Print "Constraints not yet implemented" 
     671                useConstr = 0 
     672        endif    
     673        WAVE/Z constr=constr            //will be a null reference 
     674         
     675//      // do not construct constraints for any of the coefficients that are being held 
     676//      // -- this will generate an "unknown error" from the curve fitting 
     677//      Make/O/T/N=0 constr 
     678//      if(useConstr) 
     679//              String constraintExpression 
     680//              Variable i, nPnts=DimSize(lolim, 0),nextRow=0 
     681//              for (i=0; i < nPnts; i += 1) 
     682//                      if (strlen(lolim[i]) > 0 && hold[i] == 0) 
     683//                              InsertPoints nextRow, 1, constr 
     684//                              sprintf constraintExpression, "K%d > %s", i, lolim[i] 
     685//                              constr[nextRow] = constraintExpression 
     686//                              nextRow += 1 
     687//                      endif 
     688//                      if (strlen(hilim[i]) > 0 && hold[i] == 0) 
     689//                              InsertPoints nextRow, 1, constr 
     690//                              sprintf constraintExpression, "K%d < %s", i, hilim[i] 
     691//                              constr[nextRow] = constraintExpression 
     692//                              nextRow += 1 
     693//                      endif 
     694//              endfor 
     695//      endif 
     696 
     697        if(useCursors) 
     698                Print "Cursors not yet implemented" 
     699                useCursors = 0 
     700        endif    
    607701///// NO CURSORS for 2D waves 
    608702        //if useCursors, and the data is USANS, need to feed a (reassigned) trimmed matrix to the fit 
    609 //      Variable pt1,pt2,newN 
     703        Variable pt1,pt2,newN 
     704        pt1 = 0 
     705        pt2 = numpnts(inten)-1 
    610706//      if(useCursors && (dimsize(resW,1) > 4) ) 
    611707//              if(pcsr(A) > pcsr(B)) 
     
    631727// dispatch the fit 
    632728        //      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  
     729        Variable t1=StopMSTimer(-2) 
     730        Variable tb = 0         //no textbox 
     731 
     732// /NTHR=1 means just one thread for the fit (since the function itself is threaded) 
     733// NTHR = 0 == "Auto" mode, using as many processors as are available (not appropriate here since the function itself is threaded?) 
    633734 
    634735        do 
    635                 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 
     736         
     737                        // now useCursors, useEps, and useConstr are all handled w/ /NWOK, just like FitWrapper 
     738 
     739 
     740//              if(useResol && useResiduals && useTextBox)              //do it all 
     741//                      FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /R /NWOK 
     742//                      break 
     743//              endif 
     744//               
     745//              if(useResol && useResiduals)            //res + resid 
     746//                      FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /R /NWOK 
     747//                      break 
     748//              endif 
     749// 
     750//               
     751//              if(useResol && useTextBox)              //res + text 
     752//                      FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /NWOK 
     753//                      break 
     754//              endif 
     755                 
     756                if(useResol)            //res only 
     757                        Print "useRes only" 
     758                        FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /NWOK 
    637759                        break 
    638760                endif 
    639                  
    640                 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 
     761                         
     762                                 
     763/////   same as above, but all without useResol (no /STRC flag) 
     764                if(useResiduals && useTextBox)          //resid+ text 
     765                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /R /NWOK 
    642766                        break 
    643767                endif 
    644768                 
    645                 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 
     769                if(useResiduals)                //resid 
     770                        FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /R /NWOK 
    647771                        break 
    648772                endif 
    649                  
    650                 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 
     773 
     774                 
     775                if(useTextBox)          //text 
     776                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK 
    652777                        break 
    653778                endif 
    654779                 
    655                 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 
    657                         break 
    658                 endif 
    659                  
    660                 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 
    662                         break 
    663                 endif 
    664          
    665                 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 
    667                         break 
    668                 endif 
    669                  
    670                 if(useRes)              //just res 
    671                         FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten /X={Qx,Qy} /M=mask /W=sw /I=1 /D 
    672                         break 
    673                 endif 
    674                  
    675 /////   same as above, but all without useRes (no /STRC flag) 
    676                 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 
    678                         break 
    679                 endif 
    680                  
    681                 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 
    683                         break 
    684                 endif 
    685                  
    686                  
    687                 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 
    689                         break 
    690                 endif 
    691                  
    692                 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 
    694                         break 
    695                 endif 
    696                  
    697                 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 
    699                         break 
    700                 endif 
    701                  
    702                 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 
    704                         break 
    705                 endif 
    706          
    707                 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 
    709                         break 
    710                 endif 
    711                  
    712780                //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 
    714          
     781 
     782                FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK 
     783                 
    715784        while(0) 
    716785         
     786        Print "elapsed 2D fit time  = ",(StopMSTimer(-2) - t1)/1e6," s = ",(StopMSTimer(-2) - t1)/1e6/60," min" 
     787 
    717788        // append the fit 
    718789        // need to manage duplicate copies 
     
    735806        print w_sigma 
    736807        String resultStr="" 
     808         
     809        if(waveexists(W_sigma)) 
     810                //append it to the table, if it's not already there 
     811                CheckDisplayed/W=WrapperPanel#T0 W_sigma 
     812                if(V_flag==0) 
     813                        //not there, append it 
     814                        AppendtoTable/W=wrapperPanel#T0 W_sigma 
     815                else 
     816                        //remove it, and put it back on to make sure it's the right one (do I need to do this?) 
     817                        // -- not really, since any switch of the function menu takes W_Sigma off 
     818                endif 
     819        endif 
    737820                 
    738821        //now re-write the results 
     
    754837         
    755838        if(yesReport) 
    756                 String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //this is *hopefully* one wave 
     839                String parStr = getFunctionParams(funcStr) 
     840//              String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //old way, but doesn't work in 2D folders 
    757841                String topGraph= TopGizmoWindow()       //this is the topmost Gizmo (XOP) window 
    758842                 
     
    914998         
    915999End 
     1000 
     1001// This routine assumes that the 2D data was loaded with the NCNR loader, so that the 
     1002// data is in a data folder, and the extensions are known. A more generic form could  
     1003// be made too, if needed. 
     1004// 
     1005// X- need error on I(q) 
     1006// -- need to set "proper" number of data points (delta and qMax?) 
     1007// X- need to remove points at high Q end 
     1008// 
     1009// -- like the routines in CircSectAve, start with 500 points, and trim after binning is done. 
     1010//      you'l end up with < 200 points. 
     1011// 
     1012// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed 
     1013//  
     1014//Function fDoBinning_QxQy2D(inten,qx,qy,qz) 
     1015Function fDoBinning_QxQy2D(folderStr) 
     1016        String folderStr 
     1017 
     1018//      Wave inten,qx,qy,qz 
     1019 
     1020        SetDataFolder $("root:"+folderStr) 
     1021         
     1022        WAVE inten = $(folderStr + "_i") 
     1023        WAVE qx = $(folderStr + "_qx") 
     1024        WAVE qy = $(folderStr + "_qy") 
     1025        WAVE qz = $(folderStr + "_qz") 
     1026         
     1027        Variable xDim=numpnts(qx),yDim 
     1028        Variable ii,jj,delQ 
     1029        Variable qTot,nq,var,avesq,aveisq 
     1030        Variable binIndex,val 
     1031         
     1032        nq = 500 
     1033         
     1034        yDim = XDim 
     1035        Make/O/D/N=(nq) iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy 
     1036        delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2))               //use bins of 1 pixel width  
     1037        qBin_qxqy[] =  p*       delQ     
     1038        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning 
     1039 
     1040        iBin_qxqy = 0 
     1041        iBin2_qxqy = 0 
     1042        eBin_qxqy = 0 
     1043        nBin_qxqy = 0   //number of intensities added to each bin 
     1044         
     1045        for(ii=0;ii<xDim;ii+=1) 
     1046                qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2) 
     1047                binIndex = trunc(x2pnt(qBin_qxqy, qTot)) 
     1048                val = inten[ii] 
     1049                if (numType(val)==0)            //count only the good points, ignore Nan or Inf 
     1050                        iBin_qxqy[binIndex] += val 
     1051                        iBin2_qxqy[binIndex] += val*val 
     1052                        nBin_qxqy[binIndex] += 1 
     1053                endif 
     1054        endfor 
     1055 
     1056//calculate errors, just like in CircSectAve.ipf 
     1057        for(ii=0;ii<nq;ii+=1) 
     1058                if(nBin_qxqy[ii] == 0) 
     1059                        //no pixels in annuli, data unknown 
     1060                        iBin_qxqy[ii] = 0 
     1061                        eBin_qxqy[ii] = 1 
     1062                else 
     1063                        if(nBin_qxqy[ii] <= 1) 
     1064                                //need more than one pixel to determine error 
     1065                                iBin_qxqy[ii] /= nBin_qxqy[ii] 
     1066                                eBin_qxqy[ii] = 1 
     1067                        else 
     1068                                //assume that the intensity in each pixel in annuli is normally 
     1069                                // distributed about mean... 
     1070                                iBin_qxqy[ii] /= nBin_qxqy[ii] 
     1071                                avesq = iBin_qxqy[ii]^2 
     1072                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii] 
     1073                                var = aveisq-avesq 
     1074                                if(var<=0) 
     1075                                        eBin_qxqy[ii] = 1e-6 
     1076                                else 
     1077                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1)) 
     1078                                endif 
     1079                        endif 
     1080                endif 
     1081        endfor 
     1082         
     1083        // find the last non-zero point, working backwards 
     1084        val=nq 
     1085        do 
     1086                val -= 1 
     1087        while(nBin_qxqy[val] == 0) 
     1088         
     1089//      print val, nBin_qxqy[val] 
     1090        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy 
     1091         
     1092        SetDataFolder root: 
     1093         
     1094        return(0) 
     1095End 
  • sans/Release/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf

    r698 r766  
    9696//              Variable/G root:Packages:NIST:USANS_dQv = 0.117 
    9797//      endif 
    98                  
     98        String angst = StrVarOrDefault("root:Packages:NIST:gAngstStr", "A" ) 
     99 
    99100        // if no fileStr passed in, display dialog now 
    100101        if (cmpStr(fileStr,"") == 0) 
     
    302303                        Duplicate/O $("root:"+n2), $w2 
    303304                        Duplicate/O $("root:"+n3), $w3 
    304                         Duplicate/O $("root:"+n0), $w0 // Set qb wave to nominal measured Q values 
     305                        Duplicate/O $("root:"+n0), $w4 // Set qb wave to nominal measured Q values 
    305306                        Duplicate/O $("root:"+n0), $w5 // Make wave of appropriate length 
     307                        $w4 = $w0 
    306308                        $w5 = 1                                           //  Set all shadowfactor to 1 
    307309         
     
    415417                                        ModifyGraph tickUnit(left)=1 
    416418                                        Label left "I(q)" 
    417                                         Label bottom "q (A\\S-1\\M)" 
     419                                        Label bottom "q ("+angst+"\\S-1\\M)" 
    418420                                        Legend 
    419421                                endif 
     
    426428                                ModifyGraph tickUnit(left)=1 
    427429                                Label left "I(q)" 
    428                                 Label bottom "q (A\\S-1\\M)" 
     430                                Label bottom "q ("+angst+"\\S-1\\M)" 
    429431                                Legend 
    430432                        endif 
     
    10951097 
    10961098        // SANS Reduction bits 
    1097         tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;Monte_SANS_W3;Monte_SANS_W4;Monte_SANS;FractionReachingDetector;" 
    1098         list = RemoveFromList(tmp, list  ,";") 
     1099        tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;" 
     1100        list = RemoveFromList(tmp, list  ,";") 
     1101        tmp = "NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;Monte_SANS_W3;Monte_SANS_W4;Monte_SANS;FractionReachingDetector;TwoLevel_EC;SmearedTwoLevel_EC;" 
     1102        list = RemoveFromList(tmp, list  ,";") 
     1103 
    10991104 
    11001105        // USANS Reduction bits 
     
    12481253         
    12491254        /// items for everyone 
    1250         val = NumVarOrDefault("root:Packages:NIST:gXML_Write", 1 ) 
     1255        val = NumVarOrDefault("root:Packages:NIST:gXML_Write", 0 ) 
    12511256        Variable/G root:Packages:NIST:gXML_Write = val 
    12521257         
  • sans/Release/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf

    r570 r766  
    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