/* 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 #include #include #include #include #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/: %%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; }