from random import random
from math import sqrt,pi
from Gnuplot import Gnuplot,Data
from math import pi
from Modul4 import det_pi,det_area

# Berechne pi einmal mit 100000 Punkten
Npi=1
Npoints=100000
pimc,sigpimc,inc,outc=det_pi(Npi,Npoints,1)

# Grafische Darstellung der Integration
pinc=Data(inc,with_="d lt 1")
poutc=Data(outc,with_="d lt 3")
gp=Gnuplot()
gp("set size square")
gp("f(x)=sqrt(1-x**2)")
gp("set xrange[0:1]")
gp("set yrange[0:1]")
gp("set title 'Monte-Carlo-Integration, {/Symbol= p} = %f'"%(pimc))
gp("set label '%i Punkte' at 0.5,0.5 center front textcolor lt 1"%(len(inc)))
gp("set label '%i Punkte' at 0.75,0.9 center front tc lt 3"%(len(outc)))
gp.plot(pinc,poutc,"f(x) with l lt -1 lw 1 notitle")
gp.hardcopy("pi.ps",color=1,fontsize=24)
print "pi aus %i Punkten = "%(Npoints),pimc
raw_input("Weiter")

# Berechne pi 1000 mal mit 1000 Punkten
Npi=1000
Npoints=1000
pimc,sigpimc=det_pi(Npi,Npoints)

# Gebe das Ergebnis sammt Fehler aus und vergleiche die Abweichung
print "pith = ",pi
print "Mittelwert von %i bestimmten Werten fuer pi und die Standardabweichung des Mittelwertes. Fuer jeden Einzelwert wurden %i Punkte verwendet."%(Npi,Npoints)
print "pimc = ",pimc," +- ",sigpimc
print "|pith-pimc|/sigpimc = ",abs(pi-pimc)/sigpimc
raw_input("Weiter")

# Berechne die Flaeche die von 4 Geraden eingeschlossen wird
parameter=[[0.3,0.],[0.7,0.3],[3.,-0.1],[4.,-1.]]
Npoints=100000
A,ina,outa=det_area(Npoints,parameter)

# Gebe das Ergebnis aus.
print "Flaeche (%i Punkte) = "%(Npoints),A

# Grafische Darstellung der Integration
pouta=Data(outa)
pouta.set_option_colonsep("with","d lt 3")
pina=Data(ina)
pina.set_option_colonsep("with","d lt 1")
gp2=Gnuplot()
gp2("set size square")
gp2("f1(x)=%f*x+%f"%(parameter[0][0],parameter[0][1]))
gp2("f2(x)=%f*x+%f"%(parameter[1][0],parameter[1][1]))
gp2("f3(x)=%f*x+%f"%(parameter[2][0],parameter[2][1]))
gp2("f4(x)=%f*x+%f"%(parameter[3][0],parameter[3][1]))
gp2("set xrange[0:1]")
gp2("set yrange[0:1]")
gp2("set title 'Monte-Carlo-Integration, A = %f'"%(A))
gp2("set label 'A' at 0.2,0.25 center front")
gp2("set mxtics 2")
gp2("set mytics 2")
gp2("set border front lw 3")
gp2.plot(pouta,pina,"f1(x) notitle with l lt -1","f2(x) notitle with l lt -1","f3(x) notitle with l lt -1","f4(x) notitle with l lt -1")
gp2.hardcopy("Flaeche.ps",color=1,fontsize=24)
raw_input("Ende")