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

Last change on this file since 245 was 245, checked in by ajj, 15 years ago
  • Property svn:executable set to *
File size: 8.2 KB
Line 
1/*      ResolutionSmearing.c
2
3Generalized smearing routines
4
5*/
6
7
8#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
9#include "SANSAnalysis.h"
10#include "GaussWeights.h"
11#include "ResolutionSmearing.h"
12#include "libSANSAnalysis.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        int ii;
26        double nord,va,vb,summ,yyy,zi,Pi;
27        double answer,Resoln,i_shad,i_qbar,i_sigq;
28       
29        //Get name of function into string here with
30        char func[255];
31        GetCStringFromHandle(p->funcname, func, sizeof(func));
32       
33        //Make sure function that we want to call exists
34        //GetFunctionInfo(func,&fi);
35        if (err = GetFunctionInfo(func,&fi)){
36                //      if (err = GetFunctionInfo("Cyl_PolyRadiusX",&fi)){
37                XOPNotice("Error at GetFunctionInfo\r");
38                return err;
39                }
40        requiredParameterTypes[1] = NT_FP64;
41        requiredParameterTypes[0] = WAVE_TYPE;
42        if (err = CheckFunctionForm(&fi, 2, requiredParameterTypes, &badParameterNumber, NT_FP64)){
43                XOPNotice("Error at CheckFunctionForm\r");
44                return err;
45        };
46
47
48       
49        Pi = 4.0*atan(1.0);
50       
51        i_shad = p->i_shad;
52        i_qbar = p->i_qbar;
53        i_sigq = p->i_sigq;
54       
55        //sprintf(buf, "i_shad = %g, i_qbar = %g, i_sigq = %g", i_shad, i_qbar, i_sigq);
56        //XOPNotice(buf);
57       
58        pp.waveHandle = p->waveHandle;
59        pp.x = p->x;
60        //need interpolation
61       
62        if (i_sigq >= 0){
63               
64                nord = 20;
65                va = -3*i_sigq + i_qbar;
66                if (va < 0){
67                        va = 0;
68                }
69                vb = 3*i_sigq + i_qbar;
70               
71                summ = 0.0;
72                ii=0;
73                while (ii < nord) {
74                        zi = (Gauss20Z[ii]*(vb-va) + va + vb)/2.0;
75                        pp.x = zi;
76                        Resoln = i_shad/sqrt(2*Pi*i_sigq*i_sigq);
77                        Resoln *= exp((-1*pow((zi-i_qbar),2))/(2*i_sigq*i_sigq));
78                        if (err = CallFunction(&fi, (void*)&pp, &ans)){
79                                XOPNotice("Error at CallFunction in Gaussian Smear\r");
80                                return err;
81                        }
82                        //                      CallFunction(&fi, (void*)&pp, &ans);
83                        yyy = Gauss20Wt[ii]*Resoln*ans;
84                        summ += yyy;
85                       
86                        ii++;
87                } 
88                answer = (vb-va)/2.0*summ;
89        } else {
90                //USANS so do trapezoidal integration
91                int maxiter = 20;
92                double tol = 1e-4;
93               
94                va=0;
95                vb=fabs(i_sigq);
96               
97                //sprintf(buf, "i_sigq = %g, va = %g, vb = %g\r",i_sigq, va, vb);
98                //XOPNotice("USANS Data. Great\r");
99                //XOPNotice(buf);
100               
101                answer = qtrap_USANS(fi,pp,va,vb,tol,maxiter);
102                answer /= vb;
103        }
104       
105       
106        p->result = answer;
107       
108        return 0;
109        }
110
111int Smear_Model_76_X(GenSmearParamsPtr p)
112{
113        PassParams pp;
114        FunctionInfo fi;
115        int badParameterNumber;
116        int requiredParameterTypes[2];
117        int err;
118        double ans;
119        //char buf[256];
120       
121        int ii;
122        double nord,va,vb,summ,yyy,zi,Pi;
123        double answer,Resoln,i_shad,i_qbar,i_sigq;
124
125        //Get name of function into string here with
126        char func[255];
127        GetCStringFromHandle(p->funcname, func, sizeof(func));
128       
129        //Make sure function that we want to call exists
130        //GetFunctionInfo(func,&fi);
131        if (err = GetFunctionInfo(func,&fi)){
132                //      if (err = GetFunctionInfo("Cyl_PolyRadiusX",&fi)){
133                XOPNotice("Error at GetFunctionInfo\r");
134                return err;
135                }
136        requiredParameterTypes[1] = NT_FP64;
137        requiredParameterTypes[0] = WAVE_TYPE;
138        if (err = CheckFunctionForm(&fi, 2, requiredParameterTypes, &badParameterNumber, NT_FP64)){
139                XOPNotice("Error at CheckFunctionForm\r");
140                return err;
141        }
142
143       
144        Pi = 4.0*atan(1.0);
145       
146        i_shad = p->i_shad;
147        i_qbar = p->i_qbar;
148        i_sigq = p->i_sigq;
149       
150        //sprintf(buf, "i_shad = %g, i_qbar = %g, i_sigq = %g", i_shad, i_qbar, i_sigq);
151        //XOPNotice(buf);
152       
153        pp.waveHandle = p->waveHandle;
154        pp.x = p->x;
155        //need interpolation
156       
157        if (i_sigq >= 0){
158               
159                nord = 20;
160                va = -3*i_sigq + i_qbar;
161                if (va < 0){
162                        va = 0;
163                }
164                vb = 3*i_sigq + i_qbar;
165               
166                summ = 0.0;
167                ii=0;
168                while (ii < nord) {
169                        zi = (Gauss76Z[ii]*(vb-va) + va + vb)/2.0;
170                        pp.x = zi;
171                        Resoln = i_shad/sqrt(2*Pi*i_sigq*i_sigq);
172                        Resoln *= exp((-1*pow((zi-i_qbar),2))/(2*i_sigq*i_sigq));
173                        if (err = CallFunction(&fi, (void*)&pp, &ans)){
174                                XOPNotice("Error at CallFunction in Gaussian Smear\r");
175                                return err;
176                        }
177                        //                      CallFunction(&fi, (void*)&pp, &ans);
178                        yyy = Gauss76Wt[ii]*Resoln*ans;
179                        summ += yyy;
180                       
181                        ii++;
182                } 
183                answer = (vb-va)/2.0*summ;
184        } else {
185                //USANS so do trapezoidal integration
186                int maxiter = 20;
187                double tol = 1e-4;
188               
189                va=0;
190                vb=fabs(i_sigq);
191               
192                //sprintf(buf, "i_sigq = %g, va = %g, vb = %g\r",i_sigq, va, vb);
193                //XOPNotice("USANS Data. Great\r");
194                //XOPNotice(buf);
195               
196                answer = qtrap_USANS(fi,pp,va,vb,tol,maxiter);
197                answer /= vb;
198        }
199       
200       
201        p->result = answer;
202       
203        return 0;
204        }
205
206
207double
208qtrap_USANS(FunctionInfo fi, PassParams p, double aa, double bb, double eps, int maxit){
209        double olds,ss;
210        double xx,tnm,summ,del,arg1,arg2,ans;
211        int it,jj,kk;
212        int err;
213       
214        //char buf[256];
215       
216        double qval = p.x;
217        double temp = 0;
218       
219        //XOPNotice("In qtrap_USANS\r");
220       
221        olds = -1e-30;
222        for (jj = 1; jj <= maxit; jj++){
223               
224                if (jj == 1){
225                        arg1 = pow(qval,2)+pow(aa,2);
226                        arg2 = pow(qval,2)+pow(bb,2);
227                        p.x = sqrt(arg1);
228                        if (err = CallFunction(&fi, (void*)&p, &ans)){
229                                XOPNotice("Error at 1st CallFunction in qtrap\r");
230                        }
231                        temp = ans;
232                        p.x = sqrt(arg2);
233                        CallFunction(&fi, (void*)&p, &ans);
234                        temp += ans;
235                        ss = 0.5*(bb-aa)*temp;
236                        //sprintf(buf, "jj = 1: aa = %g, bb = %g, temp = %g, ss = %g\r", aa, bb, temp, ss);
237                        //XOPNotice(buf);
238                } else {
239                        it = pow(2,jj-2);
240                        tnm = (double)it;
241                        del = (bb-aa)/tnm;
242                        xx = aa+(0.5*del);
243                        summ = 0;
244                        for (kk = 1; kk <= it; kk++){
245                                arg1 = pow(qval,2)+pow(xx,2);
246                                p.x = sqrt(arg1);
247                                if (err = CallFunction(&fi, (void*)&p, &ans)){
248                                        XOPNotice("Error at 2nd CallFunction in qtrap\r");
249                                }
250                                summ += ans;
251                                xx += del;
252                        }
253                        ss = 0.5*(ss+(bb-aa)*summ/tnm);
254                        //sprintf(buf, "jj = %d: arg1 = %g, summ = %g, tnm= %g, ss = %g\r", jj, arg1, summ, tnm, ss);
255                        //XOPNotice(buf);
256                }
257               
258                //sprintf(buf,"ss = %g\r",ss);
259                //XOPNotice(buf);
260               
261                //              ss = trapzd_USANS(fi,p,aa,bb,jj);
262                if ( fabs(ss-olds) < eps*fabs(olds))
263                        return ss;
264                if ( (ss == 0.0) && (olds == 0.0) && (jj > 6) ) // no progress?
265                        return ss;
266                olds = ss;
267        }
268       
269        XOPNotice("Maxit exceeded in qtrap. If you're here, there was a problem in qtrap");
270        return (ss); // should never get here if function is well behaved
271       
272}
273
274//static double
275//trapzd_USANS(FunctionInfo fi, PassParams p, double aa, double bb, int nn){
276//      double xx,tnm,summ,del,arg1,arg2,ans;
277//      int it,jj;
278//     
279//      static double sval;
280//     
281//      double qval = p.x;
282//      double temp = 0;
283//
284//      //XOPNotice("In trapzd_USANS\r");
285//
286//      if (nn == 1){
287//              arg1 = pow(qval,2)+pow(aa,2);
288//              arg2 = pow(qval,2)+pow(bb,2);
289//              p.x = arg1;
290//              CallFunction(&fi, (void*)&p, &ans);
291//              p.x = arg2;
292//              temp = ans;
293//              CallFunction(&fi, (void*)&p, &ans);
294//              temp += ans;
295//              sval = 0.5*(bb-aa);
296//              return sval;
297//      } else {
298//              it = pow(2,nn-2);
299//              tnm = it;
300//              del = (bb-aa)/tnm;
301//              xx = aa+0.5*del;
302//              summ = 0;
303//              for (jj = 1; jj <= it; jj++){
304//                      arg1 = pow(qval,2)+pow(xx,2);
305//                      p.x = arg1;
306//                      CallFunction(&fi, (void*)&p, &ans);
307//                      summ += ans;
308//                      xx += del;
309//              }
310//              sval = 0.5*(sval+(bb-aa)*summ/tnm);
311//              return sval;
312//      }
313//     
314//
315//}
316
317int
318SmearedCyl_PolyRadiusX(SmearParamsPtr p)
319{
320        double* dp;
321        double q;
322        double ans;
323
324        int ii;
325        double nord,va,vb,summ,yyy,zi,Pi;
326        double answer,Resoln,i_shad,i_qbar,i_sigq;
327//      char buf[256];
328       
329        Pi = 4.0*atan(1.0);
330        i_shad = p->i_shad;
331        i_qbar = p->i_qbar;
332        i_sigq = p->i_sigq;
333       
334        dp = WaveData(p->waveHandle);
335        q = p->x;
336       
337        p->result = Cyl_PolyRadius(dp,q);
338        return 0;
339       
340        if (i_sigq >= 0){
341               
342                nord = 20;
343                va = -3*i_sigq + i_qbar;
344                if (va < 0){
345                        va = 0;
346                }
347                vb = 3*i_sigq + i_qbar;
348                summ = 0.0;
349                ii=0;
350                while (ii < nord) {
351                        zi = (Gauss20Z[ii]*(vb-va) + va + vb)/2.0;
352                        q = zi;
353                        Resoln = i_shad/sqrt(2*Pi*i_sigq*i_sigq);
354                        Resoln *= exp((-1*pow((zi-i_qbar),2))/(2*i_sigq*i_sigq));
355                        ans = Cyl_PolyRadius(dp,q);
356                        yyy = Gauss20Wt[ii]*Resoln*ans;
357                        summ += yyy;
358                       
359                        ii++;
360                } 
361                answer = (vb-va)/2.0*summ;
362        } else {
363                //USANS so do trapezoidal integration
364                //int maxiter = 20;
365                //double tol = 1e-4;
366               
367                //va=0;
368                //vb=fabs(i_sigq);
369               
370                //sprintf(buf, "i_sigq = %g, va = %g, vb = %g\r",i_sigq, va, vb);
371                //XOPNotice("USANS Data. Great\r");
372                //XOPNotice(buf);
373               
374                //answer = qtrap_USANS(fi,pp,va,vb,tol,maxiter);
375                //answer /= vb;
376               
377                //Just return a value that is clearly nonsense, but is non-zero
378                //USANS Smearing code not in form suitable for this test of direct calling
379                //of C library function. AJJ May 9 2007
380                answer = -10;
381        }
382        p->result = answer;
383       
384        return 0;
385}
386
Note: See TracBrowser for help on using the repository browser.