#########################################################################
###          Stochastische Geometrie II - Übungsblatt 2               ###
###          WS 2016/17			   Dr. Jürgen Kampf       	    ###
#########################################################################


### Aufgabe 2
#############

# a)
# Wir wandeln die Funktion aus Aufgabe 3e) des 1. Übungsblatts für 3-dimensionale
# Kugeln ab.

rball3d<-function(n){
Erg=data.frame(x=rep(0,n), y=rep(0,n), z=rep(0,n) );
while(n>0){
x=runif(min=-1,max=1, 1);
y=runif(min=-1,max=1, 1);
z=runif(min=-1,max=1, 1);
if(x^2+y^2+z^2<1){
	Erg$x[n]=x;
	Erg$y[n]=y;
	Erg$z[n]=z;
	n=n-1;
	}
}#while
Erg;
}

# Test
rball3d(5)


#b)
Euklid<-function(v){
sqrt(sum(v^2));
}

# Das Kreuzprodukt des R^3 ist nicht im Standard-Umfang von R enthalten, sondern
# muss (z.B. als Teil des RSEIS-Pakets nacinstalliert werden).
# Wir ziehen es aber vor, dieses selbst zu implementieren. 
# Die Funktionen outer() oder crossprod() liefern nicht das gewünschte.
Kreuzprodukt <- function(v,w){
c( v[2]*w[3]-v[3]*w[2],  v[3]*w[1]-v[1]*w[3],  v[1]*w[2]-v[2]*w[1])
} 

Lebesgue_Parallel<-function(u,v=NULL,w=NULL){
if (is.null(w)){
if (is.null(v)){
Euklid(u)
} else{
Euklid(Kreuzprodukt(u,v))
}
}else{
abs(u %*% Kreuzprodukt(v,w))[[1,1]]
}
}

# Tests
Lebesgue_Parallel(c(2,2,1))
Lebesgue_Parallel(c(2,2,1), c(2,2,0))
Lebesgue_Parallel(c(2,2,1), c(2,2,0), c(1,1,1))


# c)
l1 = apply(rball3d(1000),1, Lebesgue_Parallel) # Die Längen von 1000 zufällig
					# erzeugten Strecken
mean(l1)
k=1
gamma(5/2)/gamma(5/2+k/2)*gamma(3/2+k/2)/gamma(3/2)

mean(l1^4)
k=4
gamma(5/2)/gamma(5/2+k/2)*gamma(3/2+k/2)/gamma(3/2)


l2 = rep(0,1000) # Die Flächeninhalte von 1000 zufällig erzeugten Strecken
for(i in (1:1000)){
	u=rball3d(1)
	v=rball3d(1)
	l2[i] = Lebesgue_Parallel(c(u$x, u$y,u$z), c(v$x, v$y, v$z))
	}

mean(l2)
k=1
gamma(5/2)^2/gamma(5/2+k/2)^2*gamma(3/2+k/2)/gamma(3/2)*gamma(1+k/2)/gamma(1)

mean(l2^4)
k=4
gamma(5/2)^2/gamma(5/2+k/2)^2*gamma(3/2+k/2)/gamma(3/2)*gamma(1+k/2)/gamma(1)


l3 = rep(0,1000) # Die Flächeninhalte von 1000 zufällig erzeugten Strecken
for(i in (1:1000)){
	u=rball3d(1)
	v=rball3d(1)
	w=rball3d(1)
	l3[i] = Lebesgue_Parallel(c(u$x, u$y,u$z), c(v$x, v$y, v$z), c(w$x, w$y, w$z))
	}

mean(l3)
k=1
gamma(5/2)^3/gamma(5/2+k/2)^3*gamma(3/2+k/2)/gamma(3/2)*
				gamma(1+k/2)/gamma(1)*gamma(1/2+k/2)/gamma(1/2)

mean(l3^4)
k=4
gamma(5/2)^3/gamma(5/2+k/2)^3*gamma(3/2+k/2)/gamma(3/2)*
				gamma(1+k/2)/gamma(1)*gamma(1/2+k/2)/gamma(1/2)

## Es zeigt sich eine sehr gute Übereinstimmung zwischen theoretischen und 
## simulierten Werten.
## Je mehr Punkte beteiligt sind, desto kleiner werden die Polytope