source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/USANS_SlitSmearing.ipf @ 153

Last change on this file since 153 was 145, checked in by ajj, 15 years ago

Changes:

  • Made USANS weighting matrix double precision
  • Changed Smear_Model_(5,10,20,76) to be AAO and added test for presence of USANS matrix. If the number of columns in the resolution wave is > 4 then do the matrix multiplication, otherwise do gaussian quadrature via Smear_Model_N.
  • Changed Sphere.ipf and SchulzSpheres?.ipf to use AAO form of Smear_Model_20
File size: 3.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=3.00
3#pragma IgorVersion=6.0
4
5//Functions for doing USANS Slit smearing by method of weight matrix
6//Routines originally from J Barker fortran code
7//Translated to IGOR by M-H Kim
8//Updated to use IGOR features and integrated into SANS Macros by A J Jackson
9//
10//Calculation of weights takes a while, but then each recalculation of the smeared model takes no longer
11//than the unsmeared calculation.
12//
13//Should be used when just USANS data is being fitted. Fitting of combined SANS/USANS must
14//still be done by integral method as per S Kline macros.
15//This can be avoided if IGOR Global fit is used - different functions can be fitted to the two data sets
16//with variables tied across the two. There is a loss of simplicity for the user however.
17//
18
19// AJJ
20// July 26 2007 - Modified functions to work with new SANS Analysis Macros
21//     - pass basestr to functions to determine wave names and avoid globals
22//     - pass dQv to functions to  avoid globals
23//     - pass N to CalcR to avoid globals
24
25Function USANS_CalcWeights(basestr, dQv)
26        String basestr
27        Variable dQv
28
29        Variable/G USANS_N=numpnts($(basestr+"_q"))
30        Variable/G USANS_dQv = dQv
31        String/G USANS_basestr = basestr
32
33        Make/D/O/N=(USANS_N,USANS_N) $(basestr+"_res")
34        //Make/O/N=(N,N) W1mat
35        //Make/O/N=(N,N) W2mat
36        //Make/O/N=(N,N) Rmat
37        Wave weights = $(basestr+"_res")
38
39        Variable/G USANS_m = EnterSlope()
40
41        Variable tref = startMSTimer
42        print "Calculating W1..."
43        weights = (p <= q ) && (q < USANS_N-1) ? CalcW1(p,q)  : 0
44        print "Calculating W2..."
45        weights += (p+1 <= q ) && (q < USANS_N) ?  CalcW2(p,q) : 0
46        print "Calculating Remainders..."
47        weights += (q == USANS_N-1) ? CalcR(p) : 0
48//      print "Summing weights..."
49//      Weights = W1mat + W2mat + Rmat
50        print "Done"
51        Variable ms = stopMSTimer(tref)
52        print "Time elapsed = ", ms/1e6, "s"
53        //return Weights
54End
55
56// This is far from satisfactory!
57Function EnterSlope()
58       
59        Variable slope=-4
60
61        Prompt slope "Enter a slope"
62        DoPrompt "Enter Slope", slope
63                If (V_Flag)
64                        return -1
65                Endif   
66        print "slope=", slope
67        return slope
68       
69End
70
71Function CalcW1(i,j)
72
73        Variable i,j
74        SVAR USANS_basestr
75        NVAR dQv = USANS_dQv
76       
77        Variable UU,UL,dqj,rU,rL,wU,wL,dqw
78        Wave Qval = $(USANS_basestr+"_q")
79       
80        UU =sqrt(Qval[j+1]^2-Qval[i]^2)
81        UL = sqrt(Qval[j]^2-Qval[i]^2)
82        dqw = Qval[j+1]-Qval[j]
83        rU = sqrt(UU^2+Qval[i]^2)
84        rL = sqrt(UL^2+Qval[i]^2)
85       
86        wU = (1.0/dQv)*(Qval[j+1]*UU/dqw - 0.5*UU*rU/dqw - 0.5*Qval[i]^2*ln(UU+rU)/dqw )
87        wL = (1.0/dQv)*(Qval[j+1]*UL/dqw - 0.5*UL*rL/dqw - 0.5*Qval[i]^2*ln(UL+rL)/dqw )
88       
89        Return wU-wL
90
91End
92
93Function CalcW2(i,j)
94
95        Variable i,j
96       
97        SVAR USANS_basestr
98        NVAR dQv = USANS_dQv
99       
100        variable UU,UL,dqw,rU,rL,wU,wL
101       
102        Wave Qval = $(USANS_basestr+"_q")
103
104        UU = sqrt(Qval[j]^2-Qval[i]^2)                 
105        UL = sqrt(Qval[j-1]^2-Qval[i]^2)               
106        dqw = Qval[j]-Qval[j-1]                 
107        rU = sqrt(UU^2+Qval[i]^2)
108        rL = sqrt(UL^2+Qval[i]^2)
109        wU = (1.0/dQv)*( -Qval[j-1]*UU/dqw + 0.5*UU*rU/dqw + 0.5*Qval[i]^2*ln(UU+rU)/dqw )
110        wL = (1.0/dQv)*( -Qval[j-1]*UL/dqw + 0.5*UL*rL/dqw + 0.5*Qval[i]^2*ln(UL+rL)/dqw )
111
112        Return wU-wL
113
114End
115
116Function CalcR(i)
117
118        Variable i
119
120        SVAR USANS_basestr
121        NVAR N = USANS_N
122        NVAR dQv = USANS_dQv
123       
124        Variable retval
125        Wave Qval = $(USANS_basestr+"_q")
126        Variable/G USANS_intQpt = Qval[i]
127       
128        Variable lower = sqrt(qval[N-1]^2-qval[i]^2)
129        Variable upper = lower +dQv
130       
131        retval = Integrate1D(Remainder,lower,upper)
132       
133        Return retval
134
135End
136
137Function Remainder(i)
138       
139        Variable i
140       
141        SVAR USANS_basestr
142        NVAR m = USANS_m
143        NVAR qi = USANS_intQpt
144        NVAR N = USANS_N
145        WAVE Qval = $(USANS_basestr+"_q")       
146        Variable retVal
147       
148        retVal=Qval[N-1]^(-m)*(i^2+qi^2)^(m/2)
149
150        return retval
151
152End
Note: See TracBrowser for help on using the repository browser.