jose-solorzano

BASE_PERIOD = 157.44;

BKJD_BASELINE = 2454833.0; # Kepler mission baseline.
EVANGELINE_WTF_DAY = 404; # Dip day from chart at wherestheflux.com.
WTF_BASELINE = 2457800; # Baseline used in wherestheflux.com charts.

D140 = 140.5437;
D216 = 216.3751;
D260 = 260.89969;
D359 = 359.0791;

jose-solorzano

"""
Runs random dip simulations in the Kepler timespan
and looks for intervals that are approximate multiples
of the 157.44-day base period.

Other utilities are:
    print_best_pairs()
    get_median_intervals_below_threshold()
"""
import numpy as np

jose-solorzano

# Example:
# runIntervalsSimulation(20000, multiplesOf = 24.2, viewingAngle = 7 * 2 * pi / 13)

BASE_GAP = 48.4;

getPeriods = function() {    
    return(c(
        p(3),
        p(4),
        p(5)

jose-solorzano

# Functions for analysis of century-long periodicity of KIC 8462852.

library(lomb)

readDaschData = function(fileName = "./dasch/KIC8462852_K8462852_0000.txt",
                    removeFlagged = F, maxError = NULL, magField = "magcal_local", maxOutlierSdFactor = NULL) {
    con = file(fileName);
    lines = readLines(con);
    close(con);
    lines = lines[c(-1, -3)];

jose-solorzano

# Input file generated with:
# awk -F ',' '{print $1 "," $2 "," $3 "," $4 "," $5}' aavsodata_58adb26376a9d.txt > aavsodata_2017_02_22.csv

generateSimulatedTrends = function(band, ar1 = 0.33, numTrends = 1000, aavsoFilePath = "./aavso/aavsodata_2017_02_22.csv") {
    aavsoFrame = read.csv(aavsoFilePath);    
    # Filter by band and remove "outliers"
    aavsoFrame = aavsoFrame[aavsoFrame$Band==band & aavsoFrame$Magnitude >= 11.724 & aavsoFrame$Magnitude <= 11.964,];
    numRows = nrow(aavsoFrame);
    

jose-solorzano

# AAVSO matching with Kepler

# AAVSO will produce a poorly-formatted comma-separated file 
# named something like aavsodata_58adb26376a9d.txt.
# Convert it to a CSV file as follows:
# awk -F ',' '{print $1 "," $2 "," $3 "," $4 "," $5}' aavsodata_58adb26376a9d.txt > aavsodata_2017_02_22.csv

# Baseline BJD of Kepler mission is 2454833.

analyzeSeriesCorrelationPartial = function(fromRow, toRow, aavsoFilePath = "./aavso/aavsodata_2017_02_22.csv", keplerFilePath = "./wtf/KIC8462852_normalized_flux.csv") {

jose-solorzano

# The kf parameter of functions below should be a frame with
# Time and n_flux variables, like data from http://wherestheflux.com/

# amplitudePhaseAnalysis
# Produces amplitude and phase series by fitting sine/cosine 
# functions to a light curve.
amplitudePhaseAnalysis = function(kf, windowSize = 700, stepSize = 21, period = 0.88) {
    sequence = seq(from = 1, to = nrow(kf) - windowSize + 1, by = stepSize);
    lenseq = length(sequence);

jose-solorzano

seriesWithNoise = function(noiseSD = 0.00004) {
    previous = seriesWithPeriodicity();
    return(data.frame(
        Time = previous$Time,
        n_flux = previous$n_flux + rnorm(n = nrow(previous)) * noiseSD
    ));        
}
    
seriesWithPeriodicity = function(period = 0.88, amplitude = 0.00008) {

jose-solorzano

# Note: A complete AAVSO file should have over 18000 rows at this point. 
# Note: A complete AAVSO file should have over 18000 rows at this point. 
# R might not read it fully.
# Had to do this to AAVSO light curve file for R to read it properly:
# awk -F ',' '{print $1 "," $2 "," $3 "," $4 "," $5}' aavsodata_20161023.csv > aavsodata_20161023_5fields.csv

produceAavsoTrend = function(band, aavsoFilePath = "./aavso/aavsodata_20161023_5fields.csv") {
    aavsoFrame = read.csv(aavsoFilePath);    
    aavsoFrameBand = aavsoFrame[aavsoFrame$Band==band,];
    

jose-solorzano

# Note: A complete AAVSO file should have over 18000 rows at this point. 
# R might not read it fully.
# Had to do this to AAVSO light curve file for R to read it properly:
# awk -F ',' '{print $1 "," $2 "," $3 "," $4 "," $5}' aavsodata_20161023.csv > aavsodata_20161023_5fields.csv

analyzeVvsR = function(aavsoFilePath = "./aavso/aavsodata_20161023_5fields.csv") {
    aavsoFrame = read.csv(aavsoFilePath);    
    aavsoFrameV = aavsoFrame[aavsoFrame$Band=='V',];
    aavsoFrameR = aavsoFrame[aavsoFrame$Band=='R',];