; PLOTORB_CAT_SUBS ; Library of subroutines for use with "plotorb_ps.pro" ; Includes the following subroutines: ; orbget orbit orbpub ; plotcirc plotsqr plothip ; plottyc plotstar resid ; plotorb_ps_subs ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO orbget,star,pos,rx,ry,x,y,nsymbol,ssize,npt,xl,yl,xcalc,ycalc,wt, $ ref,refa,wstart,wend,xwedge,ywedge,orbgrid ; Read plotorb file ; ULTRIX IDL Version - read a Hartkopf plotorb input file ; Unix conversion by Doug Gies, Bill Hartkopf and Ed Dombrowski 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 p=0. & a=0. & xi=0. xnode=0. & t0=0. & ecc=0. omega=0. & ref=' ' & refa=' ' ; read header line and orbital elements star=' ' decs=' ' & pos=' ' ref=' ' & refa=' ' & pcode=' ' acode=' ' & tcode=' ' if (orbgrid eq 'N') then begin readf,10,pos,star,rah,ram,decs,decd,decm,ipm, $ format='(a10,a13,t1,f2.0,i3,a1,2f2.0,t81,i4)' pm=float(ipm)/1000. readf,10,p,pcode,a,acode,xi,xnode,t0,tcode,ecc,omega,ref,refa, $ format='(f13.6,a1,f10.5,a1,2f9.4,f13.6,a1,f9.6,f9.4,2x,a7,a1)' if (pcode eq 'c') then p=p*100. if (pcode eq 'd') then p=p/365.2421987 if (pcode eq 'h') then p=p/365.2421987/24. if (pcode eq 'm') then p=p/365.2421987/24./60. if (acode eq 'm') then a=a/1000. if (tcode eq 'd') then t0=(t0-33282.423)/365.2421987+1950. endif if (orbgrid eq 'Y') then begin readf,10,pos,star,rah,ram,decs,decd,decm, $ format='(t56,a10,t17,a13,t56,f2.0,i3,a1,2f2.0)' readf,10,ref,format='(a1)' readf,10,p,a,xi,xnode,t0,ecc,omega,pm endif ram=float(ram) ram=ram/10.0 print,star,ref,refa,format='(a13,2x,a7,a1)' orbit,rx,ry ; run orbit procedure ra=(15.*rah+ram/4.)/!radeg ; ra, dec to radians dec=(decd+decm/60.)/!radeg if (decs eq '-') then dec=-dec corr=0. if (pm ne 0.0) then corr=(0.0056*sin(ra)/cos(dec)) + 0.00417*pm*sin(dec) corr=corr/!radeg ; read observations npt=0 & x=fltarr(2500) & y=fltarr(2500) wt=fltarr(2500) & xcalc=fltarr(2500) & ycalc=fltarr(2500) nsymbol=intarr(2500) & ssize=fltarr(2500) & equ=2000.0 code1=' ' & code2=' ' & code3=' ' apa=' ' & asep=' ' & code3b=' ' elt=' ' & eltb=' ' while (not EOF(10)) do begin rdl: readf,10,bess,pa,elt,sep,code1,code2,weight,code3,code3b,apa,asep, $ format='(t7,f9.4,t20,f6.2,t32,a1,f10.6,t50,a1,a2,f5.1,t64,a3,t66,a3,t20,a5,t35,a5)' ; check for blank line (spacer for multiple plots) if (bess eq 0.) then goto, ENDREAD ; check for bad entries if ((apa eq ' ') or (asep eq ' ') or $ (elt eq '>') or (code1 eq 'Z')) then goto, rdl ; transform to x,y pa=double(pa+180.)/!radeg ; change pa to radians pa=pa+corr*(equ-bess) x(npt)=-sin(pa)*sep y(npt)= cos(pa)*sep wt(npt)=weight resid,p,a,xi,xnode,t0,ecc,omega,bess,xcc,ycc xcalc(npt)=xcc ycalc(npt)=ycc ; determine proper symbol for data point nsymbol(npt)=2 ; speckle, etc. if ((code1 eq 'X') and (code2 eq 'VI')) then nsymbol(npt)=1 ; eyepiece interf. if (code1 eq 'V') then nsymbol(npt)=3 ; micrometry if ((code1 eq 'V') and (code2 eq 'PH')) then nsymbol(npt)=4 ; photographic if ((code1 eq 'O') and (code2 eq 'CC')) then nsymbol(npt)=5 ; occultation if ((code1 eq 'H') and (code2 eq 'IP')) then nsymbol(npt)=6 ; Hipparcos if ((code1 eq 'W') and (code2 ne 'RH')) then nsymbol(npt)=7 ; WSI if ((code1 eq 'T') and (code2 eq 'YC')) then nsymbol(npt)=8 ; Tycho2 if ((code1 eq 'X') and (code3 eq 'MKT')) then nsymbol(npt)=9 ; multi-aperture if ((code1 eq 'X') and (code3 eq 'NOI')) then nsymbol(npt)=9 ; interferometry if ((code1 eq 'X') and (code3b eq 'MkT')) then nsymbol(npt)=9 ; multi-aperture if ((code1 eq 'X') and (code3b eq 'NOI')) then nsymbol(npt)=9 ; interferometry ; determine symbol size for data point if (nsymbol(npt) eq 1) then weight2=weight/5.0 > 0.1 if (nsymbol(npt) eq 2) then weight2=weight/20.0 > 0.1 if (nsymbol(npt) eq 3) then weight2=weight/0.5 > 0.1 if (nsymbol(npt) eq 4) then weight2=weight/1.0 > 0.1 if (nsymbol(npt) eq 5) then weight2=2.0 if (nsymbol(npt) eq 6) then weight2=weight/20.0 > 0.1 if (nsymbol(npt) eq 7) then weight2=weight/5.0 > 0.1 if (nsymbol(npt) eq 8) then weight2=weight/20.0 > 0.1 if (nsymbol(npt) eq 9) then weight2=weight/20.0 > 0.1 ; ssize(npt)=(0.5*sqrt(weight2)) > 0.35 < 1.0 ssize(npt)=1.0 ; set all symbols to same size for catalog plots npt=npt+1 ; increment npt endwhile close,10 ; close file ENDREAD: ; truncate vectors npt=npt-1 & x=x(0:npt) & y=y(0:npt) nsymbol=nsymbol(0:npt) & ssize=ssize(0:npt) & wt=wt(0:npt) ; determine xl,yl for line of nodes xl=fltarr(2) & yl=fltarr(2) thetal=xnode rhol=a*(1.-ecc*ecc)/(1.+ecc*cos(omega)) theta2=xnode+!pi rho2=a*(1.-ecc*ecc)/(1.-ecc*cos(omega)) xl(0)=-sin(thetal+!pi)*rhol yl(0)=cos(thetal+!pi)*rhol xl(1)=-sin(theta2+!pi)*rho2 yl(1)=cos(theta2+!pi)*rho2 ; calculate points for wedge xwedge=fltarr(103) & ywedge=fltarr(103) xwedge(0)=0. & ywedge(0)=0. xwedge(102)=0. & ywedge(102)=0. wcenter=2008.0 wdiff=wend-wstart wstep=wdiff/100. nwedge=0 bess1=wstart-wstep for nwedge=1,101 do begin bess1=bess1+wstep resid,p,a,xi,xnode,t0,ecc,omega,bess1,xcc,ycc xwedge(nwedge)=xcc ywedge(nwedge)=ycc endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO orbit,x,y ; Hartkopf Orbit Program ; ULTRIX IDL Version 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 print, ' P, a, i, W: ',p,a,xi,xnode, format='(a14,4f13.7)' print, ' T, e, w, pm: ',t0,ecc,omega,pm,format='(a14,4f13.7)' print, ' ' degrad=57.2957795 ee=sqrt((1.+ecc)/(1.-ecc)) aee=a*(1.-ecc*ecc) omega=omega/degrad xi=xi/degrad cosi=cos(xi) xnode=xnode/degrad norb=400 ; no. of points for curve dele=(2.*!pi)/(norb-1) xnu=findgen(norb)*dele + 0.0001 r=aee/(1.+ecc*cos(xnu)) theta=xnode+atan(cosi*tan(xnu+omega)) rho=r*cos(xnu+omega)/cos(theta-xnode) count=0 low=where((rho lt 0.),count) if (count gt 0) then begin rho(low)=-rho(low) theta(low)=theta(low)+!pi endif low=where((theta lt 0.),count) if (count gt 0) then theta(low)=theta(low)+2.*!pi high=where((theta gt 2.*!pi),count) if (count gt 0) then theta(high)=theta(high)-2.*!pi x=-sin(theta+!pi)*rho y= cos(theta+!pi)*rho return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO orbpub,file,x,y ; read a published orbit openr,1,file p=fltarr(2,360) readf,1,p close,1 pa=fltarr(360) sep=fltarr(360) for i=0,359 do begin pa(i)=(p(0,i)+180.)/!radeg sep(i)=p(1,i) endfor x=-sin(pa)*sep y= cos(pa)*sep return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plotcirc, x, y, nsymbol, ssize ;routine to plot and fill circle around points. ; written by Ed Dombrowski, 12/89; modified by W. Hartkopf 1/90, 7/90 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 ;define colors (from colortable 12) ; ;c_black=0 & c_red=120 & c_green=15 ;c_dkblue=65 & c_rayin=55 & c_rayout=45 ;c_fuchsia=85 & c_white=255 ;set up x and y dimensions about the point. bx=fltarr(17) & by=fltarr(17) & x2=fltarr(17) & by2=fltarr(17) A=fltarr(17) & for i=0,16 do A(i)=i & A=A*(!pi *2/16.) bx=cos(A) & by=sin(A) ;get the total number of elements in the array x and y. n=N_ELEMENTS(x) ;add the following to avoid max value indicator being plotted. !c=0 ;for loop to plot the figure and polyfill it. for i=0,n-1 do begin bx2=bx*ssize(i) by2=by*ssize(i) if (nsymbol(i) eq 2) then polyfill,bx2+x(i),by2+y(i),color=c_dkblue if (nsymbol(i) le 2) then oplot,bx2+x(i),by2+y(i),color=c_dkblue endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plotsqr, x, y, nsymbol, ssize ; routine to plot and fill square around points (for multi-aperture interf. data). ; written W. Hartkopf 6/00 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 ;define colors (from colortable 12) ; ;c_black=0 & c_red=120 & c_green=15 ;c_dkblue=65 & c_rayin=55 & c_rayout=45 ;c_fuchsia=85 & c_white=255 ;c_black=0 & c_red=155 & c_green=15 ;c_dkblue=80 & c_rayin=167 & c_rayout=175 ;c_fuchsia=106 & c_white=255 ;set up x and y dimensions about the point. bx=fltarr(5) & by=fltarr(5) & bx2=fltarr(5) & by2=fltarr(5) bx(0) =-1.0 & bx(1) = 1.0 & bx(2) = 1.0 & bx(3) =-1.0 & bx(4) =-1.0 by(0) =-1.0 & by(1) =-1.0 & by(2) = 1.0 & by(3) = 1.0 & by(4) =-1.0 ;get the total number of elements in the array x and y. n=N_ELEMENTS(x) ;add the following to avoid max value indicator being plotted. !c=0 ;for loop to plot the figure and polyfill it. for i=0,n-1 do begin bx2=bx*ssize(i) & by2=by*ssize(i) if (nsymbol(i) eq 9) then polyfill,bx2+x(i),by2+y(i),color=c_dkblue endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plothip, x, y, nsymbol, ssize ;routine to plot and fill "H" around points (for Hipparcos data). ; written W. Hartkopf 12/97 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 ;define colors (from colortable 12) ; ;c_black=0 & c_red=120 & c_green=15 ;c_dkblue=65 & c_rayin=55 & c_rayout=45 ;c_fuchsia=85 & c_white=255 ;set up x and y dimensions about the point. bx=fltarr(13) & by=fltarr(13) & bx2=fltarr(13) & by2=fltarr(13) bx(0) = 1.0 & bx(1) = 0.5 & bx(2) = 0.5 & bx(3) =-0.5 bx(4) =-0.5 & bx(5) =-1.0 & bx(6) =-1.0 & bx(7) =-0.5 bx(8) =-0.5 & bx(9) = 0.5 & bx(10)= 0.5 & bx(11)= 1.0 bx(12)= 1.0 by(0) = 1.0 & by(1) = 1.0 & by(2) = 0.25 & by(3) = 0.25 by(4) = 1.0 & by(5) = 1.0 & by(6) =-1.0 & by(7) =-1.0 by(8) =-0.25 & by(9) =-0.25 & by(10)=-1.0 & by(11)=-1.0 by(12)= 1.0 ;get the total number of elements in the array x and y. n=N_ELEMENTS(x) ;add the following to avoid max value indicator being plotted. !c=0 ;for loop to plot the figure and polyfill it. for i=0,n-1 do begin bx2=bx*ssize(i) & by2=by*ssize(i) if (nsymbol(i) eq 6) then polyfill,bx2+x(i),by2+y(i),color=c_red endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plottyc, x, y, nsymbol, ssize ;routine to plot and fill "T" around points (for Tycho data). ; written W. Hartkopf 3/00 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 ;define colors (from colortable 12) ; ;c_black=0 & c_red=120 & c_green=15 ;c_dkblue=65 & c_rayin=55 & c_rayout=45 ;c_fuchsia=85 & c_white=255 ;set up x and y dimensions about the point. bx=fltarr(9) & by=fltarr(9) & bx2=fltarr(9) & by2=fltarr(9) bx(0) = 1.0 & bx(1) =-1.0 & bx(2) =-1.0 & bx(3) =-0.25 bx(4) =-0.25 & bx(5) = 0.25 & bx(6) = 0.25 & bx(7) = 1.0 bx(8) = 1.0 by(0) = 1.0 & by(1) = 1.0 & by(2) = 0.5 & by(3) = 0.5 by(4) =-1.0 & by(5) =-1.0 & by(6) = 0.5 & by(7) = 0.5 by(8) = 1.0 ;get the total number of elements in the array x and y. n=N_ELEMENTS(x) ;add the following to avoid max value indicator being plotted. !c=0 ;for loop to plot the figure and polyfill it. for i=0,n-1 do begin bx2=bx*ssize(i) & by2=by*ssize(i) if (nsymbol(i) eq 8) then polyfill,bx2+x(i),by2+y(i),color=c_red endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plotstar, x, y, nsymbol, ssize ;routine to plot and fill star around points (for USNO data). [W. Hartkopf 5/98] 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 ;define colors (from colortable 12) ; ;c_black=0 & c_red=120 & c_green=15 ;c_dkblue=65 & c_rayin=55 & c_rayout=45 ;c_fuchsia=85 & c_white=255 ;set up x and y dimensions about the point. bx=fltarr(5) & by=fltarr(5) & bx2=fltarr(5) & by2=fltarr(5) bx(0) = 0. & bx(1) = 0.22 & bx(2) = 0.95 & bx(3) = -0.59 bx(4) = 0.0 by(0) = 1.0 & by(1) = 0.31 & by(2) = 0.31 & by(3) = -0.81 by(4) = 1.0 ;get the total number of elements in the array x and y. n=N_ELEMENTS(x) ;add the following to avoid max value indicator being plotted. !c=0 ;for loop to plot the figure and polyfill it. for i=0,n-1 do begin bx2=bx*ssize(i)*1.5 by2=by*ssize(i)*1.5 if (nsymbol(i) eq 7) then polyfill,bx2+x(i),by2+y(i),color=c_dkblue if (nsymbol(i) eq 7) then polyfill,-bx2+x(i),by2+y(i),color=c_dkblue endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO resid, p,a,xi,xnode,t0,ecc,omega,date,xcc,ycc ; determines predicted x,y from orbital elements ee=sqrt((1.+ecc)/(1.-ecc)) aee=a*(1.-ecc*ecc) cosi=cos(xi) xmu=2.*!pi/p em=xmu*(date-t0) e0=em+ecc*sin(em)+0.5*ecc*ecc*sin(2.*em) for l=0,9 do begin em0=e0-ecc*sin(e0) e0=e0+(em-em0)/(1.-ecc*cos(e0)) endfor xnu=2.*atan(ee*tan(0.5*e0)) r=aee/(1.+ecc*cos(xnu)) theta=xnode+atan(cosi*tan(xnu+omega)) rho=r*cos(xnu+omega)/cos(theta-xnode) xcc=-sin(theta+!pi)*rho ycc= cos(theta+!pi)*rho return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plotorb_ps_subs return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>