import numpy as np
import matplotlib.pyplot as plt


#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 np.arange(-1.4, 1.4, dt):     #pohyb exoplanety po disku
    A = []
    for x in np.arange(-1., 1., d):  #x-ova suradnica
        for y in np.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 - np.sqrt(1 - x**2 - y**2))])  #intenzita hviezdy
    A = np.array(A)
    intensity.append(sum(A[:, 2]))
    plt.scatter(A[:, 0], A[:, 1], c=A[:, 2], lw = 0, s=5, vmin=0, vmax=1)   #prikaz scatter vztvori 2D mapu rozlozenia intenzity
    plt.axis('equal')
    plt.gray()
    plt.show()
    plt.pause(0.01)
t = np.arange(-1.4, 1.4, dt)
I = -2.5*np.log10(intensity)
plt.figure()
plt.plot(t, I, '*')
plt.axis([min(t), max(t), max(I), min(I)])
plt.show()


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

x1 = 0
FI = []
for y1 in np.arange(-D, D, d):
    for z1 in np.arange(-D, D, d):    #odstranujeme body z1 = 0, v ktorych je potencial -Inf
        fi = 0
        for x2 in np.arange(-R, R, d):
            for y2 in np.arange(-R, R, d):
                r = np.sqrt(x2**2 + y2**2)
                r2 = np.sqrt((x1 - x2)**2 + (y1 - y2)**2 + z1**2)    #vzdialenost bodu [x1,y1,z1] od jednotlivych elementov disku [x2,y2,z2=0]
                Sigma = (1 + r**2)**(-3.0/2)    #definicia plosnej hustoty pre bod o polohe [x2,y2]
                M = Sigma*d**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 = np.array(FI)
plt.figure()
plt.scatter(FI[:, 0], FI[:, 1], c=FI[:, 2], lw = 0, s=10)
plt.colorbar()
plt.axis('equal')