from math import *
import numpy
import astropy.io.fits

f = astropy.io.fits.open("M87_R_cal.fits")
intensity = f[0].data

# stred galaxie
i0 = 241
j0 = 351

# pozadi
bg = 3.6e5

# meritko 1 pix = 0.960 arcsec
s2 = 0.96**2

# priblizny fotonovy tok v R pro mag = 0
phflux0 = 0.5e8*150  # ph/s/m2

# expozicni doba a plocha
TA = 3600 * 0.2

bmag = - 2.5*log10((bg / TA / s2) / phflux0)
print("prumerna magnituda pozadi:",bmag)

# scitani
flux = 0
n = 0
for i in range(-100,100):
    x = i + i0
    for j in range(-100,100):
        y = j + j0
        r = sqrt(i**2 + j**2)
        if 45 < r and r < 55:
            n = n + 1
            flux = flux + (intensity[x,y] - bg)

print("pocet pixelu=",n)
print("tok v mezikruzi=",flux)

print("prumerna intensita v mezikruzi [ph/s/m2/arcsec]:",flux/n/s2/TA)


# median
flux = []
n = 0
for i in range(-100,100):
    x = i + i0
    for j in range(-100,100):
        y = j + j0
        r = sqrt(i**2 + j**2)
        if 45 < r and r < 55:
            n = n + 1
            flux.append(intensity[x,y] - bg)

flux1 = numpy.array(flux)           
            
int1 = numpy.median(flux1) / TA
mag = - 2.5*log10(int1 / phflux0 / s2)
print("prumerna intensita v mezikruzi:",int1)
print("prumerna magnituda v mezikruzi:",mag)