pro test0629

snp0=[27.3400,27.3400,27.3300,27.2800,27.2400,27.2100,27.2000,27.2100,27.2300,27.2700,27.2900,27.3000,27.3200,27.3300]
x=[-0.554224354,-0.870084408,-0.996202661,-0.900463399,-0.615854357,-0.204564225,0.264551346,0.663428967,0.927203582,0.999533773,0.861360632,0.523007220,0.078724050]
y=[0.750954032,0.435821632,0.010211775,-0.402947513,-0.730959329,-0.907937549,-0.895938146,-0.699028797,-0.359352296,0.055026800,0.458410238,0.722291855,0.903976706]
z=[0.325489267,0.188873305,0.004362431,-0.174751924,-0.316952531,-0.395015752,-0.388470635,-0.303102172,-0.155839469,0.023809651,0.198692451,0.334774805,0.391868610]
vx=[-14501990,-8461598,-390478,7330288,13365541,16625663,16323805,12717728,6473851,-1229256,-8840223,-14858646,-17428145]
vy=[-8940977,-13941116,-15838038,-14208715,-9663982,-3222311,4068588,10283760,14461160,15704543,13636065,8328533,1238079]
vz=[-3875321,-6043040,-6865462,-6159132,-4188865,-1396186,1764642,4459151,6270202,6809246,5912498,3611514,537627]
vx=vx*1.0E-9
vy=vy*1.0E-9
vz=vz*1.0E-9
rse=fltarr(14) & ve=rse & vs=rse & snpc0=rse & snpc1=rse & snpc2=rse
sdp=25.38
for i=0,12 do begin
    rse(i)=SQRT(x(i)*x(i)+y(i)*y(i)+z(i)*z(i))
    ve(i)=SQRT(vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i))
    vs(i)=rse(i)*2*!pi/sdp                                   ;in AU/day
    ;ve=2*!pi*rse/365.
    snpc0(i)=sdp+(ve(i)*sdp)/vs(i)
    snpc1(i)=sdp+(ve(i)*27.2753)/vs(i)
    snpc2(i)=sdp/(1-ve(i)/vs(i))
endfor
psout,031,'Almanac.ps'
!p.multi=[0,1,5]
plot,x,xrange=[0,12],yrange=[-1.0,1.0],line=1,ystyle=1,ytitle='x,y,z',title='Astronomic Alamanac 2008'
oplot,y,line=2
oplot,z,line=3
plot,rse,xrange=[0,12],yrange=[0.90,1.05],ystyle=1,ytitle='Rse'
oplot,[0,12],[1,1]
plot,ve,xrange=[0,12],yrange=[0.016,0.018],ystyle=1,ytitle="Earth's Ve"
plot,vs,xrange=[0,12],ytitle="Sun's sidereal rotation speed"
plot,snp0,xrange=[0,12],yrange=[27.0,27.6],ystyle=1,ytitle='Synodic period'
oplot,snpc0(0:11),line=1
oplot,snpc1(0:11),line=2
oplot,snpc2(0:11),line=3
oplot,[0,12],[27.2753,27.2753]
!p.multi=0
psdone
print,'Almanac.ps'

;snpc0=[27.3341,27.3185,27.2931,27.2658,27.2416,27.2307,27.2203,27.2275,27.2452,27.2702,27.2973,27.3060,27.3353,27.3400]
;snpc0(13)=snpc0(0)
;snpc1(13)=snpc1(0)
;snpc2(13)=snpc2(0)
;print,'rse:',rse
;print,'ve:',ve
;print,'vs:',vs
;snpc0(11)=0.5*(snpc0(10)+snpc0(12))       ; Calculated snp(11) is wrong!
;snpc1(11)=0.5*(snpc1(10)+snpc1(12))       ; Calculated snp(11) is wrong!
;snpc2(11)=0.5*(snpc2(10)+snpc2(12))       ; Calculated snp(11) is wrong!
;tek_color
;plot,snp0,yrange=[25,30],line=1
;oplot,snpc0,line=2
;oplot,snpc1,line=3
;oplot,snpc2,line=4

end