# if !defined(PS_ASTROM_H)
# define PS_ASTROM_H

/** \file psAstrom.h
 *  \brief Astrometric transformations
 *  \ingroup AstroGroup
 */

/* Astrometric transformations for a focal plane array.
 * There are five coordinate frames that we need to worry about:
 * 1. Cell:          x,y in pixels (raw coordinates; lower-left pixel center is 0.5,0.5 TBR)
 * 2. Chip:          X,Y in pixels
 * 3. FPA:           P,Q in microns 
 * 4. Tangent plane: L,M in arcsec (TBR - angles or linear?) from the telescope boresight.
 * 5. Sky frame:     R,D in degrees (TBR - radians?)
 *
 * 1 <--> 2: 2 2-D polynomials: (X,Y) = f(x,y)
 * 2 <--> 3: 2 2-D polynomials: (P.Q) = f(X,Y)
 * 3 <--> 4: 2 4-D polynomials: (L,M) = f(P,Q,m,c) (m & c are magnitude and color of object)
 * 4 <--> 5: SLALib transformation using exposure parameters (psExposure or psGrommit)
 */

/** Observatory data */
typedef struct {
    const char *name;			///< Name of observatory
    const double latitude;		///< Latitude of observatory, east positive
    const double longitude;		///< Longitude of observatory
    const double height;		///< Height of observatory
    const double tlr;			///< Tropospheric Lapse Rate
} psObservatory;

psObservatory *
psObservatoryAlloc(const char *name,	///< Name of observatory
		   double latitude,	///< Latitude of observatory, east positive
		   double longitude,	///< Longitude of observatory
		   double height,	///< Height of observatory
		   double tlr		///< Tropospheric Lapse Rate
		   );

/** Information needed (by SLALIB) to convert Apparent to Observed Position */
typedef struct {
    const double latitude;		///< geodetic latitude (radians)
    const double sinLat, cosLat;	///< sine and cosine of geodetic latitude
    const double abberationMag;		///< magnitude of diurnal aberration vector
    const double height;		///< height (HM)
    const double temperature;		///< ambient temperature (TDK)
    const double pressure;		///< pressure (PMB)
    const double humidity;		///< relative humidity (RH)
    const double wavelength;		///< wavelength (WL)
    const double lapseRate;		///< lapse rate (TLR)
    const double refractA, refractB;	///< refraction constants A and B (radians)
    const double longitudeOffset;	///< longitude + eqn of equinoxes + ``sidereal UT'' (radians)
    const double siderealTime;		///< local apparent sidereal time (radians)
} psGrommit;

/** Constructor */
psGrommit *
psGrommitAlloc(const psExposure *exp	///< Exposure details from which to derive grommit
    );

/** Exposure information from the telescope */
typedef struct {
    const double ra, dec;		///< Telescope boresight
    const double ha;			///< Hour angle
    const double zd;			///< Zenith distance
    const double az;			///< Azimuth
    const psTime *time;			///< Time of observation
    const float rotAngle;		///< Rotator position angle
    const float temp;			///< Air temperature, for estimating refraction
    const float pressure;		///< Air pressure, for calculating refraction
    const float humidity;		///< Relative humidity, for calculating refraction
    const float exptime;		///< Exposure time
    const float wavelength;		///< Wavelength of observation
    const psObservatory *observatory;	///< Observatory data: longitude, latitude, height etc.
    /* Derived quantities */
    const psTime *lst;			///< Local Sidereal Time
    const float posAngle;		///< Position angle
    const float parallactic;		///< Parallactic angle
    const float airmass;		///< Airmass, calculated from zenith distance
    const float pf;			///< Parallactic factor
    const char *cameraName;		///< name of camera which provided exposure
    const char *telescopeName;		///< name of telescope which provided exposure
} psExposure;

/** Constructor */
psExposure *
psExposureAlloc(double ra, double dec,  ///< Telescope boresight
                double ha,              ///< Hour angle
                double zd,              ///< Zenith distance
                double az,              ///< Azimuth
		const psTime *time,	///< Time of observation
                float rotAngle,         ///< Rotator position angle
                float temp,             ///< Temperature
                float pressure,         ///< Pressure
                float humidity,         ///< Relative humidity
                float exptime,		///< Exposure time
		float wavelength,	///< Wavelength of observation
		const psObservatory *observatory ///< Observatory data
		);

/** The fixed pattern residual offsets.  These are specified via a coarse grid of x and y offsets. */
typedef struct {
    int nX, nY;				///< Number of elements in x and y
    double x0, y0;			///< Position of the lower-left corner of the grid on the focal plane
    double xScale, yScale;		///< Scale of the grid
    double **x, **y;			///< The grid of offsets in x and y
} psFixedPattern;

/** Constructor */
psFixedPattern *
psFixedPatternAlloc(double x0, double y0, ///< Position of the lower-left corner of the grid on the focal plane
		    double xScale, double yScale, ///< Scale of the grid
		    const psImage *x,	///< Grid of offsets in x
		    const psImage *y	///< Grid of offsets in y
		    );


/** a Focal plane array: a collection of chips.  Not all chips in a camera need to be listed in an instance of
 *  psFPA.
 */
typedef struct {
    int nChips;                         ///< Number of Cells assigned
    int nAlloc;                         ///< Number of Cells available
    struct psChip *chips;		///< Chips in the Focal Plane Array

    psMetadata *md;			///< FPA-level metadata 
    psPlaneDistort *TPtoFP;		///< Transformation term from 
    psPlaneDistort *FPtoTP;		///< Transformation term from 
    psFixedPattern *pattern;		///< Fixed pattern residual offsets
    const psExposure *exp;		///< information about this exposure
    const psGrommit *grommit;		///< Information for conversion from TP to Sky.
    psPhotSystem colorPlus, colorMinus; ///< Colour reference
    float rmsX, rmsY;                   ///< Dispersion in astrometric solution
    float chi2;                         ///< chi^2 of astrometric solution
} psFPA;

/** Constructor: makes an empty FPA */
psFPA *
psFPAAlloc(int nAlloc,			///< Number of Cells to allocate
	   const psExposure *exp	///< information about this exposure
	   );


/** a Chip: a collection of cells.  Not all valid cells in a chip need to be listed in an
 *  instance of psChip.
 */
typedef struct {
    int nCells;                         ///< Number of Cells assigned
    struct psCell *cells;		///< Cells in the Chip

    psMetadata *md;			///< Chip-level metadata
    psPlaneTransform *chipToFPA;	///< Transformations from chip coordinates to FPA coordinates
    psPlaneTransform *FPAtoChip;	///< Transformations from FPA coordinates to chip

    struct psFPA *parentFPA;		///< FPA which contains this chip
} psChip;

/** Constructor */
psChip *
psChipAlloc(int nCells,			///< Number of Cells assigned
	    psFPA *parentFPA		///< FPA which contains this chip
	    );


/** a Cell: a collection of readouts.
 */
typedef struct {
    int nReadouts;			///< number of readouts in this cell realization; each may have its
					///< own image, objects and overscan.
    struct psReadout *readouts;		///< Readouts from the cell
    psMetadata *md;			///< Cell-level metadata
    psPlaneTransform *cellToChip;	///< Transformations from cell coordinates to chip coordinates
    psPlaneTransform *cellToFPA;	///< Transformations from cell coordinates to FPA coordinates
    psPlaneTransform *toTP;		///< Quick and Dirty transformations from cell coordinates to sky
    struct psChip *parentChip;		///< chip which contains this cell
} psCell;

/** Constructor */
psCell *
psCellAlloc(int nReadouts,		///< number of readouts in this cell
	    struct psChip *parentChip	///< chip which contains this cell
	    );


/** a Readout: a collection of pixels */
typedef struct {
    const int x0, y0;			///< Offset from the lower-left corner of the physical device
    const int xBin, yBin;		///< Image binning
    psImage *image;                     ///< imaging area of cell 
    psDlist *objects;			///< objects derived from cell
    psMetadata *md;			///< Readout-level metadata
} psReadout;
    

/** Constructor */
psReadout *
psReadoutAlloc(int x0, int y0,		///< Offset from the lower-left corner of the physical device
	       const psImage *image	///< imaging area of cell
	       );


/** Functions **************************************************************/
/** \addtogroup AstroGroup Astronomy-Specific Utilities
 *  \{
 */

/*** Calculating and applying astrometric solutions ***/

/** Return the cell in FPA which contains the given FPA coordinates */
psCell *
psCellInFPA(const psPlane *coord,	///< Coordinate in FPA
	    const psFPA *fpa)		///< FPA description
;


/** returns Chip in FPA which contains the given FPA coordinate */
psChip *
psChipInFPA(const psPlane *coord,	///< coordinate in FPA
	    const psFPA *fpa) 		///< FPA description
;

/** returns Cell in Chip which contains the given chip coordinate */
psCell *
psCellInChip(const psPlane *coord,	///< coordinate in chip
	     const psChip *chip)	///< chip description
;

/************ conversions low -> high ****************/
/** converts the specified Cell coord to the coord on the parent Chip */
psPlane *
psCoordCelltoChip (psPlane *out,	///< Coordinates to return, or NULL
		   const psPlane *in,	///< input Cell coordinate
		   const psCell *cell) 	///< Cell description
;

/** converts the specified Chip coord to the coord on the parent FPA */
psPlane *
psCoordChiptoFPA (psPlane *out,		///< Coordinates to return, or NULL
		  const psPlane *in,	///< input Chip coordinate
		  const psChip *chip) 	///< Chip description
;

/** Convert focal plane coords to tangent plane coordinates */
psPlane *
psCoordFPAToTP(psPlane *out,		///< Coordinates to return, or NULL
	       const psPlane *in,	///< input FPA coordinate
	       float color,		///< Color of source, using reference defined in the fpa
	       float mag,		///< Magnitude of source
	       const psFPA *fpa)	///< FPA description
;

/** Convert tangent plane coords to (RA,Dec) */
psSphere *
psCoordTPtoSky(psSphere *out,		///< Coordinates to return, or NULL
	       const psPlane *in,	///< input TP coordinate
	       const psGrommit *grommit) ///< Grommit for fast conversion
;

/** Convert Cell coords to FPA coordinates */
psPlane *
psCoordCellToFPA(psPlane *out,		///< Coordinates to return, or NULL
		 const psPlane *in,	///< Input cell coordinates
		 const psCell *cell)	///< Cell description
;

/** Convert cell and cell coordinate to (RA,Dec) */
psSphere *
psCoordCellToSky(psSphere *out,		///< Coordinates to return, or NULL
		 const psPlane *in,	///< cell coordinates to transform
		 float color,		///< Color of source, using reference defined in the fpa
		 float mag,		///< Magnitude of source
		 const psCell *cell)	///< Cell to get coordinates for
;

/** Convert cell and cell coordinate to (RA,Dec).  Don't need color,mag here because it's quick and dirty. */
psSphere *
psCoordCellToSkyQuick(psSphere *out,	///< Coordinates to return, or NULL
		      const psPlane *in, ///< cell coordinates to transform
		      const psCell *cell) ///< Cell to get coordinates for
;


/************ conversions low -> high ****************/
/** Convert (RA,Dec) to tangent plane coords */
psPlane *
psCoordSkyToTP(psPlane *out,		///< Coordinates to return, or NULL
	       const psSphere *in,	///< input Sky coordinate
	       const psGrommit *grommit) ///< Grommit for fast conversion
;

/** Convert tangent plane coords to focal plane coordinates */
psPlane *
psCoordTPtoFPA(psPlane *out,		///< Coordinates to return, or NULL
	       const psPlane *in,	///< input TP coordinate
	       float color,		///< Color of source, using reference defined in the fpa
	       float mag,		///< Magnitude of source
	       const psFPA *fpa)	///< FPA description
;

/** converts the specified FPA coord to the coord on the given Chip */
psPlane *
psCoordFPAtoChip (psPlane *out,		///< Coordinates to return, or NULL
		  const psPlane *in,	///< input FPA coordinate
		  const psChip *chip)	///< Chip of interest
;

/** converts the specified Chip coord to the coord on the given Cell */
psPlane *
psCoordChiptoCell (psPlane *out,	///< Coordinates to return, or NULL
		   const psPlane *in,	///< input Chip coordinate
		   const psCell *cell) 	///< Cell of interest
;

/** Convert (RA,Dec) to cell and cell coordinates */
psPlane *
psCoordSkyToCell(psPlane *out,		///< Coordinates to return
		 const psSphere *in,	///< Input coordinates
		 float color,		///< Color of source, using reference defined in the fpa
		 float mag,		///< Magnitude of source
		 const psCell *cell)	///< Cell of interest
;

/** Quick and dirty cell to (RA,Dec) --- employs cellToSky transformation */
psSphere *
psCoordSkyToCellQD(psPlane *out,	///< Coordinates to return
		   const psSphere *in,	///< Input coordinates
		   const psCell *cell)	///< Cell of interest
;

/**** other utilities ****/

/** Get the airmass for a given position and sidereal time */
float
psGetAirmass(const psSphere *coord,	///< Position on the sky
             psTime *lst,		///< Local sidereal time
	     float height)		///< Height above sea level
;

/** Get the parallactic angle for a given position and sidereal time */
float
psGetParallactic(const psSphere *coord, ///< Position on the sky
                 double siderealTime)	///< Sidereal time
;

/** Estimate atmospheric refraction, along the parallactic */
float
psGetRefraction(float colour,           ///< Colour of object
		psPhotSystem colorPlus,	///< Colour reference
		psPhotSystem colorMinus, ///< Colour reference
                const psExposure *exp)	///< Telescope pointing information, for airmass, temp and pressure
;

/** Calculate the parallax factor */
psSphere *
psGetParallaxFactor(const psExposure *exp) ///< Exposure details
;

/* \} */ // End of AstroGroup Functions

#endif

/*
  example images:

  GPC:          1 Chip per file, 1 Cell per extension 
  Megacam:      1 FPA per file, 1 Chip per extension, 2 Cell per extension (DSEC) 
  CFH12K (MEF): 1 FPA per file, 1 Chip per extension, 1 Cell per extension
  CFHT-IR:      1 FPA per file, 1 Chip per extension, 4 Cell per extension (NAXIS = 3)
  WIRCAM:       1 FPA per file, 1 Chip per extension, 1 Cell per extension (?)
   

*/
