from pylab import *
from numpy import *  #je dobre zvyknutis si aj na tento import
#uvedene dva importy v neinteraktivne spustanych .py suboroch robia nieco podobne
# ako ked sa interaktivny ipython spusti s option --pylab 

G=6.674e-11 #gravitacna konstanta
M=2.    # hmotnost gravitujuceho telesa
R=3. #polomer sfery z ktorej ratame vytok pola
rM=array([1.,0.,0.]) 	#poloha gravitujuceho telesa

vytok=0.  # vynulovany sumator na vytoku


def g(r):
	#funkcia, ktora ma vratit gravitacne zrychlenie v mieste r
	rel=r-rM  #vektor od gravitujuceho telesa k polohe r
	abs0=sqrt(rel[0]**2+rel[1]**2+rel[2]**2) #velkost vektora rel, suradnice sa cisluju 0,1,2 v Pythone
	return(-G*M*rel/abs0**3)  # vracia gravitacne zrychlenie

#teraz ide vypocet vytoku, vykachlickujeme sferu so stredom v pociatku o plomere R 
#pomocou sferickych uhlov theta, fi
dtheta = pi/100. #krok v uhle theta
dfi=2*pi/100.  #krok v uhle fi 
for i in range(0,100):
	theta=i*dtheta  #aktualna uhol theta kachlicky
	for j in range(0,100):
		fi=j*dfi  #aktualny uhol fi kachlicky
		n=array([sin(theta)*cos(fi),sin(theta)*sin(fi), ????]) #doplnte treetiu komponentu jednotkovy vektor kolmy na aktualnu kachlicku
		dS=n*R**2*sin(theta)*dtheta*dfi #vektor plochy kachlicky kolmy na kachlicku
		#vsimnite si, ze dS je array, lebo vzniklo nasobenim n co bolo array
		r=????  #tu vyjafrite polohovy vektor aktuaqlnej kachlicky pomocou vektora n a polomeru sfery R
		gr=g(r)  #pole (gravitacne zrychlenie) v mieste aktualnej kachlicky
		vytok=vytok+  #sem pridajte vytok cez kachlicku ako skalarny sucin g(r).dS
	
print(vytok)

#teraz pre kontrolu vypocet vytoku podla Gaussovej vety
spravnyVytok=???? 	#sem vlozte vypocet sptravnej hodnoty vytoku pola podla 
					#Gaussovej vety podla toho aka hmotnost je uzavreta v sfere
print(spravnyVytok)