utils.py

# -*- coding: utf-8 -*-
"""utils.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1zWLknVQuT4SXTo0-6fxRovzjK8tSY1Pr
"""

import numpy as np
# import pylab as plt
import pandas as pd

"ASTROPY"
from astropy.time import Time
from astropy.table import Table
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, get_sun
from astropy.constants import c
import astropy.units as u

# Obtains the Astropy Position Vector of given earth coordinates
pos = lambda x,h: EarthLocation(lon=x['Longitude']*u.deg,lat=x['Latitude']*u.deg,height=h*u.earthRad)

# Transforms an Astropy Vector into a 3-tuple
vec = lambda a:  (a.x.value,
                  a.y.value,
                  a.z.value)

# Obtains the difference between two Astropy Vectors
vecdiff = lambda a,b: (a.x.value-b.x.value,
                       a.y.value-b.y.value,
                       a.z.value-b.z.value)

# Obtains the magnitude of a vector
dist = lambda v: np.linalg.norm(v)

# Calculates the angle between two vectors in degrees
angle = lambda v1,v2: np.rad2deg( np.arccos( np.dot(v1,v2) / (dist(v1) * dist(v2)) ) )

# Divides a distance (with astropy units) between light speed c
light_time = lambda d: (d/c).decompose()

def resol_angle(obs):
    """
    Given an observatory object, it returns the resolution angle AKA beam width.
    """
    if "-" in obs['Frequency Range']:
        print("Antenna Frequency Range:",obs['Frequency Range'])
        freq = float(input("At what frequency will the antenna operate in MHz? "))
    else:
        freq = float(obs['Frequency Range'].split(" ")[0])
    freq *= u.MHz
    lamb = c/(freq)
    diameter = float(obs['Diameter'][:-1]) * u.m
    angle = np.arcsin(1.22 * lamb / diameter)
    return angle

def size_angle(obs,apo):
    """
    Given an observatory and an apophis ephemeris row, it returns the apparent size
    Apophis would have in the sky.
    """
    APOPHIS_SIZE = 450 * u.m # Consider this size!
    pos1 = pos(obs,0)
    posA = pos(apo,apo['delta (Rt)'])
    vec1 = vecdiff(posA,pos1)
    distance = dist(vec1)*u.earthRad
    angle = np.arctan(APOPHIS_SIZE / distance)
    return angle

def reflected_power(obs1,apo,print_info=False):
    """
    Given an emitting antenna and an apophis ephemeris, this function calculates
    the power reflected by apophis considering a radio albedo of 1.
    """
    if obs1['Power'] is not np.nan:      
        power = float(obs1['Power'][:-2]) * u.kW # Power of emitting antenna
        apparent_angle = size_angle(obs1,apo)    # Apparent angle of Apophis
        resolution_angle = resol_angle(obs1)     # Resolution of emitter
        power_at_apophis = power * (apparent_angle / resolution_angle)**2 # Power reflected by Apophis
        if print_info:
            print("\nEMITTING:",obs1['Name'],obs1['Power'])
            print("\nPower reflected by Apophis: {:.2e}".format(power_at_apophis))
        return power_at_apophis
    else:
        print("First Antenna does not emmit!")
        return None

def recieved_power(obs2,apo,reflected,print_info=False):
    """
    Given the reflected power of apophis, it calculates the flux in Janskys that
    arrives at the obs2.
    """
    APOPHIS_SIZE = 450 * u.m
    pos2 = pos(obs2,0)
    posA = pos(apo,apo['delta (Rt)'])
    vec2 = vecdiff(posA,pos2)
    distance = dist(vec2)*u.earthRad.to(u.m) * u.m
    recieved_flux = reflected / (2*np.pi* distance**2)
    diameter = float(obs2['Diameter'][:-1]) * u.m
    # recieved_power = recieved_flux * (np.pi * (diameter/2)**2) * np.cos(angle(vec(pos2),vec(posA)))
    BW = 50 * u.MHz # Assumed bandwith
    recieved_power = recieved_flux * abs(np.cos(angle(vec(pos2),vec(posA)))) / (BW)
    recieved_power = recieved_power.to(u.Jansky)
    if print_info:
        print("Power recieved by {:s}: {:.2e}".format(obs2['Name'],recieved_power))
    return recieved_power

def observatory_pair(obs1,obs2,ephemeris):
    pos1 = pos(obs1,0)
    pos2 = pos(obs2,0)

    # Print the observatories information
    print("\nFIRST OBSERVATORY:",obs1['Name'])
    print("SECOND OBSERVATORY:",obs2['Name'])

    # Initializes the pairing dictionary
    pair = {'date':[],
            'distance':[],
            'elevation1':[],
            'elevation2':[],
            }

    # Iterate over the whole transit
    for r,row in ephemeris.iterrows():
        posA = pos(row,row['delta (Rt)'])
        vec1 = vecdiff(posA,pos1) # Vector from obs1 to Apophis
        vec2 = vecdiff(posA,pos2) # Vector from obs2 to Apophis
        dist1 = dist(vec1)*u.earthRad # Magnitudes
        dist2 = dist(vec2)*u.earthRad
        distance = row['delta (Rt)'] # Length of signal's path

        # Apophis elevation
        elev1 = (90-angle( vec(pos1), vec1 )) #*u.deg
        elev2 = (90-angle( vec(pos2), vec2 )) #*u.deg

        # Flux Density

        pair['date'].append(row['datetime_str'])
        pair['distance'].append( distance )
        pair['elevation1'].append( elev1 )
        pair['elevation2'].append( elev2 )

    pair_df = pd.DataFrame(pair)
    return pair_df

def graph_elevations(obs1,obs2,ephemeris):
    pair = observatory_pair(obs1,obs2,ephemeris)
    ephemeris['hours'] = 24*(ephemeris['datetime_jd']-ephemeris['datetime_jd'][0])


    fig,ax = plt.subplots(figsize=(8,5))
    ax.set_title(f"Apophis Elevations from {obs1['Name']} and {obs2['Name']}")
    ax.plot(ephemeris['hours'],pair['elevation1'],label=f"Elevation from {obs1['Name']}",lw=3,color='b')
    ax.plot(ephemeris['hours'],pair['elevation2'],label=f"Elevation from {obs2['Name']}",lw=3,color='r')
    ax.plot([-100],[-100],'k--',label='Distance',lw=3)
    ax.set_ylabel("Apophis' Elevation (°)",size=14)

    plt.xlabel("2029-04-14 time UT",size=15)
    plt.xticks([0,5,10,15],['-5:00','0:00','5:00','10:00'],size=15)
    plt.yticks(size=12)
    plt.ylim([0,90])
    plt.xlim([0,ephemeris['hours'].iloc[-1]])
    ax.grid()
    plt.legend()
    ax2=ax.twinx()
    ax2.plot(ephemeris['hours'],ephemeris['delta (Rt)']-1,'k--',lw=3)
    ax2.set_ylabel("Distance to Apophis ($R_t$)",size=14)
    # ax.legend(loc=2)
    plt.show()

sign = lambda x: x/abs(x)

def compareSign(prevElevs, actualElevs, time, obsNames, obj):
    for i in range(2):
        prevElev = prevElevs[i]
        actualElev = actualElevs[i]
        obs = obsNames[i]
        if prevElev is not None:
            if sign(prevElev) != sign(actualElev):
                if sign(actualElev) == 1:
                    obj[i]['Rise'].append(time)
                    # print("  Rising at ",obs,time)
                else:
                    obj[i]['Set'].append(time)
                    # print("  Setting at ",obs,time)
    return obj

def table2(obs1, obs2, ephemeris):
    """
    Prints out the important information for table 2 of the paper
    """
    pos1 = pos(obs1,0)
    pos2 = pos(obs2,0)

    obj = {}

    # Print the observatories information
    obj[0] = {"Name":obs1['Name'],"Rise":[],"Set":[]}
    obj[1] = {"Name":obs2['Name'],"Rise":[],"Set":[]}
    obsNames = (obs1['Name'],obs2['Name'])

    # Set the elevations and times tuples
    prevElevs = [None,None]
    maxElevs = [0,0]
    maxTimes = [None,None]

    # Iterate over the whole transit
    for r,row in ephemeris.iterrows():
        time = row['datetime_str']

        posA = pos(row,row['delta (Rt)'])
        vec1 = vecdiff(posA,pos1) # Vector from obs1 to Apophis
        vec2 = vecdiff(posA,pos2) # Vector from obs2 to Apophis

        # Apophis elevation
        elev1 = (90-angle( vec(pos1), vec1 )) #*u.deg
        elev2 = (90-angle( vec(pos2), vec2 )) #*u.deg
        elevs = (elev1,elev2)

        # If elevation changes sign, print if it rises or sets
        obj = compareSign(prevElevs, elevs, time, obsNames, obj)

        # Check for maximum elevation
        for i in range(2):
            if (elevs[i] > maxElevs[i]):
                maxElevs[i] = elevs[i]
                maxTimes[i] = time
        prevElevs = elevs

    # print("Max Elevations:")
    for i in range(2):
        if len(obj[i]['Rise']) == 1:
            obj[i]['Rise'] = obj[i]['Rise'][0]
        if len(obj[i]['Set']) == 1:
            obj[i]['Set'] = obj[i]['Set'][0]
        obj[i]["Max Elevation"] = {"Elevation":f"{round(maxElevs[i],2)} deg", "Time":maxTimes[i]}
    
    return obj