/* 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 <sys/types.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include "Light.h"
/* #ifdef HAVE_SYS_TIME_H */
/* #include <sys/time.h> */
/* #endif */
/* #ifdef HAVE_SYS_TIMES_H */
/* #include <sys/times.h> */
/* #endif */
/* #include <time.h> */
#ifdef TM_IN_SYS_TIME
#include <sys/time.h>
#else
#include <time.h>
#endif
#include <pwd.h>
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
#ifdef _WITH_PGPLOT
#ifndef _WITH_GNUPLOT
#include "cpgplot.h"
#endif
#endif
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Initialize Output flux tables
@param (void)
@return (void)
@heading Output
****************************************************************************/
void InitOutFlux()
{
int j, band; /* Loop count */
for (j = 0; j < PhaseSteps; ++j) { /* loop over phase */
for (band = 0; band < NUM_MAG; ++band) {
FluxOut[j].Mag[band] = 0;
}
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Write Output flux tables
@param (void)
@return (void)
@heading Output
****************************************************************************/
void WriteOutput()
{
FILE *file_out; /* output filehandle */
int j, band; /* Loop count */
double RVp = 0, RVs = 0; /* radial velocities */
file_out = fopen(OUT_FILE, "w");
if (file_out == NULL)
nf_error(_(errmsg[3]));
/* ----------- print header ------------------------------------ */
OutfileHeader(file_out);
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Normalization> at <Phase> was:\n"));
for (band = 0; band < NUM_MAG; ++band) {
fprintf(file_out, _("# Band %3d: %15.6g %15.6g\n"),
band, F90[band], P90[band] );
}
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
if (Flags.monochrome == OFF)
fprintf(file_out, _("# SEQ PHASE U B V R I J H K u v b y RV(Primary) RV(Secondary) \n"));
else
fprintf(file_out, _("# SEQ PHASE %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f RV(Primary) RV(Secondary) \n"),
mono_zero[0], mono_zero[1], mono_zero[2], mono_zero[3],
mono_zero[4], mono_zero[5], mono_zero[6], mono_zero[7],
mono_zero[8], mono_zero[9], mono_zero[10], mono_zero[11]);
fprintf(file_out, "# \n");
/* ----------- print data ------------------------------------ */
for (j = 0; j < PhaseSteps; ++j) { /* loop over phase */
if (Flags.elliptic == OFF){
RVp = -(FluxOut[j].RV[Primary]) /1000.;
RVs = -(FluxOut[j].RV[Secondary])/1000.;
} else {
RVp = (FluxOut[j].RV[Primary]) /1000.;
RVs = (FluxOut[j].RV[Secondary])/1000.;
}
/* 2002-05-03 rwichmann: more significant digits
*/
fprintf(file_out, "%3d %7.4f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %11.4g %11.4g\n",
j, ((FluxOut[j].Phase/(2.0*PI)) - 0.5),
FluxOut[j].Mag[Umag],
FluxOut[j].Mag[Bmag],
FluxOut[j].Mag[Vmag],
FluxOut[j].Mag[Rmag],
FluxOut[j].Mag[Imag],
FluxOut[j].Mag[Jmag],
FluxOut[j].Mag[Hmag],
FluxOut[j].Mag[Kmag],
FluxOut[j].Mag[umag],
FluxOut[j].Mag[vmag],
FluxOut[j].Mag[bmag],
FluxOut[j].Mag[ymag],
RVp,
RVs);
}
(void) fclose(file_out);
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Plot the light curve
@param (void)
@return (void)
@heading Plotting
****************************************************************************/
void PlotOutput()
{
/* >>>>>>>>>>>>>>>>>>>> PLOT <<<<<<<<<<<<<<<<<<<<<<<<< */
#ifdef _WITH_PGPLOT
float MinMag, MaxMag; /* plot scaling */
float MinMagP, MaxMagP; /* plot scaling */
float MinMagS, MaxMagS; /* plot scaling */
float *Xplot, *Yplot, *Zplot; /* data */
float *XPplot, *YPplot; /* residuals */
int j; /* Loop count */
int jstat; /* temp variable */
int FlagBand; /* passband to plot */
char Title[64]; /* title string */
float X1=0.0, X2=1.0; /* ranges */
float Y1=-1.0, Y2=0.0; /* ranges */
float CalcMax = 0.75; /* upper limit (phase) */
char file_out[256+4]; /* ps output file */
float xtick, ytick; /* tickmarks */
sprintf(file_out, "%s%s", Out_Plot_File, "/CPS");
/* --------- allocate memory ------------------------- */
Xplot = malloc(2*PHASESTEPS*sizeof(float));
Yplot = malloc(2*PHASESTEPS*sizeof(float));
Zplot = malloc(2*PHASESTEPS*sizeof(float));
XPplot = malloc(4*MAXOBS*sizeof(float));
YPplot = malloc(4*MAXOBS*sizeof(float));
/* check allocation */
if (Xplot == NULL || Yplot == NULL || Zplot == NULL ||
XPplot == NULL || YPplot == NULL)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
if (Xplot != NULL) free(Xplot);
if (Yplot != NULL) free(Yplot);
if (Zplot != NULL) free(Zplot);
if (XPplot != NULL) free(XPplot);
if (YPplot != NULL) free(YPplot);
make_dialog(_(errmsg[0]));
return;
} else nf_error(_(errmsg[0]));
}
#else
nf_error(_(errmsg[0]));
#endif
/* if we plot a lightcurve */
if (Flags.PlotBand < NUM_MAG) {
MinMag = 100000.; MaxMag = -100000.;
MinMagP = 100000.; MaxMagP = -100000.;
MinMagS = 100000.; MaxMagS = -100000.;
FlagBand = Flags.PlotBand;
/* ---------- loop over phase ----------------------- */
for (j = 0; j < PhaseSteps; ++j) {
Xplot[j] = (FluxOut[j].Phase/(2.0*PI)) - 0.5;
Xplot[j+PhaseSteps] = Xplot[j] + 1.0;
Yplot[j] = - FluxOut[j].Mag[FlagBand];
Yplot[j+PhaseSteps] = - FluxOut[j].Mag[FlagBand];
MinMag = MIN(MinMag,Yplot[j]); MaxMag = MAX(MaxMag,Yplot[j]);
if ((Xplot[j] > 0.45) && (Xplot[j] < 0.55)) {
MinMagP = MIN(MinMagP,Yplot[j]); MaxMagP = MAX(MaxMagP,Yplot[j]);
}
if ((Xplot[j] > -0.1) && (Xplot[j] < 0.1)) {
MinMagS = MIN(MinMagS,Yplot[j]); MaxMagS = MAX(MaxMagS,Yplot[j]);
}
}
CalcMax = Xplot[PhaseSteps-1];
if (Flags.eps == OFF) {
if(cpgopen("/XSERVE") != 1)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
free(Xplot);
free(Yplot);
free(Zplot);
free(XPplot);
free(YPplot);
make_dialog (_(errmsg[2]));
return;
} else nf_error(_(errmsg[2]));
}
#else
nf_error(_(errmsg[2]));
#endif
} else {
if(cpgopen(file_out) != 1)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
free(Xplot);
free(Yplot);
free(Zplot);
free(XPplot);
free(YPplot);
make_dialog (_(errmsg[1]));
return;
} else nf_error(_(errmsg[1]));
}
#else
nf_error(_(errmsg[1]));
#endif
++Flags.eps;
}
/* ------------------ start plot ----------------------- */
#ifdef _WITH_GNUPLOT
gnu_start();
#endif
cpgscf(2); cpgslw(1); cpgsch(0.9);
xtick = 0.0;
if (Flags.plot == ON) {
X1= -0.25; X2=0.75; Y1=MinMag-0.05; Y2=MaxMag+0.05;
xtick = 0.5;
}
else if (Flags.plot == 2 ) {
X1= -0.25; X2=1.75; Y1=MinMag-0.05; Y2=MaxMag+0.05;
xtick = 0.5;
}
if (Flags.plot == 3 ) {
X1= 0.45; X2=0.55; Y1=MinMagP-0.01;Y2=MaxMagP;
xtick = 0.05;
}
if (Flags.plot == 4 ) {
X1= -0.1; X2=0.1; Y1=MinMagS-0.01;Y2=MaxMagS;
xtick = 0.05;
}
ytick = 0.01 * (int)(100 * fabs( (Y1 - Y2) / 5));
if (Flags.Passbands[FlagBand] > 0) {
cpgsvp(0.2, 0.9, 0.4, 0.9);
cpgswin(X1, X2, Y1, Y2);
#ifdef _WITH_GNUPLOT
cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
#else
cpgbox("BCST", 0, 0, "BCNST", 0, 0);
#endif
} else {
cpgsvp(0.2, 0.9, 0.2, 0.9);
cpgswin(X1, X2, Y1, Y2);
#ifdef _WITH_GNUPLOT
cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
#else
cpgbox("BCNST", 0, 0, "BCNST", 0, 0);
#endif
}
#ifdef _WITH_GNUPLOT
cpglab(_("Phase"), _("Delta mag"), "");
sprintf(Title, "%s", Filters[FlagBand]);
cpgtext(0.0, MaxMag, Title);
#endif
/* ---------------------- no data ------------------- */
if (Flags.Passbands[FlagBand] == 0) {
cpgline(2*PhaseSteps, Xplot, Yplot);
}
/* ---------------------- observed data ------------- */
else if (Flags.Passbands[FlagBand] > 0) {
Y1 = Observed[FlagBand][0].residual;
Y2 = Observed[FlagBand][0].residual;
for (j = 0; j < Flags.Passbands[FlagBand]; ++j) {
XPplot[j] = Observed[FlagBand][j].phi;
if (XPplot[j] > CalcMax) XPplot[j] = XPplot[j] -1.0;
YPplot[j] = -Observed[FlagBand][j].mag;
XPplot[j+Flags.Passbands[FlagBand]] = 1+XPplot[j];
YPplot[j+Flags.Passbands[FlagBand]] = YPplot[j];
Y1 = MIN(Y1,Observed[FlagBand][j].residual);
Y2 = MAX(Y2,Observed[FlagBand][j].residual);
}
#ifndef _WITH_GNUPLOT
sprintf(Title, "%s", Filters[FlagBand]);
cpgtext(-0.1, MaxMag, Title);
cpgline(2*PhaseSteps, Xplot, Yplot);
cpgpt(2*j, XPplot, YPplot, 17);
cpglab("", _("Delta mag"), "");
#else
cpglinept(2*PhaseSteps, 2*j, Xplot, Yplot, XPplot, YPplot);
#endif
/* ---------------------- residuals ------------------- */
#ifndef _WITH_GNUPLOT
cpgsvp(0.2, 0.9, 0.2, 0.4);
#else
/* no second plot if hardcopy - multiplot not supp. */
if (Flags.eps < ON){
cpgsvp(0.2, 0.9, 0.0, 0.4);
#endif
if (Y1 == Y2)
{
Y1 = -1.0;
Y2 = 1.0;
}
cpgswin(X1, X2, Y1, Y2);
ytick = 0.01 * (int)(100 * fabs( (Y1 - Y2) / 5));
cpgbox("ABCNST", 0, 0, "BCNST", ytick, 0);
cpglab("", _("Delta mag"), "");
for (j = 0; j < Flags.Passbands[FlagBand]; ++j) {
YPplot[j] = Observed[FlagBand][j].residual;
YPplot[j+Flags.Passbands[FlagBand]] = Observed[FlagBand][j].residual;
}
cpgpt(2*j, XPplot, YPplot, 17);
#ifdef _WITH_GNUPLOT
}
#endif
}
cpglab(_("Phase"), _("Delta mag"), "");
if (Flags.eps != OFF) my_cpgend();
else cpgend();
}
/* >>>>>>>>>>>>>> RADIAL VELOCITY <<<<<<<<<<<<<<<<<<< */
if (Flags.PlotBand == NUM_MAG || Flags.PlotBand == (NUM_MAG+1) ) {
MinMag = FluxOut[0].RV[Primary]/1000.;
MaxMag = FluxOut[0].RV[Primary]/1000.;
MinMagP = FluxOut[0].RV[Primary]/1000.;
MaxMagP = FluxOut[0].RV[Primary]/1000.;
MinMagS = FluxOut[0].RV[Secondary]/1000.;
MaxMagS = FluxOut[0].RV[Secondary]/1000.;
FlagBand = Flags.PlotBand;
/* -------------- loop over phase -------------------- */
for (j = 0; j < PhaseSteps; ++j) { /* loop over phase */
Xplot[j] = (FluxOut[j].Phase/(PI+PI)) - 0.5;
Xplot[j+PhaseSteps] = Xplot[j] + 1.0;
if (Flags.elliptic == ON) {
Yplot[j] = FluxOut[j].RV[Primary]/1000.;
Yplot[j+PhaseSteps] = FluxOut[j].RV[Primary]/1000.;
Zplot[j] = FluxOut[j].RV[Secondary]/1000.;
Zplot[j+PhaseSteps] = FluxOut[j].RV[Secondary]/1000.;
} else {
Yplot[j] = - FluxOut[j].RV[Primary]/1000.;
Yplot[j+PhaseSteps] = - FluxOut[j].RV[Primary]/1000.;
Zplot[j] = - FluxOut[j].RV[Secondary]/1000.;
Zplot[j+PhaseSteps] = - FluxOut[j].RV[Secondary]/1000.;
}
MinMagP = MIN(MinMagP,Yplot[j]); MaxMagP = MAX(MaxMagP,Yplot[j]);
MinMagS = MIN(MinMagS,Zplot[j]); MaxMagS = MAX(MaxMagS,Zplot[j]);
}
CalcMax = Xplot[PhaseSteps-1];
MinMag = MIN(MinMagP,MinMagS);
MaxMag = MAX(MaxMagP,MaxMagS);
/* -------------- open plot -------------------- */
if (Flags.eps == OFF) {
if(cpgopen("/XSERVE") != 1)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
free(Xplot);
free(Yplot);
free(Zplot);
free(XPplot);
free(YPplot);
make_dialog (_(errmsg[2]));
return;
} else nf_error(_(errmsg[2]));
}
#else
nf_error(_(errmsg[2]));
#endif
} else {
if(cpgopen(file_out) != 1)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
free(Xplot);
free(Yplot);
free(Zplot);
free(XPplot);
free(YPplot);
make_dialog (_(errmsg[1]));
return;
} else nf_error(_(errmsg[1]));
}
#else
nf_error(_(errmsg[1]));
#endif
++Flags.eps;
}
/* -------------- plot -------------------- */
#ifdef _WITH_GNUPLOT
gnu_start();
#endif
cpgscf(2); cpgslw(1); cpgsch(0.9);
xtick = 0.0;
if (Flags.plot == ON) {
X1= -0.25; X2=0.75; Y1=MinMag-0.05; Y2=MaxMag+0.05;
xtick = 0.5;
}
else if (Flags.plot == 2 ) {
X1= -0.25; X2=1.75; Y1=MinMag-0.05; Y2=MaxMag+0.05;
xtick = 0.5;
}
if (Flags.plot == 3 ) {
X1= 0.45; X2=0.55; Y1=MinMagP-0.01;Y2=MaxMagP;
xtick = 0.05;
}
if (Flags.plot == 4 ) {
X1= -0.1; X2=0.1; Y1=MinMagS-0.01;Y2=MaxMagS;
xtick = 0.05;
}
if (Flags.Passbands[NUM_MAG] > 0 || Flags.Passbands[NUM_MAG+1] > 0) {
if (Y1 == Y2) { Y1 = -20.; Y2 = 20.; } /* a dirty hack */
if (Flags.Passbands[NUM_MAG+1] == 0 && Flags.plot == ON) {
Y1=MinMagP-0.01; Y2=MaxMagP+0.01;
}
if (Flags.Passbands[NUM_MAG] == 0 && Flags.plot == ON) {
Y1=MinMagS-0.01; Y2=MaxMagS+0.01;
}
cpgsvp(0.2, 0.9, 0.4, 0.9);
cpgswin(X1, X2, Y1, Y2);
ytick = 0.1 * (int)(10 * fabs( (Y1 - Y2) / 5));
#ifdef _WITH_GNUPLOT
cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
#else
cpgbox("BCST", 0, 0, "BCNST", 0, 0);
#endif
} else {
ytick = 0.1 * (int)(10 * fabs( (Y1 - Y2) / 5));
cpgsvp(0.2, 0.9, 0.2, 0.9);
cpgswin(X1, X2, Y1, Y2);
#ifdef _WITH_GNUPLOT
cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
#else
cpgbox("BCNST", 0, 0, "BCNST", 0, 0);
#endif
}
/* -------------- get data -------------------- */
if (Flags.Passbands[NUM_MAG] > 0) {
FlagBand = NUM_MAG;
Y1 = 1e20;
Y2 = -1e20;
for (j = 0; j < Flags.Passbands[FlagBand]; ++j) {
XPplot[j] = Observed[FlagBand][j].phi;
if (XPplot[j] > CalcMax) XPplot[j] = XPplot[j] - 1.0;
YPplot[j] = Observed[FlagBand][j].mag ;
XPplot[j+Flags.Passbands[FlagBand]] = 1+XPplot[j];
YPplot[j+Flags.Passbands[FlagBand]] = YPplot[j];
Y1 = MIN(Y1,Observed[FlagBand][j].residual);
Y2 = MAX(Y2,Observed[FlagBand][j].residual);
}
}
jstat = 2*Flags.Passbands[NUM_MAG];
if (Flags.Passbands[NUM_MAG+1] > 0) {
FlagBand = NUM_MAG+1;
if (Flags.Passbands[NUM_MAG] > 0) {
Y1 = MIN(Y1,Observed[FlagBand][0].residual);
Y2 = MAX(Y2,Observed[FlagBand][0].residual);
} else {
Y1 = 1e20;
Y2 = -1e20;
}
/* count up from current state */
for (j = jstat;
j < (jstat + Flags.Passbands[FlagBand]); ++j) {
XPplot[j] = Observed[FlagBand][j-jstat].phi;
if (XPplot[j] > CalcMax) XPplot[j] = XPplot[j] -1.0;
YPplot[j] = Observed[FlagBand][j-jstat].mag;
XPplot[j+Flags.Passbands[FlagBand]] = 1+XPplot[j];
YPplot[j+Flags.Passbands[FlagBand]] = YPplot[j];
Y1 = MIN(Y1,Observed[FlagBand][j-jstat].residual);
Y2 = MAX(Y2,Observed[FlagBand][j-jstat].residual);
}
}
/* -------------- plot everything -------------------- */
#ifdef _WITH_GNUPLOT
cpglab(_("Phase"), _("Radial Velocity"), "");
if (Flags.Passbands[NUM_MAG+1] > 0 || Flags.Passbands[NUM_MAG] > 0 ) {
cpgline2pt(2*PhaseSteps, 2*PhaseSteps,
(2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]),
Xplot, Yplot, Xplot, Zplot, XPplot, YPplot);
} else {
cpgline2(2*PhaseSteps, 2*PhaseSteps, Xplot, Yplot, Xplot, Zplot);
}
#else
cpgline(2*PhaseSteps, Xplot, Yplot);
cpgline(2*PhaseSteps, Xplot, Zplot);
if (Flags.Passbands[NUM_MAG+1] > 0 || Flags.Passbands[NUM_MAG] > 0 ) {
cpgpt( (2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]) ,
XPplot, YPplot, 17);
}
#endif
cpglab("", _("Radial Velocity"), "");
/* ----------- Now plot the residuals -------------- */
if (Flags.Passbands[NUM_MAG] > 0 ) {
FlagBand = NUM_MAG;
for (j = 0; j < Flags.Passbands[FlagBand]; ++j) {
YPplot[j] = Observed[FlagBand][j].residual;
YPplot[j+Flags.Passbands[FlagBand]] = Observed[FlagBand][j].residual;
}
}
if (Flags.Passbands[NUM_MAG+1] > 0 ) {
FlagBand = NUM_MAG+1;
for (j = jstat; j < (jstat+Flags.Passbands[FlagBand]); ++j) {
YPplot[j] = Observed[FlagBand][j-jstat].residual;
YPplot[j+Flags.Passbands[FlagBand]] =
Observed[FlagBand][j-jstat].residual;
}
}
#ifndef _WITH_GNUPLOT
if (Flags.Passbands[NUM_MAG] > 0 || Flags.Passbands[NUM_MAG+1] > 0) {
cpgsvp(0.2, 0.9, 0.2, 0.4);
cpgswin(X1, X2, Y1, Y2);
cpgbox("ABCNST", 0, 0, "BCNST", 0, 0);
cpgpt( (2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]) ,
XPplot, YPplot, 17);
cpglab(_("Phase"), _("Radial Velocity"), "");
}
#else
if (Flags.Passbands[NUM_MAG] > 0 || Flags.Passbands[NUM_MAG+1] > 0) {
if (Flags.eps < ON){ /* no second plot if hardcopy */
cpgsvp(0.2, 0.9, 0.0, 0.4);
if (Y1 == Y2)
{
Y1 = -1.0;
Y2 = 1.0;
}
cpgswin(X1, X2, Y1, Y2);
ytick = 0.1 * (int)(10 * fabs( (Y1 - Y2) / 5));
cpgbox("ABCNST", 0, 0, "BCNST", ytick, 0);
cpglab("", _("Radial Velocity"), "");
cpgpt( (2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]) ,
XPplot, YPplot, 17);
}
}
#endif
if (Flags.eps != OFF) my_cpgend();
else cpgend();
}
free(Xplot);
free(Yplot);
free(Zplot);
free(XPplot);
free(YPplot);
return;
#endif
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.37
@short Write header of output file
@param (FILE*) file_out The output file
@return (void)
@heading Plotting
****************************************************************************/
void OutfileHeader(FILE *file_out)
{
long j, i; /* loop counters */
int band; /* bandpass */
int runs = 0, ul = 0, ll = 0; /* runs, upper+lower limit */
double Fillout = 0.0; /* Fillout Factor */
double ScaleVolume=1.; /* units conversion */
double ScaleMass; /* units conversion */
double ResMean[NUM_MAG+2]; /* residuals */
double ResDev[NUM_MAG+2]; /* residuals */
double fratio1 = 1.0; /* async rotation */
double fratio2 = 1.0; /* async rotation */
double RochePotential; /* roche potential */
double G_SI = 6.6726e-11; /* gravitational constant */
char WhoAmI[32] = "Unknown"; /* user name */
char *AsciiTime; /* local time */
char deftime[] = "66.6"; /* default time */
/********* UNIX stuff begin *************/
/* */
/* */
time_t time_now;
struct tm *time_ptr;
#ifdef HAVE_UNISTD_H
struct passwd *user_passwd;
user_passwd = getpwuid(geteuid());
strncpy(WhoAmI, user_passwd->pw_gecos, 31);
WhoAmI[31] = '\0';
#endif
(void) (time(&time_now));
time_ptr = localtime(&time_now);
if (time_ptr != NULL) AsciiTime = asctime(time_ptr);
else AsciiTime = &deftime[0];
/********* UNIX stuff end ***************/
if (Flags.fill == ON)
Fillout = (Binary[Primary].RochePot - Binary[Primary].RCLag1)
/ (Binary[Primary].RCLag2 - Binary[Primary].RCLag1);
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# Thou invoked Nightfall, and it answered thy plea.\n"));
fprintf(file_out, "# \n");
/* for some stupid reason, asctime() inserts a terminating EOL */
fprintf(file_out, _("# Thy Name is %s, \n"), WhoAmI );
fprintf(file_out, _("# and thy time is %s"), AsciiTime);
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
if (Flags.elliptic == OFF) {
fprintf(file_out, _("# System Parameters (circular orbit):\n"));
fprintf(file_out, "# \n");
} else {
fprintf(file_out, _("# System Parameters (elliptic orbit):\n"));
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Eccentricity> %8.3f (dimensionless)\n"),
Orbit.Excentricity);
fprintf(file_out, _("# <Omega> %8.3f (degree)\n"),
Orbit.Omega);
}
fprintf(file_out, _("# <Inclination> %8.3f (degree)\n"),
(RTOD * Orbit.Inclination));
fprintf(file_out, _("# <Mass Ratio> %8.3f (dimensionless)\n"),
Binary[Primary].Mq);
fprintf(file_out, _("# <Lagrange One> %8.3f (dimensionless)\n"),
Binary[Primary].RLag1);
/* provide details about the computational model used */
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# Model used:\n"));
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Flux> %s\n"),
Flags.blackbody == ON ? _("Blackbody") : _("Model atmophere"));
if (Flags.blackbody == OFF) {
fprintf(file_out, _("# <Log g Primary> %3.1f\n"),
Binary[Primary].log_g);
fprintf(file_out, _("# <Log g Secondary> %3.1f\n"),
Binary[Secondary].log_g);
}
fprintf(file_out, _("# <Fractional visibility> %s\n"),
Flags.fractional == ON ? _("On") : _("Off"));
fprintf(file_out, _("# <Limb darkening> %s\n"),
Flags.limb == 0 ? _("linear") :
(Flags.limb == 1 ? _("quadratic") :
(Flags.limb == 2 ? _("square root") : _("inverse"))));
fprintf(file_out, _("# <Detailed reflection> %d\n"), Flags.reflect);
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# System Size (absolute units):\n"));
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Period> %8.3e (seconds)\n"),
Orbit.TruePeriod);
fprintf(file_out, _("# <Total Mass> %8.3e (kilogramm)\n"),
Orbit.TrueMass);
fprintf(file_out, _("# <Distance> %8.3e (meter)\n"),
Orbit.TrueDistance);
fprintf(file_out, _("# <Lagrange One> %8.3e (meter)\n"),
Binary[Primary].RLag1 * Orbit.TrueDistance);
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Period> %8.3f (days)\n"),
Orbit.TruePeriod / 86400.);
fprintf(file_out, _("# <Total Mass> %8.3f (solar mass)\n"),
Orbit.TrueMass / 1.989E30);
fprintf(file_out, _("# <Distance> %8.3f (solar radius)\n"),
Orbit.TrueDistance / 6.960E8);
fprintf(file_out, _("# <Lagrange One> %8.3f (solar radius)\n"),
Binary[Primary].RLag1 * Orbit.TrueDistance / 6.960E8);
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
ScaleMass = Binary[Primary].Mq / ( 1.0 + Binary[Primary].Mq);
fprintf(file_out, _("# Component Parameters:\n"));
fprintf(file_out, _("# Primary Secondary\n"));
fprintf(file_out, _("# <Max Velocity> %8.3f %8.3f (km/sec)\n"),
Binary[Primary].K, Binary[Secondary].K);
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Mass> %8.3f %8.3f (dimensionless)\n"),
1.0 - ScaleMass, ScaleMass);
fprintf(file_out, _("# <Gravity> %8.3f %8.3f (dimensionless)\n"),
Binary[Primary].Gravity, Binary[Secondary].Gravity);
/* 2002-05-03 rwichmann: more significant digits
*/
fprintf(file_out, _("# <Polar Radius> %8.3f %8.3f (dimensionless)\n"),
Binary[Primary].Radius, Binary[Secondary].Radius);
if (Flags.fill == OFF)
{
fprintf(file_out,
_("# <Point Radius> %8.6f %8.6f (dimensionless)\n"),
Binary[Primary].Xp, Binary[Secondary].Xp);
fprintf(file_out,
_("# <Mean Radius> %8.6f %8.6f (dimensionless)\n"),
Binary[Primary].RadMean, Binary[Secondary].RadMean);
fprintf(file_out,
_("# <Volume> %8.6f %8.6f (dimensionless)\n"),
Binary[Primary].Volume, Binary[Secondary].Volume);
}
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Mass> %8.2e %8.2e (kilogramm)\n"),
Orbit.TrueMass * (1.0 - ScaleMass), Orbit.TrueMass * ScaleMass);
fprintf(file_out, _("# <Gravity> %8.3f %8.3f (m kg/s^2)\n"),
Binary[Primary].Gravity
* Orbit.TrueMass * (1.0 - ScaleMass)
* G_SI
/ (Orbit.TrueDistance * Orbit.TrueDistance),
Binary[Secondary].Gravity
* Orbit.TrueMass * (1.0 - ScaleMass)
* G_SI
/ (Orbit.TrueDistance * Orbit.TrueDistance));
fprintf(file_out, _("# <Polar Radius> %8.2e %8.2e (meter)\n"),
Binary[Primary].Radius * Orbit.TrueDistance,
Binary[Secondary].Radius * Orbit.TrueDistance);
if (Flags.fill == OFF)
{
fprintf(file_out, _("# <Point Radius> %8.2e %8.2e (meter)\n"),
Binary[Primary].Xp * Orbit.TrueDistance,
Binary[Secondary].Xp * Orbit.TrueDistance);
fprintf(file_out, _("# <Mean Radius> %8.2e %8.2e (meter)\n"),
Binary[Primary].RadMean * Orbit.TrueDistance,
Binary[Secondary].RadMean * Orbit.TrueDistance);
ScaleVolume = Orbit.TrueDistance * Orbit.TrueDistance *
Orbit.TrueDistance;
fprintf(file_out,
_("# <Volume> %8.2e %8.2e (cubic meter)\n"),
Binary[Primary].Volume * ScaleVolume,
Binary[Secondary].Volume * ScaleVolume);
}
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Mass> %8.3f %8.3f (solar mass)\n"),
Orbit.TrueMass * (1.0 - ScaleMass) / 1.989E30,
Orbit.TrueMass * ScaleMass / 1.989E30 );
fprintf(file_out, _("# <Polar Radius> %8.3f %8.3f (solar radius)\n"),
Binary[Primary].Radius * Orbit.TrueDistance / 6.960E8,
Binary[Secondary].Radius * Orbit.TrueDistance / 6.960E8);
if (Flags.fill == OFF) {
fprintf(file_out, _("# <Point Radius> %8.3f %8.3f (solar radius)\n"),
Binary[Primary].Xp * Orbit.TrueDistance / 6.960E8,
Binary[Secondary].Xp * Orbit.TrueDistance / 6.960E8);
fprintf(file_out, _("# <Mean Radius> %8.3f %8.3f (solar radius)\n"),
Binary[Primary].RadMean * Orbit.TrueDistance / 6.960E8,
Binary[Secondary].RadMean * Orbit.TrueDistance / 6.960E8);
ScaleVolume = ScaleVolume / (4.0 * PI * 6.960E8 * 6.960E8 * 6.960E8 / 3.0);
fprintf(file_out, _("# <Volume> %8.3f %8.3f (solar volumes)\n"),
Binary[Primary].Volume * ScaleVolume,
Binary[Secondary].Volume * ScaleVolume);
}
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Mean Temp.> %8.1f %8.1f (Kelvin)\n"),
Binary[Primary].TempMean, Binary[Secondary].TempMean);
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# Primary Secondary\n"));
if (Flags.fill == OFF) {
fprintf(file_out, _("# <Fill Factor> %8.3f %8.3f (dimensionless)\n"),
Binary[Primary].RocheFill, Binary[Secondary].RocheFill);
} else {
fprintf(file_out, _("# <Fill Factor> %8.3f (dimensionless)\n"),
Binary[Primary].RocheFill);
}
fprintf(file_out, _("# <Temperature> %8.1f %8.1f (Kelvin)\n"),
Binary[Primary].Temperature,
Binary[Secondary].Temperature);
if (Flags.asynchron1 == ON) fratio1 = Binary[Primary].Fratio;
if (Flags.asynchron2 == ON) fratio2 = Binary[Secondary].Fratio;
fprintf(file_out, _("# <Asynchronicity>%8.3f %8.3f (dimensionless)\n"),
fratio1, fratio2);
fprintf(file_out, _("# <Albedo> %8.3f %8.3f (dimensionless)\n"),
Binary[Primary].Albedo, Binary[Secondary].Albedo);
fprintf(file_out, _("# <Grav. Dark.> %8.3f %8.3f (dimensionless)\n"),
Binary[Primary].GravDarkCoeff ,
Binary[Secondary].GravDarkCoeff);
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out,
_("# <Surface Potential Primary> %10.4f (dimensionless)\n"),
Binary[Primary].RochePot);
fprintf(file_out,
_("# <Lagrange 1 Potential> %10.4f (dimensionless)\n"),
Binary[Primary].RCLag1);
fprintf(file_out,
_("# <Lagrange 2 Potential> %10.4f (dimensionless)\n"),
Binary[Primary].RCLag2);
if (Flags.fill == ON) {
fprintf(file_out,
_("# <Fillout Factor> %10.4f (dimensionless)\n"),
Fillout);
}
fprintf(file_out,
_("# <Surface Potential Secondary> %10.4f (dimensionless)\n"),
Binary[Secondary].RochePot);
RochePot = 0.0; /* Binary[Primary].RochePot; */
Mq = Binary[Primary].Mq;
RochePotential = RochePerpendSecond(Binary[Secondary].Radius);
/*
RochePotential = 1.0/Binary[Secondary].Radius
+ Binary[Secondary].Mq/sqrt(1 + SQR(Binary[Secondary].Radius));
*/
fprintf(file_out,
_("# <... in Primary frame> %10.4f (dimensionless)\n"),
RochePotential);
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Luminosities>:\n"));
for (band = 0; band < NUM_MAG; ++band) {
fprintf(file_out, _("# Band %3d: %15.6g %15.6g"),
band, L90[Primary][band], L90[Secondary][band] );
#ifdef HAVE_DISK
if (Flags.disk == ON)
{
fprintf(file_out, " %15.6g\n", L90[Disk][band]);
}
else
{
fprintf(file_out, "\n");
}
#else
fprintf(file_out, "\n");
#endif
}
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# Spots (if any)\n"));
fprintf(file_out, "# \n");
if (Flags.Spots1 > 0) {
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Spots on Primary>\n"));
fprintf(file_out,
_("# Longitude Latitude Radius DimFactor\n"));
for (j = 0; j < Flags.Spots1; ++j) {
fprintf(file_out, _("# <Spot %3ld> %8.3f %8.3f %8.3f %8.3f\n"),
j,
Spot[Primary][j].longitude,
Spot[Primary][j].latitude,
Spot[Primary][j].radius,
Spot[Primary][j].dimfactor );
}
fprintf(file_out, "# \n");
}
if (Flags.Spots2 > 0) {
fprintf(file_out, "# \n");
fprintf(file_out, _("# <Spots on Secondary>\n"));
fprintf(file_out,
_("# Longitude Latitude Radius DimFactor\n"));
for (j = 0; j < Flags.Spots2; ++j) {
fprintf(file_out, _("# <Spot %3ld> %8.3f %8.3f %8.3f %8.3f\n"),
j,
Spot[Secondary][j].longitude,
Spot[Secondary][j].latitude,
Spot[Secondary][j].radius,
Spot[Secondary][j].dimfactor );
}
}
fprintf(file_out, "# \n");
fprintf(file_out, "# -------------------------------------------------\n");
fprintf(file_out, "# \n");
fprintf(file_out, _("# Fit Results (if any)\n"));
fprintf(file_out, "# \n");
if (Flags.IsComputed == ON) {
for (band=0; band < (NUM_MAG+2); ++band) {
if (Flags.Passbands[band] > 0) {
Runs(Observed[band], Flags.Passbands[band], &runs, &ul, &ll);
ResMean[band] = 0.0; ResDev[band] = 0.0;
for (i = 0; i < Flags.Passbands[band]; ++i) {
ResMean[band] =
ResMean[band] + Observed[band][i].residual;
}
ResMean[band] = ResMean[band]/Flags.Passbands[band];
for (i = 0; i < Flags.Passbands[band]; ++i) {
ResDev[band] = ResDev[band] +
SQR(ResMean[band]-Observed[band][i].residual);
}
ResDev[band] = sqrt(ResDev[band])/(Flags.Passbands[band]-1);
fprintf(file_out,
_("# <%12s> Mean Residual %7.4f SDV Residuals %7.4f\n"),
Filters[band], ResMean[band], ResDev[band]);
fprintf(file_out,
_("# Runs %5d Upper Limit %5d Lower Limit %6d\n"),
runs, ul, ll);
fprintf(file_out,"# \n");
}
}
}
/*@i@*/ return;
}
syntax highlighted by Code2HTML, v. 0.9.1