# Estimates and subtracts the spatial bandpass while exluding a user-specified region
# User specifies region of cube to use for bandpass calibration plus region to 
# do the recalibration on; program finds median across spatial bandpass.

import math
import numpy
import pyfits
import os

# Changes to the directory where the script is located
abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)

infile = 'AGES5523_WGTMED_h3_small.fits'

outfile = str(infile.split('.fits')[0]+str('_SimMed.fits'))

FitsFile = pyfits.open(infile)

image= FitsFile[0].data
header=FitsFile[0].header

sizez = image.shape[0]
sizey = image.shape[1]
sizex = image.shape[2]

print sizex,sizey,sizez

for xp in range(xmin,xmax):
	for yp in range(ymin,ymax):
		for zp in range(zmin,zmax):
		
			print xp
		
			vlist = []
		
			bpass = numpy.array(image[zp,yp,:])
			
			for i in range(0,len(bpass)):
				vlist.append(bpass[i])

			# Remove nans from vlist.
			# Have to convert vlist into numpy array first
			vlist = numpy.array(vlist)
			vlist = vlist[numpy.logical_not(numpy.isnan(vlist))]

					
			med =  numpy.median(vlist)
			
			bpass = bpass - med
			
			image[zp,yp,:] = bpass
			
FitsFile = pyfits.writeto(outfile,image,header=header)