/* NIGHTFALL Light Curve Synthesis Program */
/* Copyright (C) 1998 Rainer Wichmann */
/* */
/* This program is free software; you can redistribute it */
/* and/or modify */
/* it under the terms of the GNU General Public License as */
/* published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include "Light.h"
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Compute merit function (chi square)
@tip Computed globally over all data
@param (void)
@return (double) The merit
@heading Light Curve
****************************************************************************/
double MeritFunction()
{
int band; /* passband */
long data; /* # of data in passband */
long middle;
long i; /* loop counter */
int Comp; /* stellar component */
double Parabolic, Phi, Datum, Weight; /* parabolic fit */
double Hi = 0.0, Lo = 0.0, Mid = 0.0; /* parabolic fit */
double FHi =0.0, FLo =0.0, FMid =0.0; /* parabolic fit */
double *PhaseShifted, Shift; /* phase shift */
double Add = 0.0, ChiSquare = 0.0; /* merit */
static double OneOverTwoPi = 1.0 / (PI+PI); /* convenience constant */
long dof = 0; /* degrees of freedom */
PhotoOut *FluxPtr = FluxOut; /* pointer to flux */
/* ----------- allocate memory -------------------------------------- */
PhaseShifted = malloc(MAXOBS * sizeof(double));
if (PhaseShifted == NULL)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
free(PhaseShifted);
make_dialog (_(errmsg[0]));
return(0.0);
} else nf_error(_(errmsg[0]));
}
#else
nf_error(_(errmsg[1]));
#endif
/* create an ordered array of phases starting at zero, ranging to one */
/* note that in an elliptic orbit, the starting value can change */
/* thus we have to do this each time anew */
Shift = OneOverTwoPi * FluxPtr->Phase;
for (i = 0; i < PhaseSteps; ++i) {
PhaseShifted[i] = (OneOverTwoPi * FluxPtr->Phase) - Shift;
PhaseShifted[i] = PhaseShifted[i] - floor(PhaseShifted[i]);
if (Flags.debug[CHISQR] == ON)
fprintf(stderr,
_("Mag %9.5f Phase %9.5f Shifted %9.5f \n"),
FluxPtr->Mag[Vmag], OneOverTwoPi * FluxPtr->Phase,
(*PhaseShifted));
++FluxPtr;
}
for (band = 0; band < NUM_MAG+2; ++band) {
if (Flags.Passbands[band] > 0) dof = dof + Flags.Passbands[band];
}
/* >>>>>>>>>>>>>> loop over PASSBANDS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
for (band = 0; band < NUM_MAG; ++band) {
/* we have observational data in this band */
if (Flags.Passbands[band] > 0) {
if (Flags.debug[CHISQR] == ON)
fprintf(stderr,
_("MeritFunction: %ld data found for band %s\n"),
Flags.Passbands[band], Filters[band]);
for (data = 0; data < Flags.Passbands[band]; ++data) {
/* the observed datum */
Datum = Observed[band][data].mag;
/* and its weight */
Weight = Observed[band][data].weight;
/* at the phase */
Phi=(Observed[band][data].phi+0.5)-Shift;
/* get it between zero and one */
Phi = Phi - floor(Phi);
/* the corresponding model curve index */
/* do a hash search */
middle = (long) (floor(PhaseSteps * Phi));
if (Flags.debug[CHISQR] == ON)
fprintf(stderr,
_("No %ld Index %ld X %10.6f Y %10.6f W %10.6f\n"),
data, middle, Phi, Datum, Weight);
/* a parabolic fit thru three points */
/* first get the three points */
if ( middle > 0 && middle < PhaseSteps-1) {
Lo = PhaseShifted[middle-1];
Mid = PhaseShifted[middle];
Hi = PhaseShifted[middle+1];
FLo = FluxOut[middle-1].Mag[band];
FMid = FluxOut[middle].Mag[band];
FHi = FluxOut[middle+1].Mag[band];
}
else if (middle == 0) {
/* wrap around low border */
Lo = PhaseShifted[PhaseSteps-1] - 1.0;
Mid = PhaseShifted[middle];
Hi = PhaseShifted[middle+1];
FLo = FluxOut[PhaseSteps-1].Mag[band];
FMid = FluxOut[middle].Mag[band];
FHi = FluxOut[middle+1].Mag[band];
}
else if (middle == PhaseSteps-1) {
/* wrap around high border */
Lo = PhaseShifted[middle-1];
Mid = PhaseShifted[middle];
Hi = 1.0 + PhaseShifted[0];
FLo = FluxOut[middle-1].Mag[band];
FMid = FluxOut[middle].Mag[band];
FHi = FluxOut[0].Mag[band];
}
/* interpolate by Lagrange's formula */
Parabolic = FLo * ((Phi-Mid)*(Phi-Hi )) / ((Lo -Mid)*(Lo -Hi ))
+ FMid * ((Phi-Lo )*(Phi-Hi )) / ((Mid-Lo )*(Mid-Hi ))
+ FHi * ((Phi-Lo )*(Phi-Mid)) / ((Hi -Lo )*(Hi -Mid));
/* the ChiSquare */
Observed[band][data].residual = (float) (Parabolic-Datum);
Add = SQR(Parabolic-Datum)/(Weight*Weight);
ChiSquare = ChiSquare + Add;
if (Flags.debug[CHISQR] == ON)
fprintf(stderr,
_(" ChiSquare %10.6f Add %10.6f Parabolic %10.6f\n"),
ChiSquare, Add, Parabolic);
} /* END loop over data */
}
} /* END loop over passbands */
/* >>>>>>>>>>>>>>>>>>>> loop over RADIAL VELOCITIES <<<<<<<<<<<<<<<<<<< */
for (band = NUM_MAG; band < (NUM_MAG+2); ++band) {
Comp = band - NUM_MAG;
/* we have observational data in this band */
if (Flags.Passbands[band] > 0) {
if (Flags.debug[CHISQR] == ON)
printf(_("MeritFunction: %ld data found for band %s\n"),
Flags.Passbands[band], Filters[band]);
for (data = 0; data < Flags.Passbands[band]; ++data) {
/* the observed datum */
/* sign of radial velocity ? */
Datum = Observed[band][data].mag;
/* and its weight */
Weight = Observed[band][data].weight;
/* at the phase RVTEST */
if (Flags.elliptic == ON)
Phi = (Observed[band][data].phi + 0.5) - Shift;
else
Phi = (Observed[band][data].phi) - Shift;
/* get it between zero and one */
Phi = Phi - floor(Phi);
/* the corresponding model curve index */
/* do a hash search */
middle = (long) (floor(PhaseSteps * Phi));
if (Flags.debug[CHISQR] == ON)
fprintf(stderr,
_("No %ld Index %ld X %10.6f Y %10.6f W %10.6f\n"),
data, middle, Phi+Shift, Datum, Weight);
/* a parabolic fit thru three points */
/* first get the three points */
if ( middle > 0 && middle < PhaseSteps-1) {
Lo = PhaseShifted[middle-1];
Mid = PhaseShifted[middle];
Hi = PhaseShifted[middle+1];
FLo = FluxOut[middle-1].RV[Comp];
FMid = FluxOut[middle].RV[Comp];
FHi = FluxOut[middle+1].RV[Comp];
}
else if (middle == 0) {
/* wrap around low border */
Lo = PhaseShifted[PhaseSteps-1] - 1.0;
Mid = PhaseShifted[middle];
Hi = PhaseShifted[middle+1];
FLo = FluxOut[PhaseSteps-1].RV[Comp];
FMid = FluxOut[middle].RV[Comp];
FHi = FluxOut[middle+1].RV[Comp];
}
else if (middle == PhaseSteps-1) {
/* wrap around high border */
Lo = PhaseShifted[middle-1];
Mid = PhaseShifted[middle];
Hi = 1.0 + PhaseShifted[0];
FLo = FluxOut[middle-1].RV[Comp];
FMid = FluxOut[middle].RV[Comp];
FHi = FluxOut[0].RV[Comp];
}
FLo = FLo / 1000.;
FMid = FMid / 1000.;
FHi = FHi / 1000.;
/* interpolate by Lagrange's formula */
Parabolic = FLo * ((Phi-Mid)*(Phi-Hi )) / ((Lo -Mid)*(Lo -Hi ))
+ FMid * ((Phi-Lo )*(Phi-Hi )) / ((Mid-Lo )*(Mid-Hi ))
+ FHi * ((Phi-Lo )*(Phi-Mid)) / ((Hi -Lo )*(Hi -Mid));
/* the ChiSquare */
Observed[band][data].residual = (float) (Parabolic - Datum);
Add = SQR(Observed[band][data].residual)/(Weight*Weight);
ChiSquare = ChiSquare + Add;
if (Flags.debug[CHISQR] == ON)
printf(
_(" ChiSquare %10.6f Add %10.6f Parabolic %10.6f\n"),
ChiSquare, Add, Parabolic);
} /* END loop over data */
}
} /* END loop over radial velocities */
if (Flags.debug[VERBOSE] == ON || Flags.debug[CHISQR] == ON)
printf(_("Chi Square = %9.3g\n"), ChiSquare);
free(PhaseShifted);
if (dof != 0) return( (ChiSquare/dof) );
else return( 0.0 );
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Performs the runs test
@param (Photo3*) data The data
@param (long) NumData The number of data
@param (int*) Runs The number of runs
@param (int*) UpLim Upper 5 per cent limit
@param (int*) LowLim Lower 5 per cent limit
@return (void)
@heading Light Curve
****************************************************************************/
void Runs(Photo3 *data, long NumData, int *Runs, int *UpLim, int *LowLim)
{
int CountRun; /* # of runs */
long i; /* loop counter */
int CountBelow; /* # residuals above mean */
int CountAbove; /* # residuals below mean */
double ResMean; /* mean residual */
double Mean; /* expected runs */
double Sigma; /* confidence interval for expected */
ResMean = 0.0;
for (i = 0; i < NumData; ++i) {
ResMean = ResMean + data[i].residual;
}
ResMean = ResMean/NumData;
SortPhotoPoints(NumData, data);
CountRun = 0; CountBelow=0; CountAbove=0;
for (i = 0; i < NumData; ++i) {
if ((double) data[i].residual >= (ResMean+FLT_EPSILON)) ++CountAbove;
else ++CountBelow;
}
for (i = 1; i < NumData; ++i) {
if ((double) data[i-1].residual <= (ResMean-FLT_EPSILON)
&& (double) data[i].residual >= ResMean) ++CountRun;
else if ((double) data[i-1].residual >= (ResMean+FLT_EPSILON)
&& (double) data[i].residual <= ResMean) ++CountRun;
}
Mean = 1.0 + (2.0*CountAbove*CountBelow)/NumData;
if (NumData > 1)
Sigma = sqrt(2.0*CountAbove*CountBelow*(2*CountAbove*CountBelow-NumData)
/(NumData*NumData*(NumData-1.0)));
else
Sigma = 0.0;
*Runs = CountRun;
*LowLim = (int)(Mean - 1.96*Sigma - 0.5);
*UpLim = (int)(Mean + 1.96*Sigma + 0.5);
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Copy data point b to a
@param (Photo3*) a Output
@param (Photo3*) b Input
@return (void)
@heading Light Curve
****************************************************************************/
void V2CopyPhoto(/*@out@*/ Photo3 *a, Photo3 *b)
{
a->time = b->time;
a->mag = b->mag;
a->weight = b->weight;
a->phi = b->phi;
a->residual = b->residual;
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Sort in ascending order in phi w/ Heapsort Algorithm
@long D.E. Knuth, Sorting and Searching, 1973
Contrary to Quicksort, this is guaranteed to be always of
order NlogN, even in worst case.
Steps Hn refer to description by D.E. Knuth.
@tip Sorting is done in place
@param (long) n Number of data in array
@param (Photo3*) seq Array to order
@return (void)
@heading Light Curve
****************************************************************************/
void SortPhotoPoints(long n, Photo3 *seq)
{
long i, j, l, r; /* counters */
int test; /* test variable */
Photo3 tmp; /* temporary storage */
/* assume n .ge. 2 */
if (n < 2) return;
/* H1 initialize */
l = n/2 + 1;
r = n;
while (1) {
/* H2 decrease l or r */
if (l > 1) {
--l;
V2CopyPhoto(&tmp, &seq[l-1]);
} else {
V2CopyPhoto(&tmp, &seq[r-1]);
V2CopyPhoto(&seq[r-1], &seq[0]);
--r;
if (r == 1) {
V2CopyPhoto(&seq[0], &tmp);
return;
}
}
/* H3 prepare for 'sift-up' */
j=l;
/* H4 advance downward */
i = j;
j = 2 * j;
test = OFF;
while ( (j <= r) && (test == OFF)) {
/* H5 find */
if ((j < r) && ( (seq[j-1].phi+FLT_EPSILON) <= seq[j].phi)) ++j;
/* H6 larger than K ? */
if ((tmp.phi+FLT_EPSILON) <= seq[j-1].phi) {
/* H7 move it up */
V2CopyPhoto(&seq[i-1], &seq[j-1]);
/* go back to H4 */
i=j;
j = j * 2;
} else {
test = ON;
}
}
/* H8 store */
V2CopyPhoto(&seq[i-1], &tmp);
} /* return to step H2 */
/* not reached */
return;
}
syntax highlighted by Code2HTML, v. 0.9.1