import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

#a)
#Vysledky urcime zo skriptu z predchadzajuceho cvicenia.
#Pozadie (background) je subor vsetkych pixelov bez odlahlych bodov
#Intenzita pozadia je najvhodnejsie definovana ako median pozadia - 1721
#Hodnotu sumu najlepsie popisuje robustna odchylka pozadia, ktorej hodnota je 190.
#Je mozne pouzit aj menej vhodnu standardnu odchylku pozadia s hodnotou 197.

#b)
def my_median(matrix):
    matrix = np.array(matrix)
    matrix = np.sort(matrix.flatten())
    n = len(matrix)
    if n%2 == 1:
        return matrix[(n-1)/2]
    else:
        return 0.5*(matrix[n/2]+matrix[n/2-1])

#c)
def my_modus(matrix):
    matrix = np.array(matrix)
    matrix = np.sort(matrix.flatten())
    N = 1   #najvyssi pocet vyskytov 
    i = 0
    result = matrix[0]   #aktualny modus
    while i <= len(matrix):
        n = 1
        while i+n < len(matrix) and matrix[i] == matrix[i+n]:
            n += 1
        if n > N:   #ak najdeme novy modus s castejsim vyskytom, menime N a result
            N = n
            result = matrix[i]
        i += n
    return result, N   #vracia modus a jeho pocetnost

#d)
m44 = np.loadtxt("E:\m44.dat")
m44 = np.array(m44)
m44_smooth = m44
sigma_r = 1.482*np.median(np.abs(m44-np.median(m44)))   #robustna odchylka zodpovedajuca sumu
background = 1721
for i in range(len(m44)):
    for j in range(len(m44)):
        if np.abs(m44_smooth[i][j]-background)<3.5*sigma_r:
            m44_smooth[i][j] = background
x,y = np.meshgrid(range(N),range(N))   
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x,y,m44,rstride=1,cstride=1,cmap=cm.coolwarm)
#posledne 4 riadky su skopirovane z predchadzajuceho cvicenia

#e)
def inst_mag(matrix,background,x,y,r):
    n = (2*r+1)**2.0   #pocet pixelov vyrezu z okolia hviezdy
    I = np.sum(matrix[x-r:x+r+1,y-r:y+r+1])-n*background
    dI = np.sqrt(n)*sigma_r
    i_mag = -2.5*np.log10(I)
    di_mag1 = -2.5*np.log10(I+dI)+2.5*np.log10(I)   #vypocet hornej chyby
    di_mag2 = -2.5*np.log10(I)+2.5*np.log10(I-dI)   #vypocet dolnej chyby
    return i_mag, di_mag1, -di_mag2

#f)
def find_stars(matrix,n):   #najde n najjasnejsich pixelov spolu s ich poziciou
    sequence = np.sort(matrix.flatten())   #zoradi prvky matice vzostupne
    sequence = sequence[::-1]   #zoradi vektor zostupne
    result = []
    for i in range(n):
        x,y = np.where(matrix == sequence[i])   #vrati pozicie najjasnejsich pixelov
        result.append([sequence[i],x,y])
    return np.array(result)

print find_stars(m44,20)
imag_gCnc = inst_mag(m44,1721,193,19,3)
imag_M44 = inst_mag(m44,1721,125,65,30)
imag_lim = -2.5*np.log10(5*sigma_r)
mag_corr = 4.65-imag_gCnc[0]   #korekcia medzi instrumentalnou a zdanlivou velkostou
print 'instrumentalna hviezdna velkost gamma Cnc:', imag_gCnc
print 'instrumentalna hviezdna velkost M44:', imag_M44
print 'zdanliva hviezdna velkost M44:', imag_M44[0]+mag_corr
print 'limitna hviezdna velkost:', imag_lim+mag_corr