espadonsutils.py

import matplotlib
matplotlib.use("Agg") 
from astropy.io import fits as pyfits
import numpy as np
from numpy import median,sqrt,array,exp
import scipy
from scipy import signal,special,optimize,interpolate
import scipy.special as sp
import copy
import glob
import os
from pylab import *
import sys
base = '../'
sys.path.append(base+"utils/GLOBALutils")
import GLOBALutils



def FileClassify(diri, log,mode='so',amps='a'):
	"""
	Classifies all files in a directory and writes a night log of science images
	"""
	if mode == 'so':
		md = 'star only'
	elif mode == 'ss':
		md = 'star+sky'
	else:
		md = 'Polarimetry'

	if amps == 'a':
		samps = 'a'
	elif amps == 'ab':
		samps = 'a,b'

	# define output lists
	simThAr_sci    = []
	simSky_sci     = []
	biases         = []
	bias_dates     = []
	flat_dates     = []
	flats          = []
	ThAr_ref       = []
	ThAr_ref_dates = []
	mode           = []

	f = open(log,'w')
	bad_files = []
	if os.access(diri+'bad_files.txt',os.F_OK):
		bf = open(diri+'bad_files.txt')
		linesbf = bf.readlines()
		for line in linesbf:
			bad_files.append(diri+line[:-1])
		bf.close()

	all_files = glob.glob(diri+"*fits")

	jj = 0

	for archivo in all_files:
		jj+=1
		dump = False
		for bf in bad_files:
			if archivo == bf:
				dump = True
				break

		if not dump:
			h = pyfits.open(archivo)
			#print archivo, h[0].header['EXPTYPE'], h[0].header['INSTMODE']
			if md in h[0].header['INSTMODE'] and samps == h[0].header['AMPLIST']:
				print archivo, h[0].header['EXPTYPE'], h[0].header['INSTMODE'], h[0].data.shape, h[0].header['AMPLIST']
				if h[0].header['EXPTYPE'] == 'OBJECT':
					simThAr_sci.append(archivo)
					obname = h[0].header['OBJNAME']
					ra     = h[0].header['OBJRA']
					delta  = h[0].header['OBJDEC']
					airmass= h[0].header['AIRMASS']
					texp   = h[0].header['EXPTIME']
					date   = h[0].header['DATE-OBS'] + h[0].header['UTC-OBS']
					line = "%-15s %10s %10s %8.2f %4.2f %s %8s %s\n" % (obname, ra, delta, texp, airmass, h[0].header['EXPTYPE'], date, archivo)
					f.write(line)

				elif h[0].header['EXPTYPE'] == 'BIAS':
					biases.append(archivo)
					mjd,mjd0 = mjd_fromheader(h)
					bias_dates.append(mjd)

				elif h[0].header['EXPTYPE'] == 'FLAT':
					flats.append(archivo)
					mjd,mjd0 = mjd_fromheader(h)
					flat_dates.append(mjd)

				elif h[0].header['EXPTYPE'] == 'COMPARISON':
					mjd, mjd0 = mjd_fromheader(h)
					ThAr_ref.append(archivo) 
					ThAr_ref_dates.append( mjd )

	f.close()
	biases, bias_dates = np.array(biases), np.array(bias_dates)
	flats, flat_dates  = np.array(flats), np.array(flat_dates)
	IS = np.argsort(bias_dates)
	biases, bias_dates = biases[IS], bias_dates[IS]
	IS = np.argsort(flat_dates)
	flats, flat_dates = flats[IS], flat_dates[IS]

	return biases, flats, ThAr_ref, simThAr_sci, ThAr_ref_dates

def mjd_fromheader(h):
    """
    return modified Julian date from header
    """
    
    datetu = h[0].header['DATE-OBS']
    timetu = h[0].header['UTC-OBS'].split(':')

    mjd0,mjd,i = GLOBALutils.iau_cal2jd(int(datetu[:4]),int(datetu[5:7]),int(datetu[8:10]))

    ho = int(timetu[0])
    mi = int(timetu[1])
    se = float(timetu[2])
    ut = float(ho) + float(mi)/60.0 + float(se)/3600.0
    mjd_start = mjd + ut/24.0

    secinday = 24*3600.0
    fraction = 0.5
    texp     = h[0].header['EXPTIME'] #sec

    mjd = mjd_start + (fraction * texp) / secinday

    return mjd, mjd0

def MedianCombine(ImgList, zero_bo=False, zero='MasterBias.fits'):
	"""
	Median combine a list of images
	"""
	if zero_bo:
		BIAS = pyfits.getdata(zero)

	n = len(ImgList)
	if n==0:
		raise ValueError("empty list provided!")

	h = pyfits.open(ImgList[0])[0]
	d = h.data
	d = OverscanTrim(d)
	if zero_bo:
		d -= BIAS
		d = np.round(d).astype('int')
	factor = 1.25
	if (n < 3):
		factor = 1
	ronoise = factor * h.header['RDNOISEA'] / np.sqrt(n)
	gain    = h.header['GAINA']
	if (n == 1):
		return d, ronoise, gain
	else:
		for i in range(n-1):
			h = pyfits.open(ImgList[i+1])[0]
			ot = OverscanTrim(h.data)
			if zero_bo:
				d = np.dstack((d,np.round((ot - BIAS)).astype('int')))
			else:
				d = np.dstack((d,np.round(ot).astype('int')))
	return np.median(d,axis=2), ronoise, gain

def OverscanTrim(d):
    """
    Overscan correct and Trim a refurbished FEROS image
    """
    os = d[:,2060:]
    overscan = np.median(os,axis=1)[:4608]
    newdata = np.zeros(d[:4608,:2049].shape)
    for i in range(len(overscan)):
        newdata[i,:] = d[i,:2049] - overscan[i]
    return newdata

def get_them(sc,exap,ncoef,mode=1,shift=0.):
	def gauss(params,x):
		amp1 = params[0]
		med1 = params[1]
		sig1 = params[2]
		g1 = amp1 * np.exp(-0.5*((x-med1)/sig1)**2)
		return g1
	def res_gauss(params,g,x):
		return g-gauss(params,x)

	def gauss2(params,x):
		amp1 = params[0]
		amp2 = params[1]
		med1 = params[2]
		med2 = params[3]
		sig1 = params[4]
		sig2 = params[5]
		g1 = amp1 * np.exp(-0.5*((x-med1)/sig1)**2)
		g2 = amp2 * np.exp(-0.5*((x-med2)/sig2)**2)
		return g1 + g2

	def res_gauss2(params,g,x):
		return g-gauss2(params,x)

	def fit_drift(ref,dif,oref):
		coef  = np.polyfit(ref,dif,2)

		residuals = dif - np.polyval(coef,ref)
		rms = np.sqrt(np.var(residuals))
		I = np.where(np.absolute(residuals)>3*rms)[0]
		cond = True
		if len(I)==0:
			cond = False
		while cond:
			im  = np.argmax(np.absolute(residuals))
			dif = np.delete(dif,im)
			ref = np.delete(ref,im)
			coef = np.polyfit(ref,dif,2)
			residuals = dif - np.polyval(coef,ref)
			rms = np.sqrt(np.var(residuals))
			I = np.where(np.absolute(residuals)>3*rms)[0]
			if len(I)==0:
				cond = False

		cdif = np.polyval(coef,oref)
		ref = oref + cdif
		return ref
	
	sc_or = sc.copy()  
	medc = int(.5*sc.shape[1])
	d = np.median(sc[:,medc-exap:medc+exap+1],axis=1)
	#plot(d)
	ejx = np.arange(len(d))

	ref = []
	cents = np.loadtxt('ref_ords.txt').astype('int') - shift
	#plot(d)
	#plot(ejx[cents[:,0]],d[cents[:,0]],'ro')
	#plot(ejx[cents[:,1]],d[cents[:,1]],'go')
	#show()

	posO, posC = cents[:,0],cents[:,1]
	if mode == 0:
		posO = cents[:,2]
		posC = cents[:,2]
	refO,refC = [],[]
	exap2 = exap + .2*exap
	dev   = exap2/4.
	for i in range(len(posO)):
		if posO[i]-exap2 < 0:
			x = ejx[:posC[i]+exap2+1]
			y = d[:posC[i]+exap2+1]
		elif posC[i]+exap2+1 > len(d):
		    x = ejx[posO[i]-exap2:]
		    y = d[posO[i]-exap2:]
		else:	
			x = ejx[posO[i]-exap2:posC[i]+exap2+1]
			y = d[posO[i]-exap2:posC[i]+exap2+1]
		y -= y.min()
		tx1 = np.arange(x[0]-dev,x[0],1)
		tx2 = np.arange(x[-1]+1,x[-1]+dev+1,1)
		ty1 = np.zeros(len(tx1))
		ty2 = np.zeros(len(tx2))
		x = np.hstack((tx1,x,tx2))
		y = np.hstack((ty1,y,ty2))

		midi = int(0.5*len(x))
		
		if mode == 1:
			guess = [np.max(y[:midi]),np.max(y[midi:]),x[0]+np.argmax(y[:midi]),x[0]+midi+np.argmax(y[midi:]),dev,dev]
			p, success =  scipy.optimize.leastsq(res_gauss2, guess, args=(y,x))
			refO.append(p[2])
			refC.append(p[3])
		else:
			#plot(x,y)
			#show()
			#print fcds
			guess = [np.max(y),x[0]+np.argmax(y),dev]
			p, success =  scipy.optimize.leastsq(res_gauss, guess, args=(y,x))
			refO.append(p[1])
			refC.append(p[1])

	refO = np.array(refO)
	refOt = np.around(refO).astype('int')
	matO = np.zeros((len(refO),sc.shape[1]))
	matO[:,medc] = refO
	refC = np.array(refC)
	refCt = np.around(refC).astype('int')
	matC = np.zeros((len(refC),sc.shape[1]))
	matC[:,medc] = refC

	#plot(d)
	#plot(cents[:,0],d[cents[:,2]],'ro')
	#plot(cents[:,1],d[cents[:,1]],'go')
	#plot(refOt,d[refOt],'ko')
	#plot(refCt,d[refCt],'yo')
	#show()
	#print gfd5
	i = medc -1
	while i >=0:
		d = sc[:,i]
		#if i<1000:
		#	plot(d)
		j = 0
		posO = np.around(refO).astype('int')
		posC = np.around(refC).astype('int')
		trefO,trefC = [],[]
		exap2 = exap + .2*exap
		dev   = exap2/4.

		while j < len(posO):
			#print j
			if posO[j]-exap2 < 0:
				x = ejx[:posC[j]+exap2+1]
				y = d[:posC[j]+exap2+1]
			elif posC[j]+exap2+1 > len(d):
				x = ejx[posO[j]-exap2:]
				y = d[posO[j]-exap2:]
			else:	
				x = ejx[posO[j]-exap2:posC[j]+exap2+1]
				y = d[posO[j]-exap2:posC[j]+exap2+1]
			if len(x) < 7:
				trefO.append(refO[j])
				trefC.append(refC[j])
			else:
				y -= y.min()
				tx1 = np.arange(x[0]-dev,x[0],1)
				tx2 = np.arange(x[-1]+1,x[-1]+dev+1,1)
				ty1 = np.zeros(len(tx1))
				ty2 = np.zeros(len(tx2))
				x = np.hstack((tx1,x,tx2))
				y = np.hstack((ty1,y,ty2))
				midi = int(0.5*len(x))
				if mode == 1:
					guess = [np.max(y[:midi]),np.max(y[midi:]),x[0]+np.argmax(y[:midi]),x[0]+midi+np.argmax(y[midi:]),dev,dev]
					p, success =  scipy.optimize.leastsq(res_gauss2, guess, args=(y,x))
					trefO.append(p[2])
					trefC.append(p[3])
				else:
					guess = [np.max(y),x[0]+np.argmax(y),dev]
					p, success =  scipy.optimize.leastsq(res_gauss, guess, args=(y,x))
					trefO.append(p[1])
					trefC.append(p[1])

			j+=1

		orefO = refO.copy()
		trefO = np.array(trefO)
		difO = trefO-refO
		refO = fit_drift(refO,difO,orefO)
		matO[:,i] = refO
		orefC = refC.copy()
		trefC = np.array(trefC)
		difC = trefC-refC
		refC = fit_drift(refC,difC,orefC)
		matC[:,i] = refC

		#if i<500:
		#	plot(d)
		#	I = np.where(refO<2049)[0]
		#	plot(refO[I].astype('int'),d[refO[I].astype('int')],'ko')
		#	show()
		#	print gfd

		i-=4

	refO = matO[:,medc]
	refC = matC[:,medc]

	i = medc +1
	while i < sc.shape[1]:
		d = sc[:,i]
		j = 0
		posO = np.around(refO).astype('int')
		posC = np.around(refC).astype('int')
		trefO,trefC = [],[]
		exap2 = exap + .2*exap
		dev   = exap2/4.

		while j < len(posO):
			#print j
			if posO[j]-exap2 < 0:
				x = ejx[:posC[j]+exap2+1]
				y = d[:posC[j]+exap2+1]
			elif posC[j]+exap2+1 > len(d):
				x = ejx[posO[j]-exap2:]
				y = d[posO[j]-exap2:]
			else:	
				x = ejx[posO[j]-exap2:posC[j]+exap2+1]
				y = d[posO[j]-exap2:posC[j]+exap2+1]
			if len(x) < 7:
				trefO.append(refO[j])
				trefC.append(refC[j])
			else:
				y -= y.min()
				tx1 = np.arange(x[0]-dev,x[0],1)
				tx2 = np.arange(x[-1]+1,x[-1]+dev+1,1)
				ty1 = np.zeros(len(tx1))
				ty2 = np.zeros(len(tx2))
				x = np.hstack((tx1,x,tx2))
				y = np.hstack((ty1,y,ty2))
				midi = int(0.5*len(x))
				if mode == 1:
					guess = [np.max(y[:midi]),np.max(y[midi:]),x[0]+np.argmax(y[:midi]),x[0]+midi+np.argmax(y[midi:]),dev,dev]
					p, success =  scipy.optimize.leastsq(res_gauss2, guess, args=(y,x))
					trefO.append(p[2])
					trefC.append(p[3])
				else:
					guess = [np.max(y),x[0]+np.argmax(y),dev]
					p, success =  scipy.optimize.leastsq(res_gauss, guess, args=(y,x))
					trefO.append(p[1])
					trefC.append(p[1])
			j+=1

		orefO = refO.copy()
		trefO = np.array(trefO)
		difO = trefO-refO
		refO = fit_drift(refO,difO,orefO)
		matO[:,i] = refO

		orefC = refC.copy()
		trefC = np.array(trefC)
		difC = trefC-refC
		refC = fit_drift(refC,difC,orefC)
		matC[:,i] = refC

		i+=4

	#imshow(sc_or,vmax=5000)
	for i in range(matO.shape[0]):
		y = matO[i]
		x = np.arange(len(y))
		I = np.where(y!=0)[0]
		x,y = x[I],y[I]
		coef = np.polyfit(x,y,ncoef)
		#plot(x,y,'r')
		residuals = y - np.polyval(coef,x)
		rms = np.sqrt(np.var(residuals))
		I = np.where(np.absolute(residuals)>3*rms)[0]
		cond = True
		if len(I)==0:
			cond = False
		while cond:
			im = np.argmax(np.absolute(residuals))
			x = np.delete(x,im)
			y = np.delete(y,im)
			coef = np.polyfit(x,y,ncoef)
			residuals = y - np.polyval(coef,x)
			rms = np.sqrt(np.var(residuals))
			I = np.where(np.absolute(residuals)>3*rms)[0]
			if len(I)==0:
				cond = False
		if i == 0:
			acoefsO = coef
		else:
			acoefsO = np.vstack((acoefsO,coef))
		#plot(np.polyval(coef,np.arange(len(matO[i]))),'g')
	#show()

	if mode == 1:
		for i in range(matC.shape[0]):
			y = matC[i]
			x = np.arange(len(y))
			I = np.where(y!=0)[0]
			x,y = x[I],y[I]
			coef = np.polyfit(x,y,ncoef)
			#plot(x,y,'r')
			residuals = y - np.polyval(coef,x)
			rms = np.sqrt(np.var(residuals))
			I = np.where(np.absolute(residuals)>3*rms)[0]
			cond = True
			if len(I)==0:
				cond = False
			while cond:
				im = np.argmax(np.absolute(residuals))
				x = np.delete(x,im)
				y = np.delete(y,im)
				coef = np.polyfit(x,y,ncoef)
				residuals = y - np.polyval(coef,x)
				rms = np.sqrt(np.var(residuals))
				I = np.where(np.absolute(residuals)>3*rms)[0]
				if len(I)==0:
					cond = False
			if i == 0:
				acoefsC = coef
			else:
				acoefsC = np.vstack((acoefsC,coef))
			#plot(np.polyval(coef,np.arange(len(matC[i]))),'b')
		#show()
		#print gfds
		#print acoefs.shape
		return acoefsO, acoefsC
	else:
		return acoefsO