#calculate average field for the Sun from fits files

from astropy.io import fits
from numpy import mean,sin,cos,pi,sqrt,median
from scipy.stats import mode

im1=fits.open('mrzqs081231t2354c2078_110.fits')
im2=fits.open('mrzqs191231t2354c2225_011.fits')

lat=60

data1=im1[0].data
data2=im2[0].data

print(data1)
print(data1.shape)

avg1=[] #include here all the data points above the limiting latitude
avg2=[]
avgglob1=[] #here all values
avgglob2=[]
Btot1=[] #not available with only one component of B (ASSUMED to be radial...)
Btot2=[]
sqsum1=[] #sum of squares for rms
sqsum2=[]

#print(data1[0,:])

n=180*360 #len(data1)*len(data1)[0] #number of grid points
w=[] #weights

for i in range(0,len(data1)): #WEIGHT!!!!!!!
    latitude=-90+i #CHECK DETAILS (+-1 OR SO...)
    a=sin(latitude*pi/180) #grid element size a=1deg*2pi*R*sin(latitude)/360
    wei=(pi/180)*cos(latitude*pi/180)*n/(2*360) #CHECK!!!!! w=a*n/4pi (full solid angle = 4pi)
    #print(wei)
    #print(latitude*pi/180)
    for j in range(0,len(data1[i])):
        avgglob1.append(wei*data1[i,j])
        sqsum1.append(wei*data1[i,j]**2)
        if (i+1 >= 90+lat):
            avg1.append(wei*data1[i,j])

m1=mean(avg1)
mg1=mean(avgglob1) #global mean
rms1=sqrt(sum(sqsum1)/len(sqsum1))
mode1=mode(avgglob1)[0] #global mode
median1=median(avgglob1) #global median

for i in range(0,len(data2)):
    latitude=-90+i #CHECK DETAILS (+-1 OR SO...)
    a=sin(latitude*pi/180) #grid element size a=1deg*2pi*R*sin(latitude)/360
    wei=(pi/180)*cos(latitude*pi/180)*n/(2*360) #CHECK!!!!! w=a*n/4pi (full solid angle = 4pi)
    for j in range(0,len(data2[i])):
        avgglob2.append(wei*data2[i,j])
        sqsum2.append(wei*data2[i,j]**2)
        if (i+1 >= 90+lat):
            avg2.append(wei*data2[i,j])

m2=mean(avg2)
mg2=mean(avgglob2) #global mean
rms2=sqrt(sum(sqsum2)/len(sqsum2))
mode2=mode(avgglob2)[0] #global mode
median2=median(avgglob2) #global median

print(avg1[0])

print('mean 1, lat > ',str(lat),': ',str(m1))
print('mean 2, lat > ',str(lat),': ',str(m2))
print('mean 1, global: ',str(mg1))
print('mean 2, global: ',str(mg2))
print('rms 1, global: ',str(rms1))
print('rms 2, global: ',str(rms2))
print('mode 1, global: ',str(mode1))
print('mode 2, global: ',str(mode2))
print('median 1, global: ',str(median1))
print('median 2, global: ',str(median2))