import pyfits
import numpy as np
import matplotlib.pyplot as plt
import os

cmap = plt.get_cmap('jet_r')
file_in = 'cutout_rings.v3.skycell.1631.047.stk.i.unconv.fits'
hdu_image = pyfits.open(file_in)

scale=0.25 
zeropnt=25. # see for details in Sneth et al. 2010 PASP 122, 1397 
name='chandra1'

file_outbt = 'ellipse_'+name+'.fits'
file_plt='plots_'+name+'2.png'
hdu_table = pyfits.open(file_outbt)
tab = hdu_table[1].data

exptime = hdu_image[0].header['EXPTIME']
x=tab.field('SMA')*scale
int=tab.field('INTENS')
eint=0.5*tab.field('INT_ERR')*np.sqrt(tab.field('NDATA'))
mag=(zeropnt-2.5*np.log10(int/scale**2)) + 2.5*np.log10(exptime)
mag_ep=(zeropnt-2.5*np.log10((int+eint)/scale**2))+ 2.5*np.log10(exptime)
mag_en=(zeropnt-2.5*np.log10((int-eint)/scale**2))+ 2.5*np.log10(exptime)
#mag2=tab.field('MAG')+2.5*np.log10(scale**2)
#mag_ep2=mag2+tab.field('MAG_LERR')-2.5*np.log10(scale**2)
#mag_en2=mag2-tab.field('MAG_UERR')+2.5*np.log10(scale**2)
y=mag
#plt.figure(1,figsize=(10,10))
#plt.subplot(321)
plt.plot(x,y,'.')
plt.plot(x,mag_ep,'-')
plt.plot(x,mag_en,'-')
#plt.plot(x,mag2,'.')
#plt.plot(x,mag_ep2,'-')
#plt.plot(x,mag_en2,'-')
plt.ylim(25.0,18.0)
#plt.xlim(0,85)
plt.xlabel("Radius, arcsec")
plt.ylabel("$\mu$, mag/arcsec$^2$")
'''
ellip=tab.field('ELLIP')
incl=np.rad2deg( np.arccos(1.0-ellip) )
plt.subplot(323)
plt.plot(x,incl,'.')
plt.ylim(30,80)
plt.xlabel("Radius, arcsec")
plt.ylabel("$i=arccos (1-\epsilon)$")

y=tab.field('ELLIP')
plt.subplot(322)
plt.plot(x,y,'.')
plt.ylim(0.0,1.0)
plt.xlabel("Radius, arcsec")
plt.ylabel("$\epsilon=1-b/a$")

y=tab.field('PA')
plt.subplot(324)
plt.plot(x,y,'.')
plt.ylim(-90,0)
plt.xlabel("Radius, arcsec")
plt.ylabel("$PA$")

y=tab.field('A4')*100
yy=tab.field('A4')/tab.field('SMA')/scale*100
plt.subplot(325)
plt.plot(x,y,'.')
plt.ylim(-3,3)
plt.xlabel("Radius, arcsec")
plt.ylabel("$a_4 \cdot 100$")
'''
plt.savefig(file_plt)
plt.close()
print '...done'