127 lines
3.2 KiB
C
127 lines
3.2 KiB
C
#include "common.h"
|
||
|
||
/* `vector_fprint(fp, x)` prints to file `fp`
|
||
* the components of a gsl_vector.
|
||
*/
|
||
void vector_fprint(FILE *fp, const gsl_vector *x) {
|
||
fprintf(fp, "[");
|
||
for (size_t i = 0; i < x->size; i++) {
|
||
fprintf(fp, i < x->size -1 ? "%10.4g, " : "%10.4g",
|
||
gsl_vector_get(x, i));
|
||
}
|
||
fprintf(fp, "]\n");
|
||
}
|
||
|
||
|
||
/* `matrix_fprint(fp, x)` prints to file `fp`
|
||
* the components of a gsl_vector.
|
||
*/
|
||
void matrix_fprint(FILE *fp, const gsl_matrix *x) {
|
||
gsl_vector row;
|
||
for (size_t i = 0; i < x->size1; i++) {
|
||
row = gsl_matrix_const_row(x, i).vector;
|
||
fputs(" ", fp);
|
||
vector_fprint(fp, &row);
|
||
}
|
||
}
|
||
|
||
|
||
/* `vector_map(f, v)` applies the function
|
||
* `f` to each element of the vector `v`, inplace.
|
||
* The function return a GSL exit code to either
|
||
* report success or a failure.
|
||
*/
|
||
int vector_map(double (*f)(double x), gsl_vector *v) {
|
||
double y;
|
||
for (size_t i = 0; i < v->size; i++ ) {
|
||
y = (*f)(gsl_vector_get(v, i));
|
||
gsl_vector_set(v, i, y);
|
||
}
|
||
return GSL_SUCCESS;
|
||
}
|
||
|
||
|
||
/* Angular distribution function F(θ,φ; α,β,γ)
|
||
* Takes three parameters α,β,γ and the angles θ,φ.
|
||
*/
|
||
double distr(const gsl_vector *par,
|
||
const struct event *event){
|
||
double a = gsl_vector_get(par, 0),
|
||
b = gsl_vector_get(par, 1),
|
||
c = gsl_vector_get(par, 2);
|
||
|
||
double th = event->th,
|
||
ph = event->ph;
|
||
|
||
double A = 3 /(4 * M_PI);
|
||
double B = 0.5 * (1 - a);
|
||
double C = 0.5 * (3 * a - 1) * pow(cos(th), 2);
|
||
double D = b * pow(sin(th), 2) * cos(2 * ph);
|
||
double E = sqrt(2) * c * sin(2 * th) * cos(ph);
|
||
|
||
return A * (B + C - D - E);
|
||
}
|
||
|
||
/* `grad_distr(par, event, grad)`
|
||
* computes the gradient of F(α,β,γ; α,β) as:
|
||
*
|
||
* ⎡ 3/2⋅cos²(θ) - 1/2 ⎤
|
||
* ⎢ ⎥
|
||
* ⎢ sin²(θ)⋅cos(2φ) ⎥ 3/(4π)
|
||
* ⎢ ⎥
|
||
* ⎣ √2⋅sin(2θ)⋅cos(φ) ⎦
|
||
*
|
||
* where `par` is [α,β,γ], `event` is the angle
|
||
* structure. The results is stored in the
|
||
* gsl_vector `grad`
|
||
*/
|
||
void grad_distr(const gsl_vector *par,
|
||
const struct event *event,
|
||
gsl_vector *grad) {
|
||
double th = event->th,
|
||
ph = event->ph;
|
||
|
||
gsl_vector_set(grad, 0, 1.5*pow(cos(th), 2) - 0.5);
|
||
gsl_vector_set(grad, 1, -pow(sin(th), 2) * cos(2*ph));
|
||
gsl_vector_set(grad, 2, -sqrt(2)*sin(2*th) * cos(ph));
|
||
gsl_vector_scale(grad, 3.0/(4.0 * M_PI));
|
||
}
|
||
|
||
|
||
/* Test the compatibility between the original
|
||
* parameters `def` and the estimates `est.par`
|
||
* with standard errors `est.err`.
|
||
*/
|
||
void compatibility(gsl_vector *def,
|
||
min_result est) {
|
||
|
||
char* names[] = { "α", "β", "γ" };
|
||
double a, b;
|
||
double t, sigma, alpha;
|
||
|
||
fprintf(stderr, "| par | σ | α-value |\n");
|
||
fprintf(stderr, "|-----|----------|----------|\n");
|
||
for (size_t i = 0; i < def->size; i++) {
|
||
// discrepancy
|
||
a = gsl_vector_get(def, i);
|
||
b = gsl_vector_get(est.par, i);
|
||
sigma = gsl_vector_get(est.err, i);
|
||
|
||
// normal t-statistic
|
||
t = fabs(a - b)/sigma;
|
||
// α-value / significance
|
||
alpha = 1 - erf(t/sqrt(2));
|
||
|
||
fprintf(stderr, "| %s | %-8.2g | %-8.2g |\n",
|
||
names[i], sigma, alpha);
|
||
}
|
||
}
|
||
|
||
|
||
/* Free a min_result structure.
|
||
*/
|
||
void min_result_free(min_result x) {
|
||
gsl_vector_free(x.par);
|
||
gsl_matrix_free(x.cov);
|
||
}
|