cc <- c(1, 2, 3)
max <- 3
min <- -max
step <- 0.01
t <- c(seq(min, max, step), seq(max, min, -step))

# für alle Werte in cc
for (k in 1:length(cc)) {
	x <- c(0)
	y <- c(0)
	cnt <- 1
	lt2 <- length(t)/2
	
	# obere Hälfte der Ellipse berechnen
	for (i in 1:lt2) {
		d <- 16 / 121 * t[i]^2 - 48 / 11 * (4 / 11 * t[i]^2 - cc[k])
		if (d >= 0) {
			x[cnt] <- t[i]
			y[cnt] <- (4 / 11 * x[cnt] + sqrt(d)) / (24 / 11)
			cnt <- cnt + 1
		}
	}
	
	# untere Hälfte der Ellipse berechnen
	lt2 <- lt2+1
	for (i in lt2:length(t)) {
		d <- 16 / 121 * t[i]^2 - 48 / 11 * (4 / 11 * t[i]^2 - cc[k])
		if (d >= 0) {
			x[cnt] <- t[i]
			y[cnt] <- (4 / 11 * x[cnt] - sqrt(d)) / (24 / 11)
			cnt <- cnt + 1
		}
	}
	
	# Endpunkt = Anfangspunkt (sonst Loch im Plot)
	x[length(x)+1] <- x[1]
	y[length(y)+1] <- y[1]
	
	plot(x, y, type = "l", xlim = c(min, max), ylim = c(min, max))
	par(new = TRUE)
}

# Koordinatenursprung
lines(c(-0.1, 0.1), c(-0.1, 0.1))
lines(c(-0.1, 0.1), c(0.1, -0.1))

# Eigenvektoren der Kovarianzmatrix bestimmen
s <- matrix(c(3, 0.5, 0.5, 1), 2, 2)
eig <- eigen(s)
alpha <- eig$vectors

# alpha_1 und alpha_2 einzeichnen
arrows(c(0, 0), c(0, 0), alpha[1,], alpha[2,], col = c("red", "blue"))
text(alpha[1,1]+0.5, alpha[2,1]+0.3, "a_1", col = "red")
text(alpha[1,2]+0.1, alpha[2,2]+0.6, "a_2", col = "blue")