pro plot_sw_mag_plasma,file_sw_mfi,file_sw_plasma,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,sw_type=sw_type,sw_dr_type=sw_dr_type,dst_peak=dst_peak,time_dst_peak=time_dst_peak,comment=comment,win=win ;+ ; project : CME, ICME and Goemagnetic Storm ; ; Name :plot_sw_mag_plasma ; ; Purpose :plot Dst index, solar wind magnetic field and plasma data, ; and parameterize solar wind properties ; ; Explanation : ; ; CALLING : IDL> plot_sw_mag_plasma,file_sw_mfi,file_sw_plasma ; ; INPUTS: ; file_sw_mfi: solar wind magnetic field data file in CDF format ; file_sw_plasma: solar wind plasma data file in CDF format ; ; 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 ; sw_type : type of solar wind stream, e.g., SH+ICME ; sw_dr_type : type of storm driver, e.g., MC ; dst_peak : peak or dip Dst value ; time_dst_peak : time of Dst peak ; comment : comment ; ; OUTPUT: ; ; OPTIONAL OUTPUT: ; par_file: an acsii file holding the calculated parameters ; ; KEYWORDS: ; hc : make hard copy in PS format ; ; ; COMMENTS: : In this version, input data are ACE solar wind data files ; in CDF format. In the future, input will be defaul, and only need ; to specify the time range. It should also accept WIND data ; ; :this program makes use of CDAWLIB for accessing CDF files ; ; CATAGORY : Plot ; ; WRITTEN : Watanachak Poomvises, Jie Zhang, GMU, 2006 March 18 ; ; ; ;- ; ;----------- check the calling----- ; if n_params() lt 2 then begin print,'USAGE: plot_sw_mag_plasma,file_sw_mfi,file_sw_plasma' return endif ; ;---------- setup the display ------ ; ;============================================================== editing===================================================== year=anytim2utc(timerange[0],/ecs) year=strmid(year,0,4) ;=========================================================================================================================== ;=========================================================================================================================== ;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,0,xsize=600,ysize=1000,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.0 ;the printed parameters endif else begin set_plot,'PS' ;using ps as display device if not keyword_set(ps_file) then ps_file='./test_event.eps' device,filename=ps_file device,/color,bits_per_pixel=8,/portr,/inch,xsize=6.0,ysize=9.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 thick_vline=4.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,7,0,0] ;!p.multi=[0,1,11,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 point cdf_data_plasma=read_mycdf('',file_sw_plasma+'.cdf',all=2) cdf_data_mfi=read_mycdf('',file_sw_mfi+'.cdf',all=2) num_mfi=n_elements(cdf_data_mfi.epoch.dat) num_plasma=n_elements(cdf_data_plasma.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 sec_plasma = fltarr(num_plasma) ; time in sec from first pint vx_gse = fltarr(num_plasma) ; Vx in GSE vy_gse = fltarr(num_plasma) ; Vy in GSE vz_gse = fltarr(num_plasma) ; Vz in GSE np = fltarr(num_plasma) ; proton density v_total= fltarr(num_plasma) ; velocity of the solar wind tpr = fltarr(num_plasma) ; temperature ramp_press = fltarr(num_plasma) ; ram pressure calculted ;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) ;read plasma data time_string=strarr(num_plasma) for i=0L,num_plasma-1 do time_string[i]=cdf_encode_epoch(cdf_data_plasma.epoch.dat[i],epoch=0) time_utc_plasma=str2utc(time_string) time_tai_n=utc2tai(time_utc_plasma) time_tai_t=time_tai_n time_tai_plasma=time_tai_n sec_plasma=(time_tai_n-time_tai_n[0]) vx_gse=reform(cdf_data_plasma.v_gse.dat(0,*)) vy_gse=reform(cdf_data_plasma.v_gse.dat(1,*)) vz_gse=reform(cdf_data_plasma.v_gse.dat(2,*)) arr=where(vx_gse le -1000,count) if (count ge 0) then begin for i=0,count-1 do begin if arr[i] ne 0 then begin vx_gse(arr[i])=vx_gse(arr[i]-1) endif endfor endif arr=where(vy_gse le -1000,count) if (count ge 0) then begin for i=0,count-1 do begin if arr[i] ne 0 then begin vy_gse(arr[i])=vy_gse(arr[i]-1) endif endfor endif arr=where(vz_gse le -1000,count) if (count ge 0) then begin for i=0,count-1 do begin if arr[i] ne 0 then begin vz_gse(arr[i])=vz_gse(arr[i]-1) endif endfor endif v_total = sqrt(((vx_gse)^2+(vy_gse)^2+(vz_gse)^2)) arr=where(v_total gt 3000,count) if (count ge 0) then begin for i=0,count-1 do begin if arr[i] ne 0 then begin v_total(arr[i])=v_total(arr[i]-1) endif if arr[i] eq 0 then v_total[0]=v_total[1] endfor endif np=reform(cdf_data_plasma.np.dat) ;remove bad points in proton density, using previous value to replace it arr=where(np le -1.0,count) snp=size(np) sarr=size(arr) if float(sarr(1))/float(snp(1)) ge .25 then begin aceswep=get_acedata(timerange(0),timerange(1), /daily,/swepam, /bad_data) np=aceswep.P_DENSITY arr=where(np le -1.0,count) time_tai_n=utc2tai(aceswep) endif if (count ge 0) then begin for i=0,count-1 do begin if arr[i] ne 0 then begin np(arr[i])=np(arr[i]-1) endif endfor endif tpr=reform(cdf_data_plasma.tpr.dat); arr=where(tpr le -1.0,count) stp=size(tpr) sarr=size(arr) if float(sarr(1))/float(stp(1)) ge .25 then begin aceswep=get_acedata(timerange(0),timerange(1), /daily,/swepam, /bad_data) tpr=aceswep.ION_TEMP arr=where(tpr le -1.0,count) time_tai_t=utc2tai(aceswep) endif if (count ge 0) then begin for i=0,count-1 do begin if arr[i] ne 0 then begin tpr(arr[i])=tpr(arr[i]-1) endif if arr[i] eq 0 then tpr[0]=tpr[1] endfor endif dsttest=1 if dsttest ne 1 then begin ;read Dst data dst_path='/home/phess4/solarwind/' ;the directory holding the yearly Dst data year=strmid(utc2str(str2utc(timerange[0])),0,4) dst_filename=dst_path+'dst'+year+'.dat' ;print, dst_filename read_dst,dst_filename,time=time_dst,value=dst if n_elements(dst) ne 0 then begin dst_peak=dst(0) time_dst_peak_ut=time_dst(0) A=size(dst) timeint=utc2int(timerange) ;stop for i=1,A(1)-1 do begin if (dst(i) lt dst_peak) and (time_dst(i).MJD gt timeint(0).MJD) and (time_dst(i).MJD lt timeint(1).MJD) then begin ;if (time_dst(i) gt timeint(0).MJD) then begin dst_peak=dst(i) time_dst_peak_ut=time_dst(i) endif endfor time_dst_peak=utc2str(int2utc(time_dst_peak_ut), /ECS, /TRUNCATE) endif ;stop endif ; ;-------- constants used in calculating and data smoothing ; if not keyword_set(box) then box=13; defaul smoothing size is 13 points mp = 1.6726 ; mass of proton, unit in 10^-24 grams K=1.38; Boltzmann constant, unit in 10^-16 (erg / k /s) Re= 6.37*10^3.0; radius of the Earth, unit in km lz=7.0*Re ;l value used in calculating incoming energy ; ;--------calculating derived 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_tmp=atan(tmp)*180.0/!pi ; change unit from radias to degree theta=90-theta_tmp ;calculate phi angle of magnetic field (azimuthal angle) tmp=atan(by_gse,bx_gse)*180.0/!pi phi=tmp arr=where(by_gse lt 0, count) if count gt 0 then phi(arr)=phi(arr)+360 v_total_unit = v_total*1000*100 ; change velocity unit from km to centimeter ;stop ; calculate plasma beta2 (B) ; interpolate to have the same time sequence in the data array np_mfi = interpol(np,time_tai_n,time_tai_mfi) tpr_mfi = interpol(tpr,time_tai_t,time_tai_mfi) beta2 = np_mfi*K*tpr_mfi/(b_total^2.0/8.0/!pi)/10^6.0 ;stop ; calculate tpr_exp: expected temperature tpr_exp=fltarr(num_plasma) for i=0,num_plasma-1 do begin if v_total(i) ge 500.0 then $ tpr_exp(i)= (0.51*v_total(i) - 142)*10^3.0 if v_total(i) lt 500.0 then $ tpr_exp(i)= (0.031*v_total(i) - 5.1)^2*10^3.0 endfor ;Calculate E: energy input into the magnetosphere from solar wind ;formula from Akasofu 1983, Space Science Review vel_mfi=interpol(v_total,time_tai_plasma,time_tai_mfi) E= vel_mfi*(b_total)^2.0*(sin(theta/2.0*!pi/180.0))^4.0*lz^2.0 ;in Watt ; Calculate ramp pressure and thermal temperature ;ramp_press = np*mp*v_total^2.0; in unit of 10^-24 ;vth =sqrt((K*tpr/(0.5*mp))*10^8.0)/(10^5.0) ;stop ; ;plotting ; if n_elements(dst) ne 0 then begin ;plot Dst y1=floor(min(dst)/50.0)*50 & y2=ceil(max(dst)/50.0)*50 & yticks=(y2-y1)/100 ;print, yticks utplot,time_dst,dst,timerange=timerange,yrange=[y1,y2],ytitle='Dst(nT)', xstyle=1,ystyle=1,xticklen=0.0, xtitle='',xchar= 0.01,yticks=yticks,ymargin=[0,2] endif ;indicate the plotting time range at the top xyouts,0.1,0.99,timerange[0]+' UT--'+timerange[1]+' UT',/norm,charsize=charsize ;plot B_total and Bz y1=floor(min(bz_gse)/10.0)*10 ;lower boundary is determined by Bz y2=ceil(max(b_total)/10.0)*10 ;upper bounardy is defined by B yticks=(y2-y1)/10 if n_elements(dst) ne 0 then utplot,time_utc_mfi,median(b_total,box),timerange=timerange, yrange=[y1,y2], ytitle= 'B(nT)',ystyle=1, xstyle=1, yticks=yticks,xtitle=' ', xticklen=0.0,xchar= 0.01 if n_elements(dst) eq 0 then utplot,time_utc_mfi,median(b_total,box),timerange=timerange, yrange=[y1,y2], ytitle= 'B(nT)',ystyle=1, xstyle=1, yticks=yticks,xtitle=' ', xticklen=0.0,xchar= 0.01, ymargin=[0,2] 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 ;over-plot Bz outplot,time_utc_mfi,median(bz_gse,box),color=color_red axis,yaxis=1,ytitle='Bz(nT)',color=color_red ;plot zero line of B outplot,timerange,[0,0],thick=thick_vline,color=color_red,linestyle=2 ;y1=floor(min(bz_gse)/10.0)*10 & y2=ceil(max(bz_gse)/10.0)*10 & yticks=(y2-y1)/10 ;outplot,time_utc_mfi,median(bz_gse,box),timerange=timerange,ytitle= 'Bz(nT)',xtitle=' ', ystyle=1,xstyle=1,xchar=0.01, yrange=[y1,y2],yticks=yticks ;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 vel total tmpv=where(v_total gt 1000000 or v_total lt 0) if tmpv(0) ne -1 then v_total(tmpv)=0 y1=floor(min(v_total)/100.0)*100 & y2=ceil(max(v_total)/100.0)*100 & yticks=(y2-y1)/100 ;print, v_total utplot,time_utc_plasma, median(v_total,box),timerange=timerange, xtitle='',ytitle= 'V(km/s)', xstyle=1,ystyle=1,xchar= 0.01, yrange=[y1,y2],yticks=yticks 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 proton density tmpn=where(np gt 1000000 or np lt 0) if tmpn(0) ne -1 then np(tmpn)=0 tmpn2=(max(np)-np)/max(np) while n_elements(where(tmpn2 lt .1)) lt 10 do begin np(where(np eq max(np)))=average(np) tmpn2=(max(np)-np)/max(np) endwhile y1=floor(min(np)/10.0)*10 & y2=ceil(max(np)/10.0)*10 & yticks=(y2-y1)/10 utplot,time_utc_plasma,median(np,box),ytitle= 'n(cm!E-3!N)',timerange=timerange,xtitle=' ',xstyle=1,ystyle=1,xchar= 0.01,yrange=[y1,y2],yticks=yticks 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 proton temperature tpr_p=tpr/10^5.0 tmpt=where(tpr_p gt 1000000 or tpr_p lt 0) if tmpt(0) ne -1 then tpr_p(tmpt)=0 tmpt2=(max(tpr_p)-tpr_p)/max(tpr_p) while n_elements(where(tmpt2 lt .1)) lt 10 do begin tpr_p(where(tpr_p eq max(tpr_p)))=average(tpr_p) tmpt2=(max(tpr_p)-tpr_p)/max(tpr_p) endwhile y1=floor(min(tpr_p)/1.0)*1 & y2=ceil(max(tpr_p)/1.0)*1 & yticks=(y2-y1)/1 utplot,time_utc_plasma,median(tpr_p,box),timerange=timerange, ytitle= 'Tp(10!E5!N K)',xtitle='',xstyle=1,ystyle=9,xchar= 0.01,yrange=[y1,y2],yticks=yticks 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 and overlay the expected temperature outplot,time_utc_plasma,tpr_exp/10^5.0, color=color_red axis,yaxis=1,ytitle='Tp_exp (10!E5!N K)',color=color_red ;Plot plasma bet beta2=median(beta2,box) y1=0 & y2=min([ceil(max(beta2)/2.0)*2.0,1.5]) & yticks=(y2-y1)/4.0 print, yticks xtitle=timerange[0]+' UT--'+timerange[1]+' UT' utplot,time_utc_mfi,beta2,timerange=timerange ,ytitle= 'plasma '+string('!7b!3') ,xtitle=xtitle,xstyle=1,ystyle=1,yrange=[y1,y2];,yticks=yticks outplot,[time_icme_end,time_icme_end],[y1,y2],thick=thick_vline,color=color_blue,linestyle=2 ;plot the 0.3 threshold line outplot,timerange,[0.1,0.1],thick=thick_vline,color=color_red,linestyle=2 if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[y1,y2],thick=thick_vline,color=color_red ;;stop 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 ratio of tpr and tpr_exp ;tpr_ratio=median(tpr/tpr_exp,box) ;y1=floor(min(tpr_ratio)/1.0)*1 & y2=ceil(max(tpr_ratio)/1.0)*1 & yticks=(y2-y1)/1 ;utplot,time_utc_plasma,tpr_ratio,timerange=timerange, ytitle= 'Tp/Tp_exp',xtitle='',xstyle=1,ystyle=1,xchar= 0.01, yrange=[y1,y2];,yticks=yticks ;plot the 0.5 threshold line ;outplot,timerange,[0.5,0.5],thick=thick_vline,color=color_red,linestyle=2 ;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 jmap=0 if jmap eq 1 then begin loadct, 0 obj=obj_new('IDL_Savefile', '/home/phess4/jmaps/20100403/jmaparray.sav') obj->restore, 'm' ;img=congrid(m, 300, 200) img=congrid(m,504, 300) ;tv, img, 48, 83 tv, img, 100,200 a=time_utc_mfi(2) b=0 utplot, a,b, timerange=timerange, ytitle='Elongation Angle', xtitle=xtitle, xstyle=1, ystyle=1, yrange=[0,90], position=[.06,.078,.96,.365] ; utplot, a,b, timerange=timerange, ytitle='Elongation Angle', xtitle=xtitle, xstyle=1, ystyle=1, yrange=[0,90], position=[.078,.078,.92,.365] ;outplot, [timerange(0), timerange(1)], [65,65], thick=thick_vline, linestyle=2 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 if keyword_set(time_shock_start) then outplot,[time_shock_start,time_shock_start],[0,90],thick=thick_vline,color=color_red if keyword_set(time_icme_start) then outplot,[time_icme_start,time_icme_start],[0,90],thick=thick_vline,color=color_blue if keyword_set(time_icme_end) then outplot,[time_icme_end,time_icme_end],[0,90],thick=thick_vline,color=color_blue,linestyle=2 ;write_png, 'insitujmap20100403.png', tvrd(/true) end test=1 if test ne 0 then begin ; ;--------Derive relevant solar wind parameters------- ; ;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.234 ;starting position for the text dx_v=0.18 ;seperating between variable name and value dx=0.5 ;seperating between columns dy=0.018 ;seperating between rows xyouts,x00,y00,'--------------------------------------------------------------------------------',/norm,charsize=charsize ;print Dst peal value and Dst peak time, if present x0=x00 & y0=y00-1.0*dy xyouts,x0,y0,'Dst: ',/norm,charsize=charsize if keyword_set(dst_peak) then begin xyouts,x0+dx_v,y0,dst_peak,/norm,charsize=charsize print,'DST: ',dst_peak if keyword_set(par_file) then printf,lu,'DST: ',dst_peak endif x0=x00+dx & y0=y00-1.0*dy xyouts,x0,y0,'DST_PEAK_TIME: ',/norm,charsize=charsize if keyword_set(time_dst_peak) then begin xyouts,x0+dx_v,y0,time_dst_peak,/norm,charsize=charsize print,'DST_PEAK_TIME: ',time_dst_peak if keyword_set(par_file) then printf,lu,'DST_PEAK_TIME: ',time_dst_peak endif ;print solar wind stream type ;e.g., SH+MC, SH+ICME, CIR, ICME(M) ;done externally, from the optional input x0=x00 & y0=y00-2.0*dy xyouts,x0,y0,'SW_TYPE: ',/norm,charsize=charsize if keyword_set(sw_type) then begin xyouts,x0+dx_v,y0,sw_type,/norm,charsize=charsize print,'SW_TYPE: ',sw_type if keyword_set(par_file) then printf,lu,'SW_TYPE: ',sw_type endif ;print solar wind drivers of storm, e.g., which part in the solar wind x0=x00+dx & y0=y00-2.0*dy xyouts,x0,y0,'DRIVER_TYPE: ',/norm,charsize=charsize if keyword_set(sw_dr_type) then begin xyouts,x0+dx_v,y0,sw_dr_type,/norm,charsize=charsize print,'DRIVER_TYPE: ',sw_dr_type if keyword_set(par_file) then printf,lu,'DRIER_TYPE: ',sw_dr_type endif ;print the shock start time, if present x0=x00 & y0=y00-3.0*dy xyouts,x0,y0,'SH_START_TIME: ',/norm,charsize=charsize if keyword_set(time_shock_start) then begin xyouts,x0+dx_v,y0,time_shock_start,/norm,charsize=charsize print,'SH_START_TIME: ',time_shock_start if keyword_set(par_file) then printf,lu,'SH_START_TIME: ',time_shock_start endif ;print the shock end time, if present ;usually the same as ICME start time x0=x00 & y0=y00-4.0*dy xyouts,x0,y0,'SH_END_TIME: ',/norm,charsize=charsize if not keyword_set(time_shock_end) then begin if keyword_set(time_icme_start) then time_shock_end=time_icme_start endif if keyword_set(time_shock_end) then begin xyouts,x0+dx_v,y0,time_shock_end,/norm,charsize=charsize print,'SH_END_TIME: ',time_shock_END if keyword_set(par_file) then printf,lu,'SH_END_TIME: ',time_shock_end endif ;print the icme start time, if present x0=x00+0.5 & y0=y00-3.0*dy xyouts,x0,y0,'ICME_START_TIME: ',/norm,charsize=charsize if keyword_set(time_icme_start) then begin xyouts,x0+dx_v,y0,time_icme_start,/norm,charsize=charsize print,'ICME_START_TIME: ',time_icme_start if keyword_set(par_file) then printf,lu,'ICME_START_TIME: ',time_icme_start endif ;print the icme end time, if present x0=x00+0.5 & y0=y00-4.0*dy xyouts,x0,y0,'ICME_END_TIME: ',/norm,charsize=charsize if keyword_set(time_icme_end) then begin xyouts,x0+dx_v,y0,time_icme_end,/norm,charsize=charsize print,'ICME_END_TIME: ',time_icme_END if keyword_set(par_file) then printf,lu,'ICME_END_TIME: ',time_icme_end endif ;calculate shock length in unit of hour and AU, and energy input by shock E_plasma=interpol(E,time_tai_mfi,time_tai_plasma) if keyword_set(time_shock_start) and keyword_set(time_shock_end) then begin time_tai_sh_start=utc2tai(str2utc(time_shock_start)) time_tai_sh_end=utc2tai(str2utc(time_shock_end)) time_shock_sheath=(time_tai_sh_end-time_tai_sh_start)/3600.0 ; in hr time_sh_str=strtrim(string(time_shock_sheath,format='(f5.2)'),2) ;calculate the linear size L=V*T arr_sh=where(time_tai_plasma ge time_tai_sh_start and time_tai_plasma le time_tai_sh_end,count) sh_size=total(v_total(arr_sh)) ;size in km in time resulution /sec time_step=time_tai_plasma[2]-time_tai_plasma[1] sh_size=sh_size*time_step ; in unit of km sh_size=sh_size/1.5e8 ; in unit of AU sh_size_str=strtrim(string(sh_size,format='(f4.2)'),2) sh_size_str=time_sh_str+'(hr) or '+sh_size_str+'(AU)' ;calcuate energy input energy_sh=total(E_plasma(arr_sh))*time_step ;in Joule energy_sh_str=strtrim(string(energy_sh,format='(e8.2)'),2) endif ;calculate icme length in unit of hour and AU 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)) time_icme=(time_tai_icme_end-time_tai_icme_start)/3600.0 ; in hr time_icme_str=strtrim(string(time_icme,format='(f5.2)'),2) ;calculate the linear size L=V*T arr_icme=where(time_tai_plasma ge time_tai_icme_start and time_tai_plasma le time_tai_icme_end,count) icme_size=total(v_total(arr_icme)) ;size in km in time resulution /sec time_step=time_tai_plasma[2]-time_tai_plasma[1] icme_size=icme_size*time_step ; in unit of km icme_size=icme_size/1.5e8 ; in unit of AU icme_size_str=strtrim(string(icme_size,format='(f4.2)'),2) icme_size_str=time_icme_str+'(hr) or '+icme_size_str+'(AU)' ;stop ;calcuate energy input energy_icme=total(E_plasma(arr_icme))*time_step ;in Joule energy_icme_str=strtrim(string(energy_icme,format='(e8.2)'),2) endif ;print shock size and icme size x0=x00 & y0=y00-5.0*dy xyouts,x0,y0,'SH_SIZE: ',/norm,charsize=charsize if keyword_set(sh_size_str) then begin xyouts,x0+dx_v,y0,sh_size_str,/norm,charsize=charsize print,'SH_SIZE: ',sh_size_str if keyword_set(par_file) then printf,lu,'SH_SIZE: ',sh_size_str endif x0=x00+dx & y0=y00-5.0*dy xyouts,x0,y0,'ICME_SIZE: ',/norm,charsize=charsize if keyword_set(icme_size_str) then begin xyouts,x0+dx_v,y0,icme_size_str,/norm,charsize=charsize print,'ICME_SIZE: ',icme_size_str if keyword_set(par_file) then printf,lu,'ICME_SIZE: ',icme_size_str endif ;print energy input, by SH, by ICME, and their percentage energy=0 if keyword_set(energy_sh) then energy=energy+energy_sh else energy_sh=0 if keyword_set(energy_icme) then energy=energy+energy_icme else energy_icme=0 if keyword_set(energy) then energy_str=strtrim(string(energy,format='(e8.2)'),2)+'(J)' ;calcuate percentage if keyword_set(energy) then begin e_sh_per=energy_sh/energy*100.0 ;percentage of shock energy in total energy input e_sh_per_str=strtrim(string(e_sh_per,format='(f4.1)'),2)+'%' e_icme_per=energy_icme/energy*100.0 ;percentage of icme energy in total energy input e_icme_per_str=strtrim(string(e_icme_per,format='(f4.1)'),2)+'%' endif ;print energy x0=x00 & y0=y00-7.0*dy xyouts,x0,y0,'ENERGY: ',/norm,charsize=charsize if keyword_set(energy) then begin xyouts,x0+dx_v,y0,energy_str,/norm,charsize=charsize print,'ENERGY: ',energy_str if keyword_set(par_file) then printf,lu,'ENERGY: ',energy_str endif x0=x00+dx & y0=y00-7.0*dy xyouts,x0,y0,'SIZE: ',/norm,charsize=charsize if keyword_set(time_shock_start) then begin total_size = sh_size + icme_size time_total_size = time_icme + time_shock_sheath endif else begin total_size=icme_size time_total_size = time_icme endelse total_size_str = strtrim(string(total_size,format='(f4.2)'),2) time_total_size_str = strtrim(string(time_total_size,format='(f5.2)'),2) xyouts,x0+dx_v,y0,time_total_size_str+'(hr) or '+total_size_str+'(AU)',/norm,charsize=charsize print,'SIZE: ',total_size_str x0=x00 & y0=y00-6.0*dy xyouts,x0,y0,'ENERGY_SH: ',/norm,charsize=charsize if keyword_set(energy_sh) then begin e_sh_str=energy_sh_str+'(J) ('+e_sh_per_str+')' xyouts,x0+dx_v,y0,e_sh_str,/norm,charsize=charsize print,'ENERGY_SH: ',e_sh_str if keyword_set(par_file) then printf,lu,'ENERGY_SH: ',e_sh_str endif x0=x00+dx & y0=y00-6.0*dy xyouts,x0,y0,'ENERGY_ICME: ',/norm,charsize=charsize if keyword_set(energy_ICME) then begin e_icme_str=energy_icme_str+'(J) ('+e_icme_per_str+')' xyouts,x0+dx_v,y0,e_icme_str,/norm,charsize=charsize print,'ENERGY_ICME: ',e_icme_str if keyword_set(par_file) then printf,lu,'ENERGY_ICME: ',e_icme_str endif 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) np_sh=np_mfi(arr_sh) tpr_p_sh=tpr_mfi(arr_sh) beta2_sh=beta2(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) np_icme=np_mfi(arr_icme) tpr_p_icme=tpr_mfi(arr_icme) beta2_icme=beta2(arr_icme) endif ; ;--------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.037 & y00=0.135 ;starting position for the text dx_s=0.092 ;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 ;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+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 np in SH and ICME x0=x00 & y0=y00-4.0*dy xyouts,x0,y0,'Np(cm!E-3!N): ',/norm,charsize=charsize if keyword_set(np_sh) then begin ;calculate Np in shock tmp=np_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,'Np_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' if keyword_set(par_file) then printf,lu,'Np_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' endif x0=x00+dx & y0=y00-4.0*dy if keyword_set(np_icme) then begin ;calculate Np in icme tmp=np_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,'Np_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' if keyword_set(par_file) then printf,lu,'Np_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' endif ;print np in SH and ICME x0=x00 & y0=y00-5.0*dy xyouts,x0,y0,'Tp(10!E5!N K): ',/norm,charsize=charsize if keyword_set(tpr_p_sh) then begin ;calculate Tp in shock tmp=tpr_p_sh/1e4 mom=moment(tmp,sdev=sdev) ave_str=strtrim(string(mom[0],format='(f5.2)')) sdev_str=strtrim(string(sdev,format='(f5.2)')) max_str=strtrim(string(max(tmp),format='(f5.2)')) min_str=strtrim(string(min(tmp),format='(f5.2)')) 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,'Tp_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' if keyword_set(par_file) then printf,lu,'Tp_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' endif x0=x00+dx & y0=y00-5.0*dy if keyword_set(tpr_p_icme) then begin ;calculate Tp in icme tmp=tpr_p_icme/1e4 mom=moment(tmp,sdev=sdev) ave_str=strtrim(string(mom[0],format='(f5.2)')) sdev_str=strtrim(string(sdev,format='(f5.2)')) max_str=strtrim(string(max(tmp),format='(f5.2)')) min_str=strtrim(string(min(tmp),format='(f5.2)')) 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,'Tp_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' if keyword_set(par_file) then printf,lu,'Tp_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' endif ;print beta2 in SH and ICME x0=x00 & y0=y00-6.0*dy xyouts,x0,y0,'Plas beta2: ',/norm,charsize=charsize if keyword_set(beta2_sh) then begin ;calculate Np in shock tmp=beta2_sh mom=moment(tmp,sdev=sdev) ave_str=strtrim(string(mom[0],format='(f5.2)')) sdev_str=strtrim(string(sdev,format='(f5.2)')) max_str=strtrim(string(max(tmp),format='(f5.2)')) min_str=strtrim(string(min(tmp),format='(f5.2)')) 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,'Beta_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' if keyword_set(par_file) then printf,lu,'Beta_SH: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' endif x0=x00+dx & y0=y00-6.0*dy if keyword_set(beta2_icme) then begin ;calculate beta2 in icme tmp=beta2_icme mom=moment(tmp,sdev=sdev) ave_str=strtrim(string(mom[0],format='(f5.2)')) sdev_str=strtrim(string(sdev,format='(f5.1)')) max_str=strtrim(string(max(tmp),format='(f5.2)')) min_str=strtrim(string(min(tmp),format='(f5.2)')) 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,'Beta_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' if keyword_set(par_file) then printf,lu,'Beta_ICME: ',ave_str+' '+sdev_str+' '+max_str+' '+min_str+' ' endif 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/events0716/'+fdate+'/ -d', result ;if result[0] eq '' then spawn, 'mkdir /home/phess4/solarwind/bz/events0716/'+fdate spawn,'ls '+year,result3 if result3[0] eq '' then spawn, 'mkdir '+year spawn,'ls '+year+'/plots',result4 if result4[0] eq '' then spawn, 'mkdir '+year+'/plots' write_png,year+'/plots/'+fdate+'.png',tvrd(/true) ;close the parameter file if keyword_set(lu) then close,lu ;if keyword_set(hc) then device,/close empty SET_PLOT, 'X' return end