# orbits program
# Fraley June 2001
#  a dynamic model of a system of three bodies: sun, earth, and saturn
#  model 1 = aristotle: earth at center
#  model 2 = copernicus: sun at center
#  In model 1, program models planetary motion and records
#   the angle between saturn and fixed point, with earth at a dynamic origin
#  In model 2, program models planetary motion and 
#    records the angle between saturn and fixed point, with earth at origin 0,0

par(pin=c(10,3.5))
plot(0,0,ylim=c(-2,2),xlim=c(-2,10),axes=F,xlab="",ylab="")
axis(side=4,pos=2,labels=F)
axis(side=2,pos=10,labels=F)
lines(2:10,rep(0,length(2:10)),lty=2)
text(2.4,.2,"0")
text(2.5,2,"180")



model<-2
fixed<-matrix(c(-2,2),2,1)
y<-c(-.50,-.5)     # earth
y1<-c(-1.2,-1.2)       # saturn
trials<-2000
k<-.3  # constant to make out planets move slower than inner planets

rmat1<-rmat<-matrix(c(1,0,0,1),2,2)
points(t(fixed),col=6,pch=4)

if(model==1){
	text(6,-2,"angle between a fixed star and planet")
}

if(model==2){
	text(6,-2,"angle between sun and planet")
}


rad.vec<-seq(0,trials,10)
trials.n<-length(rad.vec)
fripp<-seq(3,10,length=trials.n)


angle.vec<-rep(0,trials.n)
angle.a.vec<-rep(0,trials.n)
angle.b.vec<-rep(0,trials.n)

# begin loop

for(j in 1:trials.n){
i<-rad.vec[j]

if(model==1){
points(t(rmat%*%y),col= 0,pch=16)
points(t(rmat1%*%y1),col= 0,pch=16)
points(0,0,col=17,pch=16,cex=2)
}

if(model==2){
points(t(rmat%*%y),col= 0,cex=2,pch=16)
points(t(rmat1%*%y1),col= 0,pch=16)
points(0,0,col=1,pch=16)
}


# rotate earth
psi<- (i*pi)/180
rmat<-matrix(c(cos(psi),sin(psi),-sin(psi),cos(psi)),2,2)	
if(model==1){points(t(rmat%*%y),pch=16)}
if(model==2){points(t(rmat%*%y),pch=16,col=17,cex=2)}

# rotate saturn
psi<- (i*k*pi)/180
rmat1<-matrix(c(cos(psi),sin(psi),-sin(psi),cos(psi)),2,2)	
points(t(rmat1%*%y1),pch=16)

if(model==1){
# draw earth to saturn vector
lines(c((rmat%*%y)[1],(rmat1%*%y1)[1]),c((rmat%*%y)[2],(rmat1%*%y1)[2]),col=13)
#lines(c((rmat%*%y)[1],(rmat1%*%y1)[1]),c((rmat%*%y)[2],(rmat1%*%y1)[2]),col=0)

# draw earth to fixed point vector
lines(c((rmat%*%y)[1],fixed[1]),c((rmat%*%y)[2],fixed[2]),col=13)
#lines(c((rmat%*%y)[1],fixed[1]),c((rmat%*%y)[2],fixed[2]),col=0)
}

if(model==2){
# draw sun to saturn vector
lines(c(0,(rmat1%*%y1)[1]),c(0,(rmat1%*%y1)[2]),col=13)
#lines(c(0,(rmat1%*%y1)[1]),c(0,(rmat1%*%y1)[2]),col=0)

# draw sun to earth vector
lines(c(0,(rmat%*%y)[1]),c(0,(rmat%*%y)[2]),col=13)
#lines(c(0,(rmat%*%y)[1]),c(0,(rmat%*%y)[2]),col=0)
}

# find angle between planet and fixed-point vectors
# need to shift origin such that earth, not sun, is origin
#  subtract earth vector from all coordinates involved

earth.s<-(rmat%*%y)
fixed.s<-fixed-earth.s
saturn.s<-(rmat1%*%y1)-earth.s
sun.s<- (-earth.s)

cos.theta<-(t(fixed.s)%*%saturn.s)/(sqrt(t(saturn.s)%*%saturn.s)*sqrt(t(fixed.s)%*%fixed.s))
degree.theta<-acos(cos.theta)*180*(1/pi)

angle.vec[j]<-degree.theta


# find angle between sun to saturn vector and sun to fixed-point vectors
# need to shift origin such that earth, not sun, is origin
#  subtract earth vector from all coordinates involved

earth.s<-(rmat%*%y)
fixed.s<-fixed
saturn.s<-(rmat1%*%y1)
sun.s<- c(0,0)

cos.theta<-(t(fixed.s)%*%saturn.s)/(sqrt(t(saturn.s)%*%saturn.s)*sqrt(t(fixed.s)%*%fixed.s))
degree.theta<-acos(cos.theta)*180*(1/pi)

angle.a.vec[j]<-degree.theta


# find angle between sun to saturn vector and sun to earth vectors

earth.s<-(rmat%*%y)
fixed.s<-fixed
saturn.s<-(rmat1%*%y1)
sun.s<- c(0,0)

cos.theta<-(t(earth.s)%*%saturn.s)/(sqrt(t(saturn.s)%*%saturn.s)*sqrt(t(earth.s)%*%earth.s))
degree.theta<-acos(cos.theta)*180*(1/pi)

angle.b.vec[j]<-degree.theta
if(model==1){

    lines( fripp[1:j] ,angle.vec[1:j]*(2/180))
}

if(model==2){
    lines( fripp[1:j] ,angle.b.vec[1:j]*(2/180))
}


if(model==1){
# draw earth to saturn vector
#lines(c((rmat%*%y)[1],(rmat1%*%y1)[1]),c((rmat%*%y)[2],(rmat1%*%y1)[2]),col=1)
lines(c((rmat%*%y)[1],(rmat1%*%y1)[1]),c((rmat%*%y)[2],(rmat1%*%y1)[2]),col=0)

# draw earth to fixed point vector
#lines(c((rmat%*%y)[1],fixed[1]),c((rmat%*%y)[2],fixed[2]),col=1)
lines(c((rmat%*%y)[1],fixed[1]),c((rmat%*%y)[2],fixed[2]),col=0)
}

if(model==2){
# draw sun to saturn vector
#lines(c(0,(rmat1%*%y1)[1]),c(0,(rmat1%*%y1)[2]),col=6)
lines(c(0,(rmat1%*%y1)[1]),c(0,(rmat1%*%y1)[2]),col=0)

# draw sun to earth vector
#lines(c(0,(rmat%*%y)[1]),c(0,(rmat%*%y)[2]),col=6)
lines(c(0,(rmat%*%y)[1]),c(0,(rmat%*%y)[2]),col=0)
}



}# end loop

# end program draw lines that were erased previously


if(model==1){points(t(rmat%*%y),pch=16)}
if(model==2){points(t(rmat%*%y),pch=16,col=17,cex=2)}
if(model==1){
# draw earth to saturn vector
lines(c((rmat%*%y)[1],(rmat1%*%y1)[1]),c((rmat%*%y)[2],(rmat1%*%y1)[2]),col=13)

# draw earth to fixed point vector
lines(c((rmat%*%y)[1],fixed[1]),c((rmat%*%y)[2],fixed[2]),col=13)
}

if(model==2){
# draw sun to saturn vector
lines(c(0,(rmat1%*%y1)[1]),c(0,(rmat1%*%y1)[2]),col=13)

# draw sun to earth vector
lines(c(0,(rmat%*%y)[1]),c(0,(rmat%*%y)[2]),col=13)
}