from numpy import *
from matplotlib.pyplot import *


#1. uloha
d = 0.02   #suradnicove delenie, pre podrobnejsiu mapu zmensime delenie
dt = 0.05   #rotacne delenie
u = 0.5    #koeficient okrajoveho stemnenia
r = 0.2    #polomer exoplanety
intensity = []   #vysledna intenzita
for t in arange(-1.4,1.4,dt):     #pohyb exoplanety po disku
    A = []
    for x in arange(-1.0,1.0,d):  #x-ova suradnica
        for y in arange(-1,1,d):  #y-ova suradnica
            if (x**2+y**2) < 1:   #definicia disku hviezdy
                if ((x-t)**2+y**2) <= r**2:     #poloha exoplanety
                    A.append([x,y,0])   #intenzita exoplanety
                else:
                    A.append([x,y,1-u*(1-sqrt(1-x**2-y**2))])  #intenzita hviezdy
    A = array(A)
    intensity.append(sum(A[:,2]))
    #scatter(A[:,0],A[:,1],c=A[:,2],lw = 0)   #prikay scatter vztvori 2D mapu rozlozenia intenzity
    #axis('equal')
    #show()
t = arange(-1.4,1.4,dt)
I = -2.5*log10(intensity)
plot(t,I,'*')
axis([min(t), max(t), max(I), min(I)])
show()


#2. uloha
D = 3.0    #hranica x/y/z hodnot
R = 2.0    #polomer galaktickeho disku
N = 20    #suradnicove delenie
i = 0   #riadkovy index matice Sigma

x1 = 0
FI = []
for y1 in linspace(-D,D,N*D/R):
    for z1 in linspace(-D,D,N*D/R):    #odstranujeme body z1 = 0, v ktorych je potencial -Inf
        fi = 0
        for x2 in linspace(-R,R,N):
            for y2 in linspace(-R,R,N):
                r = sqrt(x2**2+y2**2)
                r2 = sqrt((x1-x2)**2+(y1-y2)**2+z1**2)    #vzdialenost bodu [x1,y1,z1] od jednotlivych elementov disku [x2,y2,z2=0]
                if r <= R:
                    Sigma = (1+r**2)**(-3.0/2)    #definicia plosnej hustoty pre bod o polohe [x2,y2]
                else:
                    Sigma = 0    #hustota mimo galaktickeho disku
                M = Sigma*(2*R/(N-1))**2    #hustota elementu v pozicii [x2,y2,z2=0] nasobena plosnym diferencialom
                fi -= M/r2    #scitavame ciastkove potencialy od jednotlivych bodov
        FI.append([y1,z1,fi])
FI = array(FI)
figure(figsize=(10,10))
scatter(FI[:,0],FI[:,1],c=FI[:,2],lw = 0)
axis('equal')