source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/XOP/ResolutionSmearing.c @ 102

Last change on this file since 102 was 97, checked in by ajj, 16 years ago

Now committing the code correctly - having relied on XCode's SVN interface which doesn't behave (quelle surprise).

I hope this isn't too screwed up.
  • Property svn:executable set to *
File size: 7.0 KB
Line 
1/*      ResolutionSmearing.c
2
3Generalized smearing routines
4
5*/
6
7#pragma XOP_SET_STRUCT_PACKING                  // All structures are 2-byte-aligned.
8
9#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
10#include "SANSAnalysis.h"
11#include "GaussWeights.h"
12#include "ResolutionSmearing.h"
13
14int
15Smear_Model_20_X(GenSmearParamsPtr p)
16{
17        PassParams pp;
18        FunctionInfo fi;
19        int badParameterNumber;
20        int requiredParameterTypes[2];
21        int err;
22        double ans;
23        //char buf[256];
24       
25        //Get name of function into string here with
26        char func[255];
27        GetCStringFromHandle(p->funcname, func, sizeof(func));
28       
29        //Make sure function that we want to call exists
30        //GetFunctionInfo(func,&fi);
31        if (err = GetFunctionInfo(func,&fi)){
32                //      if (err = GetFunctionInfo("Cyl_PolyRadiusX",&fi)){
33                XOPNotice("Error at GetFunctionInfo\r");
34                return err;
35                }
36        requiredParameterTypes[1] = NT_FP64;
37        requiredParameterTypes[0] = WAVE_TYPE;
38        if (err = CheckFunctionForm(&fi, 2, requiredParameterTypes, &badParameterNumber, NT_FP64)){
39                XOPNotice("Error at CheckFunctionForm\r");
40                return err;
41        }
42        int ii;
43        double nord,va,vb,summ,yyy,zi,Pi;
44        double answer,Resoln,i_shad,i_qbar,i_sigq;
45       
46        Pi = 4.0*atan(1.0);
47       
48        i_shad = p->i_shad;
49        i_qbar = p->i_qbar;
50        i_sigq = p->i_sigq;
51       
52        //sprintf(buf, "i_shad = %g, i_qbar = %g, i_sigq = %g", i_shad, i_qbar, i_sigq);
53        //XOPNotice(buf);
54       
55        pp.waveHandle = p->waveHandle;
56        pp.x = p->x;
57        //need interpolation
58       
59        if (i_sigq >= 0){
60               
61                nord = 20;
62                va = -3*i_sigq + i_qbar;
63                if (va < 0){
64                        va = 0;
65                }
66                vb = 3*i_sigq + i_qbar;
67               
68                summ = 0.0;
69                ii=0;
70                while (ii < nord) {
71                        zi = (Gauss20Z[ii]*(vb-va) + va + vb)/2.0;
72                        pp.x = zi;
73                        Resoln = i_shad/sqrt(2*Pi*i_sigq*i_sigq);
74                        Resoln *= exp((-1*pow((zi-i_qbar),2))/(2*i_sigq*i_sigq));
75                        if (err = CallFunction(&fi, (void*)&pp, &ans)){
76                                XOPNotice("Error at CallFunction in Gaussian Smear\r");
77                                return err;
78                        }
79                        //                      CallFunction(&fi, (void*)&pp, &ans);
80                        yyy = Gauss20Wt[ii]*Resoln*ans;
81                        summ += yyy;
82                       
83                        ii++;
84                } 
85                answer = (vb-va)/2.0*summ;
86        } else {
87                //USANS so do trapezoidal integration
88                int maxiter = 20;
89                double tol = 1e-4;
90               
91                va=0;
92                vb=fabs(i_sigq);
93               
94                //sprintf(buf, "i_sigq = %g, va = %g, vb = %g\r",i_sigq, va, vb);
95                //XOPNotice("USANS Data. Great\r");
96                //XOPNotice(buf);
97               
98                answer = qtrap_USANS(fi,pp,va,vb,tol,maxiter);
99                answer /= vb;
100        }
101       
102       
103        p->result = answer;
104       
105        return 0;
106        }
107
108int Smear_Model_76_X(GenSmearParamsPtr p)
109{
110        PassParams pp;
111        FunctionInfo fi;
112        int badParameterNumber;
113        int requiredParameterTypes[2];
114        int err;
115        double ans;
116        //char buf[256];
117       
118        //Get name of function into string here with
119        char func[255];
120        GetCStringFromHandle(p->funcname, func, sizeof(func));
121       
122        //Make sure function that we want to call exists
123        //GetFunctionInfo(func,&fi);
124        if (err = GetFunctionInfo(func,&fi)){
125                //      if (err = GetFunctionInfo("Cyl_PolyRadiusX",&fi)){
126                XOPNotice("Error at GetFunctionInfo\r");
127                return err;
128                }
129        requiredParameterTypes[1] = NT_FP64;
130        requiredParameterTypes[0] = WAVE_TYPE;
131        if (err = CheckFunctionForm(&fi, 2, requiredParameterTypes, &badParameterNumber, NT_FP64)){
132                XOPNotice("Error at CheckFunctionForm\r");
133                return err;
134        }
135        int ii;
136        double nord,va,vb,summ,yyy,zi,Pi;
137        double answer,Resoln,i_shad,i_qbar,i_sigq;
138       
139        Pi = 4.0*atan(1.0);
140       
141        i_shad = p->i_shad;
142        i_qbar = p->i_qbar;
143        i_sigq = p->i_sigq;
144       
145        //sprintf(buf, "i_shad = %g, i_qbar = %g, i_sigq = %g", i_shad, i_qbar, i_sigq);
146        //XOPNotice(buf);
147       
148        pp.waveHandle = p->waveHandle;
149        pp.x = p->x;
150        //need interpolation
151       
152        if (i_sigq >= 0){
153               
154                nord = 20;
155                va = -3*i_sigq + i_qbar;
156                if (va < 0){
157                        va = 0;
158                }
159                vb = 3*i_sigq + i_qbar;
160               
161                summ = 0.0;
162                ii=0;
163                while (ii < nord) {
164                        zi = (Gauss76Z[ii]*(vb-va) + va + vb)/2.0;
165                        pp.x = zi;
166                        Resoln = i_shad/sqrt(2*Pi*i_sigq*i_sigq);
167                        Resoln *= exp((-1*pow((zi-i_qbar),2))/(2*i_sigq*i_sigq));
168                        if (err = CallFunction(&fi, (void*)&pp, &ans)){
169                                XOPNotice("Error at CallFunction in Gaussian Smear\r");
170                                return err;
171                        }
172                        //                      CallFunction(&fi, (void*)&pp, &ans);
173                        yyy = Gauss76Wt[ii]*Resoln*ans;
174                        summ += yyy;
175                       
176                        ii++;
177                } 
178                answer = (vb-va)/2.0*summ;
179        } else {
180                //USANS so do trapezoidal integration
181                int maxiter = 20;
182                double tol = 1e-4;
183               
184                va=0;
185                vb=fabs(i_sigq);
186               
187                //sprintf(buf, "i_sigq = %g, va = %g, vb = %g\r",i_sigq, va, vb);
188                //XOPNotice("USANS Data. Great\r");
189                //XOPNotice(buf);
190               
191                answer = qtrap_USANS(fi,pp,va,vb,tol,maxiter);
192                answer /= vb;
193        }
194       
195       
196        p->result = answer;
197       
198        return 0;
199        }
200
201
202double
203qtrap_USANS(FunctionInfo fi, PassParams p, double aa, double bb, double eps, int maxit){
204        double olds,ss;
205        double xx,tnm,summ,del,arg1,arg2,ans;
206        int it,jj,kk;
207        int err;
208       
209        //char buf[256];
210       
211        double qval = p.x;
212        double temp = 0;
213       
214        //XOPNotice("In qtrap_USANS\r");
215       
216        olds = -1e-30;
217        for (jj = 1; jj <= maxit; jj++){
218               
219                if (jj == 1){
220                        arg1 = pow(qval,2)+pow(aa,2);
221                        arg2 = pow(qval,2)+pow(bb,2);
222                        p.x = sqrt(arg1);
223                        if (err = CallFunction(&fi, (void*)&p, &ans)){
224                                XOPNotice("Error at 1st CallFunction in qtrap\r");
225                        }
226                        temp = ans;
227                        p.x = sqrt(arg2);
228                        CallFunction(&fi, (void*)&p, &ans);
229                        temp += ans;
230                        ss = 0.5*(bb-aa)*temp;
231                        //sprintf(buf, "jj = 1: aa = %g, bb = %g, temp = %g, ss = %g\r", aa, bb, temp, ss);
232                        //XOPNotice(buf);
233                } else {
234                        it = pow(2,jj-2);
235                        tnm = (double)it;
236                        del = (bb-aa)/tnm;
237                        xx = aa+(0.5*del);
238                        summ = 0;
239                        for (kk = 1; kk <= it; kk++){
240                                arg1 = pow(qval,2)+pow(xx,2);
241                                p.x = sqrt(arg1);
242                                if (err = CallFunction(&fi, (void*)&p, &ans)){
243                                        XOPNotice("Error at 2nd CallFunction in qtrap\r");
244                                }
245                                summ += ans;
246                                xx += del;
247                        }
248                        ss = 0.5*(ss+(bb-aa)*summ/tnm);
249                        //sprintf(buf, "jj = %d: arg1 = %g, summ = %g, tnm= %g, ss = %g\r", jj, arg1, summ, tnm, ss);
250                        //XOPNotice(buf);
251                }
252               
253                //sprintf(buf,"ss = %g\r",ss);
254                //XOPNotice(buf);
255               
256                //              ss = trapzd_USANS(fi,p,aa,bb,jj);
257                if ( fabs(ss-olds) < eps*fabs(olds))
258                        return ss;
259                if ( (ss == 0.0) && (olds == 0.0) && (jj > 6) ) // no progress?
260                        return ss;
261                olds = ss;
262        }
263       
264        XOPNotice("Maxit exceeded in qtrap. If you're here, there was a problem in qtrap");
265        return (ss); // should never get here if function is well behaved
266       
267}
268
269//static double
270//trapzd_USANS(FunctionInfo fi, PassParams p, double aa, double bb, int nn){
271//      double xx,tnm,summ,del,arg1,arg2,ans;
272//      int it,jj;
273//     
274//      static double sval;
275//     
276//      double qval = p.x;
277//      double temp = 0;
278//
279//      //XOPNotice("In trapzd_USANS\r");
280//
281//      if (nn == 1){
282//              arg1 = pow(qval,2)+pow(aa,2);
283//              arg2 = pow(qval,2)+pow(bb,2);
284//              p.x = arg1;
285//              CallFunction(&fi, (void*)&p, &ans);
286//              p.x = arg2;
287//              temp = ans;
288//              CallFunction(&fi, (void*)&p, &ans);
289//              temp += ans;
290//              sval = 0.5*(bb-aa);
291//              return sval;
292//      } else {
293//              it = pow(2,nn-2);
294//              tnm = it;
295//              del = (bb-aa)/tnm;
296//              xx = aa+0.5*del;
297//              summ = 0;
298//              for (jj = 1; jj <= it; jj++){
299//                      arg1 = pow(qval,2)+pow(xx,2);
300//                      p.x = arg1;
301//                      CallFunction(&fi, (void*)&p, &ans);
302//                      summ += ans;
303//                      xx += del;
304//              }
305//              sval = 0.5*(sval+(bb-aa)*summ/tnm);
306//              return sval;
307//      }
308//     
309//
310//}
311
312
313#pragma XOP_RESET_STRUCT_PACKING                        // All structures are 2-byte-aligned.
314                                                                                        // All structures are 2-byte-aligned.
Note: See TracBrowser for help on using the repository browser.