# -*- coding: utf-8 -*-
#
#   pycheops - Tools for the analysis of data from the ESA CHEOPS mission
#
#   Copyright (C) 2018  Dr Pierre Maxted, Keele University
#
#   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 3 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, see <https://www.gnu.org/licenses/>.
#
"""
make_xml_files
==============
 Generate XML files for CHEOPS observing requests

Dictionaries
------------

 SpTypeToGminusV - valid keys are A1 to M9V
 SpTypeToTeff    - valid keys are A1 to M9V

Functions 
---------
 main() - make_xml_files


"""

from __future__ import (absolute_import, division, print_function,
                                unicode_literals)

import argparse
import textwrap
from astropy.table import Table, Row
from astropy.time import Time
from astropy.coordinates import SkyCoord, Distance
import astropy.units as u
import numpy as np
from warnings import warn
import re
from os.path import join,abspath,dirname,exists,isfile
from os import listdir, getcwd
from shutil import copy
from astroquery.utils.tap.core import TapPlus
from io import StringIO
from sys import exit
from contextlib import redirect_stderr, redirect_stdout
import pickle
from .instrument import visibility, exposure_time
from . import __version__

__all__ = ['SpTypeToGminusV', 'SpTypeToTeff', '_GaiaDR2match']

# G-V transformation v. spectral type from Table 3. of Jordi et al.
# (2010A&A...523A..48J) and table (V-I)C v. spectral type from 
# http://www.stsci.edu/~inr/intrins.html (extrapolated to A0).

SpTypeToGminusV = {
'A0':-0.019, 'A1':-0.024, 'A2':-0.029, 'A3':-0.034, 'A4':-0.040,
'A5':-0.047, 'A6':-0.054, 'A7':-0.062, 'A8':-0.067, 'A9':-0.073,
'F0':-0.078, 'F1':-0.084, 'F2':-0.090, 'F3':-0.096, 'F4':-0.102,
'F5':-0.109, 'F6':-0.116, 'F7':-0.135, 'F8':-0.155, 'F9':-0.152,
'G0':-0.153, 'G1':-0.155, 'G2':-0.157, 'G3':-0.159, 'G4':-0.170,
'G5':-0.181, 'G6':-0.182, 'G7':-0.192, 'G8':-0.202, 'G9':-0.206,
'K0':-0.229, 'K1':-0.221, 'K2':-0.319, 'K3':-0.297, 'K4':-0.415,
'K5':-0.453, 'K6':-0.522, 'K7':-0.595, 'K8':-0.720, 'K9':-0.826,
'M0':-0.997, 'M1':-1.111, 'M2':-1.238, 'M3':-1.911, 'M4':-2.335,
'M5':-2.979, 'M6':-2.916, 'M7':-2.951, 'M8':-3.060, 'M9':-3.337 }

# From http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt
SpTypeToTeff = { 
    'A0':9700, 'A1':9200, 'A2':8840, 'A3':8550, 'A4':8270, 
    'A5':8080, 'A6':8000, 'A7':7800, 'A8':7500, 'A9':7440, 
    'F0':7220, 'F1':7030, 'F2':6810, 'F3':6720, 'F4':6640, 
    'F5':6510, 'F6':6340, 'F7':6240, 'F8':6170, 'F9':6040, 
    'G0':5920, 'G1':5880, 'G2':5770, 'G3':5720, 'G4':5680, 
    'G5':5660, 'G6':5590, 'G7':5530, 'G8':5490, 'G9':5340, 
    'K0':5280, 'K1':5170, 'K2':5040, 'K3':4830, 'K4':4600, 
    'K5':4410, 'K6':4230, 'K7':4070, 'K8':4000, 'K9':3940, 
    'M0':3870, 'M1':3700, 'M2':3550, 'M3':3410, 'M4':3200, 
    'M5':3030, 'M6':2850, 'M7':2650, 'M8':2500, 'M9':2400}

# Define a TapPlus query object for Gaia DR2
_gaia = TapPlus(url="http://gea.esac.esa.int/tap-server/tap",verbose=False)
_query = """SELECT source_id, ra, dec, parallax, pmra, pmdec, phot_g_mean_mag
FROM gaiadr2.gaia_source
WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec),
 CIRCLE('ICRS',{},{},0.0666))=1 AND (phot_g_mean_mag<=16.5);
"""


# XML strings and formats for input for Feasibility checker and PHT2
_xml_time_critical_fmt = """<?xml version="1.0" encoding="UTF-8"?>
<!--                                                                          -->
<!--         Template Observation Request file v. 9.1.1                       -->
<!--                                                                          -->
<!--         This file can be ingested in the CHEOPS Feasibility Checker      -->
<!--                                                                          -->
<!--         This file contains time-critical observations                    -->
<!--                                                                          -->
<!--         Generated by pycheops.make_xml_files                             -->
<!--                                                                          -->
<Earth_Explorer_File xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="ext_app_observation_requests_schema.xsd">
  <Earth_Explorer_Header>
    <Fixed_Header>
      <File_Name>CH_TU2018-03-05T10-00-00_EXT_APP_ObservationRequests_V1234</File_Name>
      <File_Description>Observation requests file</File_Description>
      <Notes>Template file for CHEOPS observation Request : PHT2, FeasibilityChecker, CHEOPSim</Notes>
      <Mission>CHEOPS</Mission>
      <File_Class>TEST</File_Class>
      <File_Type>EXT_APP_ObservationRequests</File_Type>
      <Validity_Period>
        <Validity_Start>UTC=2017-01-01T00:00:00</Validity_Start>
        <Validity_Stop>UTC=2023-01-01T00:00:00</Validity_Stop>
      </Validity_Period>
      <File_Version>0001</File_Version>
      <Source>
       <System>PSO</System>
       <Creator>PHT2</Creator>
       <Creator_Version>000</Creator_Version>
       <Creation_Date>UTC={}</Creation_Date>
      </Source>
    </Fixed_Header>
    <Variable_Header>
      <Programme_Type>10</Programme_Type>
    </Variable_Header>
  </Earth_Explorer_Header>
  
  <!-- NOTE: The listed Request was selected from the GTO-P2 catalog  -->
  <Data_Block type="xml">

    <!-- TIME-CRITICAL REQUEST  -->
    <List_of_Time_Critical_Requests count="1">
      <Time_Critical_Request>
        <!--  CHEOPSim : set the programme_ID  -->
        <Programme_ID>{:d}</Programme_ID>

        <!--  CHEOPSim : set the observation_request_ID  -->
        <Observation_Request_ID>1</Observation_Request_ID>

        <!--  PHT2, CHEOPSim : set the observation category ("time critical" or "non time critical")  -->
        <Observation_Category>time critical</Observation_Category>

        <!--  PHT2, CHEOPSim : set the propietary period after first visit (maximum is 18 months = 547 days)  -->
        <Proprietary_Period_First_Visit unit="days">{:d}</Proprietary_Period_First_Visit>

        <!--  PHT2, CHEOPSim : set the propietary period after last visit (maximum is 12 months = 365 days)  -->
        <Proprietary_Period_Last_Visit unit="days">{:d}</Proprietary_Period_Last_Visit>

        <!--  PHT2, FC, CHEOPSim : set the target name  -->
        <Target_Name>{}</Target_Name>

        <!--  PHT2, CHEOPSim/GTO : set the target GAIA_ID (https://gea.esac.esa.int/archive/)  -->
        <Gaia_ID>GAIA DR2 {}</Gaia_ID>

        <!--  PHT2, CHEOPSim : set the target spectral type  -->
        <Spectral_Type>{}V</Spectral_Type>

        <!--  PHT2, FC, CHEOPSim : set the target V-band magnitude  -->
        <Target_Vmagnitude unit="mag">{:0.2f}</Target_Vmagnitude>

        <!--  PHT2, CHEOPSim : set the target V-band magnitude error  -->
        <Target_Vmagnitude_Error unit="mag">{:0.2f}</Target_Vmagnitude_Error>

        <!--  PHT2, CHEOPSim : set the readout mode  -->
        <Readout_Mode>{}</Readout_Mode>

        <!--  PHT2, FC, CHEOPSim : set the target R.A.  -->
        <Right_Ascension unit="deg">{:0.7f}</Right_Ascension>

        <!--  PHT2, FC, CHEOPSim : set the target declination  -->
        <Declination unit="deg">{:0.7f}</Declination>

        <!--  PHT2, CHEOPSim/GTO : set the target R.A. proper motion  -->
        <RA_Proper_Motion unit="mas/year">{:0.2f}</RA_Proper_Motion>

        <!--  PHT2, CHEOPSim/GTO : set the target declination proper motion  -->
        <DEC_Proper_Motion unit="mas/year">{:0.2f}</DEC_Proper_Motion>

        <!--  PHT2, CHEOPSim/GTO : set the target parallax  -->
        <Parallax unit="mas">{:0.2f}</Parallax>

        <!--  PHT2, CHEOPSim/GTO : set the target Effective temperature  -->
        <T_Eff unit="Kelvin">{:0.0f}</T_Eff>

        <!--  PHT2 : set the V-band extinction toward the target  -->
        <Extinction unit="mag">0.00</Extinction>

        <!--  The following 2 parameters are used to bound the period when the visits should be scheduled  -->
        <!--  e.g. 15 Jan 2019 = 2458498.5 (JD)  -->
        <!--  e.g. 15 Mar 2019 = 2458529.5 (JD)  -->
        <!--  These 2 parameters are optional and can be commented out if not relevant for this observation  -->
        <!--  PHT2, FC : set the earliest start date (in BJD) for this observation request  -->
        <Earliest_Start unit="BJD">{:0.3f}</Earliest_Start>

        <!--  PHT2, FC : set the latest start date (in BJD) for this observation request  -->
        <Latest_End unit="BJD">{:0.3f}</Latest_End>

        <!--  PHT2, CHEOPSim : set the exposure time for this observation request   -->
        <Exposure_Time unit="sec">{:0.2f}</Exposure_Time>

        <!--  PHT2, CHEOPSim : set the exposure repetition period (optional) . Use fastest repetition period if not specified -->
        <Number_Stacked_Images>{:d}</Number_Stacked_Images>
        
        <!-- PHT2, CHEOPSim : mandatory but not relevant for the FC  -->
        <Number_Stacked_Imagettes>{:d}</Number_Stacked_Imagettes>

        <!--  PHT2, FC, CHEOPSim/GTO : set the central time of the transit (in BJD) (time critical only)  -->
        <Transit_Time unit="BJD">{:0.6f}</Transit_Time>

        <!--  PHT2, FC, CHEOPSim/GTO : set the orbital period of the planet (in days) (time critical only)  -->
        <Transit_Period unit="days">{:0.6f}</Transit_Period>

        <!--  1 orbit  =  6000 seconds  -->
        <!--  PHT2, FC, CHEOPSim : set the visit duration (in seconds)  -->
        <Visit_Duration unit="sec">{:d}</Visit_Duration>

        <!--  PHT2 : set the number of requested visit for this observation request  -->
        <Number_of_Visits>{:d}</Number_of_Visits>

        <Continuous_Visits>false</Continuous_Visits>   <!--  Irrelevant for nominal science observations  -->

        <!--  PHT2 : set the priority of the observation request (1 = A-grade, 2 = B-grade, 3 = C-grade)  -->
        <Priority>{:d}</Priority>

        <!-- This parameter defines the minimum on-source time relative to the visit duration  -->
        <!-- (excluding interruptions due to the SAA and Earth Occultations)  -->
        <!--  PHT2, FC : set the minimum effective duration, ~ observing efficiency (in percent)  -->
        <Minimum_Effective_Duration unit="%">{:d}</Minimum_Effective_Duration>

        <!--  This parameter defines the flexibility of the visit start time in units of planetary orbital phase.  -->
        <!--  Two values are given to define the allowed start time of the visit.  -->
        <!--  You should set "Earliest" strictly lower than "Latest" and the equivalent time-difference between  -->
        <!--    "Earliest" and "Latest" should be longer than 30 minutes  (limitation of current version of the tool)  -->
        <!--  PHT2, FC : set the start phase between 0 and 1 (mid-transit occurs at phase = 0)  -->
        <Earliest_Observation_Start unit="phase">{:0.4f}</Earliest_Observation_Start>
        <Latest_Observation_Start unit="phase">{:0.4f}</Latest_Observation_Start>

    	<!--  Optional  -->
    	<!--  The set of parameters below is used to define specific (orbital) phase ranges  -->
    	<!--  within which the observing efficiency may be increased to a specific value  -->
    	<!--  Convention is that the transit is at phase=0 (or equivalently 1)  -->
    	<!--  This can be seen as a local requirement on the observing efficiency (e.g. egresses)  -->
    	<!--  In this version, you should input "Start" < "End" (will be improved in future version)  -->
    	<!--  PHT2, FC : set the critical phase ranges  -->
{}

        <Send_Data_Taking_During_SAA>false</Send_Data_Taking_During_SAA>
        <Send_Data_Taking_During_Earth_Constraints>false</Send_Data_Taking_During_Earth_Constraints>
        <PITL>true</PITL>
      </Time_Critical_Request>
    </List_of_Time_Critical_Requests>
  </Data_Block>
</Earth_Explorer_File>
"""


_xml_non_time_critical_fmt = """<?xml version="1.0" encoding="UTF-8"?>
<!--                                                               -->
<!--  Template Observation Request file v. 9.1.1                   -->
<!--                                                               -->
<!--  This file can be ingested in the CHEOPS Feasibility Checker  -->
<!--                                                               -->
<!--                                                               -->
<!--  This file contains non-time-critical observations            -->
<!--                                                               -->
<!--  Generated by pycheops.make_xml_files                         -->
<!--                                                               -->
<Earth_Explorer_File xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="ext_app_observation_requests_schema.xsd">
  <Earth_Explorer_Header>
    <Fixed_Header>
      <File_Name>CH_TU2018-03-05T10-00-00_EXT_APP_ObservationRequests_V1234</File_Name>
      <File_Description>Observation requests file</File_Description>
      <Notes>Template file for CHEOPS observation Request : PHT2, FeasibilityChecker, CHEOPSim</Notes>
      <Mission>CHEOPS</Mission>
      <File_Class>TEST</File_Class>
      <File_Type>EXT_APP_ObservationRequests</File_Type>
      <Validity_Period>
        <Validity_Start>UTC=2017-01-01T00:00:00</Validity_Start>
        <Validity_Stop>UTC=2023-01-01T00:00:00</Validity_Stop>
      </Validity_Period>
      <File_Version>0001</File_Version>
      <Source>
       <System>PSO</System>
       <Creator>PHT2</Creator>
       <Creator_Version>000</Creator_Version>
       <Creation_Date>UTC={}</Creation_Date>
      </Source>
    </Fixed_Header>
    <Variable_Header>
      <Programme_Type>10</Programme_Type>
    </Variable_Header>
  </Earth_Explorer_Header>
  
  <!-- NOTE: The listed Request was selected from the GTO-P2 catalog  -->
  <Data_Block type="xml">

    <!-- NON-TIME-CRITICAL REQUEST  -->
    <List_of_Non_Time_Critical_Requests count="1">
      <Non_Time_Critical_Request>
        <!--  CHEOPSim : set the programme_ID  -->
        <Programme_ID>{:d}</Programme_ID>
    
        <!--  CHEOPSim : set the observation_request_ID  -->
        <Observation_Request_ID>1</Observation_Request_ID>
    
        <!--  PHT2, CHEOPSim : set the observation category ("time critical" or "non time critical")  -->
        <Observation_Category>non time critical</Observation_Category>
    
        <!--  PHT2, CHEOPSim : set the propietary period after first visit (maximum is 18 months = 547 days)  -->
        <Proprietary_Period_First_Visit unit="days">{:d}</Proprietary_Period_First_Visit>
    
        <!--  PHT2, CHEOPSim : set the propietary period after last visit (maximum is 12 months = 365 days)  -->
        <Proprietary_Period_Last_Visit unit="days">{:d}</Proprietary_Period_Last_Visit>
    
        <!--  PHT2, FC, CHEOPSim : set the target name  -->
        <Target_Name>{}</Target_Name>
    
        <!--  PHT2, CHEOPSim/GTO : set the target GAIA_ID (https://gea.esac.esa.int/archive/)  -->
        <Gaia_ID>GAIA DR2 {}</Gaia_ID>
    
        <!--  PHT2, CHEOPSim : set the target spectral type  -->
        <Spectral_Type>{}V</Spectral_Type>
    
        <!--  PHT2, FC, CHEOPSim : set the target V-band magnitude  -->
        <Target_Vmagnitude unit="mag">{:0.2f}</Target_Vmagnitude>
    
        <!--  PHT2, CHEOPSim : set the target V-band magnitude error  -->
        <Target_Vmagnitude_Error unit="mag">{:0.2f}</Target_Vmagnitude_Error>
    
        <!--  PHT2, CHEOPSim : set the readout mode  -->
        <Readout_Mode>{}</Readout_Mode>
    
        <!--  PHT2, FC, CHEOPSim : set the target R.A.  -->
        <Right_Ascension unit="deg">{:0.5f}</Right_Ascension>
    
        <!--  PHT2, FC, CHEOPSim : set the target declination  -->
        <Declination unit="deg">{:0.5f}</Declination>
    
        <!--  PHT2, CHEOPSim/GTO : set the target R.A. proper motion  -->
        <RA_Proper_Motion unit="mas/year">{:0.2f}</RA_Proper_Motion>
    
        <!--  PHT2, CHEOPSim/GTO : set the target declination proper motion  -->
        <DEC_Proper_Motion unit="mas/year">{:0.2f}</DEC_Proper_Motion>
    
        <!--  PHT2, CHEOPSim/GTO : set the target parallax  -->
        <Parallax unit="mas">{:0.2f}</Parallax>
    
        <!--  PHT2, CHEOPSim/GTO : set the target Effective temperature  -->
        <T_Eff unit="Kelvin">{:0.0f}</T_Eff>
    
        <!--  PHT2 : set the V-band extinction toward the target  -->
        <Extinction unit="mag">0.00</Extinction>
    
        <!--  The following 2 parameters are used to bound the period when the visits should be scheduled  -->
        <!--  e.g. 15 Jan 2019 = 2458498.5 (JD)  -->
        <!--  e.g. 15 Mar 2019 = 2458529.5 (JD)  -->
        <!--  These 2 parameters are optional and can be commented out if not relevant for this observation  -->
        <!--  PHT2, FC : set the earliest start date (in BJD) for this observation request  -->
        <Earliest_Start unit="BJD">{:0.3f}</Earliest_Start>
    
        <!--  PHT2, FC : set the latest start date (in BJD) for this observation request  -->
        <Latest_End unit="BJD">{:0.3f}</Latest_End>
    
        <!--  PHT2, CHEOPSim : set the exposure time for this observation request   -->
        <Exposure_Time unit="sec">{:0.2f}</Exposure_Time>
    
        <!--  PHT2, CHEOPSim : set the exposure repetition period (optional) . Use fastest repetition period if not specified -->
        <Number_Stacked_Images>{:d}</Number_Stacked_Images>
        
        <!-- PHT2, CHEOPSim : mandatory but not relevant for the FC  -->
        <Number_Stacked_Imagettes>{:d}</Number_Stacked_Imagettes>

        <!--  1 orbit  =  6000 seconds  -->
        <!--  PHT2, FC, CHEOPSim : set the visit duration (in seconds)  -->
        <Visit_Duration unit="sec">{:d}</Visit_Duration>
    
        <!--  PHT2 : set the number of requested visit for this observation request  -->
        <Number_of_Visits>{:d}</Number_of_Visits>
    
        <Continuous_Visits>false</Continuous_Visits>   <!--  Irrelevant for nominal science observations  -->
    
        <!--  PHT2 : set the priority of the observation request (1 = A-grade, 2 = B-grade, 3 = C-grade)  -->
        <Priority>{:d}</Priority>
    
    
        <!-- This parameter defines the minimum on-source time relative to the visit duration  -->
        <!-- (excluding interruptions due to the SAA and Earth Occultations)  -->
        <!--  PHT2, FC : set the minimum effective duration, ~ observing efficiency (in percent)  -->
        <Minimum_Effective_Duration unit="%">{:d}</Minimum_Effective_Duration>
    
        <Send_Data_Taking_During_SAA>false</Send_Data_Taking_During_SAA>
        <Send_Data_Taking_During_Earth_Constraints>false</Send_Data_Taking_During_Earth_Constraints>
        <PITL>true</PITL>
        </Non_Time_Critical_Request>
    </List_of_Non_Time_Critical_Requests>
  </Data_Block>
</Earth_Explorer_File>
"""

_phase_range_format_1 = """
    <List_of_Phase_Ranges count="1">
        <Phase_Range>
            <Start unit="phase">{:0.4f}</Start>
            <End unit="phase">{:0.4f}</End>
            <Minimum_Phase_Duration unit="%">{:d}</Minimum_Phase_Duration>
        </Phase_Range>
    </List_of_Phase_Ranges>
"""


_phase_range_format_2 = """
    <List_of_Phase_Ranges count="2">
        <Phase_Range>
            <Start unit="phase">{:0.4f}</Start>
            <End unit="phase">{:0.4f}</End>
            <Minimum_Phase_Duration unit="%">{:d}</Minimum_Phase_Duration>
        </Phase_Range>
        <Phase_Range>
            <Start unit="phase">{:0.4f}</Start>
            <End unit="phase">{:0.4f}</End>
            <Minimum_Phase_Duration unit="%">{:d}</Minimum_Phase_Duration>
        </Phase_Range>
    </List_of_Phase_Ranges>
      <!--  If two critical phase ranges are defined above, this parameter is used to request that both ("true") or                             -->
      <!--         only one of the two phase ranges ("false") are observed. This can be seen as a AND / OR operator, respectively.              -->
      <!--  #############################   Set the critical phase ranges    -->
      <Fulfil_all_Phase_Ranges>{}</Fulfil_all_Phase_Ranges>
"""

def _GaiaDR2Match(row, fC, match_radius=1,  gaia_mag_tolerance=0.5, 
        id_check=True):

    flags = 0 

    coo = SkyCoord(row['_RAJ2000'],row['_DEJ2000'],
            frame='icrs',unit=(u.hourangle, u.deg))
    s = coo.to_string('decimal',precision=5).split()
    j = StringIO()
    with redirect_stderr(j), redirect_stdout(j):
        job = _gaia.launch_job(_query.format(s[0],s[1]))
    DR2Table = job.get_results()

    # Replace missing values for pmra, pmdec, parallax
    DR2Table['pmra'].fill_value = 0.0
    DR2Table['pmdec'].fill_value = 0.0
    DR2Table['parallax'].fill_value = 0.0
    DR2Table = Table(DR2Table.filled(), masked=True)
    # Avoid problems with small/negative parallaxes
    DR2Table['parallax'].mask = DR2Table['parallax'] <= 0.1
    DR2Table['parallax'].fill_value = 0.0999
    DR2Table = DR2Table.filled()
    # Fix units for proper motion columns
    DR2Table['pmra'].unit = 'mas / yr'
    DR2Table['pmdec'].unit = 'mas / yr'

    cat = SkyCoord(DR2Table['ra'],DR2Table['dec'],
            frame='icrs',
            distance=Distance(parallax=DR2Table['parallax'].quantity),
            pm_ra_cosdec=DR2Table['pmra'], pm_dec=DR2Table['pmdec'],
            obstime=Time(2015.5, format='decimalyear')
           ).apply_space_motion(new_obstime=Time('2000-01-01 00:00:00.0'))
    idx, d2d, _ = coo.match_to_catalog_sky(cat)
    if d2d > match_radius*u.arcsec:
        raise ValueError('No Gaia DR2 source within specified match radius')

    try:
        key = re.match('[AFGKM][0-9]', row['SpTy'])[0]
        GV = SpTypeToGminusV[key]
    except TypeError:
        flags += 1024
        GV = 0.0
    Gmag = row['Vmag'] + GV

    if abs(Gmag-DR2Table['phot_g_mean_mag'][idx]) > gaia_mag_tolerance:
        print("Input values: V = {:5.2f}, SpTy = {} -> G_est = {:5.2f}"
                .format(row['Vmag'], row['SpTy'], Gmag))
        print("Catalogue values: G = {:5.2f}, Source = {}"
                .format(DR2Table['phot_g_mean_mag'][idx], 
                    DR2Table['source_id'][idx] ))
        raise ValueError('Nearest Gaia source does not match estimated G mag')

    if (str(row['Old_Gaia_DR2']) != str(DR2Table['source_id'][idx])):
        if id_check:
            raise ValueError('Nearest Gaia DR2 source does not match input ID')
        flags += 32768

    gmag = np.array(DR2Table['phot_g_mean_mag'])
    sep = coo.separation(cat)
    if any((sep <= 51*u.arcsec) & (gmag < gmag[idx])):
        flags += 16384
    if any((sep > 51*u.arcsec) & (sep < 180*u.arcsec) & (gmag < gmag[idx])):
        flags += 8192

    gflx = np.ma.array(10**(-0.4*(gmag-gmag[idx])), mask=False,fill_value=0.0)
    gflx.mask[idx] = True
    contam = np.nansum(gflx.filled()*fC(cat.separation(cat[idx]).arcsec))

    if contam > 1:
        flags += 4096
    elif contam > 0.1:
        flags += 2048

    return DR2Table[idx], contam, flags, cat[idx]

def _choose_stacking(Texp):
    if Texp < 0.1:
        return 40, 4
    elif Texp < 0.15:
        return 39, 3
    elif Texp < 0.20:
        return 36, 3
    elif Texp < 0.40:
        return 33, 3
    elif Texp < 0.50:
        return 30, 3
    elif Texp < 0.55:
        return 28, 2
    elif Texp < 0.65:
        return 26, 2
    elif Texp < 0.85:
        return 24,23
    elif Texp < 1.05:
        return 22, 2
    elif Texp < 1.10:
        return 44, 4
    elif Texp < 1.20:
        return 40, 4
    elif Texp < 1.25:
        return 39, 3
    elif Texp < 1.30:
        return 36, 3
    elif Texp < 1.50:
        return 33, 3
    elif Texp < 1.60:
        return 30, 3
    elif Texp < 1.65:
        return 28, 2
    elif Texp < 1.75:
        return 26, 2
    elif Texp < 1.95:
        return 24, 2
    elif Texp < 2.15:
        return 22, 2
    elif Texp < 2.40:
        return 20, 2
    elif Texp < 2.70:
        return 18, 2
    elif Texp < 2.80:
        return 16, 2
    elif Texp < 2.90:
        return 15, 1
    elif Texp < 3.05:
        return 14, 1
    elif Texp < 3.20:
        return 13, 1
    elif Texp < 3.40:
        return 12, 1
    elif Texp < 3.65:
        return 11, 1
    elif Texp < 3.90:
        return 10, 1
    elif Texp < 4.25:
        return 9, 1
    elif Texp < 4.70:
        return 8, 1
    elif Texp < 5.25:
        return 7, 1
    elif Texp < 6.05:
        return 6, 1
    elif Texp < 7.25:
        return 5, 1
    elif Texp < 9.20:
        return 4, 1
    elif Texp < 12.5:
        return 3, 1
    elif Texp < 22.65:
        return 2, 1
    else:
        return 1, 0

def _choose_romode(t_exp):
    if t_exp < 1.05:
        return 'ultrabright'
    if t_exp < 2.226:
        return 'bright'
    if t_exp < 12:
        return 'faint fast'
    return 'faint'

def _creation_time_string():
    t = Time.now()
    t.precision = 0
    return t.isot

def _make_list_of_phase_ranges(Num_Ranges, 
        BegPh1, EndPh1, Effic1,
        BegPh2, EndPh2, Effic2):

    if Num_Ranges == 1 :
        return _phase_range_format_1.format(BegPh1, EndPh1, Effic1)

    if Num_Ranges == 2 :
        return _phase_range_format_2.format(BegPh1, EndPh1, Effic1,
                BegPh2, EndPh2, Effic2, 'true')

    if Num_Ranges == -2 :
        return _phase_range_format_2.format(BegPh1, EndPh1, Effic1,
                BegPh2, EndPh2, Effic2, 'false')

    return ""



def _parcheck_non_time_critical(Priority, MinEffDur, 
        Earliest_start_date, Latest_end_date):

    if not Priority in (1, 2, 3):
        return """
        The priority has to be set equal to 1, 2, or 3: 1 = A-grade, 2 = B-grade, 3 = C-grade
        """

    if (MinEffDur < 0) or (MinEffDur > 100):
        return """
        The minimum effective duration is in % and it has to be between 0 and 100
        """

    if ( (Earliest_start_date > 0) and (Latest_end_date > 0) and 
         Earliest_start_date >= Latest_end_date) :
        return """
        The earliest start date must be less than the latest start date
        """
    return None

def _parcheck_time_critical(Priority, MinEffDur,
        Earliest_start_date, Latest_end_date,
        Earliest_start_phase, Latest_start_phase, 
        Period, Num_Ranges, 
        BegPh1, EndPh1, Effic1,
        BegPh2, EndPh2, Effic2):

    if (((Earliest_start_phase > -50) and (Earliest_start_phase < 0)) or 
            (Earliest_start_phase > 1)) :
        return """
        The earliest start phase should be a number between 0 and 1, inclusive
        """
    if (((Latest_start_phase > -50) and (Latest_start_phase < 0)) or 
            (Latest_start_phase > 1)) :
        return """
        The latest start phase should be a number between 0 and 1, inclusive
        """
    if not Num_Ranges in (-2,0,1,2):
        return """
        The number of constrained ranges is invalid (not 0, 1, 2 or -2)
        """
    if abs(Num_Ranges) > 0:
        if (BegPh1 < 0) or (BegPh1 > 1) or (EndPh1 < 0) or (EndPh1 > 1):
            return """
            Invalid phase range for phase constraint 1
            """
        if (Effic1 < 0) or (Effic1 > 99):
            return """
            Invalid efficiency for phase constraint 1
            """

    if abs(Num_Ranges) > 1:
        if (BegPh2 < 0) or (BegPh2 > 1) or (EndPh2 < 0) or (EndPh2 > 1):
            return """
            Invalid phase range for phase constraint 2
            """
        if (Effic2 < 0) or (Effic2 > 99):
            return """
            Invalid efficiency for phase constraint 2
            """

    return _parcheck_non_time_critical(Priority, MinEffDur, Earliest_start_date,
        Latest_end_date)

def _target_table_row_to_xml(row, 
        progamme_id=0, proprietary_first=547, proprietary_last=365):

    period = row['Period'] 
    t_exp = row['T_exp']
    n_stack_image, n_stack_imagettes = _choose_stacking(t_exp)

    c = SkyCoord("{} {}".format(row['_RAJ2000'],row['_DEJ2000']),
            frame='icrs', obstime='J2000.0',
            unit=(u.hourangle, u.deg),
            pm_ra_cosdec=row['pmra']*u.mas/u.yr,
            pm_dec=row['pmdec']*u.mas/u.yr )
    radeg = float(c.to_string(precision=5).split()[0])
    dedeg = float(c.to_string(precision=5).split()[1])
    

    if period > 0:
        error = _parcheck_time_critical(
                row['Priority'], row['MinEffDur'],
                row['BJD_early'], row['BJD_late'],
                row['Ph_early'], row['Ph_late'],
                period, row['N_Ranges'],
                row["BegPh1"], row["EndPh1"], row["Effic1"],
                row["BegPh2"], row["EndPh2"], row["Effic2"])

        assert error is None, (
                "Failed to process data for observing request {}\n{}\n"
                .format(row['ObsReqName'],error))

        xml = _xml_time_critical_fmt.format(
              _creation_time_string(),
              progamme_id, proprietary_first, proprietary_last,
              row['Target'], row['Gaia_DR2'], row['SpTy'],
              row['Vmag'], row['e_Vmag'], 
              _choose_romode(t_exp),
              radeg, dedeg, row['pmra'], row['pmdec'],
              row['parallax'], row['T_eff'], 
              row['BJD_early'], row['BJD_late'], t_exp,
              n_stack_image, n_stack_imagettes,
              row['BJD_0'], period, 
              row['T_visit'], row['N_Visits'], row['Priority'],
              row['MinEffDur'],
              row['Ph_early'], row['Ph_late'],
              _make_list_of_phase_ranges(row['N_Ranges'],
                  row["BegPh1"], row["EndPh1"], row["Effic1"],
                  row["BegPh2"], row["EndPh2"], row["Effic2"])
              )
      
    else:
        error = _parcheck_non_time_critical(
                row['Priority'], row['MinEffDur'],
                row['BJD_early'], row['BJD_late'])

        assert error is None, (
                "Failed to process data for observing request {}\n{}\n"
                .format(row['ObsReqName'],error))

        xml = _xml_non_time_critical_fmt.format(
              _creation_time_string(),
              progamme_id, proprietary_first, proprietary_last,
              row['Target'], row['Gaia_DR2'], row['SpTy'],
              row['Vmag'], row['e_Vmag'], 
              _choose_romode(t_exp),
              radeg, dedeg, row['pmra'], row['pmdec'],
              row['parallax'], row['T_eff'],
              row['BJD_early'], row['BJD_late'], t_exp,
              n_stack_image, n_stack_imagettes,
              row['T_visit'], row['N_Visits'], row['Priority'],
              row['MinEffDur']
              )
    return xml


def main():

    # Set up command line switches
    parser = argparse.ArgumentParser(
        description='Create xml files for CHEOPS observing requests.',
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog = textwrap.dedent('''\

        Creates XML files suitable for CHEOPS PHT2, FC, and CHEOPSim from
        observation requests given in an input table file. 

        The target for each observation is defined by the input _RAJ2000 and
        _DEJ2000 coordinates. There must be a matching source in Gaia DR2 for
        each input coordinate. The G-band magnitude of the source must also
        match the G-band magnitude estimated from Vmag and SpTy. The following
        table is an abbreviated version of the look-up table used to estimate
        the G-band magnitude from  Vmag and SpTy.

          SpTy G-V             SpTy G-V
          -----------          -----------
          A0   -0.019          K0   -0.229
          F0   -0.078          K5   -0.453
          F5   -0.109          M0   -0.997
          G0   -0.153          M5   -2.979
          G5   -0.181          M9   -3.337

        The input table can be any format suitable for reading with the
        command astropy.table.Table.read(), e.g., CSV.

        The following columns must be defined in the table.
        
         ObsReqName - unique observing request identifier
         Target     - target name
         _RAJ2000   - right ascension, ICRS epoch J2000.0, hh:mm:ss.ss
         _DEJ2000   - declination, ICRS epoch J2000.0, +dd:mm:ss.s
         SpTy       - spectral type (any string starting [AFGKM][0-9])
         Vmag       - V-band magnitude
         e_Vmag     - error in V-band magnitude
         BJD_early  - earliest start date (BJD)
         BJD_late   - latest start date (BJD) 
         T_visit    - visit duration in seconds
         N_Visits   - number of requested visits
         Priority   - 1, 2 or 3 
         MinEffDur  - minimum on-source time, percentage of T_visit (integer)

        If the flag --ignore-gaia-id-check is not specified on the command
        line then the following column is also required.
         Gaia_DR2   - Gaia DR2 identification number (integer)

        If the flag --auto-expose is not specified on the command
        line then the following column is also required.
         T_exp      - exposure time (seconds)

        In addition, for time-critical observations the following columns must 
        also be defined. 
         BJD_0      - reference time for 0 phase (e.g., mid-transit), BJD
         Period     - period in days
         Ph_early   - earliest allowable start phase for visit
         Ph_late    - latest allowable start phase for visit

        The following columns will also be used if available. 
         N_Ranges   - number of phase ranges with extra efficiency constraints
         BegPh1     - start of phase range 1
         EndPh1     - end of phase range 1
         Effic1     - minimum observing efficiency (%), phase range 1 (integer)
         BegPh2     - start of phase range 1
         EndPh2     - end of phase range 1
         Effic2     - minimum observing efficiency (%), phase range 2 (integer)
        N.B. If you have 2 phase ranges with extra efficiency constraints but
        only require one of them to be satisified then use N_Ranges = -2
         
        The terminal output includes the following columns

         Gaia_DR2_ID - Gaia DR2 ID from Gaia data archive. This must match
             the value of Gaia_DR2 in the input file unless the flag
             --ignore-gaia-id-check is specified. 
             ** N.B. The PI is responsible to check the DR2 ID is correct **

         Gmag - The target mean G-band magnitude from Gaia DR2 catalogue.

         Contam - estimate of the contamination of a 30 arcsec photometric
             aperture by nearby stars relative to the target flux.

         Vis - estimate of the percentage of the orbit for which the target is
               observable by CHEOPS. This estimate is not a substitute for the
               detailed scheduling information provided by the CHEOPS
               Feasibility Checker.

        Texp - the exposure time used in the output XML file.

         _RAJ2000,_DEJ2000 - ICRS position of matching Gaia source in degrees

         Flags - sum of the following error/warnings flags.
             + 32768 = Gaia ID error - input/output IDs do not match
             + 16384 = Acquisition error, brighter star within 51"
             +  8192 = Acquisition warning, brighter star within 51"-180"
             +  4096 = Contamination error, Contam > 1
             +  2048 = Contamination warning, Contam = 0.1 - 1
             +  1024 = No spectral type match, assuming G-V = 0 
             +  512 = Visibility error, efficiency = 0
             +  256 = Visibility warning, efficiency < 50%
             +  128 = Exposure time error - target will be saturated
             +  64 = Exposure time warning - below recommended minimum time
             +  32 = Exposure time error - magnitude out of range, not set 
             +  16 = Exposure time warning - magnitude out of range, not checked

        N.B. Automatic exposure times only reliable for the range G=5.847 to 
        G=12.847

        See examples/make_xml_files/ReadMe.txt in the source distribution for a
        description of example input files included in the same folder.
        --
              
        '''))

    parser.add_argument('table', nargs='?',
        help='Table of observing requests to be processed into XML'
    )

    parser.add_argument('-p', '--programme_id', 
        default=1, type=int,
        help='''Programme ID 
        (default: %(default)d)
        '''
    )

    parser.add_argument('-r', '--match_radius', 
        default=1.0, type=float,
        help='''
        Tolerance in arcsec for cross-matching with Gaia DR2 (default:
        %(default)3.1f)
        '''
    )

    parser.add_argument('-g', '--gaia_mag_tolerance', 
        default=0.5, type=float,
        help= '''
        Tolerance in magnitudes for Gaia DR2 cross-match (default:
        %(default)3.1f)
        '''
    )

    parser.add_argument('--ignore-gaia-id-check',
        action='store_const',
        dest='id_check',
        const=False,
        default=True,
        help='''
        Use Gaia DR2 ID from Gaia data archive
        ** N.B. The PI is responsible to check the DR2 ID is correct **
        '''
    )

    parser.add_argument('-a', '--auto-expose', 
        action='store_const',
        dest='auto_expose',
        const=True,
        default=False,
        help='Use recommended maximum exposure time'
    )

    parser.add_argument('-e', '--example-file-copy', 
        action='store_const',
        dest='copy_examples',
        const=True,
        default=False,
        help='Get a copy of the example files - no other action is performed'
    )

    parser.add_argument('-f', '--overwrite', 
        action='store_const',
        dest='overwrite',
        const=True,
        default=False,
        help='Overwrite existing output files.'
    )

    parser.add_argument('-x', '--suffix', 
        default='.xml', type=str,
        help='''
        Output file name suffix
        (default: %(default)s)
        '''
    )

    parser.add_argument('--proprietary_last', 
        default=365, type=int,
        help='Propietary period after last visit'
    )

    parser.add_argument('--proprietary_first', 
        default=547, type=int,
        help='Propietary period after first visit'
    )

    args = parser.parse_args()

    if args.copy_examples:
        src = join(dirname(abspath(__file__)),'examples','make_xml_files')
        src_files = listdir(src)
        print("Copy examples files from {}".format(src))
        for file_name in src_files:
            full_file_name = join(src, file_name)
            if (isfile(full_file_name)):
                copy(full_file_name, getcwd())
        exit()


    if args.table is None:
        parser.print_usage()
        exit(1)

    table = Table.read(args.table)

    if len(set(table['ObsReqName'])) < len(table):
        raise ValueError("Duplicate observing request names in {}"
                .format(args.table))

    try:
        table['Old_Gaia_DR2'] = table['Gaia_DR2']
    except KeyError as e:
        if args.id_check:
            message = e.args[0]
            message += (" - use flag --ignore-gaia-id-check to insert GAIA"
                    " DR2 identifiers from Gaia data archive. ** N.B. The PI"
                    " is responsible to check the DR2 ID is correct ** " )
            e.args = (message,)
            raise
        else:
            table['Gaia_DR2'] = -1
            table['Old_Gaia_DR2'] = -1


    try:
        table['Old_T_exp'] = table['T_exp']
    except KeyError as e:
        if args.auto_expose:
            table['T_exp'] = -1
            table['Old_T_exp'] = -1
        else:
            message = e.args[0]
            message += (" - use flag --auto-expose to use recommended maximum"
                    " exposure time")
            e.args = (message,)
            raise

    for key in ('T_eff',):
        try:
            table['Old_{}'.format(key)] = table[key]
        except KeyError:
            table[key] = -1

    for key in ('pmra', 'pmdec', 'parallax'):
        try:
            table['Old_{}'.format(key)] = table[key]
        except KeyError:
            table[key] = 0.0

    # Create missing optional columns
    for key in ( 'Period', 'BJD_0', 'Ph_early', 'Ph_late', 'BegPh1', 
            'EndPh1', 'Effic1', 'BegPh2', 'EndPh2', 'Effic2'):
        if not key in table.columns:
            table[key] = 0.0

    if not 'N_Ranges' in table.columns:
        table['N_Ranges'] = 0

    # Load contamination function from pickle
    data_path = join(dirname(abspath(__file__)),'data','instrument')
    pfile = join(data_path,'Contamination_30arcsec_aperture.p')
    with open(pfile, 'rb') as fp: 
        fC= pickle.load(fp)

    rtol = args.match_radius
    gtol = args.gaia_mag_tolerance

    # Header to screen output
    print('# Output from: {} version {}'.format(parser.prog, __version__))
    print('# Run started: {}'.format(Time(Time.now(),precision=0).iso))
    print('# Input file: {}'.format(args.table))
    print('# Gaia match radius: {:0.1f} arcsec'.format(rtol))
    print('# Gmag tolerance: {:0.1f} mag '.format(gtol))
    print('# Output file suffix: {} '.format(args.suffix))
    ObsReqNameFieldWidth = max(12, len(max(table['ObsReqName'],key=len)))
    ObsReqNameFormat = "{{:{}s}}".format(ObsReqNameFieldWidth)
    ObsReqNameHeader = ObsReqNameFormat.format('ObsReqName')
    TerminalOutputFormat = (
        '{}'.format(ObsReqNameFormat)+
        '{:20d} {:5.2f} {:8.4f} {:+8.4f}  {:6.3f}  {:2d} {:4.1f} {:5d}')
    print('#')
    if not args.id_check:
        print('#')
        print('# ** WARNING: Gaia ID of target not checked against input **')
        print('# ** The PI is responsible to check the DR2 ID is correct **')
        print('#')

    print('#{}'.format(ObsReqNameHeader) +
    'Gaia_DR2_ID         Gmag  _RAJ2000 _DEJ2000  Contam Vis Texp Flags'
           .format(ObsReqNameHeader))

    # String of coordinates, Vmag and SpTy to enable re-use of DR2 data
    old_tag = None
    for row in table:

        coo = SkyCoord(row['_RAJ2000'],row['_DEJ2000'],
              frame='icrs',unit=(u.hourangle, u.deg))

        tag = "{}, {}, {}".format(coo.to_string(), row['Vmag'], row['SpTy'])
        if tag != old_tag:
            old_tag = tag
            DR2data,contam,flags,coords = _GaiaDR2Match(row, fC, rtol, gtol,
                    args.id_check)

        rastr  = coords.ra.to_string('hour', precision=2, sep=':', pad=True)
        decstr = coords.dec.to_string('deg', precision=1, sep=':', pad=True, 
                alwayssign=True)
        row['Gaia_DR2'] = DR2data['source_id']
        row['_RAJ2000'] = rastr
        row['_DEJ2000'] = decstr
        row['pmra'] = DR2data['pmra']
        row['pmdec'] = DR2data['pmdec']
        row['parallax'] = DR2data['parallax']

        tmin, tmax = exposure_time(DR2data['phot_g_mean_mag'])
        if not np.isfinite(tmax):
            if args.auto_expose:
                flags += 32
            else:
                flags += 16
            if DR2data['phot_g_mean_mag'] > 13:
                tmin, tmax = 35, 60
            else:
                tmin, tmax = 0.05, 0.5

        if args.auto_expose:
            row['T_exp'] = tmax
        else:
            if row['T_exp'] > tmax:
                flags += 128
            if row['T_exp'] < tmin:
                flags += 64

        try:
            if row['Old_T_eff'] <= 0:
                raise KeyError
            row['T_eff'] = row['Old_T_eff']
        except KeyError:
            try:
                key = re.match('[AFGKM][0-9]', row['SpTy'])[0]
                row['T_eff'] =  SpTypeToTeff[key]
            except KeyError:
                warn('# No Teff value for spectral type, using Teff=5999')
                row['T_eff'] = 5999


        xmlfile = "{}{}".format(row['ObsReqName'],args.suffix)
        if exists(xmlfile) and not args.overwrite:
            raise IOError("Output file {} exists, use -f option to overwrite"
                    .format(xmlfile))
        f = open(xmlfile,'w')
        f.write(_target_table_row_to_xml(row,
                    progamme_id=args.programme_id, 
                    proprietary_first=args.proprietary_first,
                    proprietary_last=args.proprietary_last)
            )
        f.close()

        vis = visibility(coords.ra.degree,coords.dec.degree)
        if vis < 50:
            flags += 256
        if vis == 0:
            flags += 256



        print(TerminalOutputFormat.format(
            row['ObsReqName'], DR2data['source_id'], DR2data['phot_g_mean_mag'],
            coords.ra.degree,coords.dec.degree, contam, vis, 
            row['T_exp'],flags))

