1 | #include <math.h> |
---|
2 | #include "invertor.h" |
---|
3 | |
---|
4 | double pi = 3.1416; |
---|
5 | |
---|
6 | /** |
---|
7 | * Deallocate memory |
---|
8 | */ |
---|
9 | void invertor_dealloc(Invertor_params *pars) { |
---|
10 | //free(pars->x); |
---|
11 | //free(pars->y); |
---|
12 | //free(pars->err); |
---|
13 | } |
---|
14 | |
---|
15 | void invertor_init(Invertor_params *pars) { |
---|
16 | pars->d_max = 180; |
---|
17 | } |
---|
18 | |
---|
19 | |
---|
20 | /** |
---|
21 | * P(r) of a sphere, for test purposes |
---|
22 | * |
---|
23 | * @param R: radius of the sphere |
---|
24 | * @param r: distance, in the same units as the radius |
---|
25 | * @return: P(r) |
---|
26 | */ |
---|
27 | double pr_sphere(double R, double r) { |
---|
28 | if (r <= 2.0*R) { |
---|
29 | 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 ); |
---|
30 | } else { |
---|
31 | return 0.0; |
---|
32 | } |
---|
33 | } |
---|
34 | |
---|
35 | /** |
---|
36 | * Orthogonal functions: |
---|
37 | * B(r) = 2r sin(pi*nr/d) |
---|
38 | * |
---|
39 | */ |
---|
40 | double ortho(double d_max, int n, double r) { |
---|
41 | return 2.0*r*sin(pi*n*r/d_max); |
---|
42 | |
---|
43 | } |
---|
44 | |
---|
45 | /** |
---|
46 | * Fourier transform of the nth orthogonal function |
---|
47 | * |
---|
48 | */ |
---|
49 | double ortho_transformed(double d_max, int n, double q) { |
---|
50 | return 8.0*pow(pi, 2.0)/q * d_max * n * pow(-1.0, n+1) |
---|
51 | *sin(q*d_max) / ( pow(pi*n, 2.0) - pow(q*d_max, 2.0) ); |
---|
52 | } |
---|
53 | |
---|
54 | /** |
---|
55 | * First derivative in of the orthogonal function dB(r)/dr |
---|
56 | * |
---|
57 | */ |
---|
58 | double ortho_derived(double d_max, int n, double r) { |
---|
59 | return 2.0*sin(pi*n*r/d_max) + 2.0*r*cos(pi*n*r/d_max); |
---|
60 | } |
---|
61 | |
---|
62 | /** |
---|
63 | * Scattering intensity calculated from the expansion. |
---|
64 | */ |
---|
65 | double iq(double *pars, double d_max, int n_c, double q) { |
---|
66 | double sum = 0.0; |
---|
67 | int i; |
---|
68 | for (i=0; i<n_c; i++) { |
---|
69 | sum += pars[i] * ortho_transformed(d_max, i+1, q); |
---|
70 | } |
---|
71 | return sum; |
---|
72 | } |
---|
73 | |
---|
74 | /** |
---|
75 | * P(r) calculated from the expansion. |
---|
76 | */ |
---|
77 | double pr(double *pars, double d_max, int n_c, double r) { |
---|
78 | double sum = 0.0; |
---|
79 | int i; |
---|
80 | for (i=0; i<n_c; i++) { |
---|
81 | sum += pars[i] * ortho(d_max, i+1, r); |
---|
82 | } |
---|
83 | return sum; |
---|
84 | } |
---|
85 | |
---|
86 | |
---|