[959eb01] | 1 | #include <math.h> |
---|
| 2 | #include "invertor.h" |
---|
| 3 | #include <memory.h> |
---|
| 4 | #include <stdio.h> |
---|
| 5 | #include <stdlib.h> |
---|
| 6 | |
---|
| 7 | double pi = 3.1416; |
---|
| 8 | |
---|
| 9 | /** |
---|
| 10 | * Deallocate memory |
---|
| 11 | */ |
---|
| 12 | void invertor_dealloc(Invertor_params *pars) { |
---|
| 13 | free(pars->x); |
---|
| 14 | free(pars->y); |
---|
| 15 | free(pars->err); |
---|
| 16 | } |
---|
| 17 | |
---|
| 18 | void invertor_init(Invertor_params *pars) { |
---|
| 19 | pars->d_max = 180; |
---|
| 20 | pars->q_min = -1.0; |
---|
| 21 | pars->q_max = -1.0; |
---|
[cb62bd5] | 22 | pars->est_bck = 0; |
---|
[959eb01] | 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 | */ |
---|
| 33 | double 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 | */ |
---|
| 46 | double 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 | */ |
---|
| 54 | double 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 | */ |
---|
| 63 | double 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 | */ |
---|
| 102 | double 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 | */ |
---|
| 109 | double 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 | */ |
---|
| 122 | double 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 | */ |
---|
| 134 | double 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 | */ |
---|
| 146 | void 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 | */ |
---|
| 169 | double 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 | */ |
---|
| 181 | double 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 | */ |
---|
| 197 | double 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 | */ |
---|
| 213 | double 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 | */ |
---|
| 229 | int 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 | */ |
---|
| 257 | double 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 | */ |
---|
| 277 | double 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 | */ |
---|
| 301 | double 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 | } |
---|