;shamelessly lifted from show_stokes.pro to go to a postscript rather than
; x-window...

; Filename only?

dimen=[512,512,40]
dir='/data/SOLMAG/DATA/IVM/NA-SURVEY/'
inim=dir+'/20031026/sodium/IMG0_20031026.1821_df_rf_ic_pl_sc_ds_sub_av_st_ca_c' ; AR 10486
;inim=dir+'/20031026/sodium/IMG0_20031026.1832_df_rf_ic_pl_sc_ds_sub_av_st_ca_c' ; AR 10484, flare residual emission.  Pretty cool, though.
; below: testing on palila
;inim='/export/home/leka/data/ivm/10486/NA/IMG0_20031029.1846_dpa_df_re_ic_pl_sc_ds_av_st_ca_unb_sw'
OPENR, lun, inim[0], /GET_LUN
im = FLTARR(dimen[0],dimen[1],4,dimen[2])
READU,lun,im
CLOSE,lun
FREE_LUN,lun

sz = SIZE( im )
ncol = dimen[0]
nrow = dimen[1]
ns = 4
nw = dimen[2]

; Location to plot column number
xt = (0.5 * ncol) / (ncol + 4 * 4 * nw)
yt = 0.95

; Avoid bad edges if this is full frame IVM
IF( (ncol EQ nrow) AND (ncol EQ 128 OR ncol EQ 256 OR ncol EQ 512) ) THEN BEGIN
	ex = 5 * ncol/256
	ey = 30 * nrow/256
; Get the mask
	mask = IVM_MASK(512/ncol)
   ENDIF ELSE BEGIN
	ex = 0
	ey = 0
	mask = REPLICATE(1., ncol, nrow)
ENDELSE

; Get average line profile
line = FLTARR(nw)
FOR i=0,nw-1 DO BEGIN
	line(i) = TOTAL( im(ex:ncol-ex-1,ey:nrow-ey-1,0,i) )
ENDFOR
line = SHIFT(line, 1)

; What wavelength to display?
ncore = WHERE( line EQ MIN(line) )
ncontin = WHERE( line EQ MAX(line) )
nwing = (ncore + ncontin)/2
;ndisp=ncore ; if want the other, ncontin, set here
ndisp=nwing ; if want the other, ncontin, set here
;IF( KEYWORD_SET(contin) ) THEN ndisp = ncontin ELSE ndisp = ncore
;IF( KEYWORD_SET(wing) ) THEN ndisp = nwing

; What Stokes parameter to display?
	nstokes = 0

; Gather image and mark the column to display spectra for
tim = REFORM( im(*,*,nstokes,ndisp) )
tim = tim * mask
mini = MIN( tim(ex:ncol-ex-1,ey:nrow-ey-1) )
maxi = MAX( tim(ex:ncol-ex-1,ey:nrow-ey-1) )

mag=1
xc = 246 ; good for the 10486 one
yc = 0
image = BYTSCL( tim, MIN=mini, MAX=maxi )
	nx = xc/mag
	ny = yc/mag

;image(nx,*) = BYTSCL( tim(nx,*) * 0.9, MIN=mini, MAX=maxi )
image(nx-1:nx+1,*) = 10.


set_plot,'ps'
device,file='PAPER_IVM_showstokes_eg.eps',xsi=20,ysi=10,bits=8,/encap,/times
!x.omargin=[2,1]
!x.margin=[2,1]
!y.omargin=[0,0]
!y.margin=[2,1]

image2_cont,image,image,lev=[-99],/aspect

minp=[2200,-0.1,-0.1,-0.3]
maxp=[34000.,0.1,0.1,0.3]
; Now gather the profiles
cmask = REBIN( mask(nx,*),nw,nrow )
FOR i=0,ns-1 DO BEGIN
	pim = TRANSPOSE( REFORM(im(nx,*,i,*) ) ) 
	pim = SHIFT(pim, 1, 0) * cmask
;	minp = MIN( pim(*,ey:nrow-ey-1) )
;       maxp = MAX( pim(*,ey:nrow-ey-1) )
; those are good to find good scaling, but now I know what we want for this particular image...
;	IF( i GT 0 ) THEN BEGIN
;		mp = ABS(minp) > ABS(maxp)
;		minp = -mp
;		maxp = mp
;	ENDIF 
;	TV, BYTSCL( REBIN( pim, mag*4*nw, mag*nrow ), $
	TV, BYTSCL( congrid( pim, mag*4*nw, mag*nrow ), $
;		 MIN=minp, MAX=maxp ),mag*(ncol+4*nw*i),0
		 MIN=minp[i], MAX=maxp[i] ),0.5+i*0.125,0.06,ysize=0.92,/norm
;		 print,minp,maxp
ENDFOR
xyouts,0.561,0.93,'I',/norm,col=220,size=1.5
xyouts,0.6875,0.93,'Q',/norm,col=220,size=1.5
xyouts,0.8125,0.93,'U',/norm,col=220,size=1.5
xyouts,0.9375,0.93,'V',/norm,col=220,size=1.5



device,/close
set_plot,'x'
reset

END