#!/usr/bin/python
# Load the WCS information from a fits header, and use it
# to convert pixel coordinates to world coordinates.

from __future__ import division, print_function

import numpy
from astropy import wcs
from astropy.io import fits
import sys, getopt
from astropy import units as u
from astropy.coordinates import SkyCoord

#def coors(argv):
#   ra = ''
#   dec = ''
#   try:
#      opts, args = getopt.getopt(argv,"r:d:",["right_ascension=","declination="])
#   except getopt.GetoptError:
#      print('test.py -r <right_ascension> -d <declination>')
#      sys.exit(2)
#for opt, arg in opts:
#    if opt == '-h':
#        print('test.py ima.fits -r <right_ascension> -d <declination>')
#        sys.exit()
#    elif opt in ("-r", "--right_ascension"):
#        ra = arg
#    elif opt in ("-d", "--declination"):
#        dec = arg
#        print('right ascension is "', ra)
#        print('declination is "', dec)
#if __name__ == "__coors__":
#    coors(sys.argv[1:])


def load_wcs_from_file(filename):
    # Load the FITS hdulist using astropy.io.fits
    hdulist = fits.open(sys.argv[1])

    ra = sys.argv[2]
    dec = sys.argv[3]

    # Parse the WCS keywords in the primary HDU
    w = wcs.WCS(hdulist[0].header)

    # Print out the "name" of the WCS, as defined in the FITS header
#    print(w.wcs.name)

    # Print out all of the settings that were parsed from the header
#    w.wcs.print_contents()

    # Some pixel coordinates of interest.
#    pixcrd = numpy.array([[ra, dec], [500.73034, 372.74316], [459.14685, 248.90399]], numpy.float_)
    pixcrd = numpy.array([[ra, dec]], numpy.float_)

    # Convert pixel coordinates to world coordinates
    # The second argument is "origin" -- in this case we're declaring we
    # have 1-based (Fortran-like) coordinates.
    world = w.wcs_pix2world(pixcrd, 1)
    #world =  numpy.array([[ra, dec]], numpy.float_)
    print(world)
    
    #w1 = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree, frame='fk5')
    w1 = SkyCoord(world, frame='fk5', unit='deg')
    w2 = w1.to_string('hmsdms')
    print(w2) 
    #w1=w1.fk5
    #print(w1) 

    # Convert the same coordinates back to pixel coordinates.
    #pixcrd2 = w.wcs_world2pix(world, 1)
    #print(pixcrd2)

    # These should be the same as the original pixel coordinates, modulo
    # some floating-point error.
 #   assert numpy.max(numpy.abs(pixcrd - pixcrd2)) < 1e-6


if __name__ == '__main__':
    load_wcs_from_file(sys.argv[-1])