;=======================================================================
; NAME :	magcont (procedure)
; PURPOSE :
;	display vector mag.field as contour
; CATEGORY :
;	idl/lib/flare
; CALLING SEQUENCE :
;	magcont,img,q,u,v,h=h,/spline,/order,save_ps=filename,/seq, $
;		magnify=magnify,bin=bin,/novect,/printout
;  or	magcont,img,Bx,By,Bl,h=h,/cal,h=h,/spline,/order,save_ps=filename, $
;		/seq,magnify=magnify,bin=bin,/nocont,/printout
; INPUTS :
;	img     : image (*,*) or i,q,u,v as (*,*,4)
;       q,u,v   : polarization degree*10000
;       Bx,By,Bl: long. & trans. mag fields, gauss*10
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	h       : header array
;       spline  : if set, draw spline contour
; 	order   : if set, up-down direction is opposite
; 	save_ps=filename  : save  result in post script file
; 	/printout : printout the result
; 	/cal	: filenamefor calibrated data	Bl=>q,  Bx=>u, By=>v
; 	/seq	: display seq.no
; 	magnify : magnification factor
; 	bin	: binning factor
; 	novect  : no linear pol.
; OUTPUTS :
;	none
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	none
; RESTRICTIONS :	none
; PROCEDURE :
;	if h(1).value eq 'Bl', it displays as calibrated data (/cal)
; MODIFICATION HISTORY :
;	K.I. '92/02/29
;	K.I. '92/10/20
;	K.I. '92/10/24	display cal.ver
;	K.I. '92/10/28	change order of Bl & bx,by
;	K.I. '92/11/19	change v sign
;	K.I. '92/11/29	h.dx -> bin
;	K.I. '93/05/16	lstyle, lthick changed
;	K.I. '93/11/15	cmpfact=max(...)
;	K.I. '93/12/21	rebin -> congrid
;----------------------------------------------------------------------
pro magcont,img0,b1,b2,b3,h=h,spline=spline,order=order,save_ps=save_ps, $
	cal=cal,printout=printout,seq=seq,magnify=magnify,bin=bin,	$
	novect=novect,nospot=nospot

;----------  parameters  --------------
;pix1=0.6738     ; arc sec. for 1 pixel
pix1=pix_size('ft1')	; arc sec. for 1 pixel
cmpfact=2	; compress factor for drawing contour
Blmax=1500.	; max. of long. mag.	(Gauss)
Btmin=150.	; min. of tang. mag.	(Gauss)
pcmax=0.2	; max. of circular pol. for display
plmin=0.002	; min. of linear pol. for display
nstep=10	; step of linear pol. bar plot
Bfact=10.	; ratio-fact  q,u,v -> q,u,v/Bfact (gauss)
rfact=10000.	; ratio-fact  q,u,v -> q,u,v/rfact
imgmin=0.2	; minimum level to display pol.
arrow=1		; display arrow for Bt

;------ initialise keywords ------
if keyword_set(bin) then bin=bin else bin=1
if keyword_set(magnify) then begin
	magnify=magnify 
	print,'magnification factor : ',magnify
endif else begin
	magnify=1
endelse
if keyword_set(spline) then begin
	spline=1
	c_labels=[0]
endif else begin
	spline=0
	c_labels=0
endelse
if keyword_set(novect) then novect=1 else novect=0
if keyword_set(cal) then cal=1 else cal=0
if keyword_set(h) then begin
	if h(1).value eq 'Bx' then begin 
		cal=1
		print,'It is calibrated data...'
	endif
	if strmid(h(0).dx,0,1) eq ' ' then dx=0 $
	else dx=float(h(0).dx)
	if dx ne 0. then bin=dx/pix1
endif
if keyword_set(printout) then begin
	scale=25
	print,'result will be printed out with a scale of',scale
	set_plot,'ps'
	device,filename='~/magcont.ps'
endif else begin
	scale=1
endelse
if keyword_set(save_ps) then begin
	s=size(save_ps)
	if s(s(0)+1) eq 7 then ps_file=save_ps else ps_file='~/magcont.ps'
	scale=25
	print,'result will be saved in PostScript file'
	set_plot,'ps'
	device,filename=ps_file,xsize=17.,ysize=17.
endif
;------ initialise array ------
s=size(img0)
nx=s(1)	& ny=s(2)
if s(0) eq 2 then img=img0 else img=img0(*,*,0)
if n_params() eq 2 then b3=b1
if n_params() eq 1 then begin	; case only img0 given
	if s(0) eq 2 then begin	; case img0 is 2-dimensional
		img=img0
	endif else begin	; case img0 is 3-dimensional
		img=img0(*,*,0)
		if s(3) eq 4 then begin
			b1=img0(*,*,1)
			b2=img0(*,*,2)
			b3=img0(*,*,3)
		endif else begin
			b3=img0(*,*,1)
		endelse
	endelse
endif
cmpfact=max([nx/256,1])

;------- initialize window & color table ----------
x0=50
y0=50
wwx=nx*magnify+2*x0	;expected window size
wwy=ny*magnify+y0+55
if !d.window ne -1 then wno=!d.window else wno=0
if not keyword_set(printout) and not keyword_set(save_ps) then begin
    if (!d.window eq -1) or (wwx gt !d.x_size+1) or (wwy gt !d.y_size+1) then $
	window,wno,xsize=wwx,ysize=wwy
endif
loadct,0
tvlct,0,0,0,!d.table_size-1	; color of frame
tvlct,0,0,0,!d.table_size-2	; color of contour
tvlct,255,255,255,0		; back ground
!order=1
if keyword_set(order) then !order=0
pos=[scale*x0,scale*y0,scale*(x0+nx*magnify-1),scale*(y0+ny*magnify-1)]

;---------------  display spot position  -----------------
erase
;img1=rebin(img,nx/cmpfact,ny/cmpfact)
img1=congrid(img,nx/cmpfact,ny/cmpfact,/interp)
if !order eq 1 then mirror,img1
iav=(total(img1(*,0))+total(img1(*,ny/cmpfact-1)))/(nx/cmpfact)/2
if not keyword_set(nospot) then begin
    ;s=size(img1)
    ;temp=replicate(min(img1),s(1)+2,s(2)+2)
    ;temp(1,1)=img1
    ;dpos=[-1,-1,1,1]*scale
    contour,img1/float(iav),	$
	pos=pos,/dev,xstyle=1+4,ystyle=1+4,	$
	xrange=[0.,nx/cmpfact],yrange=[0.,ny/cmpfact],	$
	level=[0.5,0.9],/noerase,path_filename='~/path.dat',back=255,color=1
    if keyword_set(printout) or keyword_set(save_ps) then $
	polycontour,'~/path.dat',color_index=[230,250,255],delete_file=1 $
    else	polycontour,'~/path.dat',color_index=[180,220,0],delete_file=1
endif

;----------------  display longitudinal fields  ---------------
    if keyword_set(cal) then begin 
	lfact=0.2	; factor for length of linear pol. bar
	;r=float(smooth(rebin(b3,nx/cmpfact,ny/cmpfact),2))/Bfact
	r=float(smooth(congrid(b3,nx/cmpfact,ny/cmpfact,/interp),2))/Bfact
	levp=[50,200,500,1000,1500]
	nn=n_elements(levp)
	levstr='level=!9+!3'
	for i=0,nn-2 do	$
	levstr=levstr+strcompress(string(levp(i),format='(i)'),/remove_all)+','
	levstr=levstr+string(levp(nn-1),format='(I4)')+' Gauss'
    endif else begin
	lfact=50.	; factor for length of linear pol. bar
	;r=float(smooth(rebin(-b3,nx/cmpfact,ny/cmpfact),2))/rfact
	r=float(smooth(congrid(-b3,nx/cmpfact,ny/cmpfact,/interp),2))/rfact
	levp=[0.01,0.03,0.05,0.1,0.15]
	nn=n_elements(levp)
	levstr='level = !9+!3'
	for i=0,nn-2 do levstr=levstr+string(levp(i),format='(f4.2)')+', '
	levstr=levstr+string(levp(nn-1),format='(f4.2)')+'  pol.degree'
    endelse

    levn=levp
    for i=0,nn-1 do levn(i)=-levp(nn-i-1)
    lev=[levn,levp]
    lstyle=intarr(n_elements(lev))
    lstyle(where(lev gt 0))=0	; solid
    ;lstyle(where(lev lt 0))=3	; dot
    lstyle(where(lev lt 0))=0	; dot
    lthick=intarr(n_elements(lev))
    lthick(where(lev gt 0))=3.0
    lthick(where(lev lt 0))=1.0
    levn=levp
    for i=0,nn-1 do levn(i)=-levp(nn-i-1)

    ;img=rebin(img,nx/cmpfact,ny/cmpfact)
    img=congrid(img,nx/cmpfact,ny/cmpfact,/interp)
    if !order eq 1 then begin
	mirror,r
	mirror,img
    endif
    imgmax=max(img)
    sky=where(img lt imgmax*imgmin, nsky)
    if(nsky ne 0l) then r(sky) = 0.
    contour,r,	$		
	pos=pos,/dev,xstyle=1+4,ystyle=1+4,	$
	xrange=[0.,nx/cmpfact],yrange=[0.,ny/cmpfact],		$
	level=lev,c_linestyle=lstyle,c_thick=lthick,/noerase,	$
	back=255,color=1,	$
	spline=spline,c_labels=c_labels	; use follow method and surpress labels

    if !order eq 1 then mirror,img

;----------------  display transversal fields  ---------------
if n_elements(b2) gt 1  then begin
    plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase,	$
	xrange=[0.,nx*magnify],yrange=[0.,ny*magnify],xstyle=1+4,ystyle=1+4

    nstep=nstep/magnify
    if not keyword_set(novect) then begin
       if keyword_set(cal) then begin 
	  for j=nstep,ny-nstep,nstep do begin
	  for i=nstep,nx-nstep,nstep do begin
	    if img(i/cmpfact,j/cmpfact) gt imgmax*imgmin then begin
		Bt1=sqrt((float(b1(i,j))^2+float(b2(i,j))^2)>1)/Bfact
		plen=sqrt((Bt1-Btmin)>1) 
		if Bt1 gt Btmin then begin
			th=atan(b2(i,j),b1(i,j))
			x=plen*cos(th)*lfact/magnify
			y=plen*sin(th)*lfact/magnify
			if !order eq 0 then begin
				x1=[i-x,i+x]
				y1=[j-y,j+y]
			endif else begin
				x1=[i-x,i+x]
				y1=[ny-j-1+y,ny-j-1-y]
			endelse
			oplot,x1*magnify,y1*magnify,color=1
			if arrow eq 1 then begin
				dth=25./180.*!pi
				len=sqrt((x1(1)-x1(0))^2+(y1(1)-y1(0))^2)
				larrow=min([len/4.,5.])
				ddx=larrow*(x1(0)-x1(1))/len
				ddy=larrow*(y1(0)-y1(1))/len
				xx1=x1(1)+ddx*cos(dth)+ddy*sin(dth)
				yy1=y1(1)-ddx*sin(dth)+ddy*cos(dth)
				xx2=x1(1)+ddx*cos(-dth)+ddy*sin(-dth)
				yy2=y1(1)-ddx*sin(-dth)+ddy*cos(-dth)
				oplot,  [xx1,x1(1),xx2,xx1]*magnify,	$
					[yy1,y1(1),yy2,yy1]*magnify,	$
					color=1
			endif
		endif
	    endif
 	  endfor
 	  endfor
	  message1='Bx & By not defined. Only Bl map is displayed!'
	endif else begin
	  for j=nstep,ny-nstep,nstep do begin
	  for i=nstep,nx-nstep,nstep do begin
	    if img(i/cmpfact,j/cmpfact) gt imgmax*imgmin then begin
		pl1=sqrt((float(b1(i,j))^2+float(b2(i,j))^2)>1)/rfact
		plen=sqrt((pl1-plmin)>0.0001)
		if pl1 gt plmin   then begin
			;print,plen,plmin
			th=atan(b2(i,j),b1(i,j))/2.
			x=cos(th)*plen*lfact/magnify
			y=-sin(th)*plen*lfact/magnify
			if !order eq 0 then begin
				x1=[i-x,i+x]
				y1=[j-y,j+y]
			endif else begin
				x1=[i-x,i+x]
				y1=[ny-j-1+y,ny-j-1-y]
			endelse
			oplot,x1*magnify,y1*magnify,color=1
		endif
	    endif
	  endfor
	  endfor
	  message1='Q & U not defined. Only V map is displayed!'
	endelse

    endif else begin
	;print,message1
    endelse
endif

;--------------------   frame display   ---------------------
if keyword_set(h) then begin
	if n_elements(h) gt 1 then begin
		timeint,h,time1,time2,date1,date2
		time_sys=h(0).time_sys
		g_x=h(0).g_x
		g_y=h(0).g_y
		seqno=h(0).seqno
		calib='Calib.'+string(h(1).correct,format='(I3)')
		hql=h(0)
	endif else begin
		date1=h.date1
		time1=h.time1
		date2=h.date2
		time2=h.time2
		time_sys=h.time_sys
		g_x=h.g_x
		g_y=h.g_y
		seqno=h.seqno
		hql=h
	endelse
	if keyword_set(save) then begin
		hql.nbyt=1
		s=size(img)
		hql.nx=s(1)	& hql.ny=s(2)
		hql.date1=date1
		hql.date2=date2
		hql.time1=time1
		hql.time2=time2
		hql.nset=1
		hql.value='map'
		wshow
		rd=tvrd(x0,y0,wx,wy)
		nkrsave,save,hql,rd,/nocheck
		print,'display saved in '+save
	endif
	if time_sys eq 'JST' then begin
		jst2ut,date1,time1
		jst2ut,date2,time2
	endif
	title1='!N!3Solar Flare Telescope (MTK) : vector magnetic field'
	title2=date1+'  '+strmid(time1,0,8)+'-'+strmid(time2,0,8)+' UT' $
		+'    '+str_gpos(g_x,g_y)
endif else begin
	title1=''
	title2='!N!3Solar Flare Telescope (MTK) : vector magnetic field'
endelse

nx=nx*magnify
ny=ny*magnify
pix1=pix1*bin/magnify
if keyword_set(printout) or keyword_set(save_ps) then begin
	xyouts,scale*x0,scale*(y0+ny+47),title1,/dev,charsize=1.3,color=1
	xyouts,scale*x0,scale*(y0+ny+27),title2,/dev,charsize=1.2,color=1
	if keyword_set(h) then $
	xyouts,scale*x0,scale*(y0+ny+7),levstr+',   '+calib,/dev,charsize=1.2
	if keyword_set(seq) then $
	xyouts,scale*(x0+nx-40),scale*(y0+ny+0),	$
	'('+string(seqno,format='(i3)')+')',/dev,charsize=1.2
endif else begin
	xyouts,scale*x0,scale*(y0+ny+35),title1,/dev,charsize=1.7,color=1
	xyouts,scale*x0,scale*(y0+ny+20),title2,/dev,charsize=1.4,color=1
	if keyword_set(h) then $
	xyouts,scale*x0,scale*(y0+ny+5),levstr+',   '+calib,/dev,charsize=1.4
	if keyword_set(seq) then $
	xyouts,scale*(x0+nx-35),scale*(y0+ny+5),	$
	'('+string(seqno,format='(i3)')+')',/dev,charsize=1.4
endelse
xrange=[0.,pix1*(nx-1)]
yrange=[0.,pix1*(ny-1)]
tickv=[0,100,200,300,400,500,600]
plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase,color=1,	$
	xrange=xrange,yrange=yrange,xstyle=1,ystyle=1,	$
        xticks=fix(nx*pix1/100.),xtickv=tickv,xminor=10,xticklen=0.03, $ 
        yticks=fix(ny*pix1/100.),ytickv=tickv,yminor=10,yticklen=0.03, $ 
	xtitle='arc sec.',ytitle='arc sec.'

if keyword_set(printout) then begin
	device,/close
	set_plot,'x'
	spawn,'lpr -h ~/magcont.ps'
	spawn,'rm ~/magcont.ps'
endif
if keyword_set(save_ps) then begin
	device,/close
	set_plot,'x'
	print,'mag. drawing is saved in '+ps_file
endif

end