Ignore:
Timestamp:
May 15, 2015 11:35:18 AM (8 years ago)
Author:
srkline
Message:

some reorganization of the r/w routines to generate HDF test files for SANS and VSANS (all are housed together for testing). also some reorganization of the detector binning routines to get functions grouped in more logical locations.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Utils.ipf

    r954 r955  
    11#pragma rtGlobals=3             // Use modern global access method and strict wave access. 
     2 
     3///////////////////////// 
     4// 
     5// Utility functions to: 
     6//              calculate Q, Qx, Qy, Qz 
     7//              fill the detector panels with simulated data (the model functions are here) 
     8//              bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width) 
     9// 
     10///////////////////////// 
     11 
     12 
    213 
    314 
     
    3142        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to  
    3243        //  keep the panel to keep people used to using it.) 
     44        // peak @ 0.1 ~ AgBeh 
     45        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}                
     46        // 
     47        // peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved 
     48        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}              
    3349        String funcStr = VCALC_getModelFunctionStr() 
    3450        strswitch(funcStr) 
     
    4561                        tmpInten = V_SphereForm(1,60,1e-6,0.001,qTot[p][q])      
    4662                        break 
     63                case "AgBeh": 
     64                        tmpInten = V_BroadPeak(1e-9,3,20,100.0,0.1,3,0.1,qTot[p][q]) 
     65                        break 
     66                case "Vycor": 
     67                        tmpInten = V_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q]) 
     68                        break    
     69                case "Empty Cell": 
     70                        tmpInten = V_EC_Empirical(2.2e-8,3.346,0.0065,9.0,0.016,qTot[p][q]) 
     71                        break 
     72                case "Blocked Beam": 
     73                        tmpInten = V_BlockedBeam(1,qTot[p][q]) 
     74                        break 
    4775                default: 
    4876                        tmpInten = V_Debye(10,300,0.1,qTot[p][q]) 
     
    148176 
    149177//function to calculate the overall q-value, given all of the necesary trig inputs 
    150 //NOTE: detector locations passed in are pixels = 0.5cm real space on the detector 
    151178//and are in detector coordinates (1,128) rather than axis values 
    152179//the pixel locations need not be integers, reals are ok inputs 
     
    177204//ALL inputs are in detector coordinates 
    178205// 
    179 //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
    180206//sdd is in meters 
    181207//wavelength is in Angstroms 
     
    208234//input/output is the same as CalcQval() 
    209235//ALL inputs are in detector coordinates 
    210 //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
    211236//sdd is in meters 
    212237//wavelength is in Angstroms 
     
    239264//input/output is the same as CalcQval() 
    240265//ALL inputs are in detector coordinates 
    241 //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector 
    242266//sdd is in meters 
    243267//wavelength is in Angstroms 
     
    366390End 
    367391 
     392// a sum of a power law and debye to approximate the scattering from a real empty cell 
     393// 
     394//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016} 
     395// 
     396Function V_EC_Empirical(aa,mm,scale,rg,bkg,x) 
     397        Variable aa,mm,scale,rg,bkg 
     398        Variable x 
     399         
     400        // variables are: 
     401        //[0] = A 
     402        //[1] = power m 
     403        //[2] scale factor 
     404        //[3] radius of gyration [A] 
     405        //[4] background        [cm-1] 
     406         
     407        Variable Iq 
     408         
     409        // calculates (scale*debye)+bkg 
     410        Variable Pq,qr2 
     411         
     412//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values 
     413//              return(scale+bkg) 
     414//      endif 
     415         
     416        Iq = aa*x^-mm 
     417         
     418        qr2=(x*rg)^2 
     419        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2 
     420         
     421        //scale 
     422        Pq *= scale 
     423        // then add the terms up 
     424        return (Iq + Pq + bkg) 
     425End 
     426 
     427// blocked beam 
     428// 
     429Function V_BlockedBeam(bkg,x) 
     430        Variable bkg 
     431        Variable x 
     432         
     433        return (bkg) 
     434End 
     435 
     436 
     437// 
     438// a broad peak to simulate silver behenate or vycor 
     439// 
     440// peak @ 0.1 ~ AgBeh 
     441//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}                
     442// 
     443// 
     444// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved 
     445//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}              
     446// 
     447// 
     448Function V_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x) 
     449        Variable aa,nn,cc,LL,Qzero,mm,bgd 
     450        Variable x 
     451         
     452        // variables are:                                                        
     453        //[0] Porod term scaling 
     454        //[1] Porod exponent 
     455        //[2] Lorentzian term scaling 
     456        //[3] Lorentzian screening length [A] 
     457        //[4] peak location [1/A] 
     458        //[5] Lorentzian exponent 
     459        //[6] background 
     460         
     461//      local variables 
     462        Variable inten, qval 
     463//      x is the q-value for the calculation 
     464        qval = x 
     465//      do the calculation and return the function value 
     466         
     467        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd 
     468 
     469        Return (inten) 
     470         
     471End 
    368472 
    369473 
     
    429533 
    430534 
     535 
     536////////// 
     537// 
     538//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel 
     539//              - each pixel QxQyQz has been calculated beforehand 
     540//              - if multiple panels are selected to be combined, it is done here during the binning 
     541//              - the setting of deltaQ step is still a little suspect (TODO) 
     542// 
     543// 
    431544// see the equivalent function in PlotUtils2D_v40.ipf 
    432545// 
     
    436549// -- type = FL or FR or...other panel identifiers 
    437550// 
    438 // TODO "iErr" is all messed up since it doesn't really apply here for data that is not 2D simulation 
    439 // 
     551// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation 
    440552// 
    441553Function V_fDoBinning_QxQy2D(folderStr,type) 
     
    592704 
    593705//TODO: properly define the errors here - I'll have this if I do the simulation 
    594         if(WaveExists(iErr)==0) 
     706        if(WaveExists(iErr)==0  && WaveExists(inten) != 0) 
    595707                Duplicate/O inten,iErr 
    596708                Wave iErr=iErr 
    597709//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...) 
    598710                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value 
     711        endif 
     712        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0) 
     713                Duplicate/O inten2,iErr2 
     714                Wave iErr2=iErr2 
     715//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...) 
     716                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value 
     717        endif 
     718        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0) 
     719                Duplicate/O inten3,iErr3 
     720                Wave iErr3=iErr3 
     721//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...) 
     722                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value 
     723        endif 
     724        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0) 
     725                Duplicate/O inten4,iErr4 
     726                Wave iErr4=iErr4 
     727//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...) 
     728                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value 
    599729        endif 
    600730 
     
    622752         
    623753//      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  
    624 // TODO: not sure if I want to so dQ in x or y direction... 
     754// TODO: not sure if I want to set dQ in x or y direction... 
    625755        // the short dimension is the 8mm tubes, use this direction as dQ? 
    626756        // but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0] 
     
    680810                                        iBin_qxqy[binIndex] += val 
    681811                                        iBin2_qxqy[binIndex] += val*val 
    682                                         eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj] 
     812                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj] 
    683813                                        nBin_qxqy[binIndex] += 1 
    684814                                endif 
     
    702832                                        iBin_qxqy[binIndex] += val 
    703833                                        iBin2_qxqy[binIndex] += val*val 
    704                                         eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj] 
     834                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj] 
    705835                                        nBin_qxqy[binIndex] += 1 
    706836                                endif 
     
    721851                                        iBin_qxqy[binIndex] += val 
    722852                                        iBin2_qxqy[binIndex] += val*val 
    723                                         eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj] 
     853                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj] 
    724854                                        nBin_qxqy[binIndex] += 1 
    725855                                endif 
Note: See TracChangeset for help on using the changeset viewer.