source: sans/XOP_Dev/SANSAnalysis/XOP/ResolutionSmearing.c @ 834

Last change on this file since 834 was 834, checked in by srkline, 11 years ago

Changes to the XOP code to upgrade to ToolKit? v6. Changes are the ones outlined in the Appendix A of the TK6 manual. SOme of the XOP support routines and the #pragma for 2-byte structures have changed. Per Howard Rodstein, there is no need to change the c files to cpp. c should work and compile just fine.

These changes work correctly on my mac. Next is to make sure that they work correctly on the two build machines.

  • 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 = (double*)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.