/* 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;

}