pro euvi_corrplot,image_pair,para,nsm,fov2,dcont,outfile,ct,fov
;+
; Project     : STEREO
;
; Name        : EUVI_CORRPLOT 
;
; Category    : Graphics 
;
; Explanation : This routine displays a selected FOV of a stereoscopic image pair
;		and the positions (3x3 pixel areas) of correlated point-like sources
;		stored in OUTFILE.
;
; Syntax      : IDL>euvi_corrplot,image_pair,para,nsm,fov2,dcont,outfile,ct,fov
;
; Inputs      : images_pair(2,*,*) : stereo image pair B,A
;		para               : index structure with image parameters 
;		nsm		   : smoothing parameter for highpass filter
;		fov2		   : [i1,j1,i2,j2] small field-of-view in image pixels 
;               dcont              : lowest contour level (DN), dcont*2^i (logarithmic higher levels)
;		outfile		    ; filename containing coordinates of point-like sources
;               ct                 : 0, ..., 38 IDL color table
;		fov		   = [i1,j1,i2,j2] context field-of-view in image pixels 
;
; Outputs     : postscript file names 'corrplot_a_col.ps'
;
; History     : 23-Feb-2007, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;________________________________________________________________
;			SUBIMAGES  				|
;________________________________________________________________
rsun_m  =0.696*1.e9
cdelt1	=para.cdelt1
cdelt2	=para.cdelt2
crpix1	=para.crpix1
crpix2	=para.crpix2
da	=para.da/rsun_m
db	=para.db/rsun_m
rsun   	=para.rsun
a_sep	=para.sep_deg*(!pi/180.)
hmax	=0.1
rad_pix	=rsun/cdelt1
r_max   =1.0+hmax
i1	=fov2(0)
j1	=fov2(1)
i2	=fov2(2)
j2	=fov2(3)
nx	=(i2-i1+1)
ny	=(j2-j1+1)
r_n	=1.0
im	=(i1+i2)/2.
jm	=(j1+j2)/2.
rm	=sqrt((im-crpix1)^2+(jm-crpix2)^2)
im_b	=im
jm_b	=jm
if (rm lt rad_pix) then euvi_stereo_rot,im,jm,para,r_n,im_b,jm_b,zm_b
i1_b	=im_b-(im-i1)
i2_b	=i1_b+(i2-i1)
xrange1 =[i1  ,i2  ]
xrange2 =[i1_b,i2_b]
images_fov=fltarr(nx,ny,2)		;2 spacecraft(A,B)  
image_a	=float(reform(image_pair(1,i1:i2,j1:j2)))
image_b	=float(reform(image_pair(0,i1_b:i2_b,j1:j2)))
for it=0,1 do begin
 if (it eq 0) then image0 =image_a
 if (it eq 1) then image0 =image_b
 if (nsm ge 3) then images_fov(*,*,it)=(image0-smooth(image0,nsm))
endfor

;________________________________________________________________
;			IMAGES  				|
;________________________________________________________________
i1_	=fov(0)
j1_	=fov(1)
i2_	=fov(2)
j2_	=fov(3)
nx_	=(i2_-i1_+1)
ny_	=(j2_-j1_+1)
im_	=(i1_+i2_)/2.
jm_	=(j1_+j2_)/2.
rm_	=sqrt((im_-crpix1)^2+(jm_-crpix2)^2)
im_b	=im_
jm_b	=jm_
if (rm_ lt rad_pix) then euvi_stereo_rot,im_,jm_,para,r_n,im_b_,jm_b_,zm_b_
i1_b_	=im_b_-(im_-i1_)
i2_b_	=i1_b_+(i2_-i1_)
xrange1_=[i1_  ,i2_ ]
xrange2_=[i1_b_,i2_b_]
images_fov_=fltarr(nx_,ny_,2)		;2 spacecraft(A,B)  
image_a_=float(reform(image_pair(1,i1_:i2_,j1_:j2_)))
image_b_=float(reform(image_pair(0,i1_b_:i2_b_,j1_:j2_)))
for it=0,1 do begin
 if (it eq 0) then image0_=image_a_
 if (it eq 1) then image0_=image_b_
 images_fov_(*,*,it)=image0_
endfor

;________________________________________________________________
;		GREYSCALE PLOT OF FULL FOV			|
;________________________________________________________________
for io=0,3,3 do begin
form    =1     ;landscape format
char    =1     ;smaller character size
plotname='corrplot'
fignr	 ='_a'
unit    =0
ref     =''
loopcol	=150
loopcol2= 50
sc	=['A','B']

fig_open,io,form,char,fignr,plotname,unit
loadct,ct
aspect_	=0.7*float(ny_)/float(nx_)
panel1_	=[0.55,0.65,0.95,0.65+0.4*aspect_]
panel2_	=[0.10,0.65,0.50,0.65+0.4*aspect_]
for it=0,1 do begin
 if (it eq 0) then panel=panel1_ 
 if (it eq 1) then panel=panel2_ 
 x1_	=panel(0)
 y1_	=panel(1)
 x2_	=panel(2)
 y2_	=panel(3)
 !p.position=[x1_,y1_,x2_,y2_]
 if (it eq 0) then !x.range=xrange1_
 if (it eq 1) then !x.range=xrange2_
 !y.range=[j1_,j2_]		
 !x.style=1
 !y.style=1
 !x.title='EW [pixels]'
 if (it eq 1) then !y.title='NS [pixels]'
 if (it eq 0) then !y.title=''
 plot,[0,0],[0,0]
 if (io eq 0) then begin
  nxw  =long(!d.x_vsize*(x2_-x1_)+0.5)
  nyw  =long(!d.y_vsize*(y2_-y1_)+0.5)
  zz   =congrid(images_fov_(*,*,it),nxw,nyw)    ;pixel scaling
 endif
 if (io ne 0) then zz=images_fov_(*,*,it)                ;high resolution
 c0     =median(zz)
 statistic,zz,avg,dev
 c1     =c0-dev*2
 c2     =c0+dev*2
 print,'z-range of greyscale = c1,c0,c2=',c1,c0,c2
 tv,bytscl(c2-zz,min=c1,max=c2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal
 xyouts_norm,0.02,0.90,sc(it),2
 if (it eq 0) then oplot,[i1,i2,i2,i1,i1],[j1,j1,j2,j2,j1],thick=3
 if (it eq 1) then oplot,[i1_b,i2_b,i2_b,i1_b,i1_b],[j1,j1,j2,j2,j1],thick=3
 dy	=(j2_-j1_)*0.05/(0.4*aspect_)
 if (it eq 0) then oplot,[i1,i1_],[j2,j1_-dy],/noclip
 if (it eq 0) then oplot,[i2,i2_],[j2,j1_-dy],/noclip
 if (it eq 1) then oplot,[i1_b,i1_b_],[j2,j1_-dy],/noclip
 if (it eq 1) then oplot,[i2_b,i2_b_],[j2,j1_-dy],/noclip
 !noeras=1
endfor

;________________________________________________________________
; 		CONTOUR PLOTS OF SUB-FOV			|
;________________________________________________________________
aspect	=0.7*float(ny)/float(nx)
panel1	=[0.55,0.6-0.4*aspect,0.95,0.6]
panel2	=[0.10,0.6-0.4*aspect,0.50,0.6]
for it=0,1 do begin
 if (it eq 0) then panel=panel1 
 if (it eq 1) then panel=panel2 
 x1_	=panel(0)
 y1_	=panel(1)
 x2_	=panel(2)
 y2_	=panel(3)
 image  =images_fov(*,*,it)
 !p.position=[x1_,y1_,x2_,y2_]
 if (it eq 0) then !x.range=xrange1
 if (it eq 1) then !x.range=xrange2
 !y.range=[j1,j2]		
 !x.style=1
 !y.style=1
 plot,[0,0],[0,0]
 xx	=!x.range(0)+findgen(nx) 
 yy	=!y.range(0)+findgen(ny) 
 zz	=dcont*2^findgen(10)
 contour,image,xx,yy,level=zz
 xyouts_norm,0.02,0.94,sc(it),2
 !noeras=1
endfor		;stereo image pair
oplot,[i1,i1],[j1,j2],thick=3

;________________________________________________________________
;		POSITIONS OF CORRELATED POINTS			|
;________________________________________________________________
file_exist=findfile(outfile,count=nfiles)
if (nfiles eq 0) then stop,'No previous outfile found';
readcol,outfile,iloop_,ix_,iy_,iz_,ih_pix,dh_pix,ih_mm,dh_mm,nsm_
print,'Statistics of height measurements:'
statistic,ih_mm
if (n_elements(iloop_) lt 1) then stop,'outfile empty: ',outfile
nloop	=long(max(iloop_))+1
print,'Datafile with loop coordinates found  LOOPFILE=',outfile 
print,'Number of loops found                NLOOP   ='+string(nloop,'(i3)')
di	=1
ni	=di*2+1

!p.position=panel1		;STEREO A projection of field lines
!x.range=xrange1
plot,[0,0],[0,0]
loadct,3
for iloop=0,nloop-1 do begin 
 ind	=where(iloop_ eq iloop,np) 
 xp	=ix_(ind)
 yp	=iy_(ind)
 zp	=iz_(ind)
 euvi_stereo_rot,xp,yp,para,1.,xb,yb,zb

 oplot,[xp,xp],[yp,yp],psym=1,symsize=2,color=loopcol
 xyouts,xp,yp,string(iloop+1,'(I3)'),size=char,color=loopcol
 ip	=long(xp+0.5)
 jp	=long(yp+0.5)
 oplot,ip(0)-di+[0.,ni-1,ni-1,0.,0.],jp(0)+[-1.,-1.,1.,1.,-1.],color=loopcol
endfor

loadct,ct
!p.position=panel2 		;STEREO B projection of field lines
!x.range=xrange2
plot,[0,0],[0,0]
loadct,3
for iloop=0,nloop-1 do begin 
 ind	=where(iloop_ eq iloop,np) 
 xp	=ix_(ind)
 yp	=iy_(ind)
 zp	=iz_(ind)
 euvi_stereo_rot,xp,yp,para,1.,xb,yb,zb

 oplot,[xb,xb],[yb,yb],psym=1,symsize=2,color=loopcol
 xyouts,xb,yb,string(iloop+1,'(I3)'),size=char,color=loopcol  
 ip2	=long(xb+0.5)
 jp2	=long(yb+0.5)
 oplot,ip2(0)-di+[0.,ni-1,ni-1,0.,0.],jp2(0)+[-1.,-1.,1.,1.,-1.],color=loopcol
endfor

if (io ge 1) then fig_close,io,fignr,ref
endfor
end