# -*- 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
# Suppress output to stdout on import of astroquery.gaia.Gaia
from contextlib import redirect_stdout, redirect_stderr
from io import StringIO
_ = StringIO()
with redirect_stdout(_):
    from astroquery.gaia import Gaia
from sys import exit
from .core import load_config
import pickle
from .instrument import visibility, exposure_time, count_rate, cadence
from . import __version__

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

# G-V and Teff  v. spectral type from 
# http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt
# version 2019.03.22
SpTypeToGminusV = {
'A0':+0.007, 'A1':+0.000, 'A2':+0.005, 'A3':-0.009, 'A4':-0.020,
'A5':-0.024, 'A6':-0.026, 'A7':-0.036, 'A8':-0.046, 'A9':-0.047,
'F0':-0.060, 'F1':-0.079, 'F2':-0.093, 'F3':-0.100, 'F4':-0.107,
'F5':-0.116, 'F6':-0.129, 'F7':-0.135, 'F8':-0.140, 'F9':-0.146,
'G0':-0.155, 'G1':-0.162, 'G2':-0.167, 'G3':-0.169, 'G4':-0.172,
'G5':-0.174, 'G6':-0.180, 'G7':-0.182, 'G8':-0.188, 'G9':-0.204,
'K0':-0.221, 'K1':-0.232, 'K2':-0.254, 'K3':-0.322, 'K4':-0.412,
'K5':-0.454, 'K6':-0.528, 'K7':-0.595, 'K8':-0.628, 'K9':-0.69 ,
'M0':-0.65 , 'M1':-0.82 , 'M2':-0.92 , 'M3':-1.09 , 'M4':-1.41 ,
'M5':-1.74 , 'M6':-2.14 , 'M7':-2.98 , 'M8':-3.08 , 'M9':-3.00 }

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':6060,
'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 query object for Gaia DR2
_query = """SELECT source_id, ra, dec, parallax, pmra, pmdec, \
phot_g_mean_mag, phot_g_mean_flux_over_error, bp_rp 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"?>
<!--                                                                          -->
<!--         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>
  
  <Data_Block type="xml">

    <!-- TIME-CRITICAL REQUEST  -->
    <List_of_Time_Critical_Requests count="1">
      <Time_Critical_Request>
        <Programme_ID>{:d}</Programme_ID>
        <Observation_Request_ID>1</Observation_Request_ID>
        <Observation_Category>time critical</Observation_Category>
        <Proprietary_Period_First_Visit unit="days">{:d}</Proprietary_Period_First_Visit>
        <Proprietary_Period_Last_Visit unit="days">{:d}</Proprietary_Period_Last_Visit>
        <Target_Name>{}</Target_Name>
        <Gaia_ID>GAIA DR2 {}</Gaia_ID>
        <Spectral_Type>{}</Spectral_Type>
        <Target_Magnitude unit="mag">{:0.2f}</Target_Magnitude>
        <Target_Magnitude_Error unit="mag">{:0.2f}</Target_Magnitude_Error>
        <Readout_Mode>{}</Readout_Mode>
        <Right_Ascension unit="deg">{:0.7f}</Right_Ascension>
        <Declination unit="deg">{:0.7f}</Declination>
        <RA_Proper_Motion unit="mas/year">{:0.2f}</RA_Proper_Motion>
        <DEC_Proper_Motion unit="mas/year">{:0.2f}</DEC_Proper_Motion>
        <Parallax unit="mas">{:0.2f}</Parallax>
        <T_Eff unit="Kelvin">{:0.0f}</T_Eff>
        <Extinction unit="mag">0.00</Extinction>
        <Earliest_Start unit="BJD">{:0.3f}</Earliest_Start>
        <Latest_End unit="BJD">{:0.3f}</Latest_End>
        <Exposure_Time unit="sec">{:0.2f}</Exposure_Time>
        <Number_Stacked_Images>{:d}</Number_Stacked_Images>
        <Number_Stacked_Imagettes>{:d}</Number_Stacked_Imagettes>
        <Transit_Time unit="BJD">{:0.6f}</Transit_Time>
        <Transit_Period unit="days">{:0.6f}</Transit_Period>
        <Visit_Duration unit="sec">{:0.1f}</Visit_Duration>
        <Number_of_Visits>{:d}</Number_of_Visits>
        <Continuous_Visits>false</Continuous_Visits>   <!--  Irrelevant for nominal science observations  -->
        <Priority>{:d}</Priority>
        <Minimum_Effective_Duration unit="%">{:d}</Minimum_Effective_Duration>
        <Earliest_Observation_Start unit="phase">{:0.4f}</Earliest_Observation_Start>
        <Latest_Observation_Start unit="phase">{:0.4f}</Latest_Observation_Start>
        <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"?>
<!--                                                               -->
<!--  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>
  
  <Data_Block type="xml">
    <!-- NON-TIME-CRITICAL REQUEST  -->
    <List_of_Non_Time_Critical_Requests count="1">
      <Non_Time_Critical_Request>
        <Programme_ID>{:d}</Programme_ID>
        <Observation_Request_ID>1</Observation_Request_ID>
        <Observation_Category>non time critical</Observation_Category>
        <Proprietary_Period_First_Visit unit="days">{:d}</Proprietary_Period_First_Visit>
        <Proprietary_Period_Last_Visit unit="days">{:d}</Proprietary_Period_Last_Visit>
        <Target_Name>{}</Target_Name>
        <Gaia_ID>GAIA DR2 {}</Gaia_ID>
        <Spectral_Type>{}</Spectral_Type>
        <Target_Magnitude unit="mag">{:0.2f}</Target_Magnitude>
        <Target_Magnitude_Error unit="mag">{:0.2f}</Target_Magnitude_Error>
        <Readout_Mode>{}</Readout_Mode>
        <Right_Ascension unit="deg">{:0.5f}</Right_Ascension>
        <Declination unit="deg">{:0.5f}</Declination>
        <RA_Proper_Motion unit="mas/year">{:0.2f}</RA_Proper_Motion>
        <DEC_Proper_Motion unit="mas/year">{:0.2f}</DEC_Proper_Motion>
        <Parallax unit="mas">{:0.2f}</Parallax>
        <T_Eff unit="Kelvin">{:0.0f}</T_Eff>
        <Extinction unit="mag">0.00</Extinction>
        <Earliest_Start unit="BJD">{:0.3f}</Earliest_Start>
        <Latest_End unit="BJD">{:0.3f}</Latest_End>
        <Exposure_Time unit="sec">{:0.2f}</Exposure_Time>
        <Number_Stacked_Images>{:d}</Number_Stacked_Images>
        <Number_Stacked_Imagettes>{:d}</Number_Stacked_Imagettes>
        <Visit_Duration unit="sec">{:0.1f}</Visit_Duration>
        <Number_of_Visits>{:d}</Number_of_Visits>
        <Continuous_Visits>false</Continuous_Visits>   <!--  Irrelevant for nominal science observations  -->
        <Priority>{:d}</Priority>
        <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()
    _ = StringIO()
    with redirect_stdout(_), redirect_stderr(_):
        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.15
    try:
        Gmag = float(row['Gmag'])
    except ValueError:
        raise ValueError('Invalid Gmag value ',row['Gmag'])
    except KeyError:
        Gmag = row['Vmag'] + GV

    if abs(Gmag-DR2Table['phot_g_mean_mag'][idx]) > gaia_mag_tolerance:
        if 'Gmag' in row.colnames:
            print("Input value: G = ", Gmag)
        else:
            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, 
        checker=False, user_g_mag=False):

    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 user_g_mag:
        gmag = row['Gmag']
        e_gmag = row['e_Gmag']
    else:
        gmag = row['dr2_g_mag']
        e_gmag = row['e_dr2_g_mag']

    if checker:
        try:
            key = re.match('[AFGKM][0-9]', row['SpTy'])[0]
            GV = SpTypeToGminusV[key]
        except TypeError:
            GV = -0.15
        mag, e_mag = gmag-GV, e_gmag
    else:
        mag, e_mag = gmag, e_gmag

    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'],
              mag, e_mag,
              _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'],
              mag, e_mag,
              _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']
              )

    if checker:
        xml = xml.replace('Target_Magnitude','Target_Vmagnitude')

    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 provided in the input table, or estimated
        from Vmag and SpTy from the same table if Gmag is not given. 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    Teff/K
          ------------------
          A0   -0.019   9700
          F5   -0.109   6510
          G5   -0.181   5660
          K0   -0.229   5280
          M0   -0.997   3870
          M9   -3.337   2400

        N.B. An estimate of the spectral type is needed in any case for
        observation requests because accurate flat-fielding requires an
        estimate of the star's effective temperature.

        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])
         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)

         In addition, the input table must specify either ...
         Gmag       - G-band magnitude
         e_Gmag     - error in g-band magnitude
         ... or ...
         Vmag       - V-band magnitude
         e_Vmag     - error in V-band magnitude

        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 **

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

         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.

        e-/s - The count rate in e-/s based on the star's Gaia G magnitude and
               spectral type. This value returned is suitable for use in the
               CHEOPS exposure time calculator using the option "Expected flux
               in CHEOPS passband". 

        duty - duty cycle (%)

        frac - maximim counts as a fraction of the full-well capacity (%)

        img  - image stacking order

        igt  - imaggete stacking order

         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.15
             +  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

        The exposure time can be calculated automatically based on the G-band
        magnitude and spectral type of the target. The default behaviour is to
        use an exposure time that will give 85% of the full-well capacity at
        the peak of PSF, up to the maximum allowed exposure time of 60s.  This
        percentage can be adjusted using the option --scaling-factor-percent. 

        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('-u', '--use_gaia_mag_from_table', 
        action='store_const',
        dest='user_g_mag',
        const=True,
        default=False,
        help='''
        Use Gaia magnitude from the input table instead of Gaia DR2 value for
        calculation and in the output XML file.
        '''
    )

    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='Calculate exposure time automatically'
    )

    parser.add_argument('-s', '--scaling-factor-percent', 
        default=85., type=float,
        help='Scaling factor for auto-expose calculation'
    )

    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('-c', '--checker', 
        action='store_const',
        dest='checker',
        const=True,
        default=False,
        help='Output XML suitable for use with Scheduling Feasibility Checker'
    )

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

    parser.add_argument('-d', '--directory', 
        default='.', type=str,
        help='''
        Output directory for xml files 
        (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)
        for file_name in src_files:
            full_file_name = join(src, file_name)
            if (isfile(full_file_name)):
                copy(full_file_name, getcwd())
        print("Copied examples files from {}".format(src))
        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.0
            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', 'dr2_g_mag', 'e_dr2_g_mag'):
        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

    # Ensure RA and Dec columns are wide enough to accept updates RA/Dec
    # values from Gaia catalogue to full precision
    table['_RAJ2000'] =  [s.rjust(11) for s in table['_RAJ2000']]
    table['_DEJ2000'] =  [s.rjust(11) for s in table['_DEJ2000']]

    # Load contamination function from pickle
    config = load_config()
    cache_path = config['DEFAULT']['data_cache_path']
    pfile = join(cache_path,'Contamination_33arcsec_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))
    if args.auto_expose:
        print('# Exposure time scaling factor: {:0.1f} %'.
                format(args.scaling_factor_percent))
    else:
        print('# Exposure time from input file')
    
    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}'+
        '{:9.2e} {:4.0f}  {:3.0f} {:3d} {:3d}')
    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('#')

    tstr = 'Gaia_DR2_ID         Gmag  _RAJ2000 _DEJ2000  Contam Vis Texp Flags'
    tstr += ' e-/s     frac duty img igt'

    print('#{}'.format(ObsReqNameHeader) + tstr.format(ObsReqNameHeader))

    # String of coordinates, Vmag/Gmag 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))

        if 'Gmag' in row.colnames:
            tag = "{}, {}, {}".format(coo.to_string(), row['Gmag'], row['SpTy'])
        else:
            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']
        row['dr2_g_mag'] = DR2data['phot_g_mean_mag']
        row['e_dr2_g_mag'] = 1.086/DR2data['phot_g_mean_flux_over_error']

        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

        _T = row['T_eff']
        if args.user_g_mag:
            _G = row['Gmag']
        else:
            _G = DR2data['phot_g_mean_mag']
        if args.auto_expose:
            row['T_exp'] = exposure_time(_G, _T, 
                    frac=args.scaling_factor_percent/100)
        if row['T_exp'] >60:
            raise ValueError("Maximum exposure time 60 s exceeded")
        img, igt, cad, duty, frac = cadence(row['T_exp'], _G, _T)

        if frac > 0.95:
            flags += 128
        if frac < 0.1:
            flags += 64

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

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

        c_tot, c_av, c_max = count_rate(_G, row['T_exp'])
        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, c_tot, 100*frac, duty, img, igt))

