# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/Mariagca/EscapeEarth/blob/main/Interns/Maria/TICID_190885165.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + colab={"base_uri": "https://localhost:8080/"} id="YPdxitvrCxQc" outputId="41ef1844-a45e-423a-8371-cd88f540a73b"
## typical imports
# %matplotlib inline
# !pip install lightkurve==1.9.0 #b/c non-standard library & we want a specific version
import lightkurve as lk
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# + colab={"base_uri": "https://localhost:8080/", "height": 99} id="ygBs5I04C4uP" outputId="91d9eaa6-6c51-488c-f1d5-516c3db779f4"
#find available data
search_result = lk.search_lightcurvefile('TIC 190885165')
search_result

# + colab={"base_uri": "https://localhost:8080/"} id="WTu7gSoWDACj" outputId="93f286fb-5d4b-4879-8e92-d28c6770a701"
#download all data available
lc_files = search_result.download_all()
lc_files

# + colab={"base_uri": "https://localhost:8080/", "height": 387} id="Dqlfa5MnDImj" outputId="f72fa8d4-a361-4936-8568-347203a5f66e"
#clean & visualize data

lc = lc_files[0].PDCSAP_FLUX.normalize() #select first file, select detrended flux, normalize #NOTE: if only one file available try removing [0] b/c it will cause an error
lcflat=lc.flatten() #flatten rotational modulations
lc_gtg = lcflat.remove_nans() #remove nans from data gaps
#NOTE: we skipped the remove outliers step to avoid cutting off transit bottoms 

lc.scatter(); #plot light curve before flattening (ie theres still rotational modulations in the data)


# + id="vYIrgNI1Da2h"
#functions from our LCAnalysis

def periods(N=1000):
    period=np.logspace(-0.523, 1.43, N, endpoint=True)
    return period

def duration_grid(N=3):
    duration=np.linspace(.01, 0.298, N)
    return duration
    
def BLS(periodgrid,lightcurve,flat_time,durationgrid):
    from astropy.timeseries import BoxLeastSquares

    '''
    Purppose
    ------------------
    A Box Least Squares function to print out the compute stats of the periodogram.
    Parameters
    -------------------
    period grid            - describes how often the transit is happening (arrays different value)
    duration grid          - describes the width of the transit (array of different values)
    lightcurve             - lightkurve class object
    Returns
    list of stats in the following order: period, duration, transit-time, power, depth
    ------------------
    Calculate several statistics of a candidate transit.
    '''
    #assigning parameters to variables
    period = periodgrid
    duration = durationgrid
    lc = lightcurve

    t = flat_time #time
    y = lc.flux #flux

    #dy is the uncertianty
    model = BoxLeastSquares(t, y, dy= lc.flux_err)
    periodogram = model.power(period,duration)
    max_power = np.argmax(periodogram.power)

    #calculates the max stats w/in the transit
    stats = [periodogram.period[max_power],
            periodogram.duration[max_power],
            periodogram.transit_time[max_power],
            max_power,
            periodogram.depth[max_power],
            periodogram.period,periodogram.power,periodogram.transit_time]

    #stats is the one peak, periodogram is the array
    return stats




# + colab={"base_uri": "https://localhost:8080/"} id="aWB63DqHCn8N" outputId="5f84e81f-e163-4136-8614-e2cd23a148ab"
#run BLS

##to use defaults
#pg = periods()
#dg = duration_grid()

## OR customize period & duration grids
pg = np.linspace(1.2,2, 10000) 
 # change 1st & 2nd values to change startpoint & endpoint of list of periods to check

dg = np.linspace(.03,0.04, 10000) # change same as for period grid

#run BLS
bls_output = BLS(pg, lc_gtg, lc_gtg.time, dg)
bls_output[0:5] #period, duration, transit_time, power, depth  <---meaning of order of output values

# + colab={"base_uri": "https://localhost:8080/", "height": 350} id="gZ_0KviYUaU9" outputId="413a9720-7164-4d64-be28-24b5557a7960"
#plot the periodogram

plt.figure(figsize=(10,5))
plt.plot(bls_output[5],bls_output[6]) #x-axis = periods checked, y-axis = strength of model fit at period checked
plt.xlabel('Period [day]')
plt.ylabel('Power') 
plt.title('BLS Periodogram');


#highest peak is the best fit orbital period, any other peaks should appear at multiples of the best peak OR could be more planets!

# + colab={"base_uri": "https://localhost:8080/", "height": 408} id="2lYbxpH1FVn4" outputId="6e70f9c1-7e45-4559-d309-896621e6d231"
        #plot folded light curve

#fold lc
lcfold = lc_gtg.fold(bls_output[0],t0=bls_output[2]) #inputs are best fit model's (period value, transit_time) 
#the transit time helps center the transit in the folded plot

lcfold.scatter() #actual plot


# + id="fJG5k-7EF7M8"
from astroquery.mast import Catalogs
catalog_data = Catalogs.query_criteria(catalog='Tic',ID=190885165)
catalog_df = catalog_data.to_pandas() 

# + id="JJLtmm6VLpZg"
catalog_df['Depth']= bls_output[4]

# + colab={"base_uri": "https://localhost:8080/", "height": 146} id="GE7GUd21LpfT" outputId="b99ca1db-dbb4-44aa-ccf8-8dff516bb340"
catalog_df


# + id="08uEtefqLx4J"
def planet_radii(depth, star_radius):
    import math as m
    R_Sun = 696340000 #m
    R_Earth = 6371000 #m
    r_star = star_radius * R_Sun #b/c in R_Sun units
    radius = r_star * m.sqrt(depth) 
    r_planet = radius / R_Earth #b/c need in R_Earth
    return r_planet


# + colab={"base_uri": "https://localhost:8080/"} id="j4Qp6cenL5JB" outputId="bbdb68f6-0744-4a30-fb51-d6a8152ab288"
PLANET = catalog_df['Depth'].to_numpy() #transit's depth
STAR = catalog_df['rad'].to_numpy() #star's radius
NAMES = catalog_df['ID'].to_numpy() #star id

radii = [] #empty list to hold results
for i,j in zip(PLANET,STAR):
  starr = planet_radii(i, j) #apply function
  radii.append(starr) #capture all results
  print("Estimared exoplanet radius:", starr, "Earth radii") #print out results

StarcoderdataJupyterScriptsDedupFiltered