; Plot Geotail trajectories in Y-Z, and X-Z planes overlaid
; with model bowshocks, magnetopauses, etc.
;
; 10/09/97 SHC adapted from orbit_hk.pro also by SHC
;
;-------subroutine aberration
;pro aberration,x,y,angle=angle
;    if(n_elements(angle) eq 0) then angle=5.0 ;degrees, assuming vsw=400 km/s
;    angle_rad=-angle*0.017453292 ;!PI/180.
;    a1=cos(angle_rad)
;    a2=sin(angle_rad)
;    x_buf= x*a1+y*a2
;    y_buf=-x*a2+y*a1
;    x=x_buf
;    y=y_buf
;    return
;end
;------ subroutine get_mp
;
pro get_mp,xmp,ymp,n,standoff_bs,pdyn_index
    path='/export/local/www/space/spdf/bowshock/idl/'
    theta1=320.d0
    theta2=320.d0
    DegtoRad=!PI/180.d0
    dtheta1=theta1/(n-1)
    dtheta2=theta2/(n-1)
    theta1=-160.d0
    theta2=-160.d0
    r_bs=0.
    x1=0.
    y1=0.
    z1=0.
    rnorm=0.
    bs_factor=1.0
    standoff_bs2=0.0
    n2=n*2
    n3=n*3

    if(pdyn_index eq 1) then bs_factor=1.0 $
    else bs_factor=(1.8/8.0)^(1.0/6.0)

    bs_factorx=(1.8/7.0)^(1.0/6.0)

    s1=call_external(path+'bs_mp.so','bowshock_peredo_idl',$
        10.,2.0,13.,0.0,0.0,standoff_bs2,y1,z1,rnorm,2L)

    for i=0,n-1 do begin
	if(pdyn_index eq 1) then begin
	    xmp(i)=11.000d0*cos(theta1*DegtoRad)
	    if(xmp(i) lt 0.) then xmp(i)=xmp(i)*2.0
	    r=217.2d0-0.14*xmp(i)*xmp(i)-18.2*xmp(i)
	    if(r lt 0.) then r=0.
	    ymp(i)=sqrt(r) 			; Sibeck's mp model
	    					; for 1.47<Pdyn<2.6nPa
	endif else begin	 
	    xmp(i)=8.815*cos(theta1*DegtoRad)
	    if(xmp(i) lt 0.) then xmp(i)=xmp(i)*2.3
	    r=139.2-0.18*xmp(i)*xmp(i)-14.2*xmp(i)
	    if(r lt 0.) then r=0.
	    ymp(i)=sqrt(r)			; Sibeck's mp model
						; for 4.49<Pdyn<9.9nPa
	endelse

	if(theta1 lt 0.) then ymp(i)=-ymp(i) 

; 	k=i+n
;	xmp(k)=10.3d0*cos(theta2*DegtoRad)	;Farris's mp mdel
;						;for Pdyn=1.8 or 7.5 nPa
;	if(xmp(k) lt 0.) then xmp(k)=xmp(k)*1.5
;	r=14.7-0.43*xmp(k)
;	r=r*r-xmp(k)*xmp(k)
;	if(r lt 0.) then r=0.
;	xmp(k)=xmp(k)*bs_factor
;	ymp(k)=sqrt(r)


	k=i+n

	xmp(k)=8.815*cos(theta1*DegtoRad)
	if(xmp(k) lt 0.) then xmp(k)=xmp(k)*2.3
	r=139.2-0.18*xmp(k)*xmp(k)-14.2*xmp(k)
	if(r lt 0.) then r=0.
	ymp(k)=sqrt(r)			; Sibeck's mp model
					; for 4.49<Pdyn<9.9nPa
	if(theta1 lt 0.) then ymp(k)=-ymp(k) 

	;xmp(k)=standoff_bs2*cos(theta2*DegtoRad)
	;if(xmp(k) lt 0.) then xmp(k)=xmp(k)*1.5
	;x_buf=float(xmp(k))
	;s2=call_external('bs_mp.so','bowshock_peredo_idl',$
	;	10.,2.,x_buf,9999.,0.0,x1,y1,z1,r_bs,1L) ; solve for y
	;if(y1 ge 9999.0) then ymp(k)=0. $
	;else ymp(k)=y1

	;if(theta2 lt 0.) then ymp(k)=-ymp(k) 

	k=i+n2

	xmp(k)=13.7*cos(theta2*DegtoRad) ; Farris' bs model
	if(xmp(k) lt 0.) then xmp(k)=xmp(k)*1.5
	x_buf=float(xmp(k))
	s2=call_external(path+'mp_bs.so','position_bs_mp_idl',$
	    2L,2.0,0.0,x_buf,r_bs,0L); 
	if(r_bs ge 9999.0) then ymp(k)=0. $
	else ymp(k)=r_bs


	if(theta2 lt 0.) then ymp(k)=-ymp(k) 

	k=i+n3

	ymp(k)=ymp(k-n)
	xmp(k)=xmp(k-n)-13.7*(1.-bs_factorx) 

	;xmp(k)=standoff_bs*cos(theta2*DegtoRad)
	;if(xmp(k) lt 0.) then xmp(k)=xmp(k)*1.5
	;x_buf=float(xmp(k))
	;if(pdyn_index eq 1) then begin
	;    s2=call_external('bs_mp.so','bowshock_peredo_idl',$
	;	5.,2.,x_buf,9999.,0.0,x1,y1,z1,r_bs,1L) ; solve for y
	;endif else begin
	;    s2=call_external('bs_mp.so','bowshock_peredo_idl',$
	;	15.,8.,x_buf,9999.,0.0,x1,y1,z1,r_bs,1L) ; solve for y
	;endelse
	;if(y1 ge 9999.0) then ymp(k)=0. $
	;else ymp(k)=y1

	;if(theta2 lt 0.) then ymp(k)=-ymp(k) 

;	print,k,standoff_bs,xmp(k),ymp(k)

	theta1=theta1+dtheta1
	theta2=theta2+dtheta2

    endfor
;
; correction for the flarring angle of Farris's bowshock at Pdyn=7 over 2 nPa
;
    for i=0,n-1 do begin
	k2=i+n2
	k3=i+n3
	ratio=(13.7-xmp(k2))/(13.7-xmp(n3-1))
;	ymp(k3)=ymp(k3)-0.46*ymp(n3-1)*ratio*ratio
	if(ymp(k3) ge 0) then ymp(k3)=ymp(k3)-0.11*ymp(n3-1)*ratio $
	else ymp(k3)=ymp(k3)+0.11*ymp(n3-1)*ratio
    endfor

    theta1=0.
    dtheta1=18.
    for i=0,20 do begin
	k=k+1
	xmp(k)=cos(theta1*DegtoRad)
	ymp(k)=sqrt(1.-xmp(k)*xmp(k))
	if(theta1 gt 180.) then ymp(k)=-ymp(k)
	theta1=theta1+dtheta1
    endfor
return
end
;
;---------- get_circle
;
pro get_circle,xmp,ymp,n,pdyn_index
DegtoRad=!PI/180.
theta=360.d0
theta2=360.d0
dtheta=theta/(n-1)
dtheta2=theta2/(n-1)
theta=0.d0
theta2=0.d0
    for i=0,n-1 do begin
	if(pdyn_index eq 1) then begin
	    xmp(i)=14.73*sin(theta*DegtoRad)
	    ymp(i)=14.73*cos(theta*DegtoRad)
	endif else begin
;	    xmp(i)=13.1*sin(theta*DegtoRad) ;at x=-2.3 re
;	    ymp(i)=13.1*cos(theta*DegtoRad)
	    xmp(i)=11.8*sin(theta*DegtoRad) ;at x=0
	    ymp(i)=11.8*cos(theta*DegtoRad)
	endelse
	k=i+n
	if(pdyn_index eq 1) then begin
	    xmp(k)=14.7*sin(theta*DegtoRad)	;Farris's mp mdel
					;for Pdyn=1.8 nPa
	    ymp(k)=14.7*cos(theta*DegtoRad)
	endif else begin
	    xmp(k)=16.1*sin(theta*DegtoRad)
	    ymp(k)=16.1*cos(theta*DegtoRad)
	endelse

	k=i+2*n
	if(pdyn_index eq 1) then begin
	    xmp(k)=24.8*sin(theta2*DegtoRad);Farris's bw mdel
					;for Pdyn=1.8 nPa
	    ymp(k)=24.8*cos(theta2*DegtoRad)
	endif else begin
	    xmp(k)=22.5*sin(theta2*DegtoRad);Farris's bw mdel
					;for Pdyn=8. nPa
	    ymp(k)=22.5*cos(theta2*DegtoRad)


	endelse

	theta=theta+dtheta
	theta2=theta2+dtheta2
    endfor

theta=0.
dtheta=18.
    for i=0,20 do begin
	k=k+1
	xmp(k)=sin(theta*DegtoRad)
	ymp(k)=cos(theta*DegtoRad)
	theta=theta+dtheta
    endfor
return
end
;
;--------subroutine get_scaling
;
pro get_scaling_roelofmp,standoff,standoff_bs,xmp,ymp,n,index
    path='/export/local/www/space/spdf/bowshock/idl/'
    bz0=0.
    ma=10.0
    pdyn0=2.0
    x1=0.
    y1=0.
    z1=0.
    r_mp=0.
    r_bs=0.
    index=0L
    theta1=160.0
    theta2=90.0
    DegtoRad=!PI/180.0
    dtheta1=theta1/(n-1)
    dtheta2=theta2/(n-1)
    theta1=0.0
    theta2=0.0
    ;print,' MP stand-off=',standoff,'   BS stand-off bs=',standoff_bs
    for i=0,n-1 do begin
	xmp(i)=standoff*cos(theta1*DegtoRad)-i*0.05
	x_buf=float(xmp(i))
	s=call_external(path+'mp_roelof.so','magnetopause_roelof_idl',$
	    bz0,pdyn0,x_buf,9999.,9999.,x1,y1,z1,r_mp,index)
	ymp(i)=y1

	k=i+n
	xmp(k)=10.3*cos(theta2*DegtoRad)	;Farris's mp mdel
						;for Pdyn=1.8 nPa
	r=14.7-0.43*xmp(k)
	r=r*r-xmp(k)*xmp(k)
	if(r lt 0.) then r=0.
	ymp(k)=sqrt(r)

	k=i+2*n
	xmp(k)=standoff_bs*cos(theta2*DegtoRad)
	x_buf=float(xmp(k))
	s2=call_external(path+'bs_mp.so','bowshock_peredo_idl',$
		ma,pdyn0,x_buf,0.0,9999.0,x1,y1,z1,r_bs,0L) ; solve for z
	;print,xmp(k),x1,y1,z1
	if(z1 ge 9999.0) then ymp(k)=0. $
	else ymp(k)=z1

	theta1=theta1+dtheta1
	theta2=theta2+dtheta2
    endfor

    theta1=0.
    dtheta1=9.
    for i=0,10 do begin
	k=k+1
	xmp(k)=cos(theta1*DegtoRad)
	ymp(k)=sqrt(1.-xmp(k)*xmp(k))
	theta1=theta1+dtheta1
    endfor
return
end
;
;- functions and procedures -------------------------------------------
;
;function fun_epoch,year,doy,tod
;yr=year
;if n_params() eq 0 then begin
;   print,'function,epoch0:',year,doy,tod
;   return,-1
;endif
;;ep0 is epoch at start of 20th cent.
;;tod - time of doy in seconds
;;doy - day of year starting at 1
;;year - ie 1974 is 74
;if(yr gt 1000) then yr=yr-1900
;ep0 = 5.99582304d13
;return,(yr*365d0 + fix((yr-1)/4)+ doy-1)*8.64d7 +tod*1000d0 + ep0 ;in millisec
;end
;------------------------------------------------------------------------
; to transfer epoch0 time in msec to yr doy hr min sec
;
;pro yrdhms,epoch0,yr,doy,hr,min,sec             ; in milliseconds
;ep0 = 5.99582304d13             ;  1900 jan 1 00:00:00 UT
;ep1=ep0+3.15576d12              ; year of 2000
;if(epoch0 lt ep0) or (epoch0 gt ep1) then begin
;        print,'func-yrdhms: out of range'
;        return
;endif
;depoch=0.001d0*(epoch0-ep0)
;yr=long(depoch/31536000.d0)     ; start at the year of 1900
;doy=long(depoch/86400.d0)-yr*365-fix((yr-1)/4)+1
;if(doy le 0) then begin
;        yr=yr-1
;        doy=365+doy
;endif
;days=yr*365+fix((yr-1)/4)+doy-1
;hr=long(depoch/3600.d0-24.d0*days)
;min=long(depoch/60.d0-1440.d0*days-60.d0*hr)
;sec=depoch-86400.d0*days-hr*3600.d0-min*60.d0
;return
;end
;---------------------------------------------------------------------
;			Main program 
;---------------------------------------------------------------------
;
;Device setup
;
;close,/all
path='/export/local/www/space/spdf/bowshock/idl/'
epoch0=0.d0
c_or_t=0
tilt=0
pdyn_index=1
npage=0
nlegend=0
draw_field=0
temp_char=' '
theta=0.0
sm=0
yr=0
doy=0
hr=0
min=0
sec=0.0
yr=0
doy=0
mon=0
day=0
hr=0
min=0
sec=0.0
filename=strarr(1)
legend=strarr(6)
tra_sheath=1
;
;
infile='/export/local/www/space/spdf/bowshock/idl/orbit_geotail.feed'
openr,lun_feed,infile,/get_lun
;
;read,'Read hk_tra or sheath_omni file?(1 - hk_tra; 2-sheath_omni) ',tra_sheath
;read,' Input file for Geotail? ',filename
;openr,lun_geo,filename(0),/get_lun
;read,' Input file for IMP 8? ',filename
;openr,lun_imp,filename(0),/get_lun

;read,' Time? (yr,mon,day,hr,min,sec) ',yr,mon,day,hr,min,sec
;read,yr,mon,day,hr,min,sec
readf,lun_feed,yr,mon,day,hr,min,sec
;
;read,' circles (0) or tangent (1) or both (2)?',c_or_t
c_or_t=0
;read,' GSM or SM coordinate ( 0 - GSM; 1 - SM)? ',sm
;read,' Conserve the anlge while proj. into x-z plane?(no 0 or yes 1)',conserve
conserve=0
;read,' Draw field lines? (0 - No; 1 - Yes) ',draw_field

;if(draw_field eq 1) then $
;   read,'0 tilt (0), 20 deg tilt (1), for 7418418 (2) ,or for 74262 (3)? ',tilt
;
;read,' dynamic pressure: 2 nPa (1)  5-10 nPa (2)? ',pdyn_index
;
pdyn_index=1
;read,' Legend: # of lines (0-5)? ',nlegend
;nlegend=0
;for i=1,nlegend do begin
;	print,'line',i,':'
;	read,temp_char
;	legend(i)=temp_char
;endfor
;read,' 0 - Draw both pages; 1 - X-Z plane | 2 - Y-Z plane? ',npage 
;
;read,' Draw every how many points? ',everynpts
;if(everynpts lt 1) then everynpts=1
everynpts=1

;read,' Output device (0 - postcript; 1 - X window; 2 - Zbuf)? ',idevice 
;read,idevice
readf,lun_feed,idevice
aberr=1 ; consider aberration correction due to the motion of the Earth
readf,lun_feed,aberr
angle=5.0

; 
if (idevice eq 0) then begin 
;*** for poscript ***************************************************
;	read,' Encapsulated (eps)? (0-yes; 1-no)',enca
	enca=1
        read,' Output file name ? ',filename
        set_plot,'PS'
	if(enca eq 0) then begin
            device,filename=filename(0),xoffset=2,xsize=18.5,yoffset=2,$
		ysize=24,/HELVETICA,/BOLD,FONT_SIZE=12,/encapsul 
	endif else begin
	    device,filename=filename(0),xoffset=2,xsize=18.5,yoffset=2,$
		ysize=24,/HELVETICA,/BOLD,FONT_SIZE=12
	endelse
	print,'The output file name will be ',filename
	charthick=4.0
        !p.charthick=charthick
        !p.charsize=1.5
        !x.thick=6.
        !y.thick=6.
	xyratio=8.5/11.
	xleft=0.12*xyratio	; In normal coordinate [0,1]
	xright=0.81*xyratio
	ybottom=0.05    ; In normal coordinate [0,1].
	ytop=0.95
	ymid=0.55;
;	ybottom=0.12*xyratio    ; In normal coordinate [0,1].
;	ytop=0.95*xyratio
	col=0	;color of the nightside earth
        charsize=1.
endif else begin
if(idevice eq 1) then begin
;
;*** for X-window ***************************************************
        set_plot,'X'
	loadct,5 ; standard gamma-II color table
	charsize=0.8
	!p.background=255
	!p.color=0
        !p.charsize=charsize
        !p.charthick=2.
        !x.thick=2.
        !y.thick=2.
	charthick=1.0
	xleft=0.15	; In normal coordinate [0,1]
	xright=0.95
	ymid=0.55;
	ybottom=0.05    ; In normal coordinate [0,1].
	ytop=0.95
	col=0
endif else begin
;        read,' Output file name ? ',filename
	;filename='delme.cgm'
	;print,'The output file name will be ',filename
	;set_plot,'CGM'
	;device,filename=filename(0)
	set_plot,'Z'
	loadct,5 ; standard gamma-II color table
		 ; 0 - black; 64 - blue; 
		 ; 96- dark red;
 		 ; 128 - light red; 160 - green; 192 - yellow 
		 ; 255 - black
	charsize=1.2
	charthick=1.0
	!p.background=255
	!p.color=0
        !p.charsize=charsize
        !p.charthick=1.0
        !x.thick=2
        !y.thick=2
	xyratio=0.8
	;xleft=0.2	; In normal coordinate [0,1]
	;xright=0.6*xyratio
	xleft=0.1
	xright=0.78*xyratio
	ymid=0.55;
	;ybottom=0.08   ; In normal coordinate [0,1].
	;ytop=0.93
	ybottom=0.2
	ytop=0.88

	;xleft=0.12	; In normal coordinate [0,1]
	;xright=0.95
	;ybottom=0.12    ; In normal coordinate [0,1].
	;ytop=0.95
endelse
endelse
;*********************************************************************
;
axis_type=1
if(axis_type eq 2) then begin
	xtype=1
	ytype=0
endif else begin
if(axis_type eq 3) then begin
	xtype=0
	ytype=1
endif else begin
if(axis_type eq 4) then begin
	xtype=1
	ytype=1
endif else begin
	xtype=0
	ytype=0
endelse
endelse
endelse
;
N_dimen=100000L
;time_geo=dblarr(N_dimen)
;xloc_geo=fltarr(N_dimen)
;yloc_geo=fltarr(N_dimen)
;zloc_geo=fltarr(N_dimen)
;index_geo=intarr(N_dimen)
;pdyn=fltarr(N_dimen)
;bz=fltarr(N_dimen)
;
;time_imp=dblarr(N_dimen)
;xloc_imp=fltarr(N_dimen)
;yloc_imp=fltarr(N_dimen)
;zloc_imp=fltarr(N_dimen)
;index_imp=intarr(N_dimen)
;
; *** Read data...
;
npoints_geo=0L
;npoints_imp=0L
nlines_geo=0L
;nlines_imp=0L
time_buf=0.d0
start0=0.d0
end0=0.d0
temp1=0.0
temp2=0.0
temp3=0.0
temp4=0.0
temp5=0.0
temp6=0.0
temp7=0.0
itemp1=0
re=6371.
xgsm_hk=0.
ygsm_hk=0.
zgsm_hk=0.
r_hk=0.
Pdyn_buf=0.
n_sw=0.
temp_sw=0.
vel_sw=0.
lon_vel_sw=0.
lat_vel_sw=0.
;bx_imp=0.
;by_imp=0.
;bz_imp=0.
;bt_imp=0.
;xgse_imp=0.
;ygse_imp=0.
;zgse_imp=0.
theta_bn=0.
r_mp_buf=0.
;ma0=10.; for 1996 2/15
ma0=10. ; for 1995 4/25
;pdyn0=2.
pdyn0=2.
r_mp=0.
x1=0.
y1=0.
z1=0.
rnorm=0.
standoff=0.
standoff_bs=0.
dawndusk=0.
dawndusk_bs=0.
;
; convert to epoch time 
;
;sec=sec+hr*3600.+min*60.
;start=fun_epoch(yr,doy,sec)
;
;s1=call_external(path+'bs_mp.so','bowshock_peredo_idl',$
;    ma0,pdyn0,13.,0.0,0.0,standoff_bs,y1,z1,rnorm,2L)

;standoff_bs=13.7 ; Farris' model

;s1=call_external(path+'bs_mp.so','bowshock_peredo_idl',$
;    ma0,pdyn0,0.,0.0,0.0,x1,y1,dawndusk_bs,rnorm,0L)
s2=call_external(path+'mp_roelof.so','magnetopause_roelof_idl',$
    0.,pdyn0,9999.,9999.,9999.,standoff,dawndusk,z1,r_mp,0L) 
;
; read geotail
;
;print,'Reading geotail data...'
outfile1=strarr(1)
outfile2=strarr(1)
;
readf,lun_feed,outfile1
readf,lun_feed,outfile2
readf,lun_feed,npoints_geo
if(npoints_geo gt 0) then begin 
    xloc_geo=fltarr(npoints_geo)
    yloc_geo=fltarr(npoints_geo)
    zloc_geo=fltarr(npoints_geo)
    pdyn=fltarr(npoints_geo)
    sctype=fltarr(npoints_geo)
    index_geo=intarr(npoints_geo)
    if(aberr eq 1) then begin
	angle_rad=-angle*0.017453292 ;!PI/180.
	a1=cos(angle_rad)
	a2=sin(angle_rad)
        for i=0,npoints_geo-1 do begin
            readf,lun_feed,temp1,temp2,temp3,temp4,itemp1;<<<<< input <<<<<<<<
	    x_buf= temp1*a1+temp2*a2
	    y_buf=-temp1*a2+temp2*a1
            xloc_geo(i)=x_buf
            yloc_geo(i)=y_buf
            zloc_geo(i)=temp3
    	    pdyn(i)=temp4
	    sctype(i)=itemp1
	    index_geo(i)=0
        endfor
    endif else begin
        for i=0,npoints_geo-1 do begin
            readf,lun_feed,temp1,temp2,temp3,temp4,itemp1;<<<<< input <<<<<<<<
            xloc_geo(i)=temp1
            yloc_geo(i)=temp2
            zloc_geo(i)=temp3
    	    pdyn(i)=temp4
	    sctype(i)=itemp1
	    index_geo(i)=0
        endfor
    endelse
endif
;
;    while(not EOF(lun_geo)) and (npoints_geo lt N_dimen) and $
;		 (time_buf le end0) do begin
;        readu,lun_geo,time_buf,temp1,temp2,temp3,temp4,temp5,temp6
;	if((fix(nlines_geo/everynpts)*everynpts-nlines_geo) eq 0) and $
;		(time_buf ge start) and (time_buf le end0) then begin
;	    ;print,time_buf,temp4,temp5,temp6
;	    time_geo(npoints_geo)=time_buf
;	    xloc_geo(npoints_geo)=temp4
;	    yloc_geo(npoints_geo)=temp5
;	    zloc_geo(npoints_geo)=temp6
;	    index_geo(npoints_geo)=0
;            npoints_geo=npoints_geo+1
;	endif else if(time_buf gt end0) then goto, Continue 
;	nlines_geo=nlines_geo+1
;    endwhile

Continue:
;
; read imp 8
;
time_buf=0.d0
everynpts=1
;npoints_imp=0L
;print,'Reading imp 8 data...'
;    while(not EOF(lun_imp)) and (npoints_imp lt N_dimen) and $
;		 (time_buf le end0) do begin
;	readu,lun_imp,time_buf,temp1,temp2,temp3,temp4,temp5,$
;            xgse_imp,ygse_imp,zgse_imp
;	;print,format='(f20.1,f8.1,f8.1,f8.1)',$
;	;	time_buf,xgse_imp,ygse_imp,zgse_imp
;	if((fix(nlines_imp/everynpts)*everynpts-nlines_imp) eq 0) and $
;		(time_buf ge start) and  (time_buf le end0) then begin
;	    time_imp(npoints_imp)=time_buf
;;	    pdyn(npoints_imp)=pdyn_buf
;;	    bz(npoints_imp)=bz_imp
;;	    if(bz_imp gt 7) then bz_imp=7. $
;;	    else if(bz_imp lt -7.) then bz_imp=-7.
;;    	    s2=call_external('mp_roelof.so','magnetopause_roelof_idl',$
;;		bz_imp,pdyn_buf,9999.,9999.,9999.,x1,y1,z1,r_mp,0L) 
;	    xloc_imp(npoints_imp)=xgse_imp
;	    yloc_imp(npoints_imp)=ygse_imp
;	    zloc_imp(npoints_imp)=zgse_imp
;	    index_imp(npoints_imp)=0
;	    npoints_imp=npoints_imp+1
;	endif else if(time_buf gt end0) then goto, Continue2
;	nlines_imp=nlines_imp+1
;    endwhile

Continue2:
;
;close,lun_geo
;close,lun_imp
;free_lun,lun_geo
;free_lun,lun_imp
;
;print,'Total number of geotail data points: ',npoints_geo
;print,'Total number of imp 8 data points: ',npoints_imp
;
; if sm =1 transfer from GSM to SM (z-axis aigned with dipole axis)
;
;if(sm eq 1) then begin
;    for i=0,npoints-1 do begin 
;	epoch0=time(i)
;	yrdhms,epoch0,yr,doy,hr,min,sec             ; in milliseconds
;	yr1=long(yr+1900)
;	doy1=long(doy)
;	hr1=long(hr)
;	min1=long(min)
;	sec1=float(sec)
;	theta=0.0
;	s0=call_external('dipole_tilt.so','dipoletilt_idl',$
;            yr1,doy1,hr1,min1,sec1,theta)
;	zenith_angle=(90.-theta)*!pi/180. + atan(zloc(i),xloc(i))
;	r=sqrt(xloc(i)*xloc(i)+zloc(i)*zloc(i))
;	xloc(i)=r*cos(zenith_angle)
;	zloc(i)=r*sin(zenith_angle)
;	if(i eq 0) or (i eq npoints-1) then print,yr,doy,hr,min,sec
;    endfor
;endif

;
npanels=4       ;number of panels.
!p.multi=[0,1,2] ; 2 plots vertically aligned
;
;
!x.type=xtype
!x.minor=4
!x.style=1
!y.type=ytype
!y.style=1
!y.minor=4
a=findgen(9)*(!PI*2./8.)
radius=0.5
usersym,radius*cos(a),radius*sin(a),/fill
;
; To determine the x positions of panels.
;
xgap=0.
;xleft=100	; In device coordinate 
;xright=700
;
; To determine the y positions of panels.
;
ygap=0       ; ratio of the height.
;ybottom=100 	; In device coordinate
;ytop=700 
;
n1=201
n2=2*n1
n3=3*n1
n4=4*n1
xseg=fltarr(2)
yseg=fltarr(2)
zeros=fltarr(2)
fullwidth=fltarr(2)
xmp=fltarr(n4+22)
ymp=fltarr(n4+22)
xpos=fltarr(npanels)
ypos=fltarr(npanels)
xtitle=strarr(npanels)
ytitle=strarr(npanels)
ymin=fltarr(npanels)
ymax=fltarr(npanels)
xmin=40
xmax=-40
ymin0=40
ymax0=-40
zeros(0)=0.
zeros(1)=0.
fullwidth(0)=ymin0
fullwidth(1)=ymax0
ymin(0)=ymin0
ymax(0)=ymax0
ymin(1)=ymin0
ymax(1)=ymax0
ymin(2)=ymin0
ymin(2)=ymax0
ymax(3)=ymin0
ymax(3)=ymax0
;
;title="GEOTAIL & IMP 8 Orbits !C96/02/15 1200 - 02/16 1200 UT"
;title="GEOTAIL & IMP 8 Orbits !C95/04/25 0000 - 04/25 1200 UT"
;title="GEOTAIL & IMP 8 Orbits !C96/11/6 2200 - 11/7 0400 UT"
;title="GEOTAIL & IMP 8 Orbits !C96/11/26 1000 - 11/28 0000 UT"
;title="GEOTAIL Orbits !C96/11/26 1000 - 11/28 0000 UT"
;title="GEOTAIL Orbits !C96/12/08 0000 - 12/09 0600 UT"
;ytitle(0)="(Y!Dgse!N!E2!N+Z!Dgse!N!E2!N)!E1/2!N (Re)"
;
if(aberr eq 1) then begin
	title="Bow Shock Crossings !C(5 deg. aberrated)"
	xtitle="Xgse_aberr (Re)"
	ytitle(0)="Ygse_aberr (Re)"
endif else begin
	title="Bow Shock Crossings"
	xtitle="Xgse (Re)"
	ytitle(0)="Ygse (Re)"
endelse

width=(xright-xleft)/2
height=width
xpos(0)=xleft
ypos(0)=ybottom+height+ygap
xpos(1)=xleft+width+xgap
ypos(1)=ypos(0)
xpos(2)=xleft
ypos(2)=ybottom
xpos(3)=xpos(1)
ypos(3)=ypos(2)
;
;for i=1,npanels-1 do ypos(i)=ypos(i-1)+height
;ygap=ygap*height        ; gaps in normal coordinate.
;
; the first panel
;
; model magnetopause
;
if(npage eq 2) then goto, TheSecond
;
; get model magnetpause, earth, etc.
;
;if(tra_sheath eq 1) then get_mp,xmp,ymp,n1,standoff_bs,pdyn_index $
;else get_scaling_roelofmp,standoff,standoff_bs,xmp,ymp,n1,0

get_scaling_roelofmp,standoff,standoff_bs,xmp,ymp,n1,0

;print,'standoff_bs: ',standoff_bs
;print,'pdyn_index : ',pdyn_index
get_mp,xmp,ymp,n1,standoff_bs,pdyn_index

if(idevice eq 1) then window,0,xsize=400,ysize=800
thick=1
if(idevice ne 2) then begin
    plot,xmp(0:n1-1),ymp(0:n1-1),$
	title=title,$
	xtitle=xtitle,xrange=[xmin,xmax],charthick=charthick,$
	ytitle=ytitle(0),yrange=[ymin(0),ymax(0)],$
	position=[xleft,ymid,xright,ytop],thick=thick
;	position=[xleft,ybottom,xright,ytop],thick=2
endif else begin
    plot,xmp(0:n1-1),ymp(0:n1-1),$
	title=title,$
	xtitle=xtitle,xrange=[xmin,xmax],charthick=charthick,$
	ytitle=ytitle(0),yrange=[ymin(0),ymax(0)],$
	position=[xleft,ybottom,xright,ytop],thick=thick
endelse
oplot,xmp(n1:n2-1),ymp(n1:n2-1),linestyle=2,thick=1
;if(pdyn_index eq 2) then linestyle=2 $
;else linestyle=0
oplot,xmp(n2:n3-1),ymp(n2:n3-1),linestyle=0,thick=2,color=96
oplot,xmp(n3:n4-1),ymp(n3:n4-1),linestyle=2,thick=2,color=96
oplot,xmp(n4:n4+20),ymp(n4:n4+20),linestyle=0,thick=2
;
; Draw field lines
;
;if(draw_field eq 1) then begin
;    itemp2=0
;    if(tilt eq 0) then openr,35,'fieldlines_ts89kp2_0t.asc' $
;    else if(tilt eq 1) then openr,35,'fieldlines_ts89kp2_20t.asc' $
;    else if(tilt eq 2) then openr,35,'fieldlines_ts89kp2_7418418.asc' $
;    else openr,35,'fieldlines_ts89kp5_7426202.asc'
;    xfieldlines=fltarr(100)
;    zfieldlines=fltarr(100)
;    while(not EOF(35)) do begin
;	readf,35,temp1,itemp2
;	for i=0,itemp2-1 do begin
;	    readf,35,temp1,temp2,temp3
;	    xfieldlines(i)=temp1
;	    zfieldlines(i)=temp3
;	endfor
;	oplot,xfieldlines(0:itemp2-1),zfieldlines(0:itemp2-1),thick=1
;    endwhile
;    close,35
;endif
;
; Draw a vertical line at x=0.
;
plots,zeros,fullwidth,thick=1
;
; Draw a horizontal line at y=0.
;
plots,[xmin,xmax],[0.,0.],thick=1

;
; Draw the earth
;
;for i=0,10 do begin
;    xmp(n4+i)=-xmp(n4+i)
;endfor
;xmp(11)=0.
;ymp(11)=0.
;polyfill,xmp(n4:n4+11),ymp(n4:n4+11),col=col,thick=2
;
;-----------------------------------------------------------
;
kkk=1
if(kkk eq 0) then goto, jump1
;
; Legend if any
;
if(nlegend ne 0) then begin
    xlegend=xleft+0.075
    ylegend=ytop-0.04
    for i=1,nlegend do begin
	ylegend=ylegend-0.028
	xyouts,xlegend,ylegend,legend(i),/normal,charsiz=1.2
    endfor
endif
plots,[35,23],[31,31],linestyle=0,thick=2,color=96
plots,[35,23],[36,36],linestyle=2,thick=2,color=96
xyouts,22,32,'Pdyn=2',/data,charsiz=charsize,charthick=charthick
xyouts,22,37,'Pdyn=7',/data,charsiz=charsize,charthick=charthick
;
c_or_t=0
vratio=0.5  ;ratio of vector to the size of plot
for i=0,npoints_geo-1 do begin
    if(index_geo(i) eq 1) then begin
	if(c_or_t eq 1) or (c_or_t eq 2) then begin
	    if(conserve eq 1) then begin 
		tvec=sqrt(zvec(i)*zvec(i)+yvec(i)*yvec(i))
		tloc=sqrt(zloc_geo(i)*zloc_geo(i)+yloc_geo(i)*yloc_geo(i))
	    endif else begin
		tvec=zvec(i)
		tloc=zloc_geo(i)
	    endelse
	    factor=vratio/sqrt(xvec(i)*xvec(i)+tvec*tvec)
;
; get tangential inst. of normal 
;
	    xseg(0)=xloc_geo(i)-tvec*factor   
	    yseg(0)=tloc+xvec(i)*factor
	    xseg(1)=xloc_geo(i)+tvec*factor
	    yseg(1)=tloc-xvec(i)*factor 
	    plots,xseg,yseg,thick=4
	endif else begin
	    if(conserve eq 0) then plots,xloc_geo(i),zloc_geo(i),psym=8 $
	    else plots,xloc_geo(i),$
		sqrt(zloc_geo(i)*zloc_geo(i)+yloc_geo(i)*yloc_geo(i)),$
			psym=8,thick=2
	endelse
    endif else begin
	if(c_or_t eq 0) or (c_or_t eq 2) then begin
	    if(sctype(i) eq 1) then $
		plots,xloc_geo(i),yloc_geo(i),psym=8,thick=1,color=64 $
	    else plots,xloc_geo(i),yloc_geo(i),psym=8,thick=1,color=128
;
;	    plots,xloc_geo(i),$
;		sqrt(zloc_geo(i)*zloc_geo(i)+yloc_geo(i)*yloc_geo(i)),$
;		psym=8,thick=3,color=0
	endif
    endelse
endfor
;usersym,0.6*cos(a),0.6*sin(a)
;xyouts,xloc_geo(0)-1.0,yloc_geo(0),'02/15 12:00',$
;xyouts,xloc_geo(0)+16.0,yloc_geo(0),'04/25 00:00',$
;xyouts,xloc_geo(0)-2,yloc_geo(0),'11/06 22:00',$
;xyouts,xloc_geo(0)+5,yloc_geo(0)+3,'12/08 00:00',$
;xyouts,xloc_geo(0)+5,yloc_geo(0)+3,'11/15 12:00',$
;	charsize=charsize,charthick=charthick
;plots,xloc_geo(0),yloc_geo(0),psym=8,thick=3,color=0
;plots,xloc_geo(npoints_geo/2),yloc_geo(npoints_geo/2),psym=8,thick=3,color=0
;xyouts,xloc_geo(npoints_geo-1),yloc_geo(npoints_geo-1)-2.0,$
;	'02/16 12:00',charsize=charsize,charthick=charthick
;xyouts,xloc_geo(npoints_geo-1)+3.0,yloc_geo(npoints_geo-1)-2.0,$
;	'04/25 12:00',charsize=charsize,charthick=charthick
;xyouts,xloc_geo(npoints_geo-1)-2.0,yloc_geo(npoints_geo-1),$
;	'11/07 04:00',charsize=charsize,charthick=charthick
;xyouts,xloc_geo(npoints_geo-1)+10,yloc_geo(npoints_geo-1)+3,$
;	'11/16 09:00',charsize=charsize,charthick=charthick
;plots,xloc_geo(npoints_geo-1),yloc_geo(npoints_geo-1),psym=8,thick=3,color=0
;xyouts,20,-20,'GEOTAIL',charsize=charsize*1.5,charthick=charthick
;
;
;usersym,radius*cos(a),radius*sin(a)
;for i=0,npoints_imp-1 do begin
;;    plots,xloc_imp(i),sqrt(zloc_imp(i)*zloc_imp(i)+yloc_imp(i)*yloc_imp(i)),$
;    plots,xloc_imp(i),yloc_imp(i),$
;	 psym=8,thick=3,color=0
;endfor
;usersym,0.6*cos(a),0.6*sin(a)
;xyouts,xloc_imp(0)-1.0,yloc_imp(0),'02/15 12:00',$
;xyouts,xloc_imp(0)+15.0,yloc_imp(0),'04/25 00:00',$
;xyouts,xloc_imp(0)-2,yloc_imp(0),'11/07 04:00',$
;	charsize=charsize,charthick=charthick
;plots,xloc_imp(0),yloc_imp(0),psym=8,thick=3,color=0
;xyouts,xloc_imp(npoints_imp-1)-1.0,yloc_imp(npoints_imp-1),$
;	'02/16 12:00',charsize=charsize,charthick=charthick
;plots,xloc_imp(npoints_imp-1),yloc_imp(npoints_imp-1),psym=8,thick=3,color=0
;xyouts,-7,36,'IMP 8',charsize=charsize*1.5,charthick=charthick
;
;
; historgram of solar wind pressure
;
xtitle="No. events"
ytitle="Pressure (nPa)"
title="Pdyn!C!95!X (nPa)  "
h=histogram(pdyn,binsize=0.5,min=-0.25,max=10.249)
xupper=fix(max(h)*1.2)
h_elem=n_elements(h)
plot,[0,0],[0,0],yrange=[0,10],xrang=[0,xupper],$
        xtitle=xtitle,ytitle=" ",title=title,thick=2,$
	xticks=1,yticks=5,xminor=1,yminor=2,yticklen=0.1,$
        position=[0.68,ybottom,0.78,ytop]
x_buf=0
y_buf=0
for i=0,h_elem-1 do begin
    yhist=i*0.5
    x_buf2=h(i)
    y_buf2=yhist-0.25
    oplot,[x_buf,x_buf2],[y_buf2,y_buf2],thick=2,color=160 
    oplot,[x_buf2,x_buf2],[y_buf2,y_buf2+0.5],thick=2,color=160
    x_buf=x_buf2
    y_buf=y_buf2
endfor
x_buf2=0
oplot,[x_buf,x_buf2],[y_buf2,y_buf2],thick=2
;
jump1:
;
if(npage eq 1) then goto, TheEnd

;if(idevice eq 1) then begin
;    print,'*** click the left bottom on the mouse for the next page ***'
;    cursor,xcursor,ycursor,/down
;endif
;
; The second panel
;
TheSecond:

if(idevice eq 2) then begin
    bytemap=tvrd(0,0,512,480)
    ;help,bytemap
    ;print,'creating GIF output: ',outfile1(0)
    write_gif,outfile1(0),bytemap
    erase,255
endif

xmin=-40
xmax=40
ymin0=-40
ymax0=40
fullwidth(0)=ymin0+1.
fullwidth(1)=ymax0
ymin(0)=ymin0
ymax(0)=ymax0
ytitle(0)="Zgse (Re)"
if(aberr eq 1) then begin
	xtitle="Ygse_aberr (Re)"
endif else begin
	xtitle="Ygse (Re)"
endelse
;
; model magnetopause and earth
;
get_circle,xmp,ymp,n1,pdyn_index	; get model magnetopause, earth, etc.
;
;
if(idevice ne 2) then begin
plot,xmp(0:n1-1),ymp(0:n1-1),$
	title=' ',$
	xtitle=xtitle,xrange=[xmin,xmax],charthick=charthick,$
	ytitle=ytitle(0),yrange=[ymin(0),ymax(0)],$
	position=[xleft,ybottom,xright,0.45],thick=thick
;	position=[xleft,ybottom,xright,ytop],thick=thick
endif else begin
    plot,xmp(0:n1-1),ymp(0:n1-1),$
	title=' ',$
	xtitle=xtitle,xrange=[xmin,xmax],charthick=charthick,$
	ytitle=ytitle(0),yrange=[ymin(0),ymax(0)],$
	position=[xleft,ybottom,xright,ytop],thick=thick
endelse
;if(pdyn_index eq 2) then oplot,xmp(n1:n2-1),ymp(n1:n2-1),linestyle=2,thick=2
if(pdyn_index eq 2) then linestyle=2 $
else linestyle=0
oplot,xmp(n2:n3-1),ymp(n2:n3-1),linestyle=0,thick=2,color=96
oplot,xmp(n3:n3+20),ymp(n3:n3+20),linestyle=0,thick=2
plots,[xmin,xmax],zeros,thick=1
plots,zeros,[ymin0,ymax0],thick=1
;
vratio=0.5  ;ratio of vector to the size of plot
;-----------------------------------------------------------
kkk=1
if(kkk eq 0) then goto, TheEnd 
usersym,radius*cos(a),radius*sin(a),/fill
for i=0,npoints_geo-1 do begin
    if(index_geo(i) eq 1) then begin
;    if(index(i) ne 9) then begin
	if(c_or_t eq 1) or (c_or_t eq 2) then begin
	    factor=vratio/sqrt(yvec(i)*yvec(i)+zvec(i)*zvec(i))
	    xseg(0)=yloc(i)-zvec(i)*factor  ; get tangential inst. of normal
	    yseg(0)=zloc(i)+yvec(i)*factor
	    xseg(1)=yloc(i)+zvec(i)*factor
	    yseg(1)=zloc(i)-yvec(i)*factor
	    plots,xseg,yseg,thick=4
	endif else begin
	    plots,yloc(i),zloc(i),psym=8,thick=1,color=64
	endelse
    endif else begin
	if(c_or_t eq 0) or (c_or_t eq 2) then begin
	    if(sctype(i) eq 1) then $
		plots,yloc_geo(i),zloc_geo(i),psym=8,thick=1,color=64 $
	    else plots,yloc_geo(i),zloc_geo(i),psym=8,thick=1,color=128
	endif
    endelse
endfor
;usersym,0.6*cos(a),0.6*sin(a)
;xyouts,yloc_geo(0)+1.0,zloc_geo(0)+0.5,'02/15 12:00',$
;xyouts,yloc_geo(0)+1.0,zloc_geo(0)-1.0,'04/25 00:00',$
;xyouts,yloc_geo(0)-10,zloc_geo(0)+1.0,'11/06 22:00',$
;xyouts,yloc_geo(0)-10,zloc_geo(0)-2,'12/08 00:00',$
;	charsize=charsize,charthick=charthick
;plots,yloc_geo(0),zloc_geo(0),psym=8,thick=3,color=0
;plots,yloc_geo(npoints_geo/2),zloc_geo(npoints_geo/2),psym=8,thick=3,color=0
;xyouts,yloc_geo(npoints_geo-1),zloc_geo(npoints_geo-1)+1.0,$
;	'02/16 12:00',charsize=charsize,charthick=charthick
;xyouts,yloc_geo(npoints_geo-1)-8.0,zloc_geo(npoints_geo-1)+1.,$
;	'04/25 12:00',charsize=charsize,charthick=charthick
;xyouts,yloc_geo(npoints_geo-1)+1,zloc_geo(npoints_geo-1)+1.,$
;	'11/07 04:00',charsize=charsize,charthick=charthick
;xyouts,yloc_geo(npoints_geo-1)+1,zloc_geo(npoints_geo-1),$
;	'12/09 06:00',charsize=charsize,charthick=charthick
;plots,yloc_geo(npoints_geo-1),zloc_geo(npoints_geo-1),psym=8,thick=3,color=0
;xyouts,-22,5,'GEOTAIL',charsize=charsize*1.5,charthick=charthick
;
;usersym,radius*cos(a),radius*sin(a)
;for i=0,npoints_imp-1 do begin
;	plots,yloc_imp(i),zloc_imp(i),psym=8,thick=3,color=0
;endfor
;usersym,0.6*cos(a),0.6*sin(a)
;xyouts,yloc_imp(0)+0.5,zloc_imp(0)-1.0,'02/15 12:00',$
;xyouts,yloc_imp(0)-4.,zloc_imp(0)+1.0,'04/25 00:00',$
;xyouts,yloc_imp(0)-4.,zloc_imp(0)+1.0,'11/06 22:00',$
;	charsize=charsize,charthick=charthick
;plots,yloc_imp(0),zloc_imp(0),psym=8,thick=3,color=0
;xyouts,yloc_imp(npoints_imp-1)+1.0,zloc_imp(npoints_imp-1),$
;	'02/16 12:00',charsize=charsize,charthick=charthick
;xyouts,yloc_imp(npoints_imp-1)+1.0,zloc_imp(npoints_imp-1),$
;	'04/25 12:00',charsize=charsize,charthick=charthick
;plots,yloc_imp(npoints_imp-1),zloc_imp(npoints_imp-1),psym=8,thick=3,color=0
;xyouts,30,5,'IMP 8',charsize=charsize*1.5,charthick=charthick

TheEnd:

if(idevice eq 0) then begin
    device,/close
    print,'The graphic output file will be ',filename(0)
    set_plot,'X'
endif else begin
    bytemap=tvrd(0,0,512,480)
    ;print,'creating GIF output: ',outfile2(0)
    write_gif,outfile2(0),bytemap
    ;print,'Done.'
endelse
close,lun_feed
end