;**********************************************************************

;How to read to IDL and then display, the IRAF .tab file, created by the 
;'ellipse'-routine:

;1) Make an ASCII-file in IRAF: tables.ttools.tdump file.tab > file.tab.txt
 
;2) Replace INDEF-values in the ASCII-file by some unrealistic value
;   (replace string in emacs)


;3) Run this program, which transforms the ASCII-file to IDL-format and
;plots PA vs. radius and INC vs. radius. The PA in the plot is fixed so that
;negative values are changed to positive values.
;Esim: IDL> pae_plot,'file.tab.txt'
;*************************************************************************

;changed to subroutine
;counting of lines not anymore necessary


pro s4g_pae_plot,filetab,restore=restore,title0=title0,xrange=xrange,$
parange=parange,ellrange=ellrange,meanpa0=meanpa0,meanell0=meanell0,pix0=pix0,$
clean=clean,a40=a40,pyraf=pyraf,savefile0=savefile0,noplot=noplot

if(n_params() le 0) then begin
print,'pro pae_plot,filetab,restore=restore'
print,'reads filetab-file contain iraf isophote-measurements'
print,'and saves them to idl-savefile filetab+".save"'
print,'KEYWORDS:'
print,' /restore -> do not read, plot from previous savefile'
print,' xrange   = [rmin,rmax] radius-range (def=all)'
print,' ellrange = [emin,emax] ellipticity-range (def=0,1)'
print,' parange  = [pamin,pamax] pa-range (def=0,180 = all)'
print,' meanpa   = [rmin,rmax] mean pa  measurement range (def=xr)'
print,' meanell  = [rmin,rmax] mean ell measurement range (def=xr)'
print,' pix=value  pixelsize in arcsecs (def=1)'
print,' title    = plot title (def=filetab)' 
print,' /clean   -> remove INDEFs'
return
endif


filein=filetab
savefile=filetab+'.save'

if(keyword_set(savefile0)) then savefile=savefile0

title=filetab
if(keyword_set(title0)) then title=title0

if(keyword_set(restore)) then goto,old

;***********************************************************************
;Make and IDL save-file from the ASCII-file: 
;**********************************************************************

if(keyword_set(clean)) then begin
    job_edit2,filetab,'elljunk','INDEF','9999',sile=3,max=1000
;    job_edit2,'elljunk','copy','INDEF','9999',sile=3
;    job_edit2,'elljunk','copy','INDEF','9999',sile=3
;    job_edit2,'elljunk','copy','INDEF','9999',sile=3
;    job_edit2,'elljunk','copy','INDEF','9999',sile=3
 
    filetab='elljunk'
endif

close,1
openr,1,filetab

;N=7 ;numerorivien maara ;CHANGE
;temp=fltarr(40,7);taulukko,missa 40 saraketta ja 7 rivia;CHANGE

pyraf=1
if(not keyword_set(pyraf)) then begin
;read the initial comment lines (41 lines)
   line=''
   for i=1,41 do begin
      readf,1,line
;print,line
   endfor
endif

;prepare for maximum of 999 data lines
temp=fltarr(55,9999)
templine=fltarr(55)

;read data-lines
count=-1
while not eof(1) do begin
    readf,1,templine
;    print,'----------------'
;    print,templine
    count=count+1
    temp(*,count)=templine
endwhile

;clean extras + store at proper variable names
temp=temp(*,0:count)

SMA=temp(0,*)
INTENS=temp(1,*)
INT_ERR=temp(2,*)
PIX_VAR=temp(3,*)
RMS=temp(4,*)
ELLIP=temp(5,*)
ELLIP_ERR=temp(6,*)
PA=temp(7,*)
PA_ERR=temp(8,*)
X0=temp(9,*)
X0_ERR=temp(10,*)
Y0=temp(11,*)
Y0_ERR=temp(12,*)
GRAD=temp(13,*)
GRAD_ERR=temp(14,*)
GRAD_R_ERR=temp(15,*)
RSMA=temp(16,*)
MAG=temp(17,*)
MAG_LERR=temp(18,*)
MAG_UERR=temp(19,*)
TFLUX_E=temp(20,*)
TFLUX_C=temp(21,*)
TMAG_E=temp(22,*)
TMAG_C=temp(23,*)
NPIX_E=temp(24,*)
NPIX_C=temp(25,*)
A3=temp(26,*)
A3_ERR=temp(27,*)
B3=temp(28,*)
B3_ERR=temp(29,*)
A4=temp(30,*)
A4_ERR=temp(31,*)
B4=temp(32,*)
B4_ERR=temp(33,*)
NDATA=temp(34,*)
NFLAG=temp(35,*)
NITER=temp(36,*)
STOP=temp(37,*)
A_BIG=temp(38,*)
SAREA=temp(39,*)

psma=temp(51-1,*)
sb_cor=temp(46-1,*)
esb_cor=temp(53-1,*)


;    *  PSMA: radius along the semimajor axis, in arcseconds (column #51 in the ASCII table).
;    * SB_COR: surface brightness in AB mag/sq.arcsec, including the aperture correction to take into account the scattered light across the frame. SB_COR is also corrected for foreground galactic extinction (column #46 in the ASCII table).
;    * ESB_COR: uncertainty in SB, including poisson noise of the source and the error in the determination of the sky level (which results from small and large-scale background variations, as measured with local_bck_public_median) (column #53 in the ASCII table). 

;The _COR values should be used in place of the default output photometry columns of ellipse, like MAG, INTENS, etc. 


save,filename=savefile,SMA,INTENS,INT_ERR,PIX_VAR,RMS,ELLIP,ELLIP_ERR,PA,PA_ERR,X0,X0_ERR,Y0,Y0_ERR,GRAD,GRAD_ERR,GRAD_R_ERR,RSMA,$
MAG,MAG_LERR,MAG_UERR,TFLUX_E,TFLUX_C,TMAG_E,TMAG_C,NPIX_E,NPIX_C,A3,A3_ERR,B3,B3_ERR,A4,A4_ERR,B4,B4_ERR,NDATA,NFLAG,NITER,STOP,A_BIG,SAREA,psma,sb_cor,esb_cor

print,'saved to ',savefile,' with ', count+1,' datalines'



old:
;****************************************************************************
;Use the save-file for plotting PA and INC radial variations
;*********************************************************************
 
restore,savefile  ;Takes the save-file to active IDL memory


if(keyword_set(noplot)) then goto,noplot
color_old=!p.color ; determines the type of plot
!p.multi=[0,1,2]

pix=1.0 ;pixel size in arcseconds
if(keyword_set(pix0)) then pix=pix0




;ranges
xr=[0,max(sma*pix)]
if(keyword_set(xrange)) then xr=xrange

meanpa=xr
meanell=xr
if(keyword_set(meanpa0)) then meanpa=meanpa0
if(keyword_set(meanell0)) then meanell=meanell0

par=[0,180]
if(keyword_set(parange)) then par=parange

;limit PA to r= 20-200: when necessary chance the lower limit
ind=where(pa gt par(1),count) 
if(count ge 1) then pa(ind)=pa(ind)-180.
ind=where(pa lt par(0),count) 
if(count ge 1) then pa(ind)=pa(ind)+180.


ellr=[0,1.]
if(keyword_set(ellrange)) then ellr=ellrange


;*******
; PA
;*******

;Calculate the mean and the standard deviation for PA in a certain radius interval
ind2=where(sma*pix gt meanpa(0) and sma*pix lt meanpa(1)) 
mean_pa=mean(pa(ind2))
sigma=stdev(pa(ind2))
addtitle='   <PA>= '+string(mean_pa,'(f6.1)')



!p.color=1
plot,sma*pix,pa,psym=4,xtitle='RAD(arcsec)',ytitle='PA',xr=xr,yr=par,title=title+addtitle ;CHANGE
errplot,sma*pix,pa+(pa_err),pa-(pa_err)
!p.color=color_old



print,'PA CALCULATED BETWEEN:',meanpa,' from ',n_elements(ind2),' datapoints'
print,'<PA>, STDEV(PA)',mean_pa,sigma

oplot,meanpa,[1,1]*mean_pa,lines=2



;*******
; ELL
;********

;Calculate the mean ellipticity and its standard deviation in a certain radius interval
ind2=where(sma*pix gt meanell(0) and sma*pix lt meanell(1)) 
mean_ell=mean(ellip(ind2))
sigma=stdev(ellip(ind2))
addtitle='   <ELL>= '+string(mean_ell,'(f6.3)')


!p.color=1
plot,sma*pix,ellip,psym=4,xtitle='RAD(arcsec)',ytitle='ellipticity',xr=xr,yr=ellr,ys=1,title=title+addtitle
errplot,sma*pix,ellip+(ellip_err),ellip-(ellip_err)
!p.color=color_old



print,'ELL CALCULATED BETWEEN:',meanell,' from ',n_elements(ind2),' datapoints'
print,'<ELL>, STDEV(ELL)',mean_ell,sigma

oplot,meanell,[1,1]*mean_ell,lines=2


!p.multi=0
if(keyword_set(a40)) then begin
nwin_hs
plot,sma(1:*)*pix,a4(1:*),xtitle='RAD(arcsec)',ytitle='a4',xr=xr,title=title,yr=[-1,1]*a40,psym=4
;oplot,sma(1:*)*pix,b4(1:*),lines=2
endif


;***********************************
;ADU
;***********************************
;;;;;nwin
;;;;;restore,'esim_ellfit.save'
;;;;;plot,sma*pix,intens,psym=4,xtitle='RAD(arcsec)',ytitle='ADU'
;;;;;errplot,sma,intens+(int_err),intens-(int_err),xr=[0,25],yr=[0,3000]



print,'pix',pix

noplot:
end