source: sasview/src/sas/sascalc/pr/c_extensions/invertor.c @ cb62bd5

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since cb62bd5 was cb62bd5, checked in by lewis, 7 years ago

Use manually inputted background level in Pr calculation

  • Property mode set to 100644
File size: 6.8 KB
Line 
1#include <math.h>
2#include "invertor.h"
3#include <memory.h>
4#include <stdio.h>
5#include <stdlib.h>
6
7double pi = 3.1416;
8
9/**
10 * Deallocate memory
11 */
12void invertor_dealloc(Invertor_params *pars) {
13        free(pars->x);
14        free(pars->y);
15        free(pars->err);
16}
17
18void invertor_init(Invertor_params *pars) {
19        pars->d_max = 180;
20        pars->q_min = -1.0;
21        pars->q_max = -1.0;
22        pars->est_bck = 0;
23}
24
25
26/**
27 * P(r) of a sphere, for test purposes
28 *
29 * @param R: radius of the sphere
30 * @param r: distance, in the same units as the radius
31 * @return: P(r)
32 */
33double pr_sphere(double R, double r) {
34    if (r <= 2.0*R) {
35        return 12.0* pow(0.5*r/R, 2.0) * pow(1.0-0.5*r/R, 2.0) * ( 2.0 + 0.5*r/R );
36    } else {
37        return 0.0;
38    }
39}
40
41/**
42 * Orthogonal functions:
43 * B(r) = 2r sin(pi*nr/d)
44 *
45 */
46double ortho(double d_max, int n, double r) {
47        return 2.0*r*sin(pi*n*r/d_max);
48}
49
50/**
51 * Fourier transform of the nth orthogonal function
52 *
53 */
54double ortho_transformed(double d_max, int n, double q) {
55        return 8.0*pow(pi, 2.0)/q * d_max * n * pow(-1.0, n+1)
56                *sin(q*d_max) / ( pow(pi*n, 2.0) - pow(q*d_max, 2.0) );
57}
58
59/**
60 * Slit-smeared Fourier transform of the nth orthogonal function.
61 * Smearing follows Lake, Acta Cryst. (1967) 23, 191.
62 */
63double ortho_transformed_smeared(double d_max, int n, double height, double width, double q, int npts) {
64        double sum, y, z;
65        int i, j, n_height, n_width;
66        double count_w;
67        double fnpts;
68        sum = 0.0;
69        fnpts = (float)npts-1.0;
70
71        // Check for zero slit size
72        n_height = (height>0) ? npts : 1;
73        n_width  = (width>0)  ? npts : 1;
74
75        count_w = 0.0;
76
77        for(j=0; j<n_height; j++) {
78                if(height>0){
79                        z = height/fnpts*(float)j;
80                } else {
81                        z = 0.0;
82                }
83
84                for(i=0; i<n_width; i++) {
85                        if(width>0){
86                                y = -width/2.0+width/fnpts*(float)i;
87                        } else {
88                                y = 0.0;
89                        }
90                        if (((q-y)*(q-y)+z*z)<=0.0) continue;
91                        count_w += 1.0;
92                        sum += ortho_transformed(d_max, n, sqrt((q-y)*(q-y)+z*z));
93                }
94        }
95        return sum/count_w;
96}
97
98/**
99 * First derivative in of the orthogonal function dB(r)/dr
100 *
101 */
102double ortho_derived(double d_max, int n, double r) {
103    return 2.0*sin(pi*n*r/d_max) + 2.0*r*cos(pi*n*r/d_max);
104}
105
106/**
107 * Scattering intensity calculated from the expansion.
108 */
109double iq(double *pars, double d_max, int n_c, double q) {
110    double sum = 0.0;
111        int i;
112    for (i=0; i<n_c; i++) {
113        sum += pars[i] * ortho_transformed(d_max, i+1, q);
114    }
115    return sum;
116}
117
118/**
119 * Scattering intensity calculated from the expansion,
120 * slit-smeared.
121 */
122double iq_smeared(double *pars, double d_max, int n_c, double height, double width, double q, int npts) {
123    double sum = 0.0;
124        int i;
125    for (i=0; i<n_c; i++) {
126        sum += pars[i] * ortho_transformed_smeared(d_max, i+1, height, width, q, npts);
127    }
128    return sum;
129}
130
131/**
132 * P(r) calculated from the expansion.
133 */
134double pr(double *pars, double d_max, int n_c, double r) {
135    double sum = 0.0;
136        int i;
137    for (i=0; i<n_c; i++) {
138        sum += pars[i] * ortho(d_max, i+1, r);
139    }
140    return sum;
141}
142
143/**
144 * P(r) calculated from the expansion, with errors
145 */
146void pr_err(double *pars, double *err, double d_max, int n_c,
147                double r, double *pr_value, double *pr_value_err) {
148    double sum = 0.0;
149    double sum_err = 0.0;
150    double func_value;
151        int i;
152    for (i=0; i<n_c; i++) {
153        func_value = ortho(d_max, i+1, r);
154        sum += pars[i] * func_value;
155        //sum_err += err[i]*err[i]*func_value*func_value;
156        sum_err += err[i*n_c+i]*func_value*func_value;
157    }
158    *pr_value = sum;
159    if (sum_err>0) {
160        *pr_value_err = sqrt(sum_err);
161    } else {
162        *pr_value_err = sum;
163    }
164}
165
166/**
167 * dP(r)/dr calculated from the expansion.
168 */
169double dprdr(double *pars, double d_max, int n_c, double r) {
170    double sum = 0.0;
171        int i;
172    for (i=0; i<n_c; i++) {
173        sum += pars[i] * 2.0*(sin(pi*(i+1)*r/d_max) + pi*(i+1)*r/d_max * cos(pi*(i+1)*r/d_max));
174    }
175    return sum;
176}
177
178/**
179 * regularization term calculated from the expansion.
180 */
181double reg_term(double *pars, double d_max, int n_c, int nslice) {
182    double sum = 0.0;
183    double r;
184    double deriv;
185        int i;
186    for (i=0; i<nslice; i++) {
187        r = d_max/(1.0*nslice)*i;
188        deriv = dprdr(pars, d_max, n_c, r);
189        sum += deriv*deriv;
190    }
191    return sum/(1.0*nslice)*d_max;
192}
193
194/**
195 * regularization term calculated from the expansion.
196 */
197double int_p2(double *pars, double d_max, int n_c, int nslice) {
198    double sum = 0.0;
199    double r;
200    double value;
201        int i;
202    for (i=0; i<nslice; i++) {
203        r = d_max/(1.0*nslice)*i;
204        value = pr(pars, d_max, n_c, r);
205        sum += value*value;
206    }
207    return sum/(1.0*nslice)*d_max;
208}
209
210/**
211 * Integral of P(r)
212 */
213double int_pr(double *pars, double d_max, int n_c, int nslice) {
214    double sum = 0.0;
215    double r;
216    double value;
217        int i;
218    for (i=0; i<nslice; i++) {
219        r = d_max/(1.0*nslice)*i;
220        value = pr(pars, d_max, n_c, r);
221        sum += value;
222    }
223    return sum/(1.0*nslice)*d_max;
224}
225
226/**
227 * Get the number of P(r) peaks.
228 */
229int npeaks(double *pars, double d_max, int n_c, int nslice) {
230    double r;
231    double value;
232        int i;
233        double previous = 0.0;
234        double slope    = 0.0;
235        int count = 0;
236    for (i=0; i<nslice; i++) {
237        r = d_max/(1.0*nslice)*i;
238        value = pr(pars, d_max, n_c, r);
239        if (previous<=value){
240                //if (slope<0) count += 1;
241                slope = 1;
242        } else {
243                //printf("slope -1");
244                if (slope>0) count += 1;
245                slope = -1;
246        }
247        previous = value;
248    }
249    return count;
250}
251
252/**
253 * Get the fraction of the integral of P(r) over the whole range
254 * of r that is above zero.
255 * A valid P(r) is define as being positive for all r.
256 */
257double positive_integral(double *pars, double d_max, int n_c, int nslice) {
258    double r;
259    double value;
260        int i;
261        double sum_pos = 0.0;
262        double sum = 0.0;
263
264    for (i=0; i<nslice; i++) {
265        r = d_max/(1.0*nslice)*i;
266        value = pr(pars, d_max, n_c, r);
267        if (value>0.0) sum_pos += value;
268        sum += fabs(value);
269    }
270    return sum_pos/sum;
271}
272
273/**
274 * Get the fraction of the integral of P(r) over the whole range
275 * of r that is at least one sigma above zero.
276 */
277double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice) {
278    double r;
279        int i;
280        double sum_pos = 0.0;
281        double sum = 0.0;
282        double pr_val;
283        double pr_val_err;
284
285    for (i=0; i<nslice; i++) {
286        r = d_max/(1.0*nslice)*i;
287        pr_err(pars, err, d_max, n_c, r, &pr_val, &pr_val_err);
288        if (pr_val>pr_val_err) sum_pos += pr_val;
289        sum += fabs(pr_val);
290
291
292    }
293    return sum_pos/sum;
294}
295
296/**
297 * R_g radius of gyration calculation
298 *
299 * R_g**2 = integral[r**2 * p(r) dr] /  (2.0 * integral[p(r) dr])
300 */
301double rg(double *pars, double d_max, int n_c, int nslice) {
302    double sum_r2 = 0.0;
303    double sum    = 0.0;
304    double r;
305    double value;
306        int i;
307    for (i=0; i<nslice; i++) {
308        r = d_max/(1.0*nslice)*i;
309        value = pr(pars, d_max, n_c, r);
310        sum += value;
311        sum_r2 += r*r*value;
312    }
313    return sqrt(sum_r2/(2.0*sum));
314}
Note: See TracBrowser for help on using the repository browser.