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


### Aufgabe 3
#############

# b)

faces_conv <- function(x,y,pl=TRUE){
### Diese Funktion zeichnet eine Punktwolke zusammen mit ihrer konvexen Hülle 
### und ermittelt die Zahl der Kanten der konvexen Hülle.
### Als Eingabe erwartet sich den Vektor der x-Koordinaten der Punkte und den 
### Vektor der y-Koordinaten.
### Als drittes Argument folgt noch ein Boolescher Wert, der angibt, ob die Plots 
### tatsächlich gezeichnet werden sollen.
if(length(x) != length(y)){
print("Fehler! Untrschiedlich Anzahlen von x- und y-Koordinaten")
return -1;
}

if(pl){plot(x,y);}

n=length(x)
f=0; 				# Zähler für die Zahl der Kanten

for(i in (2:n)){
for(j in (1:(i-1))){
				# Doppelte Schleife zum Durchlsaufen aller Punktepaare
orth= c(y[j]-y[i], x[i]-x[j])  # Ein Vektor, der auf der Gerade durch (x[i], y[i])
					# und (x[j],y[j]) senkrecht steht. 
scalarpro = orth[1]*x + orth[2]*y; # Die Skalarprodukte aller Punkte mit dem 
						# Normalvektor
#scalarpro_this =scalarpro[i];
extr = ((scalarpro[i] == min(scalarpro)) || (scalarpro[i] == max(scalarpro))
		|| (scalarpro[j] == min(scalarpro)) || (scalarpro[j] == max(scalarpro)) );
		# Die Verbindungsstrecke zwischen (x[i], y[i]) und (x[j], y[j])
		# ist genau dann eine Kante wegen die beiden beteiligten Punkte
		# mit dem Normalvektor entweder minimales oder maximales 
		# Skalarprodukt haben.
if(extr){
		# In diesem Fall: Zähler für die Kantenzahl um 1 erhöhen und
		# Kante einzeichnen. 
f=f+1;
if(pl){lines(c(x[i], x[j]), c(y[i], y[j]))}
}
}} # for i, j
f
} # function


# c)

x=runif(100)	
y=runif(100)
faces_conv(x,y)

# d)

Erg=rep(0,1000)
for(i in (1:1000)){
x=runif(100)	
y=runif(100)
Erg[i]=faces_conv(x,y,pl=FALSE)}
mean(Erg)
var(Erg)
hist(Erg)


# e)

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


#Test der Funktion
rball(6)
plot(rball(1000))


points=rball(100)
px=points$x
py=points$y
faces_conv(px,py)


Ergb=rep(0,1000)
for(i in (1:1000)){
points=rball(1000)
px=points$x
py=points$y
Ergb[i]=faces_conv(px,py,pl=FALSE)}
mean(Ergb)
var(Ergb)
hist(Ergb)

# f)

# Wir berechnen Näherungswerte aus den asymtotischen Formeln
8*log(100)		# Quadrat
0.789*100^(1/3)	# Kreis