; orbit_geotail.pro
;
; Plot Geotail IMP 8 trajectories in Y-Z, and X-Z planes with model
; magnetospheric field lines, bowshocks, magnetopauses, etc.
;
;
; 03/19/97 SHC adapted from orbit_hk.pro also by SHC
;
;------ subroutine get_mp
;
pro get_mp,xmp,ymp,n,standoff_bs,pdyn_index
    theta1=300.d0
    theta2=300.d0
    DegtoRad=!PI/180.d0
    dtheta1=theta1/(n-1)
    dtheta2=theta2/(n-1)
    theta1=-150.d0
    theta2=-150.d0
    r_bs=0.
    x1=0.
    y1=0.
    z1=0.
    rnorm=0.
    bs_factor=1.0
    standoff_bs2=0.0
    if(pdyn_index eq 1) then bs_factor=1.0 $
    else bs_factor=(1.8/8.0)^(1.0/6.0)

    ;s1=call_external('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)=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+2*n

	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('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
;
;	xmp(k)=xmp(k)-13.7*(1.-bs_factor); for 1996/2/15
;

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

	k=i+3*n

;	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

    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
;
;- 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
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
;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,' 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 - Z-buf)? ',idevice 
; 
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
col=0	;color of the nightside earth
charsize=1.

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='idl.ps',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
;	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'
	charsize=0.5
	!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
	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'
	charsize=2.0
        !p.charsize=charsize
        !p.charthick=2.0
        !x.thick=2.
        !y.thick=2.
	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
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('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

s2=call_external('mp_roelof.so','magnetopause_roelof_idl',$
    0.,pdyn0,9999.,9999.,9999.,standoff,dawndusk,z1,r_mp,0L) 
;
; read geotail
;
print,'Reading geotail data...'
read,npoints_geo
for i=1,npoints_geo do begin
	read,temp1,temp2,temp3
	xloc_geo(i)=temp1
	yloc_geo(i)=temp2
	zloc_geo(i)=temp3
endfor
;
;    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=5
!x.style=1
!y.type=ytype
!y.style=1
!y.minor=5
a=findgen(17)*(!PI*2./16.)
radius=0.2
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=-20
ymin0=10
ymax0=-50
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
;
xtitle="Xgse (Re)"
;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"
title="GEOTAIL Orbits !C96/11/15 1200 - 11/16 0900 UT"
ytitle(0)="Ygse (Re)"
;ytitle(0)="(Y!Dgse!N!E2!N+Z!Dgse!N!E2!N)!E1/2!N (Re)"
;
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.
;

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

;if(idevice ne 0) then window,0,xsize=300,ysize=600

;--------------------------------------------------------------
; Create x-y plot Objects
;
View=OBJ_NEW('IDLgrView',VIEWPLANE_RECT=[xleft,ybottom,xright,ytop])
Plotxy=OBJ_NEW('IDLgrPlot',xmp(0:n1-1),ymp(0:n1-1))
Xtitlexy=OBJ_NEW('IDLgrText',xtitle)
Ytitlexy=OBJ_NEW('IDLgrText',ytitle(0))
Xaxisxy=OBJ_NEW('IDLgrAxis',0,range=[xmin,xmax],title=Xtitlexy)
Yaxisxy=OBJ_NEW('IDLgrAxis',1,range=[ymin(0),ymax(0)],title=Ytitlexy)
Usym=OBJ_NEW('IDLgrSymbol',2)

;oplot,xmp(n1:n2-1),ymp(n1:n2-1),linestyle=2,thick=2
if(pdyn_index eq 2) then linestyle=2 $
else linestyle=0
Oplotxy1=OBJ_NEW('IDLgrPolyline',xmp(n2:n3-1),ymp(n2:n3-1),linestyle=0,thick=2)
Oplotxy2=OBJ_NEW('IDLgrPolyline',xmp(n3:n4-1),ymp(n3:n4-1),linestyle=5,thick=2)
Oplotxy3=OBJ_NEW('IDLgrPolyline',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.
;
Oplotxy5=OBJ_NEW('IDLgrPolyline',zeros,fullwidth,linestyle=0,thick=1)
;
; Draw a horizontal line at y=0.
;
Oplotxy6=OBJ_NEW('IDLgrPolyline',[xmin,xmax],zeros,linestyle=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,[20.,35.],[39.,39.],linestyle=0,thick=2
;plots,[20.,35.],[42.,42.],linestyle=5,thick=2
;plots,[20.,35.],[22.,22.],linestyle=0,thick=2
;plots,[20.,35.],[25.,25.],linestyle=5,thick=2
;plots,[20.,35.],[45.,45.],linestyle=2,thick=2
;;plots,[35.,20.],[-45.,-45.],linestyle=0,thick=2
;;plots,[35.,20.],[-42.,-42.],linestyle=5,thick=2
;xyouts,18.,39.,'Farris bs (p=2)',/data,charsiz=charsize,charthick=charthick
;xyouts,18.,22.,'Farris bs (p=2)',/data,charsiz=charsize,charthick=charthick
;;xyouts,18.,-45.,'Farris bs (p=2)',/data,charsiz=charsize,charthick=charthick
;xyouts,18.,42.,'Peredo bs (ma=15,p=8)',/data,charsiz=charsize,charthick=charthick
;xyouts,18.,42.,'Peredo bs (ma=10,p=2)',/data,charsiz=charsize,charthick=charthick
;xyouts,18.,25.,'Peredo bs (ma=10,p=2)',/data,charsiz=charsize,charthick=charthick
;;xyouts,18.,-42.,'Peredo bs (ma=10,p=2)',/data,charsiz=charsize,charthick=charthick
;xyouts,18.,45.,'Peredo bs (ma=10,p=2)',/data,charsiz=charsize,charthick=charthick
;
c_or_t=0
vratio=0.5  ;ratio of vector to the size of plot
;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(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
;;	    plots,xloc_geo(i),yloc_geo(i),psym=8,thick=3,color=0
;	
;;	    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

Oplotxy7=OBJ_NEW('IDLgrPolyline',xloc_geo(0:npoints_geo-1),$
	yloc_geo(0:npoints_geo-1),symbol=Usym,color=0,thick=3)

;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
;
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
;
;---------------------------------------------
; Create x-y View
;
Modelxy=OBJ_NEW('IDLgrModel')
Modelxy->ADD,plotxy
Modelxy->ADD,Oplotxy1
Modelxy->ADD,Oplotxy2
Modelxy->ADD,Oplotxy3
Modelxy->ADD,Oplotxy5
Modelxy->ADD,Oplotxy6
Modelxy->ADD,Oplotxy7
View->ADD, Modelxy
;
;------------------------------------------------
;
; The second panel
;
TheSecond:

xmin=-40
xmax=0
ymin0=-20
ymax0=20
fullwidth(0)=ymin0+1.
fullwidth(1)=ymax0
ymin(0)=ymin0
ymax(0)=ymax0
xtitle="Ygse (Re)"
ytitle(0)="Zgse (Re)"
;
; model magnetopause and earth
;
get_circle,xmp,ymp,n1,pdyn_index		; get model magnetpause, earth, etc.
;
;
;--------------------------------------------------------------
; Create y-z plot Objects
;
;Viewyz=OBJ_NEW('IDLgrView',VIEWPLANE_RECT=[xleft,ybottom,xright,0.45]);
Plotyz=OBJ_NEW('IDLgrPlot',xmp(0:n1-1),ymp(0:n1-1))
Xtitleyz=OBJ_NEW('IDLgrText',xtitle)
Ytitleyz=OBJ_NEW('IDLgrText',xtitle)
Xaxisyz=OBJ_NEW('IDLgrAxis',0,range=[xmin,xmax],title=Xtitleyz)
Yaxisyz=OBJ_NEW('IDLgrAxis',1,range=[ymin(0),ymax(0)],title=Ytitleyz)

;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
Oplotyz1=OBJ_NEW('IDLgrPolyline',xmp(n2:n3-1),ymp(n2:n3-1),linestyle=0,thick=2)
Oplotyz2=OBJ_NEW('IDLgrPolyline',xmp(n3:n3+20),ymp(n3:n3+20),linestyle=0,thick=2)

Oplotyz3=OBJ_NEW('IDLgrPolyline',[xmin,xmax],zeros,linestyle=0,thick=1)
Oplotyz4=OBJ_NEW('IDLgrPolyline',zeros,[ymin0,ymax0],linestyle=0,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=2
;	endelse
;    endif else begin
;	if(c_or_t eq 0) or (c_or_t eq 2) then begin
;	    plots,yloc_geo(i),zloc_geo(i),psym=8,thick=3,color=0
;	endif
;    endelse
;endfor
Oplotyz7=OBJ_NEW('IDLgrPolyline',yloc_geo(0:npoints_geo-1),$
	zloc_geo(0:npoints_geo-1),symbol=Usym,color=0,thick=3)

;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
;-------------------------------------------------------
;
; Create y-z View
;
Modelyz=OBJ_NEW('IDLgrModel')
Modelyz->ADD,plotyz
Modelyz->ADD,Oplotyz1
Modelyz->ADD,Oplotyz2
Modelyz->ADD,Oplotyz3
Modelyz->ADD,Oplotyz4
Modelyz->ADD,Oplotyz7
View->ADD, Modelyz

Scene=OBJ_NEW('IDLgrScene')
if(idevice eq 1) then Window=OBJ_NEW('IDLgrWindow') $
else Window=OBJ_NEW('IDLgrWindow',/pixmap)
set_view,View,window
;set_view,viewyz,Window
;Window->DRAW,View
;Window->DRAW,Viewyz
;Scene->Add,View
;Scene->Add,Viewyz
Window->DRAW,Scene
;-------------------------------------------------------
;
; Draw to a gif file
;
;set_plot,'z'
if(idevice ne 0) then begin
    myimage=Window->read()
    myimage->GetProperty,DATA=image
    image0=reform(image[0,*,*])
;    bytemap=tvrd()
    print,'creating GIF output: idl.gif...'
    write_gif,'idl0.gif',image0
endif

;-------------------------------------------------------
;
; Draw to a postcript file
;
;set_plot,'ps'
;device,filename='idl.ps',xoffset=2,xsize=18.5,yoffset=2,$
;	ysize=24,/HELVETICA,/BOLD,FONT_SIZE=12
;Window->DRAW,Scene,/create_instance
;print,'creating ps output: idl.ps...'
;device,/close

print,'Done'

TheEnd:

close,/all
;--------------------------------------------------------------
; desctroy objects
;
OBJ_DESTROY,view
OBJ_DESTROY,Scene
OBJ_DESTROY,Window

end