source: sasview/pr_inversion/c_extensions/invertor.c @ eb06cbe

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.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since eb06cbe was f71287f4, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

Allow user to set q_min/q_max

  • Property mode set to 100644
File size: 4.0 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}
23
24
25/**
26 * P(r) of a sphere, for test purposes
27 *
28 * @param R: radius of the sphere
29 * @param r: distance, in the same units as the radius
30 * @return: P(r)
31 */
32double pr_sphere(double R, double r) {
33    if (r <= 2.0*R) {
34        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 );
35    } else {
36        return 0.0;
37    }
38}
39
40/**
41 * Orthogonal functions:
42 * B(r) = 2r sin(pi*nr/d)
43 *
44 */
45double ortho(double d_max, int n, double r) {
46        return 2.0*r*sin(pi*n*r/d_max);
47}
48
49/**
50 * Fourier transform of the nth orthogonal function
51 *
52 */
53double ortho_transformed(double d_max, int n, double q) {
54        return 8.0*pow(pi, 2.0)/q * d_max * n * pow(-1.0, n+1)
55                *sin(q*d_max) / ( pow(pi*n, 2.0) - pow(q*d_max, 2.0) );
56}
57
58/**
59 * First derivative in of the orthogonal function dB(r)/dr
60 *
61 */
62double ortho_derived(double d_max, int n, double r) {
63    return 2.0*sin(pi*n*r/d_max) + 2.0*r*cos(pi*n*r/d_max);
64}
65
66/**
67 * Scattering intensity calculated from the expansion.
68 */
69double iq(double *pars, double d_max, int n_c, double q) {
70    double sum = 0.0;
71        int i;
72    for (i=0; i<n_c; i++) {
73        sum += pars[i] * ortho_transformed(d_max, i+1, q);
74    }
75    return sum;
76}
77
78/**
79 * P(r) calculated from the expansion.
80 */
81double pr(double *pars, double d_max, int n_c, double r) {
82    double sum = 0.0; 
83        int i;
84    for (i=0; i<n_c; i++) {
85        sum += pars[i] * ortho(d_max, i+1, r);
86    }
87    return sum;
88}
89
90/**
91 * P(r) calculated from the expansion, with errors
92 */
93void pr_err(double *pars, double *err, double d_max, int n_c, 
94                double r, double *pr_value, double *pr_value_err) {
95    double sum = 0.0; 
96    double sum_err = 0.0;
97    double func_value;
98        int i;
99    for (i=0; i<n_c; i++) {
100        func_value = ortho(d_max, i+1, r);
101        sum += pars[i] * func_value;
102        sum_err += err[i]*err[i]*func_value*func_value;
103    }
104    *pr_value = sum;
105    if (sum_err>0) {
106        *pr_value_err = sqrt(sum_err);
107    } else {
108        *pr_value_err = sum;
109    }
110} 
111
112/**
113 * dP(r)/dr calculated from the expansion.
114 */
115double dprdr(double *pars, double d_max, int n_c, double r) {
116    double sum = 0.0; 
117        int i;
118    for (i=0; i<n_c; i++) {
119        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));
120    }
121    return sum;
122}
123
124/**
125 * regularization term calculated from the expansion.
126 */
127double reg_term(double *pars, double d_max, int n_c, int nslice) {
128    double sum = 0.0; 
129    double r;
130    double deriv;
131        int i;
132    for (i=0; i<nslice; i++) {
133        r = d_max/(1.0*nslice)*i;
134        deriv = dprdr(pars, d_max, n_c, r);
135        sum += deriv*deriv;
136    }
137    return sum/(1.0*nslice)*d_max;
138}
139
140/**
141 * regularization term calculated from the expansion.
142 */
143double int_p2(double *pars, double d_max, int n_c, int nslice) {
144    double sum = 0.0; 
145    double r; 
146    double value;
147        int i;
148    for (i=0; i<nslice; i++) {
149        r = d_max/(1.0*nslice)*i;
150        value = pr(pars, d_max, n_c, r);
151        sum += value*value;
152    }
153    return sum/(1.0*nslice)*d_max;
154}
155
156/**
157 * Get the number of P(r) peaks.
158 */
159int npeaks(double *pars, double d_max, int n_c, int nslice) {
160    double r; 
161    double value;
162        int i;
163        double previous = 0.0;
164        double slope    = 0.0;
165        int count = 0;
166    for (i=0; i<nslice; i++) {
167        r = d_max/(1.0*nslice)*i;
168        value = pr(pars, d_max, n_c, r);
169        if (previous<=value){
170                //if (slope<0) count += 1;
171                slope = 1;
172        } else {
173                //printf("slope -1");
174                if (slope>0) count += 1;
175                slope = -1;
176        }
177        previous = value;
178    }
179    return count;
180}
181
Note: See TracBrowser for help on using the repository browser.