/* 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 "config.h"

/* This is actually nowhere used, but causes a 'redefined'   */
/* error due to a bug in the JPEG lib headers.               */
#undef HAVE_STDLIB_H

/* Fourth LD law is currently for experimental purposes only */
#define LD_LAW_MAX  2

/* --------------------------------------------------------- */
/*                                                           */
/* Default values that may be changed                        */
/*                                                           */
/* --------------------------------------------------------- */

/* default number of steps per pi (in Eta)                   */
/* pi = 3.14 = 180 degree                                    */
/* steps in Phi are calculated as 10+2*STEPS_PER_PI/sin(Eta) */  
/* -- and corresponding maximum number of surface elements   */
/*                                                           */

#ifdef HIGH_PRECISION
#define  STEPS_PER_PI  196
#define  MAXELEMENTS   65536
#else
#define  STEPS_PER_PI  64
#ifdef HAVE_DISK
#define  MAXELEMENTS   65536
#else
#define  MAXELEMENTS   8192
#endif
#endif


/* window size if Gnuplot                                    */
/* for pgplot, width will be used to scale display           */
/*   unless environment variable PGPLOT_XW_WIDTH is set      */
/*                                                           */
#define GNU_GEOMETRY "550x424+300+20"

/* maximum number of steps for a full orbit                  */
/*                                                           */
#define  PHASESTEPS   2048

/* maximum number of observations in each passband/filter    */
/*                                                           */
#define  MAXOBS  8192 

/* maximum number of spots per component                     */
/*                                                           */
#define  N_SPOT  10

/* output file  for lightcurve                               */
/*                                                           */
#define  OUT_FILE "NightfallCurve.dat"

/* output file  for fit results                              */
/*                                                           */
#define  OUT_SIMPLEX "NightfallFitted.dat"

/* output file  for chi square map                           */
/*                                                           */
#define  OUT_CHIMAP "NightfallChiMap.dat"

/* Grid size (square) for Chi Square Map                     */
/*                                                           */
#define _CHI_SCANS_ 16

/* default wavelength for line profile  (V band)             */
/*                                                           */
#define  LAMDAZERO_DEFAULT 550.0

/* array  for line profile  (multiple of two)                */
#define  PROFILE_ARRAY 128

/* resolution of profile (delta_lambda/lambda)               */
#define  PROFILE_RES 20000.0

/* switch model atmospheres  PHOENIX -> ATLAS                */
#define  SWITCH_TEMP 9000.0

/* --------------------------------------------------------- */
/*                                                           */
/* nothing to change below this line                         */
/*   (unless you know what you do ... )                      */
/*                                                           */
/* --------------------------------------------------------- */

/*  Maximum Number of Iterations for RootFind/MinFind        */
#define MAXITER 100

/*  Maximum linesize for config file                         */ 
#define  MAX_CFG_INLINE 255
       
/*  Floating-Point precision                                 */
#define EPS 3.0e-8         

/* number of passbands                                       */
#define  NUM_MAG 12

/* maximum iteration for downhill simplex                    */
#define  MAX_ITER_SIMPLEX 200

/* default tolerance  for downhill simplex                   */
#define  SIMPLEX_FIT_TOLERANCE 0.01

/* the maximum dimension of the parameter space              */
/* -- Careful: some calls to array indices are               */
/*    hardcoded (LightSimplex.c), decreasing SDIM            */
/*    can lead to disaster                                   */

#define SDIM 49

/************************************************************************/
/*                                                                      */
/* TYPEDEFS                                                             */
/*                                                                      */
/************************************************************************/

#define  WARN     0
#define  BUSY     1
#define  VERBOSE  2
#define  GEO      3
#define  OPTIMIZE 4
#define  MINFIND  5
#define  CHISQR   6
#define  INFILE   7
#define  SURFACE  8
#define  MAX_FVAL 9

#define  NF_NATIVE_FILE 0

/* -----------------  A 2D Point      -------------------------------   */

typedef struct Point2Struct {   /* 2D Point */
   double x, y;
   } Point2;
typedef Point2 Vector2;

/* -----------------  Output Table    -------------------------------   */

typedef struct PhotoOutStruct { 
   float  Phase;                     /* phase  of output datum          */
   double Mag[NUM_MAG];              /* magnitudes of output data       */
   float  RV[3];                     /* radial  velocities              */
                                     /* 2->3 components 2002/04/09, MPR */
   float  D_RV[3];                   /* disturbance of Kepler RV        */
                                     /* 2->3 components 2002/04/09, MPR */
} PhotoOut;

/* -----------------  Data Point      -------------------------------   */

typedef struct PhotoStruct {   
   double time;                      /* time of observation             */
   float  mag;                       /* value of observation            */
   float  weight;                    /* weight of observation           */
   float  phi;                       /* phase  of observation           */
   float  residual;                  /* residual of observation         */
   } Photo3;

/* -----------------  Spot Parameters -------------------------------   */

typedef struct SpotStruct { 
   double nowlongitude;               /* actual value                    */
   double longitude;                  /* latitude of spot                */
   double latitude;                   /* longitude of spot               */
   double radius;                     /* angular radius                  */
   double dimfactor;                  /* dimming (= T.new/T.old)         */
} SpotLight;
 
/* -----------------  Surface Element -------------------------------   */

typedef struct SurfaceStruct {      
   float  rad;                       /* radius                          */
   float  eta;                       /* angle to x-axis                 */
   float  phi;                       /* angle in y-z-plane              */
   float  lambda;                    /* holds rad*cos(eta)         = x  */
   float  mu;                        /* holds rad*sin(eta)cos(phi) = y  */
   float  nu;                        /* holds rad*sin(eta)sin(phi) = z  */
   float  l;                         /* surface normal x                */
   float  m;                         /* surface normal y                */
   float  n;                         /* surface normal z                */
   float  rho;                       /* distance to other stars centre  */
#ifdef HIGH_PRECISION
   double area;                      /* area of element                 */
#else
   float  area;                      /* area of element                 */
#endif
   float  grav;                      /* dimensionless gravity at element*/
   double temp;                      /* temperature of element          */
   double orig_temp;                 /* temperature of element          */
   double f_[NUM_MAG];               /* intensity in various passbands  */
   long   MyNum;                     /* number                          */
   int    ring;                      /* ring this element lives in      */
   struct SurfaceStruct *phi_ptr;    /* pointer to last in phi          */
   struct SurfaceStruct *phi1_ptr;   /* pointer to next in phi          */
   struct SurfaceStruct *next_ptr;   /* pointer to next in eta          */
   struct SurfaceStruct *last_ptr;   /* pointer to last in eta          */
   float  LOS_Pot;                   /* Potential min. on LOS           */
   float  visibility;                /* (fractional) visibility         */
   float  CosGamma;                  /* angle LOS-surface normal        */
   int    EclipseFlag;               /* eclipsed = ON/OFF               */
   int    SpotNum;                   /* number of overlapping spots     */
   double SumDimFactor;              /* average dimming factor          */
   float  Velocity;                  /* projected velocity              */
 
   /* Parameter needed for eclipse testing                              */
   int AreaType;                     /* AreaType: 
					 (TOP_SEGMENT/BOTTOM_SEGMENT/
					  OUT_RECTANGLE/IN_RECTANGLE)   */
   double RadiusIn;                  /* Inner radius of area element PR */
   double RadiusOut;                 /* Outer radius of area element PR */
   double AreaAngle;                 /* Angle of area element between 
					x-axis and middle of the element*/
   double AreaHeight;
   double AreaWidth;
   double VecE1[3];                  /* Element vector 1 of test layer  */
   double VecE2[3];                  /* Element vector 2 of test layer  */
} SurfaceElement;

/* -----------------  Disk Component --------------------------------   */
/* disk parameters MK 14.03.2001 */

/*  typedef struct DiskStruct { */
/*    double                 Incl;         */ /* disk inclination           */
/*                                         */ /* in degrees                 */ 
/*    double                 Warp;         */ /* warp angle of the disk     */ 
/*                                         */ /* in degrees                 */ 
/*    double                 Thick;        */ /* thickness of disk          */ 
/*                                         */ /* in units of R              */ 
/*    double                 HR;           */ /* H over R for the disk      */ 
/*    double                 Rout;         */ /* outer radius of the disk   */ 
/*                                         */ /* in units of seperation     */ 
/*    double                 Rin;          */ /* inner radius of the disk   */ 
/*                                         */ /* in units of seperation     */ 
/*  } DiskComponent; */

/* -----------------  Binary Component  -----------------------------   */

typedef struct ComponentStruct  { 
  double                 Mq;              /* mass ratio                 */
  double                 RocheFill;       /* Roche Lobe filling Factor  */
  double                 RocheStore;      /* Roche Lobe filling Factor  */
  double                 RLag1;           /* Lagrange One               */
  double                 RLag2;           /* Lagrange Two               */
  double                 RLag3;           /* Lagrange Three             */
  double                 RCLag1;          /* Lagrange One   Potential   */
  double                 RCLag2;          /* Lagrange Two   Potential   */
  double                 RCLag3;          /* Lagrange Three Potential   */
  double                 CosPhiLagrange;  /* Lim. angle for eclipse     */
  double                 RZatL1;          /* Throat Rad perpend. L1     */
  double                 RYatL1;          /* Throat Rad perpend. L1     */
  double                 LimRad;          /* Radius at Throat           */
  double                 LimEta;          /* Eta    at Throat           */
  double                 Fratio;          /* Nonsync Rotation           */
  double                 RXCrit;          /* Max. Radius on X axis      */
  double                 RochePotCrit;    /* Potential for Max Lobe     */
  double                 RochePot;        /* Roche Potential-Star       */
  double                 RCrit;           /* Radius for Max Lobe        */
  double                 Radius;          /* Polar Radius Star          */
  double                 Surface;         /* Surface Area               */
  double                 Gravity;         /* MeanGravity                */
  double                 Temperature;     /* Effective Temperature      */
  double                 Albedo;          /* Bolometric Albedo          */
  double                 GravDarkCoeff;   /* Gravitational Darkening    */
  double                 DarkeningGrav;   /* Denominator for dark. law  */
  long                   NumElem;         /* Number of Surface Elements */
  long                   NumPlus;         /* plus the Throat ...        */
  long                   N_PhiStep[STEPS_PER_PI+10];
                                          /* steps in phi at eta        */
  int                    N_Rings;         /* number of rings (in eta)   */
                                          /* steps in phi at eta        */
  double                 Xp;              /* Smallest enclosing radius  */
  double                 Volume;          /* volume of star             */
  double                 ScaledVolume;    /* scaled volume              */  
  double                 RadMean;         /* mean radius                */
  double                 TempMean;        /* mean temperature           */
  double                 K;               /* half-amplitude of radial v */
  double                 log_g;           /* log_g for model atmosphere */

  /* --------------------  disk parameters  --------------------------  */

  double                 Tilt;            /* disk inclination           */
                                          /* in degrees                 */
  double                 Warp;            /* warp angle of the disk     */
                                          /* in degrees                 */
  double                 Thick;           /* thickness of disk          */
                                          /* in units of R              */
  double                 HR;              /* H over R for the disk      */
  double                 Rout;            /* outer radius of the disk   */ 
                                          /* in units of seperation     */
  double                 Rin;             /* inner radius of the disk   */
                                          /* in units of seperation     */
  int                    RingsOut;        /* Number of rings on the     
                                             outer side of the disk     */
  int                    RingsIn;         /* Number of Rings on the     *
                                             inner side of the disk     */
  int                    Segments;        /* Number of Segments in Phi  
				             direction                  */

  /* --------------------  hot spot parameters  ----------------------  */

  double                 tempHS;          /* Hot Spot temperature       */
  double                 longitudeHS;     /* Hot Spot longitude         */
  double                 extentHS;        /* Hot Spot extent            */
  double                 depthHS;         /* Hot Spot depth             */

 } BinaryComponent;

/* -----------------  Flags   ---------------------------------------   */
/* added disk related flags MK 13.05.01 */

typedef struct FlagsStruct {
  unsigned long  debug[MAX_FVAL];  /* debug flags                       */
  int       disk;                  /* want disk geometry (ON/OFF)       */
  int       warp;                  /* warped disk (ON/OFF)              */
  int       plot;                  /* want lightcurve (ON/OFF)          */
  int       visualize;             /* want visualize (ON/OFF)           */
  int       animate;               /* want animated (ON/OFF)            */
#ifdef _WITH_OPENGL
  int       axes;                  /* want x-,y-, z-axis (ON/OFF)       */
  int       movie;                 /* want jpeg movie (ON/OFF)          */
  int       texture;               /* want textures (ON/OFF)            */
  int       textype;               /* what kind of textures (IMAGE,TEMP)*/
  int       frame;                 /* number of jpeg frame              */
  int       wireframe;             /* want wire frame or solid rendering*/
  int       points;                /* use POINTS instead of QUAD_STRIP  */
  int       labels;                /* want labels (ON/OFF)              */
  float     GLdistancez;           /* distance of observer in           */
                                   /* z-direction                       */
  float     GLdistancey;           /* distance of observer in           */
                                   /* y-direction                       */
  float     GLdistancex;           /* distance of observer in           */
                                   /* x-direction                       */
  int       use_opengl;            /* really use it ? (ON/OFF)          */
#endif
  int       interactive;           /* want interactive (ON/OFF)         */
  int       fill;                  /* overfill roche lobe (ON/OFF)      */
  int       lineprofile;           /* want  line profile (ON/OFF)       */
  int       smap;                  /* want  surface map (ON/OFF)        */
  int       monochrome;            /* want  monochrome wavelen (ON/OFF) */
  int       limb;                  /* limb darkening method             */
  int       reflect;               /* reflection treatment              */
  int       fractional;            /* fractional visibility (ON/OFF)    */
  int       blackbody;             /* blackbody approximation (ON/OFF)  */
  int       InEclipse1;            /* we are in eclipse w/ Primary      */
  int       InEclipse2;            /* we are in eclipse w/ Secondary    */
  int       asynchron1;            /* Primary rotates asynchron         */
  int       asynchron2;            /* Secondary rotates asynchron       */
  int       Spots1;                /* num. of spots on primary          */
  int       Spots2;                /* num. of spots on secondary        */
  int       elliptic;              /* is eccentric orbit                */
  int       first_pass;            /* first pass in main_loop()         */
  long      Passbands[NUM_MAG+2];  /* number of data read in passband N */
  int       PlotBand;              /* which passband to plot ?          */
  int       simplex[64];           /* parameter set to fit              */
  int       WantFit;               /* want to fit something  (ON/OFF)   */
  float     SimplexTol;            /* how good is the fit ?             */
  float     Step[2];               /* Stepsizes for Map                 */
  int       WantMap;               /* we want to map chi square         */
  int       anneal;                /* do simulated annealing            */
  int       eps;                   /* EPS plot                          */
  int       IsComputed;            /* Lightcurve is computed            */
  int       ProComputed;           /* Profile is computed               */
  int       parseCL;               /* do we parse the command line      */
  int       plotOpened;            /* Gnuplot window opened             */
  int       GN_exitconfirm;        /* whether confirm exit              */   
  int       GN_rempos;             /* whether remember position         */   
  int       GN_tooltips;           /* whether have tooltips             */
  int       GL;                    /* does the display support OpenGL   */
  char      ConfFile[256];         /* configuration file                */
  int       tdisk;                 /* disk temperature distribution     */
} FlagsHandle;

/* -----------------  The Orbit -------------------------------------   */

typedef struct OrbitStruct {
  char     Name[MAX_CFG_INLINE+1]; /* system name                       */
  double   Inclination;            /* orbit inclination                 */
  double   Omega;                  /* length of periastron              */
  double   OmegaZero;              /* passage of periastron             */
  double   MQuad;                  /* start (quadrature)                */ 
  double   Excentricity;           /* excentricity of orbit             */
  double   Nu;                     /* true anomaly                      */
  double   M;                      /* mean anomaly                      */
  double   E;                      /* excentric anomaly                 */
  double   MinDist;                /* distance at periastron            */
  double   Phase;                  /* the phase angle                   */
  double   Dist;                   /* the distance                      */
  double   TruePeriod;             /* Period in seconds                 */
  double   TrueMass;               /* Total mass in kg                  */
  double   TrueDistance;           /* Periastron Dist in meter          */
  double   LambdaZero;             /* Wavelength for Profile            */
  int      PhaseIndex;             /* Where do we stand ?               */
  double   Third[NUM_MAG];         /* Third Light                       */
} OrbitType;

/* -----------------  Data Files ------------------------------------   */

typedef struct FileStruct {
  char                DataFile[1024+MAX_CFG_INLINE+1];
  int                 DataFormat;
  /*@null@*/ struct  FileStruct *nextFile;
} FileType;

/*********************************/
/* one-argument macros           */
/*********************************/


/* round a to nearest integer */
#define ROUND(a)      (((a) >= FLT_EPSILON) ? (int)((a)+0.5) : -(int)(0.5-(a)))

/* square a */
extern double GArgSquare;
# define SQR(a) \
( fabs(GArgSquare = (a)) <= FLT_EPSILON ? 0.0 : (GArgSquare*GArgSquare) )

/*********************************/
/* error handling                */
/*********************************/

/* shout if something goes wrong */
#define WARNING(msg) \
{if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON) \
fprintf(stderr, "**Warning**: %s \n", msg);}

/* input error                   */
#define INERR(msg, input) \
{fprintf(stderr,"** Input Error **: %s \n (offending item: %s)\n .. exit to system\n",\
msg,input);\
exit(EXIT_FAILURE);}

/* input error                   */
#define WARNERR(msg, input) \
{fprintf(stderr, "** Input Warning **: %s \n (offending item: %s)\n", msg, input);}


/*********************************/
/* two-argument macros           */
/*********************************/

/* find minimum of a and b */
#undef  MIN
#define MIN(a,b)      ( ( (a)-(b) <= FLT_EPSILON ) ? (a) : (b))

/* find minimum of a and b */
#undef  IMIN
#define IMIN(a,b)      ( ( (a)<(b) ) ? (a) : (b))

/* find maximum of a and b */
#undef  MAX
#define MAX(a,b)      ( ( (a)-(b) >= FLT_EPSILON ) ? (a) : (b))

/* find maximum of a and b */
#undef  IMAX
#define IMAX(a,b)      ( ( (a)>(b) ) ? (a) : (b))

/* Ensures that x is between the limits set by low and high */
#undef CLAMP
#define CLAMP(x, low, high) \
 (((x) > (high)) ? (high) : (((x) < (low)) ? (low)  : (x)))


/*********************************/
/* useful constants              */
/*********************************/

#ifndef PI                      
#define PI        3.14159265358979323846
#endif
#ifndef LN10
#define LN10      2.30258509299404568402      /* ln 10; 10eX = exp(ln10*X)  */
#endif
#ifndef DTOR
#define DTOR      0.01745329251994329547      /* radian = degree * DTOR     */
#endif
#ifndef RTOD
#define RTOD      57.2957795130823228646      /* degree = radian * RTOD     */
#endif

#define ON  1
#define OFF 0

/* Surface identifier for the disk                                          */

#define TOP_SEGMENT    1
#define OUT_RECTANGLE  2
#define BOTTOM_SEGMENT 3  
#define IN_RECTANGLE   4

#define EclipsedByPrimary   1
#define EclipsedBySecondary 2
#define EclipsedByDisk      3

/* illuminated/non-illuminated part of stellar surface                      */

#define BACKSIDE        1
#define FRONTSIDE     (-1)


/*****************************************************/
/*                                                   */
/* PROTOTYPES                                        */
/*                                                   */
/*****************************************************/


/* ---------------- the gui          --------------- */

#ifdef _WITH_GTK
/* main GTK program                                  */
int the_gui(int argc, char *argv[]);
#endif

/* -------  temperature distribution --------------- */

/* compute detailed reflection                       */
void  LightReflect();

/* simple reflection treatment                       */
void SimpleReflect(int Comp);

/* Compute Gravitation Darkening, Albedo             */
void ComputeGravDark();

/* put spots on surface                              */
void MakeSpots(int Component, int phasestep);

/* ------- flux computation  ----------------------- */

/* Get a set of monochromatic wavelenghts from the   */
/* environment                                       */
int DefineMonoWavelengths ();

/* compute temperature-dependent  effective          */
/* wavelengths for blackbody                         */
void EffectiveWavelength(int Comp);

/* Compute limb darkening coeff - filter dependedent */
void LightLimbDark(int Comp);

/* flux in blackbody approximation                   */
void BlackbodyFlux(int Comp);

/* flux with atmospheric model                       */
void ModelFlux(int Comp);

/* limits for Kurucz models                          */
float ModelLimits(float log_g);

/* flux computing - filter dependent                 */
void LightFlux(int Comp, int Phase_index);

/* normalize the output flux                         */
void NormalizeFlux();

/* Initialize Output flux tables - filter dependent  */
void InitOutFlux();

/* normalize the output flux                         */
void SetNormalPhase (int band, float phase);

/* reset normalization point                         */
void ResetNormalPhaseOne (int band);

/* reset all normalization points                    */
void ResetNormalPhase ();

/* copy flux to L90 array                            */
void LightLuminosity(int component, int phaseindex);

/* ------ simplex fit algorithm -------------------  */

/* the downhill simplex code                         */
int Simplex();

/* morph Simplex                                     */
double SimplexFlow(double Factor, double M[][SDIM], 
  double Y[], double M_Sum[], int Worst, 
  int FitDim, int *ErrCode);

/* Initialize the ranges                             */
void SimplexInitRanges();

/* Initialize the Simplex                            */
int SimplexInit(double Y[], double M[][SDIM]);

/* set the global variables                          */
void SimplexSetVariables(double x[]);

/* check the Vertex                                  */
int SimplexCheckVertex(double X[]);

/* print fit result the global variables             */
void SimplexPrintVariables(double Chi, 
  double VarChiSquare);

/* simulated annealing                               */ 
float SaAnneal(double tscale, double t0,  
  double X[], double *Y_Best, 
  double BestEver[], double *BestChi);


/* Map Chi Square                                    */
int ChiMap();

/* ------   main ----------------------------------- */

/* initialize flags and variables                    */
void InitFlags();

/* Print Usage Info                                  */
void UsageInfo();

/*  Get The Command Line Arguments                   */ 
void GetArguments(int *argc, char *argv[]);

/* Here we set up the Stars, do the Lightcurve,      */
/* and return the Merit                              */
int MainLoop(/*@out@*/ double *Merit);

/* radial velocity curve                             */
void RadVel(int Phase_Index);

/* ------   statistics ----------------------------- */

/* compute the merit function - we use ChiSquare     */
/* globally over all input data                      */
double MeritFunction ();

/* perform runs test on fit residuals, return number */
/* of runs,  upper and lower 5 percent limit         */
void Runs(Photo3 *data, long NumData, int *Runs, 
int *UpLim, int *LowLim);

/* -------- Input/Output --------------------------- */

/* wrapper function to do printing                   */
void my_cpgend();

/* Get data file path                                */
char * data_data_fls();

/* Get cfg file path                                 */
char * data_cfg_fls();

/* Get doc file path                                 */
char * data_doc_fls();

/* Get pixmap file path                              */
char * data_pix_fls();

/* parse the config file                             */
void ParseConfig(const char *InputFile, int *numarg);

/* write header of output file                       */
void OutfileHeader(FILE *file_out);

/* read a single line from an open file              */
int LireLigne(char *s, int lim, FILE *fpe);

/* read from input data file                         */
void Read(const char InputFile[]);

/* visualize geometry                                */
void PlotGeometry();

/* Write Output flux tables - filter dependent       */
void WriteOutput();

/* Plot Output flux tables - filter dependent        */
void PlotOutput();

/* animated plot                                     */
void Animate(int j);

/* set PGPLOT window size                            */
void SetPgplotWindow();

/* ---------  Eclipse verification, Visibility ----- */

/* attach the throat of star a to star b             */
void LightCopyThroat();

/* copy the throat back                              */
void LightCopyThroatBack();

/* Calculate various quantities for simple           */
/* geometric tests                                   */
int LightSetupTests();

/* Eclipse Verification for stars                    */
void LightCurve(int Phase_index);

/* Eclipse Verification for disk                     */
void LightCurveDisk(int j);
void eclipsedbydisk(double lzero2, double mzero2, double nzero2, long Diskelements, int Comp, SurfaceElement *SurfPtrD, int *eclipsed);

/* Fractional Visibility                             */
void LightFractionalPot(int Component);

/* --------  Functions ----------------------------- */

/* find minimum of potential along LOS               */
int MinFinder(double *q, double *f, 
 double *t1, double *t3, 
 double *l0, double *m0, double *n0,
 double *x0, double *y0, double *z0,
 double *tmin, /*@out@*/ double *Cmin);

/* find root of Function in [low,high] w/ tolerance  */
double RootFind(double (*Function)(double), 
 const double Low, const double Up, 
 const double Tolerance, const char *caller, int *Err);

/* Zero at LagrangeOne                               */  
double LagrangeOne(double x);

/* Zero at LagrangeTwo                               */  
double LagrangeTwo(double x);

/* Zero at LagrangeThree                             */  
double LagrangeThree(double x);

/* Potential  in X-Y Plane                           */
double RocheXYPrimary(double x, double y, double q, 
 double f);

/* Potential                                         */
double RocheSurface(double x);

/* Potential                                         */
double RocheYPerpendL1(double y);

/* Potential                                         */
double RocheZPerpendL1(double z);

/* Potential                                         */
double RochePerpendSecond(double z);

/* find the phase of periastron                      */
void FindPeriastron();

/* solve kepler equation                             */
void SolveKepler(int j);

/* shift phase to correct interval                   */
void PhaseShift();

/* analytic polar radius                             */
double PolRad(double q, double r);

/* returns (volume-scaledvolume) as function of r    */
/* zero if r = RadMean(ScaledVolume)                 */
double VolOfR(double r);

/* copies PhotoPoint b to a                          */
void V2CopyPhoto(/*@out@*/ Photo3 *a, Photo3 *b);

/* Sort in ascending order in phi w/                 */
/* Heapsort Algorithm                                */
/* D.E. Knuth, Sorting and Searching, 1973           */ 
/*  contrary to Quicksort, this is guaranteed to be  */
/*  always of order NlogN, even in worst case        */
void SortPhotoPoints(long n, Photo3 *seq);

/* ---------  Surface Definition ------------------- */

/* Define Dimensionless System Parameters            */
/* -- Overcontact                                    */ 
int DefineParamOver();

/* Define Dimensionless System Parameters            */ 
int DefineParam(int n);

/* Define Dimensionless System Parameters for disk    */ 
int DefineParamDisk();

/* Surface Division  - filter dependent              */
int DivideSurface(int Comp);

/* Disk Surface Division  - filter dependent         */
/* MK 14.03.01 */
int DivideDisk(int Comp);

/* update radius and potential                       */
int UpdateGeometry(int n);

/* ---------  Surface Map        ------------------- */

/* set base path                                     */
int SetMapPath ();

/* set passband                                      */
void SetMapBand ();

/* print map                                         */
void PrintSurfaceMap (int thisPhaseStep);

/* ---------  Utilities          ------------------- */

/* error handler                                     */
void /*@exits@*/ nf_error(char error_text[]);

/* allocate a matrix                                 */
/*@only@*/ /*@null@*/ /*@out@*/ float **matrix(long nrh, long nch);

/* convert to lowercase                              */
void nf_strlwr(char * instring);

/* clear data file list                              */
void ClearList();

/* add to data file list                             */
void AddToList(const char * item, int Format);

/* determine the locale                              */
const char *guess_the_locale ();

/* Set user-defined albedo value(s)                  */
void  SetAlbedo (int Comp, const char * str);


/**************************************************************************/
/*                                                                        */
/* GLOBAL VARIABLES                                                       */
/*                                                                        */
/**************************************************************************/

#if 0
extern char * program_invocation_name;
#endif

/* -----------   Error Messages ----------------------------------------  */

extern char *errmsg[23];

/* -----------   Flags -------------------------------------------------  */

extern  FlagsHandle Flags;

/* -----------   Flux  -------------------------------------------------  */

extern  /*@out@*/    PhotoOut     *FluxOut;
extern  double       L90[3][NUM_MAG];           /* Luminosities           */
extern  double       F90[NUM_MAG];              /* Normalization Variable */
extern  double       P90[NUM_MAG];              /* Normalization Phase    */
extern  int          Umag, Bmag, Vmag;
extern  int          umag, bmag, vmag, ymag;
extern  int          Rmag, Imag;
extern  int          Jmag, Hmag, Kmag;
extern  int          mag1, mag2;
extern  char         Filters[NUM_MAG+2][24];    /* the names              */

/* -----------limb darkening coefficients ------------------------------  */

extern  double Limb[3][NUM_MAG][2];  /* component, passband, number       */
extern  double WaveL[3][NUM_MAG];    /* wavelengths                       */
extern  float mono_zero[NUM_MAG];    /* monochromatic wavelengths         */

/* ----------- Input File(s) -------------------------------------------- */

extern  Photo3       *Observed[NUM_MAG+2];
extern  /*@null@*/ FileType     *FileList;

/* ----------- Components ----------------------------------------------- */

extern  const int Primary;                  /* literal definition         */
extern  const int Secondary;                /* literal definition         */
extern  const int Disk;                     /* literal definition         */
#ifdef HAVE_DISK   
#define NUM_COMP 3
#else
#define NUM_COMP 2
#endif
extern  BinaryComponent Binary[NUM_COMP];
extern  SurfaceElement  *Surface[NUM_COMP];        /* stars+disk          */


extern  SpotLight  *Spot[2];                /* surface spots              */
extern  OrbitType  Orbit;                   /* orbit details              */

/* ---------- line profile ---------------------------------------------  */

extern  /*@out@*/ float  /* **ProData; */ ProData[PHASESTEPS][PROFILE_ARRAY+1];

/* ---------- misc  ----------------------------------------------------  */

extern  char   Out_Plot_File[256];          /* output plot file           */
extern  int    StepsPerPi;                  /* No of steps per Pi         */
extern  double MaximumElements;             /* max. no of surface elements*/ 
extern  int    PhaseSteps;                  /* phasesteps for lightcurve  */

/* --- These are Global Variables that are used (mainly) to pass -------  */
/* --- Arguments to Functions that are called by RootFind/MinFind ------  */

extern  double   Mq;                 /* mass ratio                        */
extern  double   F;                  /* nonsync rotation                  */
extern  double   RochePot;           /* Roche Potential                   */
extern  double   PhaseScaledVolume;  /* scaled volume                     */
extern  double   lambda, mu, nu;     /* Coordinates                       */


/**************************************************************************/
/*                                                                        */
/* INCLUDES                                                               */
/*                                                                        */
/**************************************************************************/

#ifndef HAVE_MKFIFO
#ifdef _WITH_GNUPLOT
#undef _WITH_GNUPLOT
#undef _WITH_PGPLOT
#endif
#endif

#ifndef HAVE_UNISTD_H

#ifdef _WITH_GNUPLOT
#undef _WITH_GNUPLOT
#undef _WITH_PGPLOT
#endif

#endif

  /* include limits for parameters       */

#include "LightLimits.h"


/* ---------- MPI   ----------------------------------------------------  */

  /* include header for MPI */
#if defined(_WITH_MPI_COARSE)
void NF_Mpi_distribute ();
void NF_Mpi_collect ();
#endif

#if defined(_WITH_MPI) || defined(_WITH_MPI_COARSE)

#include  "mpi.h"
#define NF_MPI_DEBUG
#define NF_MSG_TAG      42  

#else

#define MPI_MAX_PROCESSOR_NAME 256
#define MPI_MAX_ERROR_STRING 512
#define MPI_COMM_WORLD 0
/* dummy functions */
#define MPI_Abort(a, b)
#define MPI_Finalize()

#endif

#define NF_TAG_END        1
#define NF_TAG_MAINLOOP   4
extern int  myrank;
extern int  numprocs;
extern char processor_name[MPI_MAX_PROCESSOR_NAME];
extern char error_name[MPI_MAX_ERROR_STRING];
void NF_Mpi_call(int tag);
void NF_Mpi_endcall();
void NF_MPI_Abort();

/* ---------- MPI END  -------------------------------------------------  */

  /* include header for PGPLOT emulation */

#ifdef _WITH_GNUPLOT
#include  "LightPGemulate.h"
#endif

  /* include header for GUI, also included by gnome.h */

#ifdef _WITH_GTK

#ifdef HAVE_GNOME
#include <gnome.h>
#else
#include <gtk/gtk.h>
#endif

#include "LightGtk.h"

#ifdef _WITH_OPENGL
/* need to #undef this because /usr/include/jconfig.h
 * defines HAVE_STDDEF_H (which is plain stupid IMHO)
 */
#undef HAVE_STDDEF_H
#include "LightGL.h"
#endif 
#endif

/* end include header for GUI */

/* i18n, this stuff is already included by gnome.h  */

#ifndef HAVE_GNOME

#include "gettext.h"

#ifdef ENABLE_NLS

/* #include <libintl.h> */
#define _(string)  gettext (string)

#ifdef gettext_noop
#define N_(string) gettext_noop (string)
#else
#define N_(String) (String)
#endif

#else

#define _(string)  string
#define N_(string) string

#endif

/* #ifndef HAVE_GNOME */
#endif

/* end i18n */












syntax highlighted by Code2HTML, v. 0.9.1