OrbitalResonanceModelTester.R

# 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)
    ));  
}

getInitialTimeDelays = function(periods) {
    return(c(
        0,
        0.5 * BASE_GAP * 17,
        0.5 * BASE_GAP * 31
    ));
}

p = function(n) {
    return(n * 13 * BASE_GAP / 2);
}

runOrbitSimulation = function(timeSpan, viewingAngle = 0) {
    periods = getPeriods();
    len = length(periods);
    initialTimeDelays = getInitialTimeDelays(periods);
    if(length(initialTimeDelays) != len) {
        stop(paste("Requires", len, " initial time delays."));
    }
    results = data.frame(Name = character(0), Time = numeric(0));
    for(i in 1:len) {
        pname = intToUtf8(i + 64);        
        period = periods[i];
        occurrences = getDipOccurrences(timeSpan, period, initialTimeDelays[i], viewingAngle = viewingAngle);       
        results = mergeResultsFrame(results, pname, occurrences);
    }        
    return(results);
}

runIntervalsSimulation = function(timeSpan, multiplesOf = BASE_GAP, viewingAngle = 0) {
    dips = runOrbitSimulation(timeSpan, viewingAngle = viewingAngle);
    times = dips[["Time"]];
    names = dips[["Name"]];
    results = data.frame(From = character(0), To = character(0), NumIntervals = numeric(0));
    for(i in 2:nrow(dips)) {
        k = i - 1;
        ni = (times[i] - times[k]) / multiplesOf;
        from = names[k];
        to = names[i];
        row = data.frame(From = from, To = to, NumIntervals = ni);
        results = rbind(results, row);
    }    
    return(results);        
}

mergeResultsFrame = function(baseFrame, pname, occurrences) {
    nameVector = baseFrame[["Name"]];
    timeVector = baseFrame[["Time"]];
    timeLen = length(timeVector);
    occLen = length(occurrences);
    timeIdx = 1;
    occIdx = 1;
    newNames = character(0);
    newTimes = numeric(0);
    repeat {
        if(timeIdx > timeLen && occIdx > occLen) {
            break;
        }
        if(timeIdx > timeLen) {
            newNames = c(newNames, rep(pname, occLen - occIdx + 1));
            newTimes = c(newTimes, occurrences[occIdx:occLen]);
            break;
        }
        else if(occIdx > occLen) {
            newNames = c(newNames, as.character(nameVector[timeIdx:timeLen]));
            newTimes = c(newTimes, timeVector[timeIdx:timeLen]);
            break;
        }
        time = timeVector[timeIdx];
        occurrence = occurrences[occIdx];
        if(approxEqual(time, occurrence)) {
            combinedName = paste0(pname, as.character(nameVector[timeIdx]));
            newNames = c(newNames, combinedName);
            newTimes = c(newTimes, time);
            timeIdx = timeIdx + 1;
            occIdx = occIdx + 1;
        }        
        else if(time > occurrence) {
            newNames = c(newNames, pname);
            newTimes = c(newTimes, occurrence);            
            occIdx = occIdx + 1;
        }
        else {
            nn = as.character(nameVector[timeIdx]);
            newNames = c(newNames, nn);
            newTimes = c(newTimes, time);                    
            timeIdx = timeIdx + 1;
        }
    }
    return(data.frame(Name = newNames, Time = newTimes));
}

approxEqual = function(a, b) {
    return(abs(a - b) < 1E-5);
}

getDipOccurrences = function(timeSpan, period, timeDelay, viewingAngle = 0) {
    dipTimes = seq(0, timeSpan, by = period);
    dipTimes = dipTimes + timeDelay;
    viewingTimeOffset = angleToTime(viewingAngle, period);
    dipTimes = dipTimes - viewingTimeOffset;
    return(dipTimes[dipTimes >= 0]);
}

angleToTime = function(angle, period) {
    return(period * angle / (2 * pi));
}