Repository navigation
Expand file tree
/
Copy pathwrite_data.cpp
More file actions
159 lines (135 loc) · 6.27 KB
/
Copy pathwrite_data.cpp
File metadata and controls
159 lines (135 loc) · 6.27 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "var_struct.h"
#include "math.h"
#include "cmpfit-1.2/mpfit.h"
#include "misc.h"
#include "reddening.h"
#include "fitujeme.h"
extern double global_covar[1000][1000];
void write_covar_snpar(vars_struct data, double *covar, long n_coef) {
FILE *cov = fopen("covar_snpar.dat", "w");
int star_start, gal_start;
int n_index = 11;
int my_index[11] = { 0, 1, 2, 3, 4, 5,6,7,8,9, 12};
for (int i=0;i<data.n_stars;i++) {
star_start = data.n_per_filter*data.n_phot_bands+data.n_velo_lines+data.n_gal+data.n_per_star*i;
gal_start = data.n_per_filter*data.n_phot_bands+data.n_velo_lines;
if (data.fit_star[i] == 0) {continue;}
for (int j=0;j<n_index;j++) {
for (int k=0;k<n_index;k++) fprintf(cov, "%15.10le ", global_covar[star_start+my_index[j]][star_start+my_index[k]]);
fprintf(cov,"%15.10le ", global_covar[star_start+my_index[j]][gal_start+data.name_gal[i]]);
fprintf(cov, "\n");
}
for (int j=0;j<n_index;j++) fprintf(cov, "%15.10le ", global_covar[gal_start+data.name_gal[i]][star_start+my_index[j]]);
fprintf(cov, "%15.10le\n", global_covar[gal_start+data.name_gal[i]][gal_start+data.name_gal[i]]);
}
fclose(cov);
}
void write_sn_pars_err_full(double *a, vars_struct data, double *uncert, char *filename) {
FILE *sn_pars = fopen(filename, "w");
int star_start, nphot, nvel;
double chi2;
for (int i=0;i<data.n_stars;i++) {
if (data.fit_star[i] == 0) {continue;}
chi2 = calculate_single_chi2(i, a, data, nphot, nvel);
fprintf(sn_pars, "%-2i %-2i", i, data.name_gal[i]);
star_start = data.n_per_filter*data.n_phot_bands+data.n_velo_lines+data.n_gal+i*data.n_per_star;
for (int j=0;j<data.n_per_star;j++) fprintf(sn_pars, " %10f %10f", a[star_start+j], uncert[star_start+j]);
fprintf(sn_pars, " %5i %10f\n", nphot+nvel, chi2);
}
fclose(sn_pars);
}
void write_sn_pars_err(double *a, vars_struct data, double *uncert, double *uncert_rescaled) {
FILE *sn_pars = fopen("sn_pars_err.dat", "w");
int star_start, gal_start;
for (int i=0;i<data.n_stars;i++) {
if (data.fit_star[i] == 0) {continue;}
fprintf(sn_pars, "%-2i %-2i", i, data.name_gal[i]);
star_start = data.n_per_filter*data.n_phot_bands+data.n_velo_lines+data.n_gal+i*data.n_per_star;
gal_start = data.n_per_filter*data.n_phot_bands+data.n_velo_lines;
// t0
fprintf(sn_pars, " %8.2f %5.2f", a[star_start+0]+data.min_time[i], uncert[star_start+0]);
// plateau duration
fprintf(sn_pars, " %6.2f %5.2f", a[star_start+1], uncert[star_start+1]);
// plateau width
fprintf(sn_pars, " %7.2f %6.2f", a[star_start+2], uncert[star_start+2]);
// distance modulus
if (data.fit_galaxy) {
fprintf(sn_pars, " %6.3f %5.3f %5.3f", a[gal_start+data.name_gal[i]], uncert[gal_start+data.name_gal[i]], uncert_rescaled[gal_start+data.name_gal[i]]);
} else {
fprintf(sn_pars, " %6.3f %5.3f %5.3f", a[star_start+11], uncert[star_start+11], uncert_rescaled[star_start+11]);
}
// reddening
fprintf(sn_pars, " %6.3f %5.3f", a[star_start+12], uncert[star_start+12]);
fprintf(sn_pars, "\n");
}
fclose(sn_pars);
}
void write_sn_pars(double *a, vars_struct data) {
FILE *sn_pars = fopen("sn_pars.dat", "w");
fprintf(sn_pars, "%i\n", data.n_stars);
for (int i=0;i<data.n_stars;i++) {
fprintf(sn_pars, "%i", i);
fprintf(sn_pars, " %15.10le", a[data.n_per_filter*data.n_phot_bands+data.n_velo_lines+data.n_gal+i*data.n_per_star+0]+data.min_time[i]);
for (int j=1;j<data.n_per_star;j++) fprintf(sn_pars, " %15.10le", a[data.n_per_filter*data.n_phot_bands+data.n_velo_lines+data.n_gal+i*data.n_per_star+j]);
fprintf(sn_pars, "\n");
}
fclose(sn_pars);
}
void write_data(double *a, vars_struct data) {
FILE *phot = fopen("phot.dat", "w");
FILE *velo = fopen("velo.dat", "w");
for (int i=0;i<data.n_data; i++) {
int star_start = data.n_per_filter*data.n_phot_bands+data.n_velo_lines+data.n_gal+data.n_per_star*data.name[i];
if (data.fit_star[data.name[i]] == 0) {continue;}
if ((data.hjd[i]-a[star_start+0]) > 500) {continue;}
double theo_mag, theo_vel, wt_exp, rho, temp;
calculate_mag_vel(data, a, data.name[i], data.hjd[i], data.flt[i], data.flt[i], data.dset[i], theo_mag, theo_vel, wt_exp, rho, temp);
theo_mag += add_extinction(data, a, data.name[i], data.flt[i])+add_distance(data, a, data.name[i]);
if (data.dset[i] == 0) {
// skip if velocity obtained after transition
//if (data.hjd[i] > (a[star_start+1]+a[star_start+0])) {continue;}
fprintf(velo, "%f %f %f %i %i %f %f %f\n", data.hjd[i]+data.min_time[data.name[i]], data.val[i], data.err[i], data.flt[i], data.name[i], (data.val[i]-theo_vel)/(wt_exp*data.err[i]), rho, wt_exp);
}
if (data.dset[i] == 1) {
fprintf(phot, "%f %f %f %i %i %f %f %f %f\n", data.hjd[i]+data.min_time[data.name[i]], data.val[i], data.err[i], data.flt[i], data.name[i], (data.val[i]-theo_mag)/data.err[i], rho, temp, wt_exp);
}
}
fclose(velo);
fclose(phot);
}
void write_model(double *a, vars_struct data) {
FILE *model_phot = fopen("model_phot.dat", "w");
FILE *model_velo = fopen("model_velo.dat", "w");
int n_points = 500, star_start;
double model_time;
int flt_list[20] = {0,1,2,3,4,5,6,7,9,10,11,12,13,14,15,16,17,18,19,20};
int n_flt_list = 20;
fprintf(model_phot, "# Columns: JD");
for (int i=0;i<n_flt_list;i++) fprintf(model_phot, ", band_%-2i", flt_list[i]);
fprintf(model_phot, "\n");
for (int i=0;i<data.n_stars;i++) {
star_start = data.n_per_filter*data.n_phot_bands+data.n_velo_lines+data.n_gal+data.n_per_star*i;
for (int j=0;j<n_points;j++) {
double theo_mag, theo_vel, wt_exp, rho, temp;
model_time = double(j)/(double(n_points))*550;//(data.max_time[i]-data.min_time[i]) + data.min_time[i];
fprintf(model_phot, "%15f ", model_time+data.min_time[i]);
fprintf(model_velo, "%15f ", model_time+data.min_time[i]);
//for (int k=0;k<data.n_velo_lines;k++) {
calculate_mag_vel(data, a, i, model_time, 0, 0, 0, theo_mag, theo_vel, wt_exp, rho, temp);
fprintf(model_velo, "%15f ", theo_vel);
//}
for (int k=0;k<n_flt_list;k++) {
calculate_mag_vel(data, a, i, model_time, flt_list[k], 0, 1, theo_mag, theo_vel, wt_exp, rho, temp);
theo_mag += add_extinction(data, a, i, flt_list[k])+add_distance(data, a, i);
fprintf(model_phot, "%15f ", theo_mag);
}
fprintf(model_phot, "\n");
fprintf(model_velo, "\n");
}
}
fclose(model_phot);
fclose(model_velo);
}