read_rlrdata.c

最大似然估计算法

#include "timeseries.h"

time_series read_rlrdata(char *filename, options op, double scale_factor) {

	FILE *fpin;
	time_series ts;

	int j, n;

	int n_offset;
	int n_data;

	int n_alloc   = 512;
	int off_alloc = 16;


	double tt, up; 

	char line[512], code[64];

	if ((fpin = fopen(filename, "r")) == NULL) {
		fprintf(stderr, " read_cats_series : Cannot open file %s\n", filename);
		exit(EXIT_FAILURE);
	}


	ts.n_series   = 1;
	ts.got_errors = 0;

	ts.t            = (double *) calloc((size_t)n_alloc     , sizeof(double));
	ts.data	        = (double *) calloc((size_t)n_alloc     , sizeof(double));
	ts.offsets      = (double *) calloc((size_t)(off_alloc) , sizeof(double));
	ts.off_code     = (int *)    calloc((size_t)(off_alloc) , sizeof(int));

	n_data   = 0;
	n_offset = 0;

	while (fgets (line, 512, fpin)) {
		if (strncmp(line,"# offs",6) == 0) {
			n = sscanf(line, "# offset %lf %d", &ts.offsets[n_offset], &ts.off_code[n_offset]);
	       		if (n == 1) ts.off_code[n_offset] = 1;
			ts.off_code[n_offset] = ts.off_code[n_offset] > 1 ? 1 : ts.off_code[n_offset];
			ts.off_code[n_offset] = ts.off_code[n_offset] < 0 ? 0 : ts.off_code[n_offset];
			n_offset++;
			if (n_offset == off_alloc) {
				off_alloc += 10;
				ts.offsets  = (double *) realloc( (void *)ts.offsets,  off_alloc * sizeof(double) );
				ts.off_code = (int *)    realloc( (void *)ts.off_code, off_alloc * sizeof(int) );
			}
		} else if (strncmp(line,"#",1) == 0) continue;
		else {
			n = sscanf(line, "%lf %lf %s", &tt, &up, &code);
			if (n == 2) {
				ts.t[n_data]    = tt;
				ts.data[n_data] = up * scale_factor / 1000.0;
				n_data++;
			}
		}
		if (n_data == n_alloc) {
			n_alloc += 512;
			ts.t    = (double *) realloc( (void *)ts.t,     n_alloc * sizeof(double) );
			ts.data = (double *) realloc( (void *)ts.data,  n_alloc * sizeof(double) );
		}
	}

	ts.n_data = n_data;
	ts.n_offsets = n_offset;

	ts.index       = (int *)    calloc((size_t)(n_data)   , sizeof(int));
	ts.c_id        = (int *)    calloc((size_t)(1)        , sizeof(int));

	ts.t        = (double *) realloc( (void *)ts.t,        n_data *   sizeof(double));
	ts.data     = (double *) realloc( (void *)ts.data,     n_data *   sizeof(double));
	ts.offsets  = (double *) realloc( (void *)ts.offsets,  n_offset * sizeof(double));
	ts.off_code = (int *)    realloc( (void *)ts.off_code, n_offset * sizeof(int));

	ts.formal_error = (double *) calloc((size_t)ts.n_data     , sizeof(double));
	for (j = 0; j < ts.n_data; j++) ts.formal_error[j] = 1.0;

	ts.c_id[0] = 4;

	ts.fs = calculate_fs2(ts.t, ts.n_data, ts.index, &ts.count);

	ts.n_full = ts.index[ts.n_data-1]+1;
	ts.time_span = ts.t[n_data-1] - ts.t[0];

	if (op.verbose) fprintf(op.fpout, " Sampling frequency %.6g (Hz), %.6f days\n", ts.fs, 1.0 / ts.fs / 24.0 / 3600.0);
	if (op.verbose) fprintf(op.fpout, " Number of samples 1 period apart = %d of %d\n", ts.count, ts.n_data-1);
	if (op.verbose) fprintf(op.fpout, " Number of points in full series  = %d\n", ts.n_full);

	return(ts);

}