#! python
""" Code to compute astrometric solutions for REM frames

Context : SRP
Module  : SRPREMAstrometrySearch
Author  : Stefano Covino
Date    : 18/05/2017
E-mail  : stefano.covino@brera.inaf.it
URL:    : http://www.merate.mi.astro.it/utenti/covino
Purpose : Compute REM frame astrometry.

Usage   : SRPREMAstrometry [-h] [-b box] -i file [-m rms] [-n nsrc ncat] -o file [-s step]
            [-v] [--version]    
            -b Box size (deg)
            -i Input FITS file
            -m Max rms (arcsec) for an acceptable solution
            -n Number of objects to analyze (source catalog)
            -o Output FITS file
            -s Steps in box scanning
    
History : (25/08/2013) First version.
        : (29/08/2013) More process information.
        : (04/09/2013) Minor improvements.
        : (09/01/2014) Max residual.
        : (05/08/2014) Center pixels.
        : (23/07/2015) Bug correction.
        : (07/10/2015) Better check of input parameters.
        : (10/05/2016) Better porting.
        : (04/03/2017) Better pipes.
        : (18/05/2017) Minor update.
"""

__version__ = '1.1.4ì5'


import argparse, os, shutil
import numpy
from SRPFITS.Fits.GetHeaderValue import GetHeaderValue
import SRPFITS.Fits.FitsConstants as SSF
from SRP.SRPSystem.Pipe import Pipe
from SRP.SRPSystem.Which import Which
from SRPFITS.Fits.IsFits import IsFits
import SRPREM
from SRP.SRPMath.AngularDistance import AngularDistance


cmdlong = "SRPAstrometry -i %s -o %s -x %f %f -n %d %d -m %1.f -p %.2f %.2f"
cmdslong = "SRPAstrometry -i %s -o %s -x %f %f -n %d %d -c %.5f %.5f -P %.5f %.5f -m %.1f -p %.2f %.2f"


parser = argparse.ArgumentParser()
parser.add_argument("-b", "--boxsize", action="store", type=float, default=0.5, help="Box size (deg)", metavar="box")
parser.add_argument("-i", "--inpfile", action="store", help="Input FITS file", metavar="file", required=True)
parser.add_argument("-m", "--maxrms", action="store", type=float, default=3.0, help="Max rms (arcsec) for an acceptable solution", metavar="rms")
parser.add_argument("-n", "--nobjs", action="store", type=int, nargs=2, default=(10,10), help="Number of objects to analyze (source catalog)", metavar=("nsrc", "ncat"))
parser.add_argument("-o", "--outfile", action="store", help="Output FITS file", metavar="file", required= True)
parser.add_argument("-s", "--step", action="store", type=int, default=6, help="Steps in box scanning", metavar="step")
parser.add_argument("-v", "--verbose", action="store_true", help="Fully describe operations")
parser.add_argument("--version", action="version", version=__version__)
options = parser.parse_args()


# Component safety check
if Which(SRPREM.astro) == None:
    parser.error("%s not found." % SRPREM.astro)

if not IsFits(options.inpfile):
    parser.error("Input FITS file %s does not exist." % options.inpfile)

if options.nobjs[0] < 5 or options.nobjs[1] < 5:
    parser.error("Number of considered objects must be at least 5.")

if options.boxsize < 0.15:
    parser.error("Box size must be larger than 0.15")

if options.step < 1:
    parser.error("Number of steps must be positive.")

if options.maxrms <= 0.0:
    parser.error("Max acceptable residual must be positive.")

if options.verbose:
    print("Input FITS file : %s" % options.inpfile)
    print("Output FITS file: %s" % options.outfile)
    print("Number of source/catalogue objects to consider: %d/%d" % (options.nobjs[0], options.nobjs[1]))
    print("Box size: %.2f" % options.boxsize)
    print("Number of steps: %d" % options.step)
    print("Max rms: %.1f" % options.maxrms)

# Read header
cam = GetHeaderValue (options.inpfile, SRPREM.CAMERA)
RA = GetHeaderValue(options.inpfile, SRPREM.RA)[0]
DEC = GetHeaderValue(options.inpfile, SRPREM.DEC)[0]
if cam == None or RA == None or DEC == None:
    parser.error("Problem in reading input file %s" % options.inpfile)

#
ROSSflg = False
ROS2flg = False
REMIRflg = False
#
if cam[1] == SSF.FitsHeaderFound:
    # ROSS
    if cam[0] == SRPREM.ROSS:
        ROSSflg = True
    elif cam[0] == SRPREM.ROS2:
        ROS2flg = True
    # REMIR
    elif cam[0] == SRPREM.REMIR:
        REMIRflg == True
    else:
        parser.error("REM camera not recognized.")
elif cam[1] == SSF.FitsFileNotFound:
    parser.error("Input FITS file not found.")
elif cam[1] == SSF.FitsHeaderNotFound:
    parser.error("Input FITS file does not seem to be a REM frame.")
else:
    parser.error("Unexpcted error.")

# filename processing
fileout = os.path.join(os.path.dirname(options.outfile),'__'+os.path.basename(options.outfile))

if options.verbose:
    print("Trying at pointing coordinates...")
Astrflg = False
if ROSSflg or ROS2flg:
    cmd = cmdlong % (options.inpfile, fileout, -SRPREM.ROSSPixSize, SRPREM.ROSSPixSize, options.nobjs[0], options.nobjs[1], options.maxrms, SRPREM.ROSS2XCenter, SRPREM.ROSS2YCenter)
else:
    cmd = cmdlong % (options.inpfile, fileout, -SRPREM.REMIRPixSize, SRPREM.REMIRPixSize, options.nobjs[0], options.nobjs[1], options.maxrms, SRPREM.REMIRXCenter, SRPREM.REMIRYCenter)
res = Pipe(cmd)
if res != None:
    if res.decode().split()[0] == '1':
        Astrflg = True
    else:
        ra = numpy.linspace(RA-options.boxsize/(2.*numpy.cos(numpy.radians(DEC))),RA+options.boxsize/(2.*numpy.cos(numpy.radians(DEC))),options.step)
        dec = numpy.linspace(DEC-options.boxsize/2.,DEC+options.boxsize/2.,options.step)
        steps = []
        for r in ra:
            for d in dec:
                steps.append((r,d))
        #
        steps.sort(key=lambda e: AngularDistance((e[0],e[1]),(RA,DEC)))
        #
        istep = 1
        nstep = len(steps)
        for el in steps:
            if options.verbose:
                print("Trying at RA: %.5f, DEC: %.5f (%d of %d in box)..." % (el[0],el[1], istep, nstep))
            istep = istep + 1
            #
            if ROSSflg or ROS2flg:
                cmd = cmdslong % (options.inpfile, fileout, -SRPREM.ROSSPixSize, SRPREM.ROSSPixSize, options.nobjs[0], options.nobjs[1], el[0], el[1], el[0], el[1], options.maxrms, SRPREM.ROSS2XCenter, SRPREM.ROSS2YCenter)
            else:
                cmd = cmdslong % (options.inpfile, fileout, -SRPREM.REMIRPixSize, SRPREM.REMIRPixSize, options.nobjs[0], options.nobjs[1], el[0], el[1], el[0], el[1], options.maxrms, SRPREM.REMIRXCenter, SRPREM.REMIRYCenter)
            res = Pipe(cmd)
            if res != None:
                if res.decode().split()[0] == '1':
                    Astrflg = True
                    break
                if ROSSflg:
                    cmd = cmdslong % (options.inpfile, fileout, SRPREM.ROSSPixSize, SRPREM.ROSSPixSize, options.nobjs[0], options.nobjs[1], el[0], el[1], el[0], el[1], options.maxrms, SRPREM.ROSS2XCenter, SRPREM.ROSS2YCenter)
                    res = Pipe(cmd)
                    if res != None:
                        if res.decode().split()[0] == '1':
                            Astrflg = True
                            break
                    else:
                        parser.error("Problem with %s." % SRPREM.astro)
            else:
                parser.error("Problem with %s." % SRPREM.astro)
else:
    parser.error("Problem with %s." % SRPREM.astro)
#
#
if Astrflg:
    if options.verbose:
        print("Astrometry for frame %s computed." % options.inpfile)
        print("Astrometrized file %s generated." % options.outfile)
        print("Average residual: %s arcsec for %s stars." % (res.decode().split()[1], res.decode().split()[2])) 
        print("RA [x cos(DEC)] and DEC shift wrt the pointing coordinates: %s %s arcsec" % (res.decode().split()[3], res.decode().split()[4]))
    else:
        print(res.decode())
    #
else:
    if options.verbose:
        print("Astrometry for frame %s cannot be computed." % options.inpfile)
        print("Average residual: %s arcsec for %s stars." % (res.decode().split()[1], res.decode().split()[2]))
    else:
        print(res.decode())
#
shutil.move(fileout,options.outfile)
#
