#!/usr/local/bin/python3

# AIJ apertures for Gaia field around TIC targets
# Uses MAST for the Gaia DR2 database
# 2020-07-20: Version 1.11
# 2020-09-30: Version 1.12 corrected PPT to milli-magnitude conversion
# 2021-05-28: Version 2.0  without astropy dependence
# 2021-06-14: Version 2.2  synchronized to tk desktop version
# 2021-06-16: Version 2.3  revised to  match tk desktop version 2.3


# Optional: for debugging uncomment the following line

import cgitb
cgitb.enable() 

import cgi                               # for CGI
import os                                # for system environment
import sys	                         # for system connections
import numpy as np                       # for managing the data
from time import gmtime, strftime, time  # for utc
import re                                # for manipulating strings

import json                        # MAST API uses json data
import requests                    # used to get data from MAST
from urllib.parse import quote as urlencode    #used for request parsing

# Keep the useful quantities global
# Note that here at the top level so it is obvious later
# Dependency on astropy query is localized in one routine 

global date_obs_str 
global ra_str
global dec_str
global src_str
global magnitude_str
global depth_str
global ra_center
global dec_center
global ref_center
global tess_magnitude
global tess_depth
global tess_dmag

date_obs_str = ""
ra_str =  ""
dec_str = ""
src_str = "tic" 
magnitude_str = ""
depth_str =  ""

ra_center = 0.
dec_center = 0.
tess_magnitude = 0.
tess_depth = 0.
tess_dmag = 0.
ref_center = "tic"
search_radius_degrees = 2.5/60.
diagnostics_flag = False


# Check for validity of RA and Dec strings
# Parse coordinate strings into degrees for center of search

def parse_coords():
  
  global ra_str
  global dec_str
  global ra_center
  global dec_center
  
  rahr = 0.
  ramin = 0.
  rasec = 0.
  decdeg = 0.
  decmin = 0.
  decsec = 0.

  if (ra_str == ""):
    print ("Content-Type: text/html")
    print ("")

    print ("<br />")
    print ("<b>Missing a right ascension</b>")
    print ("<br />")
    print ("<br />")
    exit()

  if (dec_str == ""):
    print ("Content-Type: text/html")
    print ("")

    print ("<br />")
    print ("<b>Missing a declination</b>")
    print ("<br />")
    print ("<br />")
    exit()



  if ':' in ra_str:
    rahrstr, raminstr, rasecstr = ra_str.split(":")
    rahr = abs(float(rahrstr))
    ramin = abs(float(raminstr))
    rasec = abs(float(rasecstr))
    if float(rahrstr) < 0:
      rasign = -1.
    else:
      rasign = +1.
    ra_center = rasign*(rahr + ramin/60. + rasec/3600.)*15.
  else:  
    try:
      ra_center = 15.*float(ra_str)
    except:
      print ("Content-Type: text/html")
      print ("")

      print ("<br />")
      print ("<b>Cannot parse right ascension entry.</b>")
      print ("<br />")
      print ("<br />")
      exit()    
       

  if ':' in dec_str:
    decdegstr, decminstr, decsecstr = dec_str.split(":")
    decdeg = abs(float(decdegstr))
    decmin = abs(float(decminstr))
    decsec = abs(float(decsecstr))
    if float(decdegstr) < 0:
      decsign = -1.
    elif decdegstr[0] == "-":
      decsign = -1.  
    else:
      decsign = +1.
    dec_center = decsign*(decdeg + decmin/60. + decsec/3600.)
  else:  
    try:
      dec_center = float(dec_str) 
    except:
      print ("Content-Type: text/html")
      print ("")

      print ("<br />")
      print ("<b>Cannot parse declination entry.</b>")
      print ("<br />")
      print ("<br />")
      exit()
  return()          


# Parse the source of the coordinates that are submitted in the request

def parse_source():
  global src_str
  global ref_center

# Source of coordinates may be tic (default), gaia, or eod  

  if src_str == "tic":
    ref_center = "tic"
  elif src_str == "gaia":
    ref_center = "gaia"
  elif src_str == "eod":
    ref_center = "eod"
  else:
    ref_center = "tic"
  return()
          

# Parse and assign search magnitudes

def parse_magnitudes():

  global tess_magnitude
  global tess_depth
  global tess_dmag
  
  # Parse and clip the magnitude of the target
  tess_magnitude = np.clip(float(magnitude_str), 0., 20.)
    
  # Parse and clip the depth string in PPT to relative flux
  tess_depth = np.clip(float(depth_str)/1000., 0.0001, 0.9999)
  
  # Convert the depth in PPT to search delta magnitude
  # Convert ppm to ppt use 2.51188643150958*np.log10
  tess_dmag = -2.5*np.log10(tess_depth)
  
  return()

# JD routines are in ERFA (IAU), Astropy, and Skyfield
# See for example
# https://github.com/liberfa/erfa/blob/master/src/jd2cal.c
# https://github.com/liberfa/erfa/blob/master/src/cal2jd.c
# The following is adapted from Meeus Astronomical Algorithms

# Return the Julian date for a given year, month, day, and UTC

def jd(y0,m0,d0,u0):
  k = int(y0)
  m = int(m0)
  i = int(d0)
  utc = u0
  j1 = float(367*k)
  j2 = -1.*float( int( ( 7*(k + int((m+9)/12) ) )/4 ) )
  j3 = float(i + int(275*m/9))
  j4 = 1721013.5 + (utc/24.0)
  julian = j1 + j2 + j3 + j4
  return julian


# Return the Julian date for a given floating point epoch in years
# Accuracy:
#   The routine takes the difference between two large numbers
#   Python floating point has 18 digits of precision
#   JD day count uses 7 digits
#   JD UTC differencing for 1 second would be 1 part out of 86400
#   There is a 5 digit buffer at this level of time precision
#     when calling jd() twice

def jd_for_epoch(epoch):  
  int_year = int(epoch)
  fraction_of_year = epoch - float(int_year)
  jd_start_of_year = jd(int_year, 1, 1, 0.) 
  jd_start_of_next_year = jd(int_year + 1, 1, 1, 0.) 
  days_this_year = jd_start_of_next_year - jd_start_of_year
  julian_day_for_epoch  = jd_start_of_year + fraction_of_year*days_this_year
  return julian_day_for_epoch



# Format degrees or hours into a hexagesimal  +/-dd:mm:ss.sss string
# This routine is used to format the output for both RA and Dec
# If the source is RA in degrees then divide by 15 first and the output
#   will be RA in hh:mm:ss.sss 
# Or use ra_to_hms_str(degrees) copied below

def float_to_dms_str(degrees):
  if degrees > 0:
    dsign = '+'
  elif degrees == 0.0:
    dsign = ' '
  else:
    dsign = '-'
  least_count = (1.0/3600.0)/1000.0
  degrees = abs(degrees) + least_count   
  dd = int(degrees)
  subdegrees = abs( degrees - float(dd) )
  minutes = subdegrees * 60.0
  mm = int(minutes)
  subminutes = abs( minutes - float(mm))
  seconds = subminutes * 60.0
  ss = int(seconds)
  subseconds = abs( seconds - float(ss) )
  milliseconds = subseconds*1000.0
  ms = int(milliseconds)

  anglestr = "%s%02d:%02d:%02d.%03d" % (dsign, dd, mm, ss, ms) 
  return anglestr



# Format floating point Dec in degrees into a +/-dd:mm:ss.sss Dec string
# Included for reference but not used

def dec_to_dms_str(degrees):
  if degrees > 0:
    dsign = '+'
  elif degrees == 0.0:
    dsign = ' '
  else:
    dsign = '-'
  least_count = (1.0/3600.0)/1000.0
  degrees = abs(degrees) + least_count   
  dd = int(degrees)
  subdegrees = abs( degrees - float(dd) )
  minutes = subdegrees * 60.0
  mm = int(minutes)
  subminutes = abs( minutes - float(mm))
  seconds = subminutes * 60.0
  ss = int(seconds)
  subseconds = abs( seconds - float(ss) )
  milliseconds = subseconds*1000.0
  ms = int(milliseconds)

  anglestr = "%s%02d:%02d:%02d.%03d" % (dsign, dd, mm, ss, ms) 
  return anglestr


# Format floating point RA in degrees into a +/-hh:mm:ss.sss RA string
# Included for reference but not used

def ra_to_hms_str(degrees):
  hours = degrees/15.0
  if hours > 0:
    hsign = '+'
  elif hours == 0.0:
    hsign = ' '
  else:
    hsign = '-'
  least_count = (1.0/3600.0)/10000.0
  hours = abs(hours) + least_count   
  hh = int(hours)
  subhours = abs( hours - float(hh) )
  minutes = subhours * 60.0
  mm = int(minutes)
  subminutes = abs( minutes - float(mm))
  seconds = subminutes * 60.0
  ss = int(seconds)
  subseconds = abs( seconds - float(ss) )
  milliseconds = subseconds*1000.0
  ms = int(milliseconds)

  hourstr = "%s%02d:%02d:%02d.%03d" % (hsign, hh, mm, ss, ms) 
  return hourstr



# Send a request to the MAST server and return the response

def mast_query(request):
  """
    Perform a MAST query.
  
    Parameters
    ----------
    request (dictionary): The MAST request json object
    
    Returns head,content where head is the response HTTP headers, 
    and content is the returned data
  """
  
  # Base API url
  request_url='https://mast.stsci.edu/api/v0/invoke'    
  
  # Grab Python Version 
  version = ".".join(map(str, sys.version_info[:3]))
  
  # Create Http Header Variables
  headers = {"Content-type": "application/x-www-form-urlencoded",
    "Accept": "text/plain",
    "User-agent":"python-requests/"+version}
  
  # Encoding the request as a json string
  req_string = json.dumps(request)
  req_string = urlencode(req_string)
  
  # Perform the HTTP request
  resp = requests.post(request_url, data="request="+req_string, headers=headers)
  
  # Pull out the headers and response content
  head = resp.headers
  content = resp.content.decode('utf-8')
  
  return head, content


# Perform a get request to download a specified file from the MAST server

def download_request(payload, filename, download_type="file"):
  request_url='https://mast.stsci.edu/api/v0.1/Download/' + download_type
  resp = requests.post(request_url, data=payload)

  with open(filename,'wb') as FLE:
    FLE.write(resp.content)

  return filename


# Query the Gaia catalog through MAST

def query_mast_api_gaia(query_ra, query_dec, query_radius):

  # Use MAST tools for the inquiry 
  # We search around the input coordinates assuming ICRS frame
  # Search is without reqard to proper motion since we do not know it yet
  # Choosing a pagesize large enough to accommodate all of the results
  #   or choose a smaller pagesize and 
  #   page through them using the page property
  # Return the search results as a json object
   
  request = {
    'service':'Mast.Catalogs.GaiaDR3.Cone',
    'params':{'ra':query_ra, 'dec':query_dec, 'radius':query_radius},
    'format':'json',
    'pagesize':30000,
    'page':1}   
   
  headers, out_string = mast_query(request)
  gaia_json = json.loads(out_string)

  return gaia_json


# Query the TIC through MAST
# Function retained here for reference

def query_mast_api_tic(query_ra, query_dec, query_radius):

  # Use MAST tools for the inquiry 
  # We search around the input coordinates assuming ICRS frame
  # Search is without reqard to proper motion since we do not know it yet
  # Choosing a pagesize large enough to accommodate all of the results
  #   or choose a smaller pagesize and 
  #   page through them using the page property
  # Return the search results as a json object
   
  request = {
    'service':'Mast.Catalogs.TIC.Cone',
    'params':{'ra':query_ra, 'dec':query_dec, 'radius':query_radius},
    'format':'json',
    'pagesize':30000,
    'page':1}   
   
  headers, out_string = mast_query(request)
  tic_json = json.loads(out_string)

  return tic_json


# Procedure to query Gaia and generate the star list
# Write a temporary file
# Return the file name to use in responding to the request

def make_star_list():
  
  if diagnostics_flag:
    print("Request inputs -- \n <br><br>")
    print("Source of target coordinates: ", src_str, "<br>")
    print("Magnitude requested: ", magnitude_str, "<br>")
    print("Transit depth requested: ", depth_str, "<br>")
  
  
  date_obs_parts = date_obs_str.split('-')
  yr = date_obs_parts[0]
  mo = date_obs_parts[1]
  da = date_obs_parts[2]
  jd_obs = jd(yr, mo, da, 0.)

  # Convert input text to values for processing
  
  # Parse source of input coordinates
  parse_source()
  
  # Parse coordinates based on the style of the input
  parse_coords()
  
  # Parse magnitude from input
  parse_magnitudes()

  # Set the elapsed time in years from the epoch of the Gaia data to the epoch of observation
     
  delta_t = (jd_obs - jd(2016,1,1,0.5)) / 365.25

  # Set the elapsed time in years from the epoch of the Gaia data to the epoch of center
  # This will be negative if the center epoch is before 2016 and positive after 2016

  if ref_center  == "gaia":
    # Gaia EDR3 is in BCRS at 2000.0 with proper motion included to to 2016.0
    delta_t_center = 0.0 
  elif ref_center == "tic":
    # TIC is in ICRS coordinates for 2000.0
    delta_t_center = (jd(2000,1,1,0.5) - jd(2016,1,1,0.5)) / 365.25    
  elif ref_center == "eod":  
    delta_t_center = (jd_obs - jd(2016,1,1,0.5)) / 365.25       
  else:
    delta_t_center =  (jd(2000,1,1,0.5) - jd(2016,1,1,0.5)) / 365.25    

 
  # Create a sky objects counter and an empty buffer (for sorting)

  nsky = 0
  sky = []

  # Make the request and return a json file

  gaia_json = query_mast_api_gaia(ra_center, dec_center, search_radius_degrees)
  
  # Separate the nested parts
  
  nstars =  gaia_json['paging']['rowsTotal']
  labels = gaia_json['fields']
  stars = gaia_json['data']

  if nstars < 1:
    sys.exit('No objects found at  %s  %s ' % (ra_str,dec_str,))
  
  # Set up proper motion corrections from Gaia values in milli-arcseconds per year
  
  pm_scale = 1. / (1000.)

  # Add gaia stars to the sky list
  #
  # Do this star by star to make a new list that will be sorted after all stars are processed
  # Stars in the list include proper motion set to the date of observation
  #
  # The list will be sorted by how far each star would be from the field
  #   center at the time for which the field center is determined

  for i in range(nstars):

    # Coordinates are ICRS epoch 2000.0

    if isinstance(stars[i]['ra'], type(None)):
      ra_gaia = 0.0
    else:
      ra_gaia = float(stars[i]['ra'])

    if isinstance(stars[i]['dec'], type(None)):
      dec_gaia = 0.0
    else:
      dec_gaia = float(stars[i]['dec'])

    if isinstance(stars[i]['source_id'], type(None)):
      src_id = ' '
    else:
      src_id = stars[i]['source_id']  

    if isinstance(stars[i]['phot_bp_mean_mag'], type(None)):
      mag_b = 99.0
    else:
      mag_b = float(stars[i]['phot_bp_mean_mag'])
       
    if isinstance(stars[i]['phot_g_mean_mag'], type(None)):
      mag_g = 99.0
    else:
      mag_g = float(stars[i]['phot_g_mean_mag'])

    if isinstance(stars[i]['phot_rp_mean_mag'], type(None)):
      mag_r = 99.0
    else:
      mag_r = float(stars[i]['phot_rp_mean_mag'])

    if isinstance(stars[i]['bp_rp'], type(None)):
      color_bmr = 99.0
    else:
      color_bmr = float(stars[i]['bp_rp'])

    if isinstance(stars[i]['pmra'], type(None)):
      pm_ra = 0.0
    else:
      pm_ra = float(stars[i]['pmra'])

    if isinstance(stars[i]['pmdec'], type(None)):
      pm_dec = 0.0
    else:
      pm_dec = float(stars[i]['pmdec'])
      
    if isinstance(stars[i]['parallax'], type(None)):
      parallax = 0.0
    else:
      parallax = float(stars[i]['parallax'])
       
    if isinstance(stars[i]['distance'], type(None)):
      distance = 0.0
    else:
      distance = float(stars[i]['distance'])

    if isinstance(stars[i]['ref_epoch'], type(None)):
      ref_epoch = 2016.00
    else:
      ref_epoch = stars[i]['ref_epoch']  
               
    # Assign an effective TESS magnitude based on Gaia DR2 calibration in the TIC    
    # Check if the color b-r exists first

    if (color_bmr < 99.0):
      c1 = color_bmr
      c2 = c1*c1
      c3 = c2*c1
      star_magnitude = mag_g - 0.00522555*c3 + 0.0891337*c2 - 0.633923*c1 + 0.0324473
    else:   
      star_magnitude = mag_g - 0.430          

    # Select targets based on magnitude
    # Reject stars that are too faint to create the observed transit depth

    selected_flag = (star_magnitude < (tess_magnitude + tess_dmag + 0.5))
    
    # Other selection criteria would go here
 

    # Build the ordered output database for the selected stars          

    if selected_flag: 
                  
      # Use the Gaia entry to find the coordinates at the epoch of observation.       
      # PMs are in milliarcseconds per year for 
      # the angle on the sky in the direction of Dec of RA
    
      delta_ra = pm_ra * delta_t / ( 3600.0*1000.0*np.cos(np.radians(dec_gaia)) )
      delta_dec = pm_dec * delta_t / ( 3600.0*1000.0 )
      ra = ra_gaia + delta_ra
      dec = dec_gaia + delta_dec

      # Where was this star at the epoch of the center coordinates?
      
      delta_ra_center = pm_ra * delta_t_center / ( 3600.0*1000.0*np.cos(np.radians(dec_gaia)) )
      delta_dec_center = pm_dec * delta_t_center / ( 3600.0*1000.0 )
      ra_at_center_epoch = ra_gaia + delta_ra_center
      dec_at_center_epoch = dec_gaia + delta_dec_center

      # What was the angular separation of this star from the center at the center epoch?

      delta_theta_rad = np.arccos(np.sin(np.radians(dec_center))*np.sin(np.radians(dec_at_center_epoch)) + np.cos(np.radians(dec_center))*np.cos(np.radians(dec_at_center_epoch))*np.cos(np.radians(ra_center - ra_at_center_epoch)))              
      delta_theta = abs(3600.*np.degrees(delta_theta_rad))         
              
      # Find magnitude difference (flux ratio) compared to the target
      # Take smallest magnitude ratio to associate this star with the target and assign it as T1
      
      delta_mag_gaia = abs(star_magnitude - tess_magnitude)
                   
      
      # Create a selection ranking for ordering the aperture list from 0.0
      # Within uncertainty_theta it will rank most highly the stars closest to the target 
      # Outside that region it will rank by magnitude

      uncertainty_theta = 8.
      uncertainty_mag_gaia = 1.
      if delta_theta > uncertainty_theta:
        rank = delta_theta   
      else:
        rank = delta_mag_gaia/100.
                       
      # Add this selected star to the sky list with ra, dec, magnitude,  delta mag_tic,  delta_theta, and rank
      
      sky.append((ra, dec, star_magnitude, delta_theta, delta_mag_gaia, rank))
      nsky = nsky + 1
      
           
    
  # Create a numpy array to make sorting easier
  
  sky_stars = np.array(sky)
  
  # Sort by rank to make the  target star first
  
  sky_stars = sky_stars[sky_stars[:,5].argsort()]

  if diagnostics_flag:
    print("<br><br><hr><br><br>")
    print("Found ", nsky, " out of ",  nstars,  " after culling the Gaia field.")
    print("\n")
    print("<br>")
    print("   #    &nbsp;   delta_theta  &nbsp;   delta_mag_gaia  &nbsp;     rank\n")
    print("<br>")

    # List the stars in rank order from which the AIJ apertures will be made
    for i in range(nsky):  
      print("%4d    &nbsp;   %8.3f   &nbsp;    %8.3f   &nbsp;   %8.3f" % (i, sky_stars[i,3], sky_stars[i,4], sky_stars[i,5]))
      print("<br \>")
  
  # Copy the first two columns with ra and dec from columns "0" and "1" of the sorted array
  # Create a new array for the ra, dec, and magnitude list
  
  sky_coord = sky_stars[:,0:3]  
  
  # Create a unique name for temporary storage
  gaia_full_path = "/data1/www/htdocs/gaia_to_aij/tmp/"
  gaia_web_path = "/gaia_to_aij/tmp/"
  gaia_stars_file = "gaia_stars_"+str(int(time()))+".radec"
    
  # Open the file
  gaia_stars_fp = open((gaia_full_path+gaia_stars_file), 'w')
  
  # Write data line by line
  # RA in sexagesimal hours
  # Dec in sexagesimal degrees
  # Gaia magnitude based on selection 
  
  # Formatted for AstroImageJ
  # RA in decimal or sexagesimal HOURS
  # Dec in decimal or sexagesimal DEGREES
  # Ref Star=0,1,missing (0=target star, 1=ref star, missing->first ap=target, others=ref)
  # Centroid=0,1,missing (0=do not centroid, 1=centroid, missing=centroid)
  # Apparent Magnitude or missing (value = apparent magnitude, or value > 99 or missing = no mag info)
  # One comma separated line per aperture in the following format:
  #   RA, Dec, Ref Star, Centroid, Magnitude

  for i in range(nsky):
        
    ra_star = (sky_coord[i,0])/15.
    dec_star = sky_coord[i,1]
    mag_star = sky_coord[i,2]
    ra_star_str = float_to_dms_str(ra_star)
    dec_star_str = float_to_dms_str(dec_star)
    if i == 0:
      skyline = "%s,  %s, 0, 1, %6.3f  \n" % (ra_star_str, dec_star_str, mag_star) 
    else:
      skyline = "%s,  %s, 0, 0, %6.3f  \n" % (ra_star_str, dec_star_str, mag_star) 
    gaia_stars_fp.write(skyline)
  gaia_stars_fp.close()
  return((gaia_web_path+gaia_stars_file))


# Parse the cgi fields

request = cgi.FieldStorage()

date_obs_str = request.getfirst("date_obs", "")
ra_str = request.getfirst("ra", "")
dec_str = request.getfirst("dec", "")
magnitude_str = request.getfirst("mag", "")
depth_str  = request.getfirst("depth", "")
src_str = request.getfirst("src", "tic")   

# Remove all but numbers and dashes from date-obs
# Use default date if this field is empty

date_obs_str = re.sub(r'([^-0-9])+', '', date_obs_str)
if date_obs_str =="":
  date_obs_str = "2021-01-01"

# Remove non-numeric characters and spaces from magnitude and depth

magnitude_str = float(re.sub(r'([^-+\.0-9])+', '', magnitude_str))
depth_str = float(re.sub(r'([^-+\.0-9])+', '', depth_str))

# Remove unacceptable characters and trailing spaces ra and dec

ra_str = re.sub(r'([^-+\.\:0-9])+', '', ra_str).strip()
dec_str = re.sub(r'([^-+\.\:0-9])+', '', dec_str).strip()

# Check and respond to missing necessary fields
 
if (depth_str == ""):
  print ("Content-Type: text/html")
  print ("")

  print ("<br />")
  print ("<b>Missing event depth</b>")
  print ("<br />")
  print ("<br />")
  exit()

if (magnitude_str == ""):
  print ("Content-Type: text/html")
  print ("")

  print ("<br />")
  print ("<b>Missing target magnitude </b>")
  print ("<br />")
  print ("<br />")
  exit()


# Now acknowledge the request

print ("Content-Type: text/html")
print ("")

print ("<br />")
print ("<b>Request sent for Gaia stars in the field:</b>")
print ("<br />")
print ("Date = ", date_obs_str)
print ("<br />")
print ("RA = ", ra_str)
print ("<br />")
print ("Dec = ", dec_str)
print ("<br />")
print ("Magnitude = ", magnitude_str)
print ("<br />")
print ("Depth = ", depth_str)
print ("<br />")
print ("<br />")
print ("Processing ... ")
print ("<br />")
print ("<br />")



# Query Gaia and generate the list

output_file = make_star_list()

print ("<br />")
print ("<br />")
outline = "Download: <a href=\""+output_file+"\" download> Gaia stars </a>"
print (outline)


print ("<br />")

exit()