pro plot_sw_mag,file_sw_mfi,timerange=timerange,time_shock_start=time_shock_start, time_shock_end=time_shock_end, time_icme_start=time_icme_start,time_icme_end=time_icme_end,time_cir=time_cir,ps_file=ps_file,par_file=par_file,box=box,hc=hc,comment=comment,win=win

;+
; project	: CME, ICME and Goemagnetic Storm
;
; Name		:plot_sw_mag
;
; Purpose	:plot solar wind magnetic field and calculate parameters
;
; Explanation	:
;
; CALLING	: IDL> plot_sw_mag,file_sw_mfi
;
; INPUTS:
;	file_sw_mfi:	solar wind magnetic field data file in CDF form;
;
; OPTIONAL INPUTS: 
;	timerange: data time plotting range in the format of 
;			['yyyy/mm/dd hh:mm','yyyy/mm/dd hh:mm']
;	ps_file	: file name for the output ps file, default "tmp.ps"
;	time_shock_start: start time of shock, 'yyyy/mm/dd hh:mm'; 
;				vertical solid red line
;	time_shock_end	: end time of shock sheath, 'yyyy/mm/dd hh:mm'; 
;				default value is start time of icme
;				vertical dotted red line
;	time_icme_start	: start time of ICME, 'yyyy/mm/dd hh:mm'
;				vertical solid blue line
;	time_icme_end	: end time of ICME, 'yyyy/mm/dd hh:mm'
;				vertical dotted blue line
;	time_cir	: interface of corotating interaction region
;				vertical dotted yellow line
;	box		: the size of smoothing box 
;				when use smooth or median filter on data
;				default smooting size is 13 points
;	comment		: comment
;
; OUTPUT:
;	
; OPTIONAL OUTPUT:
;	par_file: an acsii file holding the calculated parameters
;
; KEYWORDS:
;	hc	: make hard copy in PS format
;
;
; COMMENTS: 
;
;		:this program makes use of CDAWLIB for accessing CDF files
; 
; CATAGORY	: Plot
;
; WRITTEN	: Watanachak Poomvises, Jie Zhang, GMU, 2006 March 18
;
; Modifed	:
;
;-

;
;----------- check the calling-----
;
if n_params() lt 1 then begin
	print,'USAGE:  plot_sw_mag,file_sw_mfi'
	return
endif

;
;---------- setup the display ------
;

;load and define colors
ncolors=1.0*!d.table_size 
loadct,15 
color_black=0 
color_white=255 
color_plot=color_black  & color_text=color_black
color_red=15 & color_blue=127 & color_yellow=192

;check display device
if not keyword_set(hc) then begin
	set_plot,'X'
	window,1,xsize=600,ysize=800,retain=2 ; using window as display
	!p.charsize=2.0
	!p.symsize=2.9
	thick_vline=2.0 ; thickness of the vertical lines
	charsize=1.2    ;the printed parameters
endif	else begin
	set_plot,'PS' ;using ps as display device
	if not keyword_set(ps_file) then ps_file='./tmp.eps'
	device,filename=ps_file
  	device,/color,bits_per_pixel=8,/portr,/inch,xsize=6.0,ysize=8.0,$
         xoff=1.,yoff=1.0,/encap
   	!x.thick=2.0 & !y.thick=2.0 & !p.thick=2.0 & !p.charthick=2.0
  	!p.charsize=1.5
  	!p.symsize=1.0	
        thich_vline=2.0 ; thickness of the vertical lines
	charsize=.75
endelse

;plot panels and colors                                                                                
!p.color=color_black
!p.background=color_white
!p.multi=[0,1,8,0,0]
!x.margin=[8,8]
!y.margin=[0,0]

;
;-------- read-in the data --------
;

if not keyword_set(box) then box=13; defaul smoothing size is 13 points

cdf_data_mfi=read_mycdf('',file_sw_mfi,all=2)
num_mfi=n_elements(cdf_data_mfi.epoch.dat)

;initialize the data variables
sec_mfi    = fltarr(num_mfi) ;time in sec from first point
bx_gse = fltarr(num_mfi)     	; Bx in GSE
by_gse = fltarr(num_mfi)	; By in GSE
bz_gse = fltarr(num_mfi)	; Bz in GSE
b_total=fltarr(num_mfi)		; Magnetic field intensity

;read mfi data
time_string=strarr(num_mfi)
for i=0L,num_mfi-1 do time_string[i]=cdf_encode_epoch(cdf_data_mfi.epoch.dat[i],epoch=0)
time_utc_mfi=str2utc(time_string)
time_tai_mfi=utc2tai(time_utc_mfi)

if not keyword_set(timerange) then begin  ;use MFI data time range as plot time
	time_plot_start=time_string[0]
	time_plot_end=time_string[num_mfi-1]
	timerange=[time_plot_start,time_plot_end]
endif

sec_mfi=(time_tai_mfi-time_tai_mfi[0])
bx_gse=reform(cdf_data_mfi.bgsec.dat(0,*))
by_gse=reform(cdf_data_mfi.bgsec.dat(1,*))
bz_gse=reform(cdf_data_mfi.bgsec.dat(2,*))
bx_gse=median(bx_gse,box)
by_gse=median(by_gse,box)
bz_gse=median(bz_gse,box)

arr=where(bx_gse le -100,count)
if (count ge 0) then begin
    	for i=0,count-1 do bx_gse(arr[i])=bx_gse(arr[i]-1)
endif
 
arr=where(by_gse le -100,count)
if (count ge 0) then begin
 for i=0,count-1 do by_gse(arr[i])=by_gse(arr[i]-1)
endif
 
arr=where(bz_gse le -100,count)
if (count ge 0) then begin
 for i=0,count-1 do bz_gse(arr[i])=bz_gse(arr[i]-1)
endif

b_total =sqrt((bx_gse)^2.0+(by_gse)^2.0+(bz_gse)^2.0)

;
;--------calculating parameters used in plotting
;

;calculate theta angle of magnetic field in GSE (elevation angle)
tmp=bz_gse/sqrt(bx_gse^2.0+by_gse^2.0)
theta=atan(tmp)*180.0/!pi
;theta=90-theta_tmp

;calculate phi angle of magnetic field (azimuthal angle)
tmp=atan(median(by_gse,box),median(bx_gse,box))*180.0/!pi
phi=tmp
arr=where(phi lt 0, count)
if count gt 0 then phi(arr)=phi(arr)+360
;remove jumping points between close to 360 to close to 0
for i=1L,n_elements(phi)-1 do begin
        if (abs(phi[i]-phi[i-1]) gt 300)  then phi[i]=phi[i-1]
endfor

;
;plotting
;

;indicate the plotting time range at the top
;xyouts,0.1,0.98,timerange[0]+' UT--'+timerange[1]+' UT',/norm,charsize=charsize

;plot B_total
y1=floor(min(b_total)/10.0)*10  ;lower boundary is determined by B
y2=ceil(max(b_total)/10.0)*10  ;upper bounardy is defined by B
yticks=(y2-y1)/10
utplot,time_utc_mfi,b_total,timerange=timerange, yrange=[y1,y2], ytitle= 'B(nT)',ystyle=1, xstyle=1, yticks=yticks,xtitle=' ', xticklen=0.0,xchar= 0.01,ymargin=[0,1]
if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[y1,y2],thick=thick_vline,color=color_red
if keyword_set(time_icme_start) then outplot,[time_icme_start,time_icme_start],[y1,y2],thick=thick_vline,color=color_blue
if keyword_set(time_icme_end) then outplot,[time_icme_end,time_icme_end],[y1,y2],thick=thick_vline,color=color_blue,linestyle=2

;plot Bx
y1=floor(min(bx_gse)) 
y2=ceil(max(bx_gse))
yticks=(y2-y1)/10
utplot,time_utc_mfi,bx_gse,timerange=timerange, yrange=[y1,y2], ytitle= 'Bx(nT)',ystyle=1, xstyle=1, yticks=yticks,xtitle=' ', xticklen=0.0,xchar= 0.01
if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[y1,y2],thick=thick_vline,color=color_red
if keyword_set(time_icme_start) then outplot,[time_icme_start,time_icme_start],[y1,y2],thick=thick_vline,color=color_blue
if keyword_set(time_icme_end) then outplot,[time_icme_end,time_icme_end],[y1,y2],thick=thick_vline,color=color_blue,linestyle=2

;plot By 
y1=floor(min(by_gse))  
y2=ceil(max(by_gse))  
yticks=(y2-y1)/10 
utplot,time_utc_mfi,by_gse,timerange=timerange, yrange=[y1,y2], ytitle= 'By(nT)',ystyle=1, xstyle=1, yticks=yticks,xtitle=' ', xticklen=0.0,xchar= 0.01
if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[y1,y2],thick=thick_vline,color=color_red
if keyword_set(time_icme_start) then outplot,[time_icme_start,time_icme_start],[y1,y2],thick=thick_vline,color=color_blue
if keyword_set(time_icme_end) then outplot,[time_icme_end,time_icme_end],[y1,y2],thick=thick_vline,color=color_blue,linestyle=2

;plot Bz 
y1=floor(min(bz_gse))
y2=ceil(max(bz_gse))  
yticks=(y2-y1)/10 
utplot,time_utc_mfi,bz_gse,timerange=timerange, yrange=[y1,y2], ytitle= 'Bz(nT)',ystyle=1, xstyle=1, yticks=yticks,xtitle=' ', xticklen=0.0,xchar= 0.01
if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[y1,y2],thick=thick_vline,color=color_red
if keyword_set(time_icme_start) then outplot,[time_icme_start,time_icme_start],[y1,y2],thick=thick_vline,color=color_blue
if keyword_set(time_icme_end) then outplot,[time_icme_end,time_icme_end],[y1,y2],thick=thick_vline,color=color_blue,linestyle=2

 
;plot theta 
y1=-90   
y2=90 
yticks=(y2-y1)/30
utplot,time_utc_mfi,theta,timerange=timerange, yrange=[y1,y2], ytitle= '!7h!3',ystyle=1, xstyle=1, yticks=yticks,xtitle=' ', xticklen=0.0,xchar= 0.01
if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[y1,y2],thick=thick_vline,color=color_red
if keyword_set(time_icme_start) then outplot,[time_icme_start,time_icme_start],[y1,y2],thick=thick_vline,color=color_blue
if keyword_set(time_icme_end) then outplot,[time_icme_end,time_icme_end],[y1,y2],thick=thick_vline,color=color_blue,linestyle=2
 
;plot phi 
xtitle=timerange[0]+' UT--'+timerange[1]+' UT'
y1=0                     
y2=360                       
yticks=(y2-y1)/60 
utplot,time_utc_mfi,phi,timerange=timerange, yrange=[y1,y2],xtitle=xtitle,ytitle= '!7u!3',ystyle=1, xstyle=1, yticks=yticks, xticklen=0.0
if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[y1,y2],thick=thick_vline,color=color_red
if keyword_set(time_icme_start) then outplot,[time_icme_start,time_icme_start],[y1,y2],thick=thick_vline,color=color_blue
if keyword_set(time_icme_end) then outplot,[time_icme_end,time_icme_end],[y1,y2],thick=thick_vline,color=color_blue,linestyle=2




;
;--------Derive relevant solar wind parameters-------
;for each variable, calculate minimum, maximum, average, and chi_square

;open output data file for derived parameters, if set 
if keyword_set(par_file) then OPENW, lu, par_file,/get_lun

x00=0.05 & y00=0.20   ;starting position for the text
dx_s=0.08            ;seperating between sub-columns
dx=0.4               ;seperating between major columns
dy=0.02             ;seperating between rows

xyouts,x00,y00,'------------------------------------------------------',/norm,charsize=charsize

;print the two columns indicating SH and ICME
x0=x00 & y0=y00-1.0*dy
xyouts,x0+0.20,y0,'SH ',/norm,charsize=charsize
x0=x00+dx & y0=y00-1.0*dy
xyouts,x0+0.20,y0,'ICME ',/norm,charsize=charsize

;print the sub-columns
x0=x00 & y0=y00-2.0*dy 
xyouts,x0,y0,'          ave    sdev    max    min',/norm,charsize=charsize
x0=x00+dx & y0=y00-2.0*dy  
xyouts,x0,y0,'          ave    sdev    max    min',/norm,charsize=charsize
 


;find the shock array
;find the shock end time, if possible
;usually the same as ICME start time
if keyword_set(time_shock_start) then begin
	if not keyword_set(time_shock_end) then begin
        	if keyword_set(time_icme_start) then time_shock_end=time_icme_start
	endif
endif


if keyword_set(time_shock_start) and keyword_set(time_shock_end) then begin
        time_tai_shock_start=utc2tai(str2utc(time_shock_start))
        time_tai_shock_end=utc2tai(str2utc(time_shock_end))
        arr_sh=where(time_tai_mfi ge time_tai_shock_start and time_tai_mfi le time_tai_shock_end,count)
	b_sh=b_total(arr_sh)
	bx_sh=bx_gse(arr_sh)
	by_sh=by_gse(arr_sh)
	bz_sh=bz_gse(arr_sh)
endif 

if keyword_set(time_icme_start) and keyword_set(time_icme_end) then begin
        time_tai_icme_start=utc2tai(str2utc(time_icme_start))
        time_tai_icme_end=utc2tai(str2utc(time_icme_end))
        arr_icme=where(time_tai_mfi ge time_tai_icme_start and time_tai_mfi le time_tai_icme_end,count) 
        b_icme=b_total(arr_icme)
        bx_icme=bx_gse(arr_icme)
        by_icme=by_gse(arr_icme)
        bz_icme=bz_gse(arr_icme)
endif  

;print B in SH and ICME
x0=x00 & y0=y00-3.0*dy 
xyouts,x0,y0,'B(nT): ',/norm,charsize=charsize 
if keyword_set(b_sh) then begin
;calculate B in shock
	tmp=b_sh
	mom=moment(tmp,sdev=sdev)
	ave_str=strtrim(string(mom[0],format='(f5.1)'))
	sdev_str=strtrim(string(sdev,format='(f5.1)'))
	max_str=strtrim(string(max(tmp),format='(f5.1)'))
	min_str=strtrim(string(min(tmp),format='(f5.1)'))

	xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
	xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
	xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
	xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
	print,'B_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'B_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif

x0=x00+dx & y0=y00-3.0*dy
if keyword_set(b_icme) then begin
;calculate B in icme
	tmp=b_icme
        mom=moment(tmp,sdev=sdev)
        ave_str=strtrim(string(mom[0],format='(f5.1)'))
        sdev_str=strtrim(string(sdev,format='(f5.1)'))
        max_str=strtrim(string(max(tmp),format='(f5.1)'))
        min_str=strtrim(string(min(tmp),format='(f5.1)'))

        xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
        xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
        xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
        xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
        print,'B_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'B_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif

;print Bx in SH and ICME 
x0=x00 & y0=y00-4.0*dy  
xyouts,x0,y0,'Bx(nT): ',/norm,charsize=charsize
if keyword_set(b_sh) then begin
;calculate B in shock 
        tmp=bx_sh 
        mom=moment(tmp,sdev=sdev)
        ave_str=strtrim(string(mom[0],format='(f5.1)'))
        sdev_str=strtrim(string(sdev,format='(f5.1)'))
        max_str=strtrim(string(max(tmp),format='(f5.1)'))
        min_str=strtrim(string(min(tmp),format='(f5.1)'))
 
        xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
        xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
        xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
        xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
        print,'Bx_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'Bx_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif 
 
x0=x00+dx & y0=y00-4.0*dy 
if keyword_set(bx_icme) then begin
;calculate Bx in icme 
        tmp=bx_icme 
        mom=moment(tmp,sdev=sdev)
        ave_str=strtrim(string(mom[0],format='(f5.1)'))
        sdev_str=strtrim(string(sdev,format='(f5.1)'))
        max_str=strtrim(string(max(tmp),format='(f5.1)'))
        min_str=strtrim(string(min(tmp),format='(f5.1)'))
 
        xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
        xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
        xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
        xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
        print,'Bx_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'Bx_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif 
  

 
;print By in SH and ICME  
x0=x00 & y0=y00-5.0*dy   
xyouts,x0,y0,'By(nT): ',/norm,charsize=charsize
if keyword_set(by_sh) then begin
;calculate B in shock  
        tmp=by_sh  
        mom=moment(tmp,sdev=sdev)
        ave_str=strtrim(string(mom[0],format='(f5.1)'))
        sdev_str=strtrim(string(sdev,format='(f5.1)'))
        max_str=strtrim(string(max(tmp),format='(f5.1)'))
        min_str=strtrim(string(min(tmp),format='(f5.1)'))
  
        xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
        xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
        xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
        xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
        print,'By_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'By_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif  
  
x0=x00+dx & y0=y00-5.0*dy  
if keyword_set(by_icme) then begin
;calculate Bx in icme  
        tmp=by_icme  
        mom=moment(tmp,sdev=sdev)
        ave_str=strtrim(string(mom[0],format='(f5.1)'))
        sdev_str=strtrim(string(sdev,format='(f5.1)'))
        max_str=strtrim(string(max(tmp),format='(f5.1)'))
        min_str=strtrim(string(min(tmp),format='(f5.1)'))
  
        xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
        xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
        xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
        xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
        print,'By_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'By_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif  



;print Bz in SH and ICME  
x0=x00 & y0=y00-6.0*dy
xyouts,x0,y0,'Bz(nT): ',/norm,charsize=charsize
if keyword_set(bz_sh) then begin
;calculate Bz in shock  
        tmp=bz_sh 
        mom=moment(tmp,sdev=sdev) 
        ave_str=strtrim(string(mom[0],format='(f5.1)'))
        sdev_str=strtrim(string(sdev,format='(f5.1)'))
        max_str=strtrim(string(max(tmp),format='(f5.1)'))
        min_str=strtrim(string(min(tmp),format='(f5.1)'))
   
        xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
        xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
        xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
        xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
        print,'Bz_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'Bz_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif   
   
x0=x00+dx & y0=y00-6.0*dy   
if keyword_set(bz_icme) then begin
;calculate Bx in icme  
        tmp=bz_icme
        mom=moment(tmp,sdev=sdev)
        ave_str=strtrim(string(mom[0],format='(f5.1)'))
        sdev_str=strtrim(string(sdev,format='(f5.1)'))
        max_str=strtrim(string(max(tmp),format='(f5.1)'))
        min_str=strtrim(string(min(tmp),format='(f5.1)'))
   
        xyouts,x0+1.0*dx_s,y0,ave_str,/norm,charsize=charsize
        xyouts,x0+2.0*dx_s,y0,sdev_str,/norm,charsize=charsize
        xyouts,x0+3.0*dx_s,y0,max_str,/norm,charsize=charsize
        xyouts,x0+4.0*dx_s,y0,min_str,/norm,charsize=charsize
        print,'Bz_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
        if keyword_set(par_file) then printf,lu,'Bz_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' '
endif
 
fdate=strmid(timerange(0),0,4)+strmid(timerange(0),5,2)+strmid(timerange(0),8,2)+strmid(timerange(0),11,2)
;spawn, 'ls /home/phess4/solarwind/bz/events0318/'+fdate+'/ -d', result

;if result[0] eq '' then	spawn, 'mkdir /ome/phess4/solarwind/bz/events0613/'+fdate

 
write_png,'/swd4/Users/projects/sun_to_earth/2012/idl/2012/plots/'+fdate+'_mag.png',tvrd(/true)



;close the parameter file
if keyword_set(lu) then close,lu

empty
SET_PLOT, 'X'

return
end