PRO plotorb_ps, dum ; plotorb plots the orbital ".in" file 
; ULTRIX IDL version, written by Edmund G. Dombrowski, 12/89, 
; modified by W. Hartkopf 1/90, 7/90. 
; Based on VMS fortran programs written by W. Hartkopf and adapted to IDL by D. Gies.
; menu versions by WIH 8/99, 9/99, 3/00, 6/00
; 2005/08/24 - add black+white option for journal figures, etc.
; 2007/10/21 - version for remote operation (ps only, no menu)

plotorb_ps_subs										; compile subroutines
!p.background=255

common orbpar,pm,p,a,xi,xnode,t0,ecc,omega,c_black,c_red,c_green, $
	c_dkblue,c_rayin,c_rayout,c_fuchsia,c_white

set_plot,'ps'
device, /inches, xoffset=0.75, yoffset=2.0      				; offsets for bottom left corner
device, /inches, xsize=7.0, ysize=7.0           				; x and y size for 8 1/2 x 11
!p.font=0                                       				; font selector to hardware option

loadct,12                                 					; load color table
tvlct,r,g,b,/get                                				; define colors for ct=12
c_black=0		&  c_red=200		&  c_green=30
c_dkblue=108	&  c_wedge=230		&  c_fuchsia=135	
c_white=255		&  c_rayin=215  		&  c_rayout=230	

MAIN:

rx=[0.]		&  ry=[0.]			
xl=[0.]		&  yl=[0.]

close,10		&  infile=' '		&  print,' '
read,' Enter plot file name: ',infile							; READ PLOT FILE
openr,10,infile
ref=' '		&  refa=' '

orbgrid='N'
read,' Is this an orbgrid plot file? (Y or N): ',orbgrid				; SOURCE OF FILE?
if (orbgrid eq 'n') then orbgrid='N'
if (orbgrid eq 'y') then orbgrid='Y'

psout=' '		
read,' Enter output ps file name: ',psout							; POSTSCRIPT FILENAME
device,/close_file
device, filename=psout  
device,/color

puborb='N'		&  pubfile=' '
read,' Do you wish to overplot a published orbit? (Y or N): ',puborb		; PLOT PUBLISHED ORBIT?
if (puborb eq 'y') then puborb='Y'
if (puborb eq 'Y') then read,' Enter name for published orbit: ',pubfile
px1=[0.]		&  py1=[0.]
if (pubfile ne ' ') then orbpub,pubfile,px1,py1

qfill='N'		&  aper=999.
read,' Do you wish to plot a Rayleigh limit circle? (Y or N): ',qfill		; PLOT RAYLEIGH LIMIT?
if (qfill eq 'y') then qfill='Y'
if (qfill eq 'Y') then read,'    Enter aperture in meters: ',aper

wedge='N'		&  wstart=2000.		&  wend=2020.
read,' Do you wish to plot a wedge? (Y or N): ',wedge					; PLOT WEDGE?
if (wedge eq 'y') then wedge='Y'
if (wedge eq 'Y') then read,'    Enter starting date: ',wstart
if (wedge eq 'Y') then read,'    Enter ending date: ',wend
print,' '

orbget,star,pos,rx,ry,x,y,nsymbol,ssize,npt,xl,yl,xcalc,ycalc,wt, $		; read in data
	ref,refa,wstart,wend,xwedge,ywedge,orbgrid

;...................................................................

PLOTIT:

xmax=max([rx,x,px1])			&  xmin=min([rx,x,px1])				; find maximum data range 
ymax=max([ry,y,py1])			&  ymin=min([ry,y,py1])
ylen=ymax-ymin				&  xlen=xmax-xmin
yxratio=ylen/xlen											; calculate y to x axis ratio 
													; to determine an even range

; Check the ratio for cases >=1 and <1 and add an appropriate width to
; xmin, xmax, ymin, and ymax to center the plot in the window. This width
; (= xymargin) is set smaller if the orbit is largely horizontal or
; vertical, since these plots are more likely to be crowded. Conversely,
; it's set somewhat larger for near circular orbits to avoid the NE 
; arrows and the semi-major axis legend.

xymargin=0.2
if((yxratio lt 0.5) or  (yxratio gt 2.0)) then xymargin=0.1
if((yxratio gt 0.9) and (yxratio lt 1.1)) then xymargin=0.25

if(yxratio ge 1.0000) then begin
  	yfactor=xymargin*ylen
  	xfactor=yfactor+(ylen-xlen)/2.0
  	endif

if(yxratio lt 1.0000) then begin
  	xfactor=xymargin*xlen
  	yfactor=xfactor+(xlen-ylen)/2.0
  	endif

newxmin=xmin-xfactor			&  newxmax=xmax+xfactor
newymin=ymin-yfactor 			&  newymax=ymax+yfactor

occxmax=xmax+xfactor*0.3		&  occxmin=xmin-xfactor*0.3
occymax=ymax+yfactor*0.3		&  occymin=ymin-yfactor*0.3

; set up other plotting variables.

!x.minor=-1					&  !y.minor=-1    			; suppress minor ticks.

STARTPLOT:

!x.margin=[0,0]				&  !y.margin=[0,0]			; set up x and y range
!x.style=1					&  !y.style=1
!x.range=[newxmin,newxmax]		&  !y.range=[newymin,newymax]

plus=0.025*(newxmax-newxmin)								; plot origin as large plus sign,
												; based on size of new axes.

plot, [-plus,plus], [0,0], thick=2,charsize=1.25,xcharsize=0.75,ycharsize=0.75,ticklen=0.01, $
xtickname=REPLICATE(' ',10),ytickname=REPLICATE(' ',10),/isotropic,/nodata
axis,0.0,yaxis=1,ticklen=-0.01,/normal
axis,0.0,xaxis=1,ticklen=-0.01,/normal
	
if (qfill eq 'Y') then begin								; shade in Rayleigh limit region 
	ray1=0.1258/aper
	xlim=fltarr(361)	&  ylim=fltarr(361)
	A=fltarr(361)
	for i=0,360 do A(i)=i
	A=A*(!pi *2/360.)
	xlim=cos(A)*ray1 > newxmin < newxmax
	ylim=sin(A)*ray1 > newymin < newymax
	polyfill,xlim,ylim,color=c_rayout
	endif

diff1=(newxmax-newxmin)			&  diff2=(newymax-newymin)
xt1=newxmin+0.08*diff1			&  xt2=newxmax-0.08*diff1
yt1=newymin+0.930*diff2			&  yt2=newymin+0.975*diff2

oplot, [-plus,plus], [0,0], thick=2
oplot, [0,0], [-plus,plus], thick=2
polyfill,[xt1,xt2,xt2,xt1,xt1],[yt1,yt1,yt2,yt2,yt1],color=255
xyouts,0.06,0.94,'WDS '+pos, charsize=1.0,alignment=0.0,/normal
xyouts,0.55,0.94,star,       charsize=1.0,alignment=0.5,/normal
if (refa eq ' ') then xyouts,0.94,0.94,'('+ref+')',charsize=1.0,alignment=1.0,color=c_black,/normal
if (refa ne ' ') then xyouts,0.94,0.94,'('+ref+refa+')',charsize=1.0,alignment=1.0,color=c_black,/normal

;plot the north-east bars in lower right corner.

nesetx1=newxmax-(0.15*diff1)		&  nesetx2=newxmax-(0.10*diff1)
nesety2=newymin+(0.15*diff2)		&  nesety1=newymin+(0.10*diff2)

nletx=nesetx1-(0.006*diff1)		&  nlety=nesety1-(0.025*diff2)
eletx=nesetx2+(0.0125*diff1)		&  elety=nesety2-(0.006*diff2)

oplot, [nesetx1,nesetx2], [nesety2,nesety2],color=c_black
oplot, [nesetx1,nesetx1], [nesety2,nesety1],color=c_black
xyouts, nletx, nlety, 'N', size=0.9,color=c_black
xyouts, eletx, elety, 'E', size=0.9,color=c_black

;  plot a bold frame around axis

oplot, [newxmin*0.999,newxmin*0.999], [newymin,newymax],color=c_black, thick=5  
oplot, [newxmin,newxmax], [newymin*0.999,newymin*0.999],color=c_black, thick=5
oplot, [newxmax*0.999,newxmax*0.999], [newymin,newymax],color=c_black, thick=5
oplot, [newxmin,newxmax], [newymax*0.999,newymax*0.999],color=c_black, thick=5


bx=fltarr(17)				&  by=fltarr(17)				; plot CW or CCW arc
A=fltarr(17)				&  for i=0,16 do A(i)=i
A=A*(!pi *1/32.) + !pi*1.5
xarc=newxmax-(0.124*diff1)		&  yarc=newymin+(0.125*diff2)
bx=cos(A)*0.043*diff1			&  by=sin(A)*0.043*diff2

oplot,bx+xarc,by+yarc,thick=1.5
if (xi le !pi/2.) then oplot, $
	[xarc+0.038*diff1,xarc+0.043*diff1,xarc+0.048*diff1,xarc+0.043*diff1,xarc+0.043*diff1], $
	[yarc,yarc+0.010*diff2,yarc,yarc+0.010*diff2,yarc],thick=1.5,color=c_black
if (xi le !pi/2.) then oplot,[xarc,xarc-0.010*diff1],[yarc-0.043*diff1,yarc-0.043*diff1],  $
	thick=1.5,color=c_black
if (xi gt !pi/2.) then oplot,  $
	[xarc,xarc-0.010*diff1,xarc,xarc-0.010*diff1,xarc], $
	[yarc-0.048*diff2,yarc-0.043*diff1,yarc-0.038*diff2,yarc-0.043*diff1,yarc-0.043*diff1],  $
	thick=1.5,color=c_black
if (xi gt !pi/2.) then oplot,[xarc+0.043*diff1,xarc+0.043*diff1],[yarc,yarc+0.010*diff2],  $
	thick=1.5,color=c_black

;................................................................

; new plotting routine.

sym=intarr(10)
sym(1)=6      ;open circle   	(eyepiece interferometry - plotted in routine plotcirc)
sym(2)=3      ;filled circle 	(speckle, etc - plotted in routine plotcirc)
sym(3)=1      ;plus sign     	(micrometry)
sym(4)=2      ;asterisk      	(photography)
sym(5)=5      ;triangle      	(occultation)
sym(6)=4      ;filled "H"    	(Hipparcos - plotted in routine plothip)
sym(7)=7      ;filled star   	(USNO - plotted in routine plotstar)
sym(8)=8      ;filled "T"    	(Tycho - plotted in routine plottyc)
sym(9)=9      ;filled square  (multi-aper. int. - plotted in routine plotsqr)

if (wedge eq 'Y') then begin    						; shade in wedge if desired
	polyfill,xwedge,ywedge,color=c_rayin
      oplot, [-plus,plus], [0,0], thick=2,color=c_black
      oplot, [0,0], [-plus,plus], thick=2,color=c_black
      endif

xx=fltarr(2)			&  yy=fltarr(2)				;plot O-C lines 
for i=0,npt do begin
	xx(0)=x(i)			&  xx(1)=xcalc(i)
	yy(0)=y(i)			&  yy(1)=ycalc(i)
	if ((wt(i) ne 0.) and (nsymbol(i) ne 5)) then begin
		if (nsymbol(i) eq 1) then oplot,xx,yy,linestyle=0,color=c_dkblue
		if (nsymbol(i) eq 2) then oplot,xx,yy,linestyle=0,color=c_dkblue
		if (nsymbol(i) eq 3) then oplot,xx,yy,linestyle=0,color=c_green
		if (nsymbol(i) eq 4) then oplot,xx,yy,linestyle=0,color=c_fuchsia
		if (nsymbol(i) eq 6) then oplot,xx,yy,linestyle=0,color=c_red
		if (nsymbol(i) eq 7) then oplot,xx,yy,linestyle=0,color=c_dkblue
		if (nsymbol(i) eq 8) then oplot,xx,yy,linestyle=0,color=c_red
		if (nsymbol(i) eq 9) then oplot,xx,yy,linestyle=0,color=c_dkblue
		endif
	if ((wt(i) eq 0.) and (nsymbol(i) ne 5)) then begin
		if (nsymbol(i) eq 1) then oplot,xx,yy,linestyle=1,color=c_dkblue
		if (nsymbol(i) eq 2) then oplot,xx,yy,linestyle=1,color=c_dkblue
		if (nsymbol(i) eq 3) then oplot,xx,yy,linestyle=1,color=c_green
		if (nsymbol(i) eq 4) then oplot,xx,yy,linestyle=1,color=c_fuchsia
		if (nsymbol(i) eq 6) then oplot,xx,yy,linestyle=1,color=c_red
		if (nsymbol(i) eq 7) then oplot,xx,yy,linestyle=1,color=c_dkblue
		if (nsymbol(i) eq 8) then oplot,xx,yy,linestyle=1,color=c_red
		if (nsymbol(i) eq 9) then oplot,xx,yy,linestyle=1,color=c_dkblue
		endif
	endfor

for i=0,npt do begin   								;plot observations 
	temp=nsymbol(i)    
	if (temp eq 3) then oplot,[x(i)],[y(i)],psym=sym(temp), $
	symsize=ssize(i),color=c_green
	if (temp eq 4) then oplot,[x(i)],[y(i)],psym=sym(temp), $
	symsize=ssize(i),color=c_fuchsia
	endfor

sssize=ssize*0.006*(newxmax-newxmin)  					; scale symbol size for plotcirc, 
											; plothip, plottyc, plotstar

plotcirc, x, y, nsymbol, sssize						; plot open and filled circles
plothip,  x, y, nsymbol, sssize						; plot filled "H"
plottyc,  x, y, nsymbol, sssize						; plot filled "T"
plotstar, x, y, nsymbol, sssize						; plot stars
plotsqr,  x, y, nsymbol, sssize						; plot filled squares

oplot,xl,yl,linestyle=3,thick=2,color=c_black				; plot line of nodes 

xx=fltarr(2)				&  yy=fltarr(2)			;plot occultation vectors
for i=0,npt do begin
	xx(0)=x(i)				&  xx(1)=xcalc(i)
	yy(0)=y(i)				&  yy(1)=ycalc(i)

	if (nsymbol(i) eq 5) then begin
		aaa1=y(i)/x(i)		&  bbb1=ycalc(i)-aaa1*xcalc(i)
		aaa2=-x(i)/y(i)		&  bbb2=y(i)-aaa2*x(i)
		xx(1)=x(i)			&  xx(0)=(bbb2-bbb1)/(aaa1-aaa2)
		yy(1)=y(i)			&  yy(0)=aaa1*xx(0)+bbb1
		oplot,xx,yy,linestyle=0,color=85
		xx(1)=xcalc(i)		&  yy(1)=ycalc(i)
		if (wt(i) ne 0) then oplot,xx,yy,linestyle=0,color=85
		if (wt(i) eq 0) then oplot,xx,yy,linestyle=1,color=85
		endif
	endfor

											; plot published orbit
if (pubfile ne ' ') then begin			
	pxx=fltarr(2)			&  pyy=fltarr(2)
	for i=0,179 do begin
		pxx(0)=px1(i+i)		&  pxx(1)=px1(i+i+1)
		pyy(0)=py1(i+i)		&  pyy(1)=py1(i+i+1)
		oplot,pxx,pyy,color=c_black
		endfor
	endif
	
	oplot, rx,ry,thick=5,color=c_black   				; plot calculated orbit over data

device,/close_file
return
end