from pyraf import iraf
import pyfits
import numpy as np
import matplotlib.pyplot as plt
import os
iraf.stsdas()
iraf.analysis()
iraf.isophote()
#parameters========================================
#--------------------------------------------------
name='n1381'
filt=''
dir_data='/storage/katkov/obs/photometry/cfht/n1381/'
dir_output='./'
file_in =dir_data+'n1381_I2_megapipe.fits'
#file_msk=dir_data+'NGC2841.1.finmask_nonan.fits'
#output file names---------------------------------
file_out=dir_output+'out_'+name+filt+'.tab'
file_outbt=dir_output+'ellipse_'+name+filt+'.fits'
file_mod=dir_output+'model_'+name+filt+'.fits'
file_res=dir_output+'resid_'+name+filt+'.fits'
file_plt=dir_output+'plots_'+name+filt+'.png'
#--------------------------------------------------
# scale and zero point are fixed for all galaxies
scale=0.18580220
zeropnt=30.0 # see for details in Sneth et al. 2010 PASP 122, 1397 
#--------------------------------------------------
params={}
params['input']=file_in
#params['dqf']=file_msk
params['output']=file_out
#params['x0']=507.0
#iraf.ellipse.geompar.x0=507.0
iraf.ellipse.x0=1656.0
iraf.ellipse.y0=1683.0
iraf.ellipse.pa=-40.0
iraf.ellipse.minsma=0.0
iraf.ellipse.maxsma=500
iraf.ellipse.step=0.08#2.0
iraf.ellipse.ellip0=0.8
iraf.ellipse.sma0=50
iraf.ellipse.linear='no'
iraf.ellipse.recenter='yes'
iraf.controlpar.hcenter='no'
iraf.controlpar.hellip='no'
iraf.controlpar.hpa='no'
iraf.controlpar.maxit=400
iraf.samplepar.usclip=3
iraf.samplepar.lsclip=3
iraf.samplepar.nclip=50
#--------------------------------------------------
#==================================================
flist=[file_out,file_outbt,file_mod,file_res]
for file in flist:
    if os.path.exists(dir_output+file): os.remove(dir_output+file)
iraf.ellipse(**params)
iraf.stwfits(file_out,file_outbt)
iraf.bmodel(file_out,file_mod)
hdu_image = pyfits.open(file_in)
hdu_model = pyfits.open(file_mod)
image = hdu_image[0].data
model = hdu_model[0].data
hdu_resid = pyfits.PrimaryHDU(image-model)
hdu_resid.header = hdu_image[0].header
hdu_resid.writeto(file_res)

#plots isophote===================================
hdu_table = pyfits.open(file_outbt)
tab = hdu_table[1].data

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))
mag_ep=(zeropnt-2.5*np.log10((int+eint)/scale**2))
mag_en=(zeropnt-2.5*np.log10((int-eint)/scale**2))
#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(27.5,14.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(-45,-35)
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'