scargle.c

最大似然估计算法

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

void scargle(double *x, double *y, int n, double ofac, double hifac, double *px, double *py, int np, int *nout, int *jmax, double *prob) {

int i, j, kout;

double ave, c, cc, cwtau, effm, expy, pnow, pymax, s, ss, sumc, sumcy, sums, sumsh;
double sumsy, swtau, var, wtau, xave, xdif, xmax, xmin, yy;

double arg, wtemp, *wi, *wpi, *wpr, *wr;

wi  = (double *) calloc((size_t)n, sizeof(double) );
wpi = (double *) calloc((size_t)n, sizeof(double) );
wpr = (double *) calloc((size_t)n, sizeof(double) );
wr  = (double *) calloc((size_t)n, sizeof(double) );

*nout = (int)(ofac) * (int)(hifac) * n / 2;
kout  = *nout;
if (kout > np) {
	fprintf(stderr, "Output arrays too short in scargle %d > %d\n", kout, np);
	exit(EXIT_FAILURE);
}

avevar(y,n,&ave,&var);

xmax = xmin = x[1];
for (j = 0; j <n; j++) {
	if (x[j] > xmax) xmax = x[j];
	if (x[j] < xmin) xmin = x[j];
}
xdif  = xmax - xmin;
xave  = 0.5 * (xmax + xmin);
pymax = 0.0;

pnow = 1.0 / (xdif * ofac);
for (j = 0; j <n; j++) {
	arg    = 2.0 * M_PI * ((x[j]-xave)*pnow);
	wpr[j] = -2.0 * sin(0.5 * arg) * sin(0.5 * arg);
	wpi[j] = sin(arg);
	wr[j]  = cos(arg);
	wi[j]  = wpi[j];
}
for (i=0; i < kout; i++) {
	px[i] = pnow;
	sumsh = sumc = 0.0;
	for (j = 0; j < n; j++) {
		c = wr[j];
		s = wi[j];
		sumsh += s * c;
		sumc  += (c-s)*(c+s);
	}
	wtau = 0.5 * atan2(2.0*sumsh,sumc);
	swtau = sin(wtau);
	cwtau = cos(wtau);
	sums = sumc = sumsy = sumcy = 0.0;
	for (j = 0; j < n; j++) {
		s      = wi[j];
		c      = wr[j];
		ss     = s*cwtau-c*swtau;
		cc     = c*cwtau+s*swtau;
		sums  += ss*ss;
		sumc  += cc*cc;
		yy     = y[j] - ave;
		sumsy += yy*ss;
		sumcy += yy*cc;
		wr[j]  = ((wtemp=wr[j]) * wpr[j]-wi[j]*wpi[j]) + wr[j];
		wi[j]  = (wi[j]*wpr[j]+wtemp*wpi[j])+wi[j];
	}
	py[i] = 0.5 * (sumcy*sumcy/sumc+sumsy*sumsy/sums) ;
	if (py[i] >= pymax) pymax = py[(*jmax = i)];
	pnow += 1.0/(ofac*xdif);
}
expy  = exp(-pymax);
effm  = 2.0 * (double) kout / ofac;
*prob = effm*expy;
if(*prob > 0.01) *prob = 1.0-pow(1.0-expy,effm);

free(wr);
free(wpr);
free(wpi);
free(wi);


}