;+
; NAME :    HR_Mag.pro (procedure)
; PURPOSE :
;   Make a vector magnetogram from calibrated Huairou magnetic field data saved in FITS
;	format file.
; CATEGORY :
;   idl/huairou
; CALLING SEQUENCE :
;   hr_mag,fname,back=back,saveps=saveps,dir=dir,ext=ext,magnify=magnify,$
;		nstep=nstep,level=level,xrange=xrange,yrange=yrange,xtitle=xtitle,$
;		ytitle=ytitle,pix=pix,/btonly,/blonly,/btcont,/blackgnd
; INPUTS :
;   fname: FITS file name of calibrated data without extension
; OPTIONAL INPUT PARAMETERS :
;   none
; OUTPUTS :
;   Show the magnetogram on the screen or save as a .PS file if set /saveps
; KEYWORD PARAMETERS :
;	back = : a file name in Huairou data format OR a variable containing the background
;	         image (e,g, photospheric spot image); OR just set this keyword to show the
;	         spot image in the read-in FITS file as the background.
;	dir=dir: directory where the observed data files are stored
;	ext=ext: extension of the file names. default: '.dat'
;	saveps : if set, save the aftermath in a PostScript file
;	btonly : if set, show transversal field component only
;	blonly : if set, show longitudinal field component only
;	nstep  : step increase to plot transversal field. default: nstep=10
;	magnify: magnifying factor when showing the magnetogram. default: =3
;	btcont : if set, show the transversal field as contours
;	level  : contour level. default: [200,500,800,1200,1600,2000,2500,3000] for Bt contour
;				and [40,80,160,320,640,960,1280,1600,1920,2240,2560,2880] for Bl contour
;	blackgnd: if set, set the background to black
;	xrange : Range for the axis X when plotting. default: auto determined
;	yrange : Range for the axis Y when plotting. default: auto determined
;	xtitle : title for axis X. default: 'Distance (arcsec)'
;	ytitle : title for axis Y. default: 'Distance (arcsec)'
;	pix=pix: pixel size in arcsec. default: =1.830
; COMMON BLOCKS :   none
; SIDE EFFECTS :
; RESTRICTIONS :
; PROCEDURE :
;
; MODIFICATION HISTORY :
;   Hui LI, write on 2001/02/25
;	Hui LI, 2001/06/04 add /xrange,/yrange,/xtitle,/ytitle,/pix
;------------------------------------------------------------------------
;-
pro hr_mag,fname,back=back,saveps=saveps,dir=dir,ext=ext,btonly=btonly,$
	blonly=blonly,nstep=nstep,magnify=magnify,btcont=btcont,nocaption=nocaption,$
	level=level,blackgnd=blackgnd,xrange=xrange,yrange=yrange,xtitle=xtitle,$
	ytitle=ytitle,pix=pix

on_error,2
!except=2
if !version.os_family eq 'Windows' then olddev='win' $
else if !version.os_family eq 'unix' then olddev='x' else olddev=!d.name
loadct,3
scale0=1 & flag=0
yp0=60 & x0=100 & y0=140 & imgmin=0.2
bfact=10. & plmin=0.002
c_labels=0

; ----------------- handle keywords ----------------------

if not keyword_set(blackgnd) then begin
	tvlct,0,0,0,!d.table_size-2     ; color of contour
	tvlct,0,0,0,!d.table_size-1     ; color of frame
	tvlct,255,255,255,0             ; back ground
endif
if keyword_set(pix) then pix=pix else pix=1.830
if keyword_set(nstep) then nstep=nstep else nstep=10
if keyword_set(ext) then ext=ext else ext='.fit'
if keyword_set(magnify) then magnify=magnify else $
	if !d.x_size gt 400 then magnify=3 else magnify=2.5
if keyword_set(dir) then begin
    if !version.os_family eq 'Windows' then begin
        if strmid(dir,strlen(dir)-1,1) eq '\' then dir=dir else dir=dir+'\'
    endif else begin
        if strmid(dir,strlen(dir)-1,1) eq '/' then dir=dir else dir=dir+'/'
    endelse
endif else begin
    dir=''
endelse
if keyword_set(saveps) then begin
	scale0=25 & flag=1
    ps_file=dir+'idl.ps'
    ss=size(saveps)
    if ss(ss(0)+1) eq 7 then begin
        sl=strlen(saveps)-2
        if (sl gt 1) and (strupcase(strmid(saveps,sl,2)) eq 'PS') then $
            ps_file=dir+saveps else ps_file=dir+saveps+'.ps'
    endif
endif

; ------------ Determine and read files ----------------------

fname=dir+fname+ext

ff=findfile(fname)
!error=0
if ff(0) eq '' then begin
    !error=1
    if keyword_set(dir) then $
        print,'File '+fname+' not found ! Procedure stopped.' else $
        print,'File '+fname+' in current directory not found ! Procedure stopped.'
    return
endif

array=readfits(fname,head,/sil)

bx=array(*,*,0)
by=array(*,*,1)
bz=array(*,*,2)
bt=sqrt(bx^2+by^2)
nx=sxpar(head,'NAXIS1')
ny=sxpar(head,'NAXIS2')

if keyword_set(back) then begin
	ss=size(back)
	if ss(ss(0)+1) eq 7 then begin
		sname=dir+back+'.dat'
		ff=findfile(sname)
		!error=0
		if ff(0) eq '' then begin
   			!error=1
   			if keyword_set(dir) then $
	       		print,'File '+sname+' not found ! Procedure stopped.' else $
       			print,'File '+sname+' in current directory not found ! Procedure stopped.'
   			return
		endif
		hr_read,sname,img,hd4
		img=congrid(float(img),nx,ny,/inter,/min)
	endif else begin
		if ss(0) eq 0 then img=array(*,*,3) $
		else img=congrid(float(back),nx,ny,/inter,/min)
	endelse
endif

bzmax=max(bz) & bzmin=min(bz)
btmax=max(bt)
print,'Bt_max = ',btmax,' (Gauss)'
print,'Bz_max = ',bzmax,' (Gauss)   Bz_min = ',bzmin,' (Gauss)'

;-------------- Open window -------------------

y0=yp0
wwx=nx*magnify+2*x0     ;expected window size
wwy=ny*magnify+y0+170
if !d.window ne -1 then wno=!d.window else wno=0
if not keyword_set(printout) and not keyword_set(saveps) 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
	pos=[scale0*x0,scale0*y0,scale0*(x0+nx*magnify-1),scale0*(y0+ny*magnify-1)]
endif else begin
    if keyword_set(xsize) then xsize=float(xsize) else xsize=18.0
   	ysize=xsize*ny/nx
   	pos=[0,0,fix(xsize*1000),fix(ysize*1000)]
   	xoff=fix((21-xsize)/2.+0.5) & yoff=fix((29.7-ysize)/2+0.5)
   	set_plot,'ps'
   	device,file=ps_file,xs=xsize,ys=ysize,xoff=xoff,yoff=yoff
endelse

erase

; ------------------- Display Spot --------------------------

if keyword_set(spotname) or keyword_set(back) then begin
   	img=10.+230.*(img-min(img))/(max(img)-min(img))
   	img=sqrt(img^3)
   	img=10.+230.*(img-min(img))/(max(img)-min(img))
   	if keyword_set(saveps) then begin
       	img1=congrid(img,nx,ny,/interp)
   	   	tv,img1,pos(0),pos(1),/dev,xsize=pos(2)-pos(0),ysize=pos(3)-pos(1)
   	endif else begin
    	img1=congrid(img,pos(2)-pos(0),pos(3)-pos(1),/interp)
    	tv,img1,pos(0),pos(1),/dev
   	endelse
endif

;------  display contour of longitudinal or transversal fields  -------

if not keyword_set(btonly) then begin
   	if keyword_set(btcont) then begin
       	r=smooth(congrid(bt,nx,ny,/interp),2)
       	if keyword_set(level) then levp=level else $
       		levp=[200,500,800,1200,1600,2000,2500,3000]
       	levstr='level='
   	endif else begin
       	r=float(smooth(congrid(bz,nx,ny,/interp),2))
       	if keyword_set(level) then levp=level else $
        	levp=[40,80,160,320,640,960,1280,1600,1920,2240,2560,2880]
   	    levstr='level=!9+!3'
   	endelse
   	nn=n_elements(levp)
    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'
    	if not keyword_set(btcont) then begin
        levn=levp
   	    for i=0,nn-1 do levn(i)=-levp(nn-i-1)
       	lev0=[levn,levp]
        lnum=n_elements(lev0)
   	    lstyle=intarr(lnum)
       	lthick=intarr(lnum)
        lcolor=replicate(!d.n_colors/1.8,lnum)
   	    lstyle(where(lev0 lt 0))=1   ; dot
       	lthick(where(lev0 lt 0))=2
        lcolor(where(lev0 lt 0))=!d.n_colors/1.6
   	endif else begin
       	lev0=levp
        lstyle=intarr(n_elements(lev0))
   	    lthick=intarr(n_elements(lev0))
    endelse
   	lstyle(where(lev0 gt 0))=0   ; solid line
    lthick(where(lev0 gt 0))=1
    	if not keyword_set(nolong) then begin
       	if flag eq 1 then $
           	contour,r,pos=pos,/dev,xstyle=1+4,ystyle=1+4,xrange=[0.,nx],$
            yrange=[0.,ny],level=lev0,c_linestyle=lstyle,c_thick=lthick,$
   	        /noerase,spline=spline,c_labels=c_labels,back=255,color=1, $
       	    xtitle=xtitle,ytitle=ytitle $
        else $
   	        contour,r,pos=pos,/dev,xstyle=1+4,ystyle=1+4,xrange=[0.,nx],$
       	    yrange=[0.,ny],level=lev0,c_linestyle=lstyle,c_thick=lthick,$
           	xtitle=xtitle,ytitle=ytitle, $
            /noerase,spline=spline,c_labels=c_labels,back=255,c_color=lcolor
   	                    ; use follow method and surpress labels
    endif
endif

;----------------  display transversal fields  --------------

tmp=strcompress(sxpar(head,'METHOD'),/remove)
if tmp eq 'OLD' then btmin=10 else btmin=30
if keyword_set(saveps) then clr=10 else $
	if keyword_set(blackgnd) then clr=!d.n_colors-10 else $
		if !d.n_colors lt 255 then clr=10 else clr=!d.n_colors-10
if not keyword_set(blonly) then begin
    if flag eq 1 then plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase,color=clr, $
       	xrange=[0.,nx*magnify],yrange=[0.,ny*magnify],xstyle=1+4,ystyle=1+4 else $
        plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase,color=150, $
   	    xrange=[0.,nx*magnify],yrange=[0.,ny*magnify],xstyle=1+4,ystyle=1+4
    nstep0=max([nstep/magnify,1])        ;initially set nstep=10
   	for j=nstep0,ny-nstep0,nstep0 do begin
       	for i=nstep0,nx-nstep0,nstep0 do begin
            Bt1=(bt(i,j)>1)/bfact
   	        plen=sqrt((Bt1-Btmin)>1)           ;Bfact=10,Btmin=100
       	    if Bt1 gt Btmin then begin
                th=atan(by(i,j),bx(i,j))
               	x=plen*cos(th)/magnify
                y=plen*sin(th)/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=clr
       	        dth=20./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=clr
            endif
       	endfor
    endfor
endif

;--------------- Plot frame -------------------

title='SOLAR VECTOR MAGNETOGRAM, BEIJING'
if keyword_set(xtitle) then xtitle=xtitle else xtitle='Distance (arcsec)'
if keyword_set(ytitle) then ytitle=ytitle else ytitle='Distance (arcsec)'

xmax=nx*pix & ymax=ny*pix    ;Field of view: 5.23 by 3.63 arc min
if keyword_set(xrange) then xrange=xrange else xrange=[0.,xmax]
if keyword_set(yrange) then yrange=yrange else yrange=[0.,ymax]

delx=xrange(1)-xrange(0) & dely=yrange(1)-yrange(0)
lgx=fix(alog10(delx))
lgy=fix(alog10(dely))
stepx=10.^lgx & stepy=10.^lgy
ntx=delx/stepx & nty=dely/stepy
if ntx lt 5 then if ntx lt 3 then stepx=stepx/4 else stepx=stepx/2
if nty lt 5 then if nty lt 3 then stepy=stepy/4 else stepy=stepy/2
ntx=delx/stepx & nty=dely/stepy
tickvx=fltarr(10) & tickvy=tickvx
for i=0,9 do tickvx(i)=stepx*i+xrange(0)
for i=0,9 do tickvy(i)=stepy*i+yrange(0)

datime='Date & Time: '+strcompress(sxpar(head,'DATE-OBS'),/remove)+' '+$
	strcompress(sxpar(head,'TIME-OBS'),/remove)
if keyword_set(blonly) then begin
	files='Files Name: '+sxpar(head,'LNAME')
	cont='Contents : Longitudinal Field'
endif else begin
	if keyword_set(btonly) then begin
		files='File Names: '+sxpar(head,'QNAME')+' '+sxpar(head,'UNAME')
		cont='Contents : Transversal Field'
	endif else begin
		files='File Names: '+sxpar(head,'LNAME')+' '+$
			sxpar(head,'QNAME')+' '+sxpar(head,'UNAME')
		cont='Contents : Vector Magnetogram'
	endelse
endelse

tmp=sxpar(head,'DISK-COY')
if float(tmp gt 0) then sn='N'+string(tmp,format='(f4.1)') $
	else sn='S'+string(abs(tmp),format='(f4.1)')
tmp=sxpar(head,'DISK-COX')
if float(tmp gt 0) then ew='W'+string(tmp,format='(f4.1)') $
	else ew='E'+string(abs(tmp),format='(f4.1)')
sunc='Sun Disk Coord. : '+sn+'   '+ew
orth='Orthogonal Coor. : '+string(sxpar(head,'ORTH-COX'),format='(f6.1)')+$
	'   '+string(sxpar(head,'ORTH-COY'),format='(f6.1)')
carr='Carrington Coor. : '+string(sxpar(head,'CARR-COX'),format='(f6.1)')+$
	'   '+string(sxpar(head,'CARR-COY'),format='(f6.1)')
wave='Wave Length : '+strcompress(sxpar(head,'WAVELENT'),/remove)
see='Seeing : '+strcompress(sxpar(head,'SEEING'),/remove)
fram='Frame Number : '+string(sxpar(head,'FRAMES'),format='(i4)')

if keyword_set(saveps) then begin
   	plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase, $
       	xrange=xrange,yrange=yrange,xstyle=1,ystyle=1,  $
        xticks=fix(ntx),xtickv=tickvx,xminor=5,xticklen=0.03, $
   	    yticks=fix(nty),ytickv=tickvy,yminor=5,yticklen=0.03, $
       	xtitle=xtitle,ytitle=ytitle,charsize=charsize,color=clr
	if not keyword_set(nocaption) then begin
	    xyouts,pos(0)*0.5,pos(3)+4500,/dev,title,chars=2.5,color=clr
    	xyouts,pos(0)+1000,pos(3)+3500,/dev,datime,chars=1.,color=clr
	    xyouts,pos(0)+1000,pos(3)+3000,/dev,files,chars=1,color=clr
    	xyouts,pos(0)+1000,pos(3)+2500,/dev,cont,chars=1,color=clr
    	xyouts,pos(0)+1000,pos(3)+2000,/dev,sunc,chars=1,color=clr
	    xyouts,pos(0)+1000,pos(3)+1500,/dev,orth,chars=1,color=clr
    	xyouts,pos(0)+1000,pos(3)+1000,/dev,carr,chars=1,color=clr
	    xyouts,pos(2)*0.6,pos(3)+3500,/dev,wave,chars=1,color=clr
    	xyouts,pos(2)*0.6,pos(3)+3000,/dev,see,chars=1,color=clr
	    xyouts,pos(2)*0.6,pos(3)+2500,/dev,fram,chars=1,color=clr
	    if not keyword_set(btonly) then $
	    	xyouts,pos(0)+1000,pos(3)+500,/dev,levstr,chars=1,color=clr
	    xyouts,1000,-yoff*850,/dev,'Created by IDL at '+string(systime()),chars=1,color=clr
	endif
   	device,/close
    print,'Result has been saved in PostScript file: '+ps_file
   	set_plot,olddev
endif else begin
   	plot,[0,0],[0,0],/nodata,pos=pos,/dev,/noerase, $
       	xrange=xrange,yrange=yrange,xstyle=1,ystyle=1,  $
        xticks=fix(ntx),xtickv=tickvx,xminor=5,xticklen=0.03, $
   	    yticks=fix(nty),ytickv=tickvy,yminor=5,yticklen=0.03, $
       	xtitle=xtitle,ytitle=ytitle,color=clr,charsize=charsize
	if not keyword_set(nocaption) then begin
	    xyouts,pos(0)*0.5,pos(3)+130,/dev,title,chars=2.,color=clr
    	xyouts,pos(0),pos(3)+105,/dev,datime,chars=0.8,color=clr
	    xyouts,pos(0),pos(3)+90,/dev,files,chars=0.8,color=clr
    	xyouts,pos(0),pos(3)+75,/dev,cont,chars=0.8,color=clr
	    xyouts,pos(0),pos(3)+60,/dev,sunc,chars=0.8,color=clr
    	xyouts,pos(0),pos(3)+45,/dev,orth,chars=0.8,color=clr
	    xyouts,pos(0),pos(3)+30,/dev,carr,chars=0.8,color=clr
    	xyouts,pos(2)*0.7,pos(3)+105,/dev,wave,chars=0.8,color=clr
	    xyouts,pos(2)*0.7,pos(3)+90,/dev,see,chars=0.8,color=clr
    	xyouts,pos(2)*0.7,pos(3)+75,/dev,fram,chars=0.8,color=clr
	    if not keyword_set(btonly) then $
		    xyouts,pos(0),pos(3)+15,/dev,levstr,chars=0.8,color=clr
	endif
endelse

; -----------

end