%1. uloha
d = 0.04;   %suradnicove delenie, pre podrobnejsiu mapu zmensime delenie
dt = 0.04;   %rotacne delenie
u = 0.5;    %koeficient okrajoveho stemnenia
r = 0.2;    %polomer exoplanety
int = [];   %vysledna intenzita
for t = -1.4:dt:1.4     %pohyb exoplanety po disku
    A = [];
    for x2 = -1:d:1  %x-ova suradnica
        for y2 = -1:d:1  %y-ova suradnica
            if (x2^2+y2^2) <= 1   %definicia disku hviezdy
                if ((x2-t)^2+y2^2) <= r^2     %poloha exoplanety
                    A = [A;x2,y2,0];  %intenzita exoplanety
                else
                    A = [A;x2,y2,1-u*(1-sqrt(1-x2^2-y2^2))];  %intenzita hviezdy
                end
            end
        end
    end
    scatter(A(:,1),A(:,2),10*ones(size(A,1),1),A(:,3),'s','filled')
    axis equal  %obe osi zobrazuje v rovnakej mierke
    colormap(gray)  %miesto automatickej farebnej mapy zvoli nazornejsiu ciernobielu
    set(gca, 'Clim', [0 1]) %zvoli rozsah farebnej mapy od 0 do 1
    int = [int;sum(A(:,3))];    %scitava celkovu intenzitu hviezdy
    pause(0.01)
end
int = int/max(int);
figure
plot(-1.4:dt:1.4,-2.5*log10(int))    %intenzitu premenime na magnitudy
xlabel('t')
ylabel('standardized intensity [mag]')
axis ij     %prevracia y-ovu os


%2. uloha
D = 3;    %hranica y/z hodnot
R = 2;    %polomer galaktickeho disku
d = 0.05;    %suradnicove delenie
x2 = 1;    %riadkovy index matice Sigma
Sigma = [];
for x1 = -R:d:R
    y2 = 1;    %stlpcovy index matice Sigma
    for y1 = -R:d:R
        r = sqrt(x1^2+y1^2);    %definicia vzdialenosti od stredu disku
        if r <= R
            Sigma(x2,y2) = (1+r^2)^(-3/2);    %definicia plosnej hustoty pre ...
            ...bod o polohe [x,y]
        else
            Sigma(x2,y2) = 0;    %hustota mimo galaktickeho disku
        end
        y2 = y2+1;
    end
    x2 = x2+1;
end

Fi = []; Y = []; Z = []; 
x1 = 0;
for y1 = -D:d:D
    for z1 = [-D:d:-0.001,0.001:d:D]    %odstranujeme body z1 = 0, v ktorych ...
        ...je potencial -Inf
        fi = 0;
        for x2 = -R:d:R
            for y2 = -R:d:R
                r = sqrt((x1-x2)^2+(y1-y2)^2+z1^2);    %vzdialenost bodu [x1,y1,z1]...
                ...od jednotlivych elementov disku [x2,y2,z2=0]
                M = Sigma(round((x2+R)/d+1),round((y2+R)/d+1))*d^2;    %hustota...
                ...elementu v pozicii [x2,y2,z2=0] nasobena plosnym diferencialom
                fi = fi-M/r;    %scitavame ciastkove potencialy od jednotlivych bodov
            end
        end
        Fi = [Fi;fi]; Y = [Y;y1]; Z = [Z;z1];
    end
end
figure
scatter(Y,Z,10,Fi,'*')    %velkost bodov = 10, Fi je reprezentovane farbou
set(gca,'CLim',[3*median(Fi),max(Fi)]);    %upravi farebne skalovanie
colorbar    %vytvori farebnu legendu
xlabel('y')
ylabel('z')