<jupyter_start><jupyter_text>Automated Transient Determination<jupyter_code># This code is intended to be tacked onto the existing reduction pipeline (after differencing).
# At this point we would ideally have an image where the only point sources remaining are transients such as 
# supernovae or asteroids. 
# This code determines which asteroids are present in the FOV at the time of each image ins a series, and logs them.
from __future__ import division, print_function, unicode_literals
import time
import glob
import numpy as np
import math
import os
import astropy.io.fits as fits
import sys
import datetime
import ephem
import matplotlib.pyplot as plt
%matplotlib inline
# paths were developed for my student account, will need updating!
paths = {'reduced':"/Users/student_a8/reduced", \
         'MPCORB': "/Users/student_a8/Huntsman_transients/MPCORB", \
         'difference': "/Users/student_a8/reduced/difference", \
         'temp': "/Users/student_a8/temp", \
         'regions': "/Users/student_a8/regions", \
         'home':"/Users/student_a8"}<jupyter_output><empty_output><jupyter_text>Define Functions:<jupyter_code># Define a function to extract the packed Epoch date format given in the MPCORB data.
# see  <http://www.minorplanetcenter.net/iau/info/PackedDates.html>  for more info.

def epoch_convert(date):
    
    "Converts the MPCORB Epoch from packed form to a regular YYYYMMDD.DDDD string to be split up and used later."
    
    packeddates = {'1':'1', '2':'2', '3':'3', '4':'4', '5':'5', '6':'6', '7':'7', '8':'8', '9':'9', \
                   '0':'0', 'A':'10', 'B':'11', 'C':'12', 'D':'13', 'E':'14', 'F':'15', 'G':'16', \
                   'H':'17', 'I':'18', 'J':'19', 'K':'20', 'L':'21', 'M':'22', 'N':'23', 'O':'24', \
                   'P':'25', 'Q':'26', 'R':'27', 'S':'28','T':'29', 'U':'30', 'V':'31' \
              }
    datestring = list(date[0:5])
    datestring2 = ""
    nums = "0123456789"
    
    # Conditionals to distinguish between packed dates, eg  avoid confusion between 1-Nov (11 1) /
    # and 11-Jan (1 11) when compiled back into a string (111). 
    # We convert all single-digit numeric month/day values to 2-digits (eg 6 --> 06):
    
    if date[3] in nums and date[4] in nums:
        datestring.insert(3,'0')
        datestring.insert(5,'0')
        for i in range(len(datestring)):
            char = packeddates['%s' % datestring[i]]
            datestring2 += str(char)
        if len(date) > 5:
            datestring2 += "." + date[5:]
        return datestring2 
    elif date[3] in nums and date[4] not in nums:
        datestring.insert(3,'0')
        for i in range(len(datestring)):
            char = packeddates['%s' % datestring[i]]
            datestring2 += str(char)
        if len(date) > 5:
            datestring2 += "." + date[5:]
        return datestring2
    elif date[3] not in nums and date[4] in nums:
        datestring.insert(4,'0')
        for i in range(len(datestring)):
            char = packeddates['%s' % datestring[i]]
            datestring2 += str(char)
        if len(date) > 5:
            datestring2 += "." + date[5:]
        return datestring2
    elif date[3] not in nums and date[4] not in nums:
        for i in range(len(datestring)):
            char = packeddates['%s' % datestring[i]]
            datestring2 += str(char)
        if len(date) > 5:
            datestring2 += "." + date[5:]
        return datestring2
# Estimate apparent magnitude under ideal conditions

def asteroid_rough_appmag(H, a, e):
    
    "A function that takes orbital parameters of asteroids and returns a rough apparent magnitude, \
    under ideal conditions."
    
       
    perihelion = float(a)*(1.0-float(e))
    appmag = float(H) + 2.5*math.log10(perihelion**2 * (perihelion - 1.0)**2)
    return float(appmag)
# phase integral approximation. Calculation of phase angle occurs within.

def phase_int(asteroid, date):
    
    "Returns the phase integral approximation for an asteroid at a specific date and time"
    
    #calculate phase angle:
    asteroid.compute(date)
    sun = ephem.Sun(date)
    phi=ephem.separation(asteroid, sun)
    phase_ang = np.arcsin((sun.earth_distance*np.sin(phi))/asteroid.sun_distance)
    
    #feed phase angle into phase integral equation:
    phase_int = ((2/3)*((1-(phase_ang/np.pi))*np.cos(phase_ang)+1/np.pi*(np.sin(phase_ang))))
    return phase_int
# Calculate the apparent magnitude more accurately of an asteroid at a spoecified time.

def asteroid_fine_appmag(body, H, phase_int):
    
    "Returns the apparent magnitude of an asteroid as calculated on a specific date. \
    Takes 3 arguments: pyephem asteroid, that asteroids H magnitude, and that asteroids phase intergral."
    
    asteroid_appmag = H+2.5*np.log10(((body.sun_distance**2)*(body.earth_distance**2))/phase_int)
    return asteroid_appmag<jupyter_output><empty_output><jupyter_text>Import MPCORB database and split into a table:<jupyter_code># Download latest MPCORB catalogue from:
# www.minorplanetcenter.org/iau/MPCORB/MPCORB.DAT

# Full list of paramters and how to read the data:
# http://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html
os.chdir(paths['MPCORB'])

# 'skip_header=41'  skip headers info, delimiter splits  lines into elements of specific character lengths.
MPCORB = np.genfromtxt('MPCORB.DAT', autostrip=True, skip_header=41, dtype=str,  delimiter = [8,6,6,6,10,11,11,11,11,12,12,3,10,6,4,10,5,4,4,11,5,28,8])

designation = MPCORB[:,0]
abs_mag = MPCORB[:,1]
slope_param = MPCORB[:,2]
epoch = MPCORB[:,3]
mean_epoch_anomaly = MPCORB[:,4]
perihelion_arg = MPCORB[:,5]
longitude = MPCORB[:,6]
inclination = MPCORB[:,7]
eccentricity = MPCORB[:,8]
mean_daily_motion = MPCORB[:,9]
semimajor_axis = MPCORB[:,10]
uncertainty = MPCORB[:,11]
reference = MPCORB[:,12]
num_obs = MPCORB[:,13]
num_opp = MPCORB[:,14]
obs_years_arc_length = MPCORB[:,15]
rms = MPCORB[:,16]
coarse_perturb = MPCORB[:,17]
precise_perturb = MPCORB[:,18]
comp_name = MPCORB[:,19]
hex_flag =  MPCORB[:,20]
read_des = MPCORB[:,21]
last_obs = MPCORB[:,22]

print('Number of minor planets:',len(MPCORB))<jupyter_output>Number of minor planets: 706783<jupyter_text>Check series of images against MPCORB:<jupyter_code>os.chdir(paths['difference'])
lightlist = glob.glob('*.resamp.diff.fits')
#lightlist = glob.glob('2014-09-20_83F011167_6*_light.bdfw.resamp.diff.fits')
mag_thresh=20.0
# specify light images to be checked, and the apparent magnitude threshold. Asteroids with ideal app_mag dimmer than
# this threshold will be ignored in ephemeris calculation.
a = time.time()
lights = []
for image in lightlist:
    lights.append([str(image), fits.getheader(image)['DATE'], '00:54:53.5', '-37:41:04'])
    # parse fits header for date/time of exposure. 
    # currently feeding the coordinates of NGC300 manually, since the fits header doesn't store telescope targets.
    
huntsman = ephem.Observer()
huntsman.lon = 149.067307
huntsman.lat = -31.274587
huntsman.elevation = 510    
# details of Coonabarabran. Needs updating when Huntsman installed at SSO.

visiblecount = 0
notviscount = 0
totcount = 0

field = 'NGC300'
# can switch to reading fits header for 'TARGET' name if appropriate.

# Create the text file to store all asteroid information:
image = "Huntsman_asteroids"
f = open("%s.txt" %(image), 'w+')

for image in lights:
    print(image[0])
    transcoords = ephem.FixedBody()   # create fixed object to hold coords of target
    #for future adaptation:   #transcoords = ephem.FixedBody(float(fits.getheader(image)['OBJECTRA']), 
    #   float(fits.getheader(image)['OBJECTDEC']))
    transcoords._ra = image[2]
    transcoords._dec = image[3]
    transcoords.compute()
    datestr = image[1]
    huntsman.date = ephem.Date((int(datestr[0:4]), int(datestr[5:7]), int(datestr[8:10]), \
                                int(datestr[11:13]), int(datestr[14:16]), int(datestr[17:19])))
    print(huntsman.date)       # ^ get the date/time from the fits header
    
    f.write("Image: %s \n" %(image[0]))
    f.write("Time: %s \n" %(huntsman.date))

    identcount = 0 # initialize count of identified objects
    totcount = 0   # initialize count of total number checked
    visible = 0    # initialize count of number of potentially visible asteroids
    identlist = []
    
    # create this images' DS9 region file:    
    os.chdir(paths['regions'])
    fout2=open(field+'_%s'%(image[0][:-28])+'_asteroids.reg','w')
    fout2.write('fk5\n')
    os.chdir(paths['difference'])

    for MP in MPCORB:
        totcount += 1
        if MP[0] == "":  # ignore break (empty) lines in the catalogue
            pass
        
        elif MP[1] == "":  # if asteroid magnitude H missing, apply H=14.0
            if asteroid_rough_appmag(14.0, MP[10], MP[8]) < mag_thresh:
                visible += 1
                newbody = ephem.EllipticalBody()
                # Longitude of ascending node:
                newbody._Om = float(MP[6])
                # Inclination:
                newbody._inc = float(MP[7])
                # Arg of perihelion
                newbody._om = float(MP[5])
                # Mean distance from Sun: 
                newbody._a = float(MP[10])
                # Eccentricity:
                newbody._e = float(MP[8])
                # Mean anomoly from perihelion:
                newbody._M = float(MP[4])
                # Epoch = J2000:
                newbody._epoch = ephem.J2000
                # Epoch for _M:
                aster_epoch_str = epoch_convert(MP[3])  # onvert the epoch of the asteroid from packed form and apply
                aster_year = int(aster_epoch_str[:4])
                aster_month = int(aster_epoch_str[4:6])
                aster_day = float(aster_epoch_str[6:])
                aster_epoch = ephem.Date((aster_year, aster_month, aster_day))  
                newbody._epoch_M = aster_epoch
                
                newbody.compute(huntsman)   # compute asteroid ephemeris as per observer location, time of image
                sep = ephem.separation(newbody, transcoords)
                if sep < ephem.degrees(0.0298):  # if separation less than 100 arcmin threshold from NGC300
                    identcount += 1
                    phase_integral = phase_int(newbody, huntsman.date)
                    asterappmag = asteroid_fine_appmag(newbody, float(params[1]), phase_integral)
                    identlist.append([MP[21], 'a_RA: %s' %(newbody.a_ra), 'a_Dec: %s' %(newbody.a_dec)])
                    identlist.append([MP[21], 'Apparent Magnitude: %s' %(asterappmag)])
                    fout2.write('circle('+str(newbody.a_ra)+','+str(newbody.a_dec)+',15") # color=cyan text={'+str(MP[21])+'_a}\n')

        else:  # asteroid magnitude H present, proceed as above
            if asteroid_rough_appmag(MP[1], MP[10], MP[8]) < mag_thresh:
                visible += 1
                newbody = ephem.EllipticalBody()
                # Longitude of ascending node:
                newbody._Om = float(MP[6])
                # Inclination:
                newbody._inc = float(MP[7])
                # Arg of perihelion
                newbody._om = float(MP[5])
                # Mean distance from Sun: 
                newbody._a = float(MP[10])
                # Eccentricity:
                newbody._e = float(MP[8])
                # Mean anomoly from perihelion:
                newbody._M = float(MP[4])
                # Epoch = J2000:
                newbody._epoch = ephem.J2000
                # Epoch for _M:
                aster_epoch_str = epoch_convert(MP[3])
                aster_year = int(aster_epoch_str[:4])
                aster_month = int(aster_epoch_str[4:6])
                aster_day = float(aster_epoch_str[6:])
                aster_epoch = ephem.Date((aster_year, aster_month, aster_day))
                newbody._epoch_M = aster_epoch
                
                newbody.compute(huntsman)
                sep = ephem.separation(newbody, transcoords)
                if sep < ephem.degrees(0.0298):
                    identcount += 1
                    phase_integral = phase_int(newbody, huntsman.date)
                    asterappmag = asteroid_fine_appmag(newbody, float(MP[1]), phase_integral)
                    identlist.append([MP[21], 'a_RA: %s' %(newbody.a_ra), 'a_Dec: %s' %(newbody.a_dec)])
                    identlist.append([MP[21], 'Apparent Magnitude: %s' %(asterappmag)])
                    fout2.write('circle('+str(newbody.a_ra)+','+str(newbody.a_dec)+',15") # color=cyan text={'+str(MP[21])+'_a}\n')
                    
                                     
    fout2.close()
    print('Difference Image:', image[0])
    
    # if no asteroids in FOV...
    if identcount == 0:
        print('No bodies located within %s degrees of position RA: %s Dec: %s at time %s' %(ephem.degrees(0.0298), transcoords._ra, transcoords._dec, huntsman.date))
        f.write('No bodies located within %s degrees of position RA: %s Dec: %s \n' %(ephem.degrees(0.0298), transcoords._ra, transcoords._dec))
    
    # if some asteroids in FOV...  
    else:
        print('%s bodies located within %s degrees of position RA: %s  Dec: %s at time %s:' %(identcount, ephem.degrees(0.0298), transcoords._ra, transcoords._dec, huntsman.date))
        f.write('%s bodies located within %s degrees of position RA: %s  Dec: %s: \n' %(identcount, ephem.degrees(0.0298), transcoords._ra, transcoords._dec))
        for ident in identlist:
            print(ident)
            f.write("%s \n" %(ident))
    f.write("\n \n")
               
f.close()            
b = time.time()
print('checked %s objects (%s brighter than H=20.0) in %f seconds' %(totcount, visible, b-a))<jupyter_output>2014-09-20_83F011167_100_light.bdfw.resamp.diff.fits
2014/9/20 18:09:13
Difference Image: 2014-09-20_83F011167_100_light.bdfw.resamp.diff.fits
10 bodies located within 1:42:26.7 degrees of position RA: 0:54:53.50  Dec: -37:41:04.0 at time 2014/9/20 18:09:13:
['(142549) 2002 TJ53', u'a_RA: 0:59:15.97', u'a_Dec: -36:51:19.0']
['(142549) 2002 TJ53', u'Apparent Magnitude: 18.5362464823']
['(216692) 2004 HW50', u'a_RA: 0:56:08.00', u'a_Dec: -38:19:19.0']
['(216692) 2004 HW50', u'Apparent Magnitude: 19.211108872']
['(256155) 2006 VO44', u'a_RA: 0:56:27.34', u'a_Dec: -37:09:00.2']
['(256155) 2006 VO44', u'Apparent Magnitude: 16.6036993932']
['2004 BY26', u'a_RA: 0:53:01.47', u'a_Dec: -37:08:58.0']
['2004 BY26', u'Apparent Magnitude: 21.9795685225']
['2013 BH16', u'a_RA: 0:50:50.54', u'a_Dec: -36:30:56.2']
['2013 BH16', u'Apparent Magnitude: 18.778988229']
['2014 QD', u'a_RA: 0:48:25.75', u'a_Dec: -37:39:18.5']
['2014 QD', u'Apparent Magnitude: 18.355927032']
['2010 NP67', u'a_RA: 0:50:12.26',[...]<jupyter_text>(Not part of main code) Use this to calculate day fractionals<jupyter_code>d = ephem.Date((2014, 9, 20, 14, 00, 20))  # Input (YYYY, MM, DD, HH, MM, SS) here
print(d.triple())<jupyter_output>(2014, 9, 20.58356481482042)
{'fraction_non_alphanumeric': 0.13831563480932346, 'fraction_numerical': 0.056368129299405594, 'mean_word_length': 2.7956907477820025, 'pattern_counts': {'":': 2, '<': 23, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '197578', 'license': 'MIT', 'n_tokens_mistral': 5708, 'n_tokens_neox': 5103, 'n_words': 1255, 'path': 'examples/Huntsman_asteroid_detection.ipynb', 'repo_name': 'AnthonyHorton/asteroid-marker', 'size': '14973'}