; PLOTORB_CAT_SUBS ; Library of subroutines for use with "plotorb_cat.pro" ; Includes the following subroutines: ; orbget orbit orbpub ; plotcirc plotsqr plottri ; plothip plottyc plotstar ; resid plotorb_rv plotorb_cat_subs ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO orbget,star,pos,rx,ry,x,y,nsymbol,ssize,npt,xl,yl,xcalc,ycalc,wt,orb,htype, $ multi,ref,refa,orbflag,xwedge1,ywedge1,xwedge2,ywedge2 ; 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 star=' ' decs=' ' & pos=' ' if (htype eq 'WDS ') 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. endif if (htype eq 'SPECKLE') then begin readf,10,star,rah,ram,decs,decd,decm,pos, $ format='(t17,a13,t56,f2.0,I3,a1,2f2.0,t56,a10)' endif ram=float(ram) ram=ram/10.0 ; read orbital elements orb='y' & ref=' ' & refa=' ' & pcode=' ' & acode=' ' & tcode=' ' if (orbflag eq 'NO') then orb='n' if (orbflag eq 'MAYBE') then readf,10,orb,format='(a1)' if ((orb eq 'y') or (orb eq 'Y')) then begin if (htype eq 'WDS ') then readf,10,p,pcode,a,acode,xi,xnode,t0,tcode,ecc,omega,ref,refa, $ format='(f14.6,a1,f10.5,a1,2f9.4,f13.6,a1,f9.6,f9.4,2x,a7,a1)' if (htype eq 'SPECKLE') then readf,10,p,a,xi,xnode,t0,ecc,omega,pm 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. 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 endif ; read from unit 10 - observational data 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 apa=' ' & asep=' ' & code3=' ' while (not EOF(10)) do begin rdl: readf,10,bess,pa,sep,code3,weight,apa,asep, $ format='(t7,f9.4,t20,f6.2,t33,f10.6,t50,a3,f5.1,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 ' ')) then goto, rdl ; transform (rho,theta) to (x,y) pa=double(pa+180.)/!radeg ; change pa to radians if (orb eq 'Y') then pa=pa+corr*(equ-bess) x(npt)=-sin(pa)*sep y(npt)= cos(pa)*sep wt(npt)=weight if ((orb eq 'Y') or (orb eq 'y')) then begin resid,p,a,xi,xnode,t0,ecc,omega,bess,xcc,ycc xcalc(npt)=xcc ycalc(npt)=ycc endif ; determine proper symbol for data point nsymbol(npt)=3 ; visual (micrometry, etc.) if ((code3 eq 'XVI') or (code3 eq 'EYE')) then nsymbol(npt)= 1 ; eyepiece interf. if ((code3 eq 'SPE') or (code3 eq 'CSP')) then nsymbol(npt)= 2 ; speckle/AO if ((code3 eq 'XSP') or (code3 eq 'WSI')) then nsymbol(npt)= 2 ; speckle/AO if ((code3 eq 'VPH') or (code3 eq 'PTG')) then nsymbol(npt)= 4 ; photographic if (code3 eq 'OCC') then nsymbol(npt)= 5 ; occultation if (code3 eq 'HIP') then nsymbol(npt)= 6 ; Hipparcos if (code3 eq '***') then nsymbol(npt)= 7 ; highlighted measure if (code3 eq 'TYC') then nsymbol(npt)= 8 ; Tycho2 if (code3 eq 'LBI') then nsymbol(npt)= 9 ; multi-aperture if (code3 eq 'CCD') then nsymbol(npt)=10 ; CCD ; 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 ; if (nsymbol(npt) eq 10) 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) if ((orb eq 'y') or (orb eq 'Y')) then begin 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 endif ; calculate points for wedges xwedge1=fltarr(103) & ywedge1=fltarr(103) xwedge1(0)=0. & ywedge1(0)=0. xwedge1(102)=0. & ywedge1(102)=0. xwedge2=fltarr(103) & ywedge2=fltarr(103) xwedge2(0)=0. & ywedge2(0)=0. xwedge2(102)=0. & ywedge2(102)=0. ; wcenter=2007.0 ; wdiff1 = 2.0 & wdiff2 = 0.0 ; wstart1=wcenter-wdiff1 & wend1=wcenter+wdiff1 ; wstart2=wcenter-wdiff2 & wend2=wcenter+wdiff2 wstart1=2010.967 & wend1=2011.035 wstart2=2010.967 & wend2=2011.035 wdiff1 = wend1-wstart1 & wdiff2 = wend2-wstart2 wstep1=wdiff1/50. & wstep2=wdiff2/50. nwedge=0 bess1=wstart1-wstep1 bess2=wstart2-wstep2 for nwedge=1,101 do begin bess1=bess1+wstep1 resid,p,a,xi,xnode,t0,ecc,omega,bess1,xcc,ycc xwedge1(nwedge)=xcc ywedge1(nwedge)=ycc bess2=bess2+wstep2 resid,p,a,xi,xnode,t0,ecc,omega,bess2,xcc,ycc xwedge2(nwedge)=xcc ywedge2(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 [Ed Dombrowski, 12/1989; W. Hartkopf 01/1990, 07/1990] 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 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) set_plot,'x' 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 set_plot,'ps' 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) [W. Hartkopf 06/2000] 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 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) set_plot,'x' & if (nsymbol(i) eq 9) then polyfill,bx2+x(i),by2+y(i),color=c_dkblue set_plot,'ps' & if (nsymbol(i) eq 9) then polyfill,bx2+x(i),by2+y(i),color=c_dkblue endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plottri, x, y, nsymbol, ssize ; routine to plot and fill triangle around points (for CCD data) [W. Hartkopf 07/2011] 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 up x and y dimensions about the point. bx=fltarr(4) & by=fltarr(4) & bx2=fltarr(4) & by2=fltarr(4) bx(0) =-1.0 & bx(1) = 1.0 & bx(2) = 0.0 & bx(3) =-1.0 by(0) =-0.577 & by(1) =-0.577 & by(2) = 1.155 & by(3) =-0.577 ;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*1.2*ssize(i) & by2=by*1.2*ssize(i) set_plot,'x' & if (nsymbol(i) eq 10) then polyfill,bx2+x(i),by2+y(i),color=c_fuchsia set_plot,'ps' & if (nsymbol(i) eq 10) then polyfill,bx2+x(i),by2+y(i),color=c_fuchsia endfor return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plothip, x, y, nsymbol, ssize ; routine to plot and fill "H" around points (for Hipparcos data) [W. Hartkopf 12/1997] 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 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) set_plot,'x' & if (nsymbol(i) eq 6) then polyfill,bx2+x(i),by2+y(i),color=c_red set_plot,'ps' & 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) [W. Hartkopf 03/2000] 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 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) set_plot,'x' & if (nsymbol(i) eq 8) then polyfill,bx2+x(i),by2+y(i),color=c_red set_plot,'ps' & 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 highlighted data) [W. Hartkopf 05/1998] 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 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 set_plot,'x' 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 set_plot,'ps' 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_rv,star,pos ; Program to plot estimated radial velocity curves for a binary ; with known visual elements and estimated mass 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,'x' & wset,2 & !x.style=0 !y.style=0 & !x.minor=0 & !y.minor=0 !p.font=-1 ; distance estimated as 100 pc, masses as equal, date range as 1990 - 2010 dist=100. & mass1=1. & mass2=1. date1=(1995.-p*0.2 < 1992.) > 1990. date2=(1995.+p*1.2 > 2000.) < 2020. au=a*dist a1=mass2*au/(mass1+mass2) a2=mass1*au/(mass1+mass2) ee=sqrt((1.0+ecc)/(1.0-ecc)) k1= 29.786 * sin(xi)*a1/(p*sqrt(1.0-ecc*ecc)) k2=-29.786 * sin(xi)*a2/(p*sqrt(1.0-ecc*ecc)) t=fltarr(3001) & v1=fltarr(3001) & v2=fltarr(3001) & dv=fltarr(3001) if (t0 gt 1990.) then t0=t0-p t0min=1990.-p while (t0 lt t0min) do t0=t0+p e=0.000 for j=0,500 do begin e=e+0.025132741 nu=2.0*atan(ee*tan(e/2)) factor=ecc*cos(omega)+cos(nu+omega) v1(j)=k1*factor & v2(j)=k2*factor dv(j)=abs(v2(j)-v1(j)) & t(j)=t0+0.159155*p*(e-ecc*sin(e)) t(j+500)=t(j)+p & v1(j+500)=v1(j) v2(j+500)=v2(j) & dv(j+500)=dv(j) t(j+1000)=t(j)+2.*p & v1(j+1000)=v1(j) v2(j+1000)=v2(j) & dv(j+1000)=dv(j) t(j+1500)=t(j)+4.*p & v1(j+1500)=v1(j) v2(j+1500)=v2(j) & dv(j+1500)=dv(j) t(j+2000)=t(j)+6.*p & v1(j+2000)=v1(j) v2(j+2000)=v2(j) & dv(j+2000)=dv(j) t(j+2500)=t(j)+8.*p & v1(j+2500)=v1(j) v2(j+2500)=v2(j) & dv(j+2500)=dv(j) endfor vmax=max(dv) ivmax=fix(vmax/25.) +1 v1=v1-5.*float(ivmax) v2=v2-5.*float(ivmax) vmin=min(v1) < min(v2) plot,t,dv,xrange=[date1,date2],yrange=[vmin-1.,vmax+1.],xtitle='Date', $ ytitle='RV in km/sec (for M1=M2, d=100pc)',psym=1, $ symsize=0.35,title=star+' WDS '+pos,charsize=1.25 oplot,t,v1,psym=3,symsize=0.25 oplot,t,v2,psym=3,symsize=0.25 set_plot,'ps' device, filename='plotorb_rv.ps' device, /inches, xoffset=0.5, yoffset=4.0 device, /inches, xsize=7.5, ysize=5.0 !p.font=0 plot,t,dv,xrange=[date1,date2],yrange=[vmin,vmax],xtitle='Date', $ ytitle='RV in km/sec (for M1=M2, d=100pc)',psym=1, $ symsize=0.25,title=star+' WDS '+pos,charsize=1.25 oplot,t,v1,psym=3,symsize=0.25 oplot,t,v2,psym=3,symsize=0.25 device,/close_file return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> PRO plotorb_cat_subs return end ; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>