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

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

Add combined XOP code. Currently only contains XCode project file to build Universal binary suitable for Igor 6.

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