options(echo = FALSE)
####################################################################

##Execute following on dos prompt (omit the "##"):
##"S:\R\rw2001\bin\Rterm.exe" -q --no-restore --no-save < ellipse_test.r > R.out

####################################################################

#--------------------------------------------------------------------------
# functions:

ellipse <- function(hlaxa = 1, hlaxb = 1, theta = 0, xc = 0, yc = 0, npoints = 100)
{
        xp <- NULL
        yp <- NULL
        for(i in 0:npoints) {
                a <- (2 * pi * i)/npoints
                x <- hlaxa * cos(a)
                y <- hlaxb * sin(a)
                if(theta != 0) {
                        alpha <- angle(x, y)
                        rad <- sqrt(x^2 + y^2)
                        x <- rad * cos(alpha + theta)
                        y <- rad * sin(alpha + theta)
                }
                xp <- c(xp, x + xc)
                yp <- c(yp, y + yc)
        }
        lines(xp, yp, type = "l")
        invisible()
}

#-----------
angle <- function(x, y)
{
        if(x > 0) {
                atan(y/x)
        }
        else {
                if(x < 0 & y != 0) {
                        atan(y/x) + sign(y) * pi
                }
                else {
                        if(x < 0 & y == 0) {
                                pi
                        }
                        else {
                                if(y != 0) {
                                        (sign(y) * pi)/2
                                }
                                else {
                                        NA
                                }
                        }
                }
        }
}

#-----------
ellipsem <- function (mu, amat, c2, npoints = 100) 
{
    if (dim(amat) == c(2, 2)) {
        eamat <- eigen(amat)
        hlen <- sqrt(c2/eamat$val)
        theta <- angle(eamat$vec[1, 1], eamat$vec[2, 1])
        ellipse(hlen[1], hlen[2], theta, mu[1], mu[2], npoints = npoints)
    }
    invisible()
}

#--------------------------------------------------------------------------
# Read in data file:
xydata <- scan("P:\\My Documents\\FrankMasci\\LinearRegress_underR\\Mvt1_Mvt2.txt", list(0,0));

# Assign variables to columns:
x <- xydata[[1]];
y <- xydata[[2]];

#--------------------------------------------------------------------------
# Simple ellipse:
# hlaxa    half-length of the major axis (parallel to the X-axis
#          when theta = 0) 
# hlaxb    half-length of the minor axis
# theta    angle of rotation, measured counter-clockwise from the
#          positive X-axis in radians
# xc       the X-coordinate of the centre
# yc       the Y-coordinate of the centre
# npoints  the number of line segments used to approximate the ellipse

# Following from my notes: define parameters for error ellipse using
# actual covariance matrix for the data.

covmat <- cov(cbind(x, y), use = "all.obs", method = "pearson");
varx <- covmat[1,1];
vary <- covmat[2,2];
cov <- covmat[1,2];
detC <- det(covmat);
nsig <- 2; # Numer of sigma for the error ellipse (joint radius measure).
alpha <- (vary / detC) / nsig^2;
beta <- (varx / detC) / nsig^2;
gamma <- (-cov / detC) / nsig^2;
theta <- 0.5 * atan( (2 * gamma) / (alpha - beta));
b <- sqrt( ((sin(theta)^4) - (cos(theta)^4))/(alpha*(sin(theta)^2) - beta*(cos(theta)^2)) );
a <- sqrt( (cos(theta)^2)/(alpha - ((sin(theta)^2)/(b^2))) );

png(file="ellipse_test.png", width = 680, height = 450, pointsize = 12, bg="transparent", res = 200);
plot(x, y)
points(mean(x), mean(y), pch="+")
ellipse(a, b, theta);

#--------------------------------------------------------------------------
# Confidence level error-ellipse for bivariate dataset.
# This draws a scatterplot, marks the mean with a "+", and draws the 95%
# confidence ellipse for the population mean. ***I don't believe the below***

#plot(x, y)
#points(mean(x), mean(y), pch="+")
#n <- length(x)
#p <- 2 # dof.
#ellipsem(c(mean(x), mean(y)), solve(var(cbind(x,y))), ((n-1)*p/(n-p))*qf(0.95,p,n-p)/n);

#------------------------------END-----------------------------------------