# include <ohana.h>
# include <dvo.h>
# include <kapa.h>
# include <signal.h>

typedef enum {
  MODE_SIMPLE,
  MODE_CHIP,
  MODE_MOSAIC,
} CoordMode;

typedef enum {ERROR_MODE_RA, ERROR_MODE_DEC, ERROR_MODE_POS} ErrorMode;

typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode;

typedef enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget;

typedef struct {
  double R, D;  /* Sky Coords    - degrees */
  double P, Q;  /* Tangent Plane - pixels  */
  double L, M;  /* Focal Plane   - pixels  */
  double X, Y;  /* Chip Coords   - pixels  */
  double Mag;
  double dMag;
  double dPos;
  int mask;
} StarData;

// structure to hold coordinate fitting terms
typedef struct {
    int Npts;
    int Nterms;
    int Norder;
    int Nsums;
    int Nelems;
    double **sum;
    double **xsum;
    double **ysum;
    double **xfit;
    double **yfit;
} CoordFit;

typedef struct {
  double Ro, dRo;
  double Do, dDo;

  double uR, duR;
  double uD, duD;

  double p, dp;

  double chisq;
  int Nfit;
} PMFit;

typedef struct {
  unsigned int start;
  unsigned int stop;
  float Mcal;
  float dMcal;
  short Xm;
  float secz;
  char flags;
  Coords coords;
} Mosaic;

typedef struct {
  double median;
  double mean;
  double sigma;
  double error;
  double chisq;
  double min;
  double max;
  double total;
  int    Nmeas;
} StatType;

/* global variables set in parameter file */
char   ImageCat[256];
char   GSCFILE[256];
char   CATDIR[256];
char   CATMODE[16];    /* raw, mef, split, mysql */
char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
char   SKY_TABLE[256];
int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */

double SIGMA_LIM;
int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
double MIN_ERROR;

int    VERBOSE;
int    VERBOSE2;

int    RESET;
int    NLOOP;
int    UPDATE;
int    PLOTSTUFF;
int    SAVEPLOT;
int    SHOW_PARAMS;
char   STATMODE[32];
int    POS_TOOFEW;
int    PM_TOOFEW;
double PM_DT_MIN;
int    PLOTDELAY;
int    CHIPORDER;

int UserCatalog;
double UserCatalogRA, UserCatalogDEC;

char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
int           NphotcodesKeep,      NphotcodesSkip;
PhotCode     **photcodesKeep,     **photcodesSkip;

int AreaSelect;
double AreaXmin, AreaXmax, AreaYmin, AreaYmax;

int ImagSelect, ImagMin, ImagMax;

double  PlotMmin, PlotMmax, PlotdMmin, PlotdMmax;

int DophotSelect, DophotValue;
int PhotFlagSelect, PhotFlagPoor, PhotFlagBad;

int TimeSelect;
time_t TSTART, TSTOP;

int FlagOutlier;
int    CLIP_THRESH;

FitMode FIT_MODE;

FitTarget FIT_TARGET;

SkyRegion UserPatch;

int DoUpdateObjects;
int DoUpdateSimple;
int DoUpdateChips;
int DoUpdateMosaics;

/*** relphot prototypes ***/
void          ConfigInit          PROTO((int *argc, char **argv));
void          GetConfig           PROTO((char *config, char *field, char *format, int N, void *ptr));
char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
off_t          *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
int           args                PROTO((int argc, char **argv));
int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
void          clean_images        PROTO((void));
void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));
void          clean_mosaics       PROTO((void));
void          clean_stars         PROTO((Catalog *catalog, int Ncatalog));
int           corner_check        PROTO((double *x1, double *y1, double *x2, double *y2));
void          dumpGrid            PROTO((void));
void          dump_grid           PROTO((void));
int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
void          findImages          PROTO((Catalog *catalog, int Ncatalog, int MATCHCAT));
int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum));
void set_db (FITS_DB *in);
int Shutdown (char *format, ...) OHANA_FORMAT(printf, 1, 2);
void TrapSignal (int sig);
void SetProtect (int mode);
int SetSignals (void);

GSCRegion    *find_regions        PROTO((Image *image, off_t Nimage, int *Nregions, GSCRegion *fullregion));
void          freeGridBins        PROTO((int Ncatalog));
void          freeImageBins       PROTO((int Ncatalog));
void          freeMosaicBins      PROTO((int Ncatalog));
void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
int           gcatalog            PROTO((Catalog *catalog, int FINAL));
// Coords       *getCoords           PROTO((off_t meas, int cat));
float         getMcal             PROTO((off_t meas, int cat));
float         getMgrid            PROTO((off_t meas, int cat));
float         getMmos             PROTO((off_t meas, int cat));
float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, off_t *Nregions));
void          getfullregion       PROTO((Image *image, off_t Nimage, GSCRegion *fullregion));
Image        *getimage            PROTO((off_t N));
Image        *getimages           PROTO((off_t *N));
void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
void          initGrid            PROTO((int dX, int dY));
void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int FULLINIT));
void          initImages          PROTO((Image *input, off_t N));
void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
void          initMosaics         PROTO((Image *image, off_t Nimage));
void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
void          initialize          PROTO((int argc, char **argv));
void          initstats           PROTO((char *mode));
int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));

void check_permissions (char *basefile);
void lock_image_db (FITS_DB *db, char *filename);
void unlock_image_db (FITS_DB *db);
void create_image_db (FITS_DB *db);
void save_catalogs (Catalog *catalog, int Ncatalog);

int           main                PROTO((int argc, char **argv));
void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat, int MATCHCAT));
void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat));
GSCRegion    *name_region         PROTO((char *name, off_t *Nregions));
double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
void          plot_defaults       PROTO((Graphdata *graphdata));
void          plot_grid           PROTO((Catalog *catalog));
void          plot_images         PROTO((void));
void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file));
void          plot_mosaic_fields  PROTO((Catalog *catalog));
void          plot_mosaics        PROTO((void));
void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog));
void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
void          reload_catalogs     PROTO((SkyList *skylist));
int           setExclusions       PROTO((Catalog *catalog, int Ncatalog));
void          setMcal             PROTO((Catalog *catalog, int Poor));
void          setMcalFinal        PROTO((void));
int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog));
void          setMgrid            PROTO((Catalog *catalog));
int           setMmos             PROTO((Catalog *catalog, int Poor));
int           setMrel             PROTO((Catalog *catalog, int Ncatalog));
void          setMrelFinal        PROTO((Catalog *catalog));
int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
void          set_ZP              PROTO((double ZERO));
int           setrefcode          PROTO((Image *image, off_t Nimage));
void          skip_measurements   PROTO((Catalog *catalog, int pass));
void          sortA               PROTO((double *X, int N));
void          sortB               PROTO((double *X, double *Y, int N));
void          sortC               PROTO((double *X, double *Y, double *F1, double *F2, int N));
void          sortD               PROTO((double *X, double *Y, double *Z, int N));
StatType      statsImageM         PROTO((Catalog *catalog));
StatType      statsImageN         PROTO((Catalog *catalog));
StatType      statsImageX         PROTO((Catalog *catalog));
StatType      statsImagedM        PROTO((Catalog *catalog));
StatType      statsMosaicM        PROTO((Catalog *catalog));
StatType      statsMosaicN        PROTO((Catalog *catalog));
StatType      statsMosaicX        PROTO((Catalog *catalog));
StatType      statsMosaicdM       PROTO((Catalog *catalog));
StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog));
StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog));
StatType      statsStarX          PROTO((Catalog *catalog, int Ncatalog));
void          wcatalog            PROTO((Catalog *catalog));
void          wimages             PROTO((void));
void          write_coords        PROTO((Header *header, Coords *coords));

double **array_init (int Nx, int Ny);
void array_free (double **array, int Nx);
CoordFit *fit_init (int order);
void fit_free (CoordFit *fit);
void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
int fit_eval (CoordFit *fit);
void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1);
double **poly2d_dx (double **poly, int Nx, int Ny);
double **poly2d_dy (double **poly, int Nx, int Ny);
double **poly2d_copy (double **poly, int Nx, int Ny);
double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
int fit_apply_coords (CoordFit *fit, Coords *coords);
int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo);
CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo);
int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
void initObjectData (Catalog *catalog, int Ncatalog);
int UpdateObjects (Catalog *catalog, int Ncatalog);
int UpdateSimple (Catalog *catalog, int Ncatalog);
int UpdateChips (Catalog *catalog, int Ncatalog);
int UpdateMosaic (Catalog *catalog, int Ncatalog);
int UpdateMeasures (Catalog *catalog, int Ncatalog);
void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im);
void FlagOutliers(Catalog *catalog);
int MeasFilterTest(Measure *measure);

int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean);
int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts);
int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);
int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);

Mosaic *getMosaicForImage (off_t N);

StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);

Mosaic *getmosaics (off_t *N);
void initMosaics (Image *image, off_t Nimage);
StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
Mosaic *getMosaicForImage (off_t im);

double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
double getMeanD (Measure *measure, Average *average, SecFilt *secfilt);
int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);

float GetAstromError (Measure *measure, int mode);
int relastro_objects (void);
int UpdateObjectOffsets (SkyList *skylist);

int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
int relastroVisualPlotScatter(double values[], double thresh, int npts);
int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, 
			       StatType statsR, StatType statsD, double thresh);



int FixProblemImages (SkyList *skylist);
int *getCatlist (int *N, off_t im);

void initCoords (void);
void getOffsets (double *dPos, off_t *nPos, off_t N);
void saveOffsets (double dPos, off_t nPos, off_t N);
void setBadCoords (off_t N);
int badCoords (off_t N);
Coords *getCoords (off_t N);
int saveCoords (Coords *coords, off_t N);
void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im);
