pro test0701,rse,ve,rse366,ve366

pv=[[-0.236014486,0.880327681,0.381576410,-16994363,-3847358,-1667831],[-0.401394416,0.828404431,0.359062331,-15991627,-6507444,-2821516],[-0.554224354,0.750954032,0.325489267,-14501990,-8940977,-3875321],[-0.702472054,0.639247556,0.277064278,-12365596,-11307895,-4902750],[-0.814628866,0.516959543,0.224045151,-10009763,-13080140,-5670388],[-0.901788006,0.379050686,0.164264377,-7384226,-14433388,-6256747],[-0.956904262,0.245032539,0.106162899,-4837685,-15287397,-6628125],[-0.990534827,0.089294225,0.038643397,-1877496,-15775976,-6838853],[-0.994449629,-0.068900404,-0.029932074,1090519,-15788339,-6844415],[-0.964701280,-0.240285823,-0.104236354,4302781,-15278833,-6624329],[-0.907662275,-0.388679470,-0.168567079,7067051,-14326740,-6210290],[-0.824182331,-0.525521892,-0.227888579,9587034,-12980922,-5627869],[-0.716825624,-0.647092157,-0.280597434,11832835,-11272666,-4887144],[-0.588778745,-0.749860968,-0.325145531,13707376,-9236147,-4003500],[-0.444024285,-0.831061459,-0.360348047,15177133,-6969805,-3022151],[-0.444024285,-0.831061459,-0.360348047,15177133,-6969805,-3022151],[-0.104293055,-0.922596976,-0.400027773,16823787,-1662210,-720444],[0.064601937,-0.926152832,-0.401572976,16880495,952535,412427],[0.231757404,-0.903565151,-0.391778799,16466916,3557489,1543092],[0.392314074,-0.855391652,-0.370889833,15570479,6050506,2622858],[0.541784919,-0.783131751,-0.339567675,14256570,8371454,3629079],[0.688479885,-0.678065564,-0.294012444,12330733,10677315,4629618],[0.801383993,-0.562169009,-0.243767754,10199069,12446340,5395299],[0.891530565,-0.430250676,-0.186579731,7784468,13880995,6018289],[0.961143877,-0.270848004,-0.117466657,4829017,15013299,6508913],[0.995296391,-0.117524377,-0.050999689,1986086,15578922,6753454],[1.000615416,0.039316120,0.016997746,-934064,15716001,6814205],[0.976544788,0.195204462,0.084584731,-3870941,15380065,6667547],[0.923579290,0.345364742,0.149679772,-6696278,14579637,6320615],[0.843161997,0.485405544,0.210397368,-9356828,13356250,5791099],[0.725410796,0.622637164,0.269892416,-11986315,11513090,4990793],[0.595183828,0.727777696,0.315472602,-13995633,9464435,4103693],[0.446776108,0.810938349,0.351532271,-15614476,7118294,3086288],[0.284515970,0.869300180,0.376831343,-16750367,4521892,1959933],[0.113520427,0.900944658,0.390552329,-17364848,1793410,778382],[-0.061089105,0.904918175,0.392280244,-17465924,-1006624,-436525]]
help,pv
sz=SIZE(pv)
x=FLTARR(sz(2)) & y=x & z=x & rse=x
vx=x & vy=vx & vz=vx & ve=x & vs=x 
for i=0,sz(2)-1 do begin
    x(i)=pv(0,i)
    y(i)=pv(1,i)
    z(i)=pv(2,i)
    vx(i)=pv(3,i)
    vy(i)=pv(4,i)
    vz(i)=pv(5,i)
endfor
;print,x
;print,'y:',y
;print,'z:',z
;print,vx
vx=vx*1.0E-9
vy=vy*1.0E-9
vz=vz*1.0E-9
snpc0=rse & snpc1=rse & snpc2=rse
sdp=25.38
for i=0,sz(2)-1 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
print,'snpc2:',snpc2
print,'snpc1:',snpc1
psout,031,'Almanac36.ps'
!p.multi=[0,1,5]
plot,x,xrange=[0,sz(2)],xsty=1,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,sz(2)],xsty=1,yrange=[0.96,1.04],ystyle=1,ytitle='Rse in AU'
oplot,[0,sz(2)],[1,1]
plot,vs,xrange=[0,sz(2)],xsty=1,yrange=[0.23,0.26],ytitle="Sun's sidereal rotation speed in AU/day"
plot,ve,xrange=[0,sz(2)],xsty=1,yrange=[0.016,0.018],ystyle=1,ytitle="Earth's Ve in AU/day"
;plot,snp0,xrange=[0,sz(2)],yrange=[27.0,27.6],ystyle=1,ytitle='Synodic period'
plot,snpc0,xrange=[0,sz(2)],xsty=1,yrange=[27.15,27.36],ystyle=1,ytitle='Synodic period in days',line=1
;oplot,snpc0,line=1
oplot,snpc1,line=2
oplot,snpc2,line=3
oplot,[0,sz(2)],[27.2753,27.2753]
oplot,[0,sz(2)],[27.2000,27.2000]
oplot,[0,sz(2)],[27.3400,27.3400]
!p.multi=0
psdone
print,'Almanac36.ps'
psout,031,'test0701_2.ps'
!p.multi=[0,1,3]
x=findgen(38)
y2=findgen(38)
y2(0)=!VALUES.F_NAN & y2(37)=!VALUES.F_NAN
y2(1:36)=snpc2
plot,x,y2,xrange=[0.,37.],yrange=[27.15,27.35],xsty=1,ysty=1,title='01:05,01:15,...,12:15,12:25'
rse38=findgen(38)
rse38(0)=rse(35) & rse38(37)=rse(0)
rse38(1:36)=rse
ve38=findgen(38)
ve38(0)=ve(35) & ve38(37)=ve(0)
ve38(1:36)=ve
snp38=findgen(38)
snp38(0)=snpc2(35) & snp38(37)=snpc2(0)
snp38(1:36)=snpc2
plot,x,snp38,xrange=[0.,37.],yrange=[27.15,27.35],xsty=1,ysty=1,title='12:25,01:05,01:15,...,12:25,01:05'
t=findgen(375)/10.0
rse375=SPLINE(x,rse38,t)
ve375=SPLINE(x,ve38,t)
snp375=SPLINE(x,snp38,t)
snp375(0:4)=!VALUES.F_NAN & snp375(370:374)=!VALUES.F_NAN
plot,t,snp375,xrange=[0.,37.5],yrange=[27.15,27.35],xsty=1,ysty=1,title='01:01,01:02,...,12:30,12:31'
rse366=rse375(5:370)
ve366=ve375(5:370)
snp366=snp375(5:370)
help,rse366,ve366,snp366
!p.multi=0
psdone
print,'test0701_2.ps'

end