/* 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 <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "Light.h"
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Read input data stream, return number of data read
@tip Number of data goes to Flags.Passbands[Band]
@param (char*) InputFile The name of the input file
@return (void)
@heading Data Fitting
****************************************************************************/
long DataRead(FILE *in_file, int *newstart)
{
long in_test; /* # of data on line (2 or 3) */
long in_count; /* # of data lines */
long i; /* loop variable */
long in_countall; /* # of lines */
int Band; /* bandpass */
double NormPoint; /* The normalisation phase */
double Period; /* The period */
double ZeroP; /* The Zeropoint */
double SysVel; /* The system velocity */
double CommWeight; /* Weight for bandpass */
double MagShift; /* Normalization correction */
char in_line[256]; /* single input line */
char junk[256]; /* temporary string */
char BandChar; /* bandpass */
char ErrorMessage[512]; /* error message */
double x, y, z; /* date, value, error */
double MinTime, MaxTime, MinMag, MaxMag; /* for normalization */
int HeadP = 0, HeadZ = 0, HeadB = 0; /* Flags */
double *thisTime, *thisMag, *thisErr; /* input data */
long prevCount; /* # of data already available*/
/* set default values for bandpass, period and phase_zero */
Band = Vmag;
BandChar = 'V';
prevCount = Flags.Passbands[Vmag];
Period = 1.0;
ZeroP = 0.0;
SysVel = 0.0;
CommWeight = 0.01;
MagShift = 0.0;
NormPoint = -666.0;
/* initial values for normalization */
MinTime = 0; MaxTime = 1; MinMag = 0; MaxMag = 1;
*newstart = 0;
/* --------------------------- allocate memory -------------------------- */
thisTime = malloc(MAXOBS * sizeof(double));
thisMag = malloc(MAXOBS * sizeof(double));
thisErr = malloc(MAXOBS * sizeof(double));
if (thisTime == NULL || thisMag == NULL || thisErr == NULL) {
#ifdef _WITH_GTK
if ( Flags.interactive == ON && Flags.parseCL == OFF ) {
make_dialog(_(errmsg[0]));
if (thisTime != NULL) free(thisTime);
if (thisMag != NULL) free(thisMag);
if (thisErr != NULL) free(thisErr);
return (-1);
}
else nf_error(_(errmsg[0]));
#else
nf_error(_(errmsg[0]));
#endif
}
/* --------------------------- read lines ------------------------------- */
in_count = 0; in_countall = 0; in_line[0] = '#';
while ( fgets(in_line, (int) sizeof(in_line), in_file) != NULL) {
/* line number */
++in_countall;
/* this is not a header line */
if (in_line[0] != '#') {
in_test = sscanf(in_line, "%lf%lf%lf", &x, &y, &z);
/* input error */
if (in_test < 2) {
sprintf(ErrorMessage,
_("Invalid line (%ld) in input file, must be Phase/Mag/<Error>: %%f %%f <%%f>\n"),
in_countall);
WARNING(ErrorMessage);
} else {
/* read data */
if (HeadP == 0)
WARNING(_("Data without period in header, use default (1)"));
if (HeadZ == 0)
WARNING(_("Data without zero phase in header, use default (0)"));
if (HeadB == 0)
WARNING(_("Data without passband in header, use default (V)"));
/* if no weight, set weight to some default */
if (in_test == 2 && Band < NUM_MAG)
{
z = CommWeight;
}
else {
if (in_test == 2)
z = 100.*CommWeight;
}
/* store data in array */
if ( (in_count + prevCount) <= (MAXOBS-1)) {
thisErr[in_count] = z;
thisMag[in_count] = y;
thisTime[in_count] = x;
if (in_count == 0) {
MinTime = x;
} else { MinTime = MIN(x,MinTime);}
if (in_count == 0) {
MaxTime = x;
} else { MaxTime = MAX(x,MaxTime);}
if (in_count == 0) {
MinMag = y;
} else { MinMag = MIN(y,MinMag);}
if (in_count == 0) {
MaxMag = y;
} else { MaxMag = MAX(y,MaxMag);}
++in_count;
} else {
sprintf(ErrorMessage,
_("%s: Data Maximum (%d) exceeded, Data discarded"),
Filters[Band], MAXOBS);
WARNING(ErrorMessage);
}
}
/* this is a header line */
} else {
/* --------------- the period --------------------------- */
if (in_line[1] == 'P') {
in_test = sscanf(in_line, "%s%lf", junk, &Period);
if (in_test < 2)
WARNING(_("Wrong Period Specification in Input File"));
HeadP = 1;
}
/* --------------- phase zero --------------------------- */
if (in_line[1] == 'Z') {
in_test = sscanf(in_line, "%s%lf", junk, &ZeroP);
if (in_test < 2)
WARNING(_("Wrong Zero_Phase Specification in Input File"));
HeadZ = 1;
}
/* --------------- system velocity ---------------------- */
if (in_line[1] == 'V') {
in_test = sscanf(in_line, "%s%lf", junk, &SysVel);
if (in_test < 2)
WARNING(_("Wrong System Velocity Specification in Input File"));
}
/* -------------- magnitude shift ----------------------- */
if (in_line[1] == 'S') {
in_test = sscanf(in_line, "%s%lf", junk, &MagShift);
if (in_test < 2)
WARNING(_("Wrong Magnitude Shift Specification in Input File"));
}
/* -------------- normalisation phase ------------------- */
if (in_line[1] == 'N') {
in_test = sscanf(in_line, "%s%lf", junk, &NormPoint);
if (in_test < 2)
WARNING(_("Wrong Normalisation Point Specification in Input File"));
}
/* ------------- common weight factor ------------------- */
if (in_line[1] == 'W') {
in_test = sscanf(in_line, "%s%lf", junk, &CommWeight);
if (in_test < 2)
WARNING(_("Wrong Common Error Specification in Input File"));
}
/* --------------- passband ----------------------------- */
if (in_line[1] == 'B') {
in_test = sscanf(in_line, "%s %c", junk, &BandChar);
if (in_test < 2)
WARNING(_("Wrong Bandpass Specification in Input File"));
HeadB = 1;
if (BandChar == 'U') Band = Umag;
else if (BandChar == 'B') Band = Bmag;
else if (BandChar == 'V') Band = Vmag;
else if (BandChar == 'R') Band = Rmag;
else if (BandChar == 'I') Band = Imag;
else if (BandChar == 'J') Band = Jmag;
else if (BandChar == 'H') Band = Hmag;
else if (BandChar == 'K') Band = Kmag;
else if (BandChar == 'u') Band = umag;
else if (BandChar == 'b') Band = bmag;
else if (BandChar == 'v') Band = vmag;
else if (BandChar == 'y') Band = ymag;
else if (BandChar == '1') Band = NUM_MAG; /* radial velocity */
else if (BandChar == '2') Band = NUM_MAG+1; /* radial velocity */
/* ------ set band to plot equal to inputfile band ---------------- */
if (Band < NUM_MAG && Flags.PlotBand == -1 ) Flags.PlotBand = Band;
prevCount = Flags.Passbands[Band];
}
}
} /* END file read */
/* store number of data */
if (in_count > 0) {
for (i = 0; i < in_count; ++i) {
Observed[Band][i + prevCount].weight = (float) thisErr[i];
Observed[Band][i + prevCount].mag = (float) thisMag[i];
Observed[Band][i + prevCount].time = thisTime[i];
Observed[Band][i + prevCount].phi = 0;
}
Flags.Passbands[Band] = in_count + prevCount;
} else {
#ifdef _WITH_GTK
if ( Flags.interactive == ON && Flags.parseCL == OFF ) {
make_dialog(_(errmsg[8]));
free(thisTime);
free(thisMag);
free(thisErr);
return (-1);
} else nf_error(_(errmsg[8]));
#else
nf_error(_(errmsg[8]));
#endif
}
/* set the normalisation point for this band */
if (NormPoint > -666.0)
SetNormalPhase (Band, NormPoint);
else
ResetNormalPhaseOne (Band);
/* Adjust Phase and Magnitudes */
for (i = 0; i < in_count; ++i) {
/* shift magnitude */
if (Band < NUM_MAG)
Observed[Band][i + prevCount].mag =
Observed[Band][i + prevCount].mag-(MinMag+MagShift);
/* shift velocity */
if (Band >= NUM_MAG)
Observed[Band][i + prevCount].mag =
Observed[Band][i + prevCount].mag-SysVel;
/* get phi in the range [0,1) */
Observed[Band][i + prevCount].phi =
(float) ((Observed[Band][i + prevCount].time - ZeroP)/Period
- floor((Observed[Band][i + prevCount].time - ZeroP)/Period));
}
free(thisTime);
free(thisMag);
free(thisErr);
return (in_count);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Read input data file
@tip Number of data goes to Flags.Passbands[Band]
@param (char*) InputFile The name of the input file
@return (void)
@heading Data Fitting
****************************************************************************/
void Read(const char * InputFile)
{
char fullPath[1024]; /* full path to file */
FILE *in_file; /* input filehandle */
long read_count; /* number of data read from stream */
int dummy = 0; /* unused here */
/* ------------------------- open file ---------------------------------- */
strncpy(fullPath, InputFile, sizeof(fullPath)-1);
in_file = fopen(InputFile, "r");
if (in_file == NULL) {
strcpy(fullPath, "./data/");
strncat(fullPath, InputFile, sizeof(fullPath)-strlen(InputFile)-8);
in_file = fopen(fullPath, "r");
if (in_file == NULL) {
strncpy(fullPath, data_data_fls(), sizeof(fullPath)-strlen(InputFile)-7);
/* strcat(fullPath, "/data/"); */
strcat(fullPath, "/");
strcat(fullPath, InputFile);
in_file = fopen(fullPath, "r");
}
}
if (in_file == NULL) {
#ifdef _WITH_GTK
if ( Flags.interactive == ON && Flags.parseCL == OFF ) {
make_dialog(_(errmsg[9]));
return;
}
else nf_error(_(errmsg[9]));
#else
nf_error(_(errmsg[9]));
#endif
}
#ifdef _WITH_GTK
if ( Flags.interactive == ON && Flags.parseCL == OFF )
my_appbar_push(_("Reading ...") );
#endif
/* ------------------------- call reader -------------------------------- */
/* setlocale (LC_NUMERIC, "C"); */
read_count = DataRead(in_file, &dummy);
/* setlocale (LC_NUMERIC, ""); */
if (read_count > 0) AddToList(InputFile, NF_NATIVE_FILE);
/* ------------------------- close file -------------------------------- */
(void) fclose(in_file);
#ifdef _WITH_GTK
if ( Flags.interactive == ON && Flags.parseCL == OFF )
my_appbar_push(_("Reading ... finished"));
#endif
return;
}
syntax highlighted by Code2HTML, v. 0.9.1