TransitAlignment.R

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;
D426 = 426.3455; 
D502 = 502.4427;
D694 = 694.8865;
D792 = 792.7199;
D1144 = 1144.607;
D1205 = 1205.888;
D1495 = 1495.902;
D1519 = 1519.523;
D1540 = 1540.385;
D1568 = 1568.482;
D1540G = (D1519 + D1568) / 2;
EVANGELINE = EVANGELINE_WTF_DAY + WTF_BASELINE - BKJD_BASELINE;

# Orbital period function for hypothesized chained
# orbital resonance system with base period of 157.44 days.
op = function(n) {
    return(BASE_PERIOD * n);
}

errorAtObserverAngle = function(dipTime1, dipTime2, period1, period2, observerAngleDeg, basePeriod = BASE_PERIOD) {
    interval = intervalAtObserverAngle(dipTime1, dipTime2, period1, period2, observerAngleDeg);
    ratio = interval / basePeriod;
    expectedInterval = round(ratio) * basePeriod;
    return(interval - expectedInterval);
}

runRandomTests = function(numSimulations = 10000, observerAngleDeg = 138) {
    errors = numeric(0)
    for(s in 1:numSimulations) {
        dipTime1 = runif(1, 100, 1600)
        dipTime2 = runif(1, 100, 1600)
        period1 = sample(20) + 5
        period2 = sample(20) + 5
        error = errorAtObserverAngle(dipTime1, dipTime2, period1, period2, observerAngleDeg)
        errors = c(errors, error)
    }
    print(paste('Mean absolute error: ', mean(abs(errors))))
}

intervalAtObserverAngle = function(dipTime1, dipTime2, period1, period2, observerAngleDeg) {
    circFraction = normalizeAngleDeg(observerAngleDeg) / 360;
    newDipTime1 = dipTime1 +  ifelse(period1 > 0, period1 * circFraction, period1 * (1 - circFraction)); 
    newDipTime2 = dipTime2 +  ifelse(period2 > 0, period2 * circFraction, period2 * (1 - circFraction));
    return(abs(newDipTime1 - newDipTime2));
}

normalizeAngleDeg = function(angle) {
    return(normalizeAngle(angle * pi / 180) * 180 / pi);
}

normalizeAngle = function(angle) {
    while(angle < 0) {
        angle = angle + 2 * pi;
    }
    while(angle >= 2 * pi) {
        angle = angle - 2 * pi;
    }
    return(angle);
}