# R source file for Sun-Jupiter-Io-Earth system.
# Source it then run the precoded plots
# Planetarry data and algorithms by Paul Shlyter
# R source by Vincent Goffin
# Fri Jan 25 00:32:56 EST 2002

# Main orbital elements
#    N = longitude of the ascending node
#    i = inclination to the ecliptic (plane of the Earth's orbit)
#    w = argument of perihelion
#    a = semi-major axis, or mean distance from Sun
#    e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
#    M = mean anomaly (0 at perihelion; increases uniformly with time)
#
# Related orbital elements
#    w1 = N + w = longitude of perihelion
#    L  = M + w1 = mean longitude
#    q  = a*(1-e) = perihelion distance
#    Q  = a*(1+e) = aphelion distance
#    P  = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)
#    T  = Epoch_of_M - (M(deg)/360) / P  = time of perihelion
#    v  = true anomaly (angle between position and perihelion)
#    E  = eccentric anomaly

earthPosition <- function(d)
{
	pi <- atan(1)*4
	d2r <- pi/180 

	# Earth orbital elements (angles in degrees)

	N <- 0.0
	i <- 0.0
	w <- 282.9404 + 4.70935E-5 * d
	a <- 1.000000					  # AU
	e <- 0.016709 - 1.151E-9 * d
	M <- 356.0470 + 0.9856002585 * d - 180

	N <- N*d2r
	i <- i*d2r
	w <- w*d2r
	M <- M*d2r

	# Excentric anomaly from the mean anomaly

	E <- eccentricAnomaly(M, e)

	# Compute Earth's distance r and its true anomaly v using
	# xv = r * cos(v) = a * (cos(E) - e)
	# yv = r * sin(v) = a * sqrt(1.0 - e*e) * sin(E)

	xv <- cos(E) - e
	yv <- sqrt(1.0 - e*e) * sin(E)
	v <- atan2(yv, xv)
	r <- sqrt(xv*xv + yv*yv)

	# Compute the planet's position in 3-dimensional space,
	# in heliocentric coordinates

	# xh <- r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) )
	# yh <- r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) )
	# zh <- r * ( sin(v+w) * sin(i) )

	xh <- r * cos(v+w)
	yh <- r * sin(v+w)
	zh <- xh * 0

	matrix(c(xh, yh, zh), ncol=3, byrow=F)
}

jupiterPosition <- function(d)
{
	pi <- atan(1)*4
	d2r <- pi/180 

	# Jupiter orbital elements (angles in degrees)

	N <- 100.4542 + 2.76854E-5 * d
	i <- 1.3030 - 1.557E-7 * d
	w <- 273.8777 + 1.64505E-5 * d
	a <- 5.20256						# AU
	e <- 0.048498 + 4.469E-9 * d
	M <-  19.8950 + 0.0830853001 * d

	N <- N*d2r
	i <- i*d2r
	w <- w*d2r
	M <- M*d2r

	#	First, compute the eccentric anomaly, E, from M, the mean anomaly,
	#	and e, the eccentricity.

	E <- eccentricAnomaly(M, e)

	#	Now compute the planet's distance and true anomaly, using
	#	xv = r * cos(v) = a * ( cos(E) - e )
	#	yv = r * sin(v) = a * ( sqrt(1.0 - e*e) * sin(E) )

	xv <- a * ( cos(E) - e )
	yv <- a * ( sqrt(1.0 - e*e) * sin(E) )
	v <- atan2( yv, xv )
	r <- sqrt( xv*xv + yv*yv )

	# Compute the planet's position in 3-dimensional space,
	# in heliocentric coordinates

	xh <- r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) )
	yh <- r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) )
	zh <- r * ( sin(v+w) * sin(i) )

	matrix(c(xh, yh, zh), ncol=3, byrow=F)
}

epochTime <- function(year, month, day)
{
	# The time scale used here is a "day number" from 2000 Jan 0.0 TDT, which
	# is the same as 1999 Dec 31.0 TDT, i.e. precisely at midnight TDT at the
	# start of the last day of this century.  With the modest accuracy we
	# strive for here, one can usually disregard the difference between
	# TDT (formerly canned ET) and UT.
	#   d  =  JD - 2451543.5  =  MJD - 51543.0
	# We can also compute d directly from the calendar date like this:
	#   d = 367*Y - (7*(Y + ((M+9)/12)))/4 + (275*M)/9 + D - 730530

	# days since jan 01 2000
	367*year - (7*(year + (month + 9)%/%12))%/%4 + (275*month)%/%9 + day - 730530
}

eccentricAnomaly <- function(M, e)
{
	E1 <- M + e * sin(M) * ( 1.0 + e * cos(M) )
	if (e > 0.04) {
		E0 <- E1
		E1 <- E0 - ( E0 - e * sin(E0) - M ) / ( 1 - e * cos(E0) )
	}
	E1
}

##############################################################################
##############################################################################

distance <- function(planet1, planet2)
{
	# planet args are 3 column matrices of x,y, z coords
	dx <- planet1 - planet2
	dx2 <- dx*dx
	sqrt(apply(dx2, 1, sum))
}

plot1 <- function()
{
	# plot Jupiter Earth distance, in AU, over the course of 2 synodic years
	# beginning 2002-01-01

	d0 <- epochTime(2002, 1, 0)
	ysyn <- 399 # days in one synodic Earth-Jupiter year
	dmax <- 2*ysyn
	d <- d0 + 1:dmax

	Epos <- earthPosition(d)
	Jpos <- jupiterPosition(d)
	EJdist <- distance(Epos, Jpos)

	plot(d - d0, EJdist, xlab="Earth days", ylab="Earth Jupiter distance [AU]",
				main="Varying Earth Jupiter Distance", sub="Two Synodic Years")
}

plot2 <- function()
{
	# plot the orbits of Earth an Jupiter, projected on the ecliptic
	# for the length of 1 synodic year beginning 2002-1-1 (opposition)

	d0 <- epochTime(2002, 1, 0)
	ysyn <- 399 # days in one synodic Earth-Jupiter year
	dmax <- 1*ysyn
	d <- d0 + 1:dmax

	Epos <- earthPosition(d)
	Jpos <- jupiterPosition(d)
	x <- c(Epos[,1], Jpos[,1])
	y <- c(Epos[,2], Jpos[,2])
	plot(x, y, pch=".", main="Orbits of Earth and Jupiter", sub="One Synodic Year")

	# mark begin and end of orbit
	x <- c(Epos[1,1], Jpos[1,1])
	y <- c(Epos[1,2], Jpos[1,2])
	points(x, y, col="blue", pch="o")

	n <- length(Epos[,1])
	x <- c(Epos[n,1], Jpos[n,1])
	y <- c(Epos[n,2], Jpos[n,2])
	points(x, y, col="red", pch="o")

	# indicate Sun
	x <- 0
	y <- 0
	points(x, y, col="black", pch="o")

	# draw lines linking begin and end between the Earth and Jupiter orbits
	# since the begin date is opposition, lines should pass through the Sun
	x <- c(Epos[1,1], Jpos[1,1])
	y <- c(Epos[1,2], Jpos[1,2])
	lines(x, y, col="black")
	x <- c(Epos[n,1], Jpos[n,1])
	y <- c(Epos[n,2], Jpos[n,2])
	lines(x, y, col="black")

}

plot3 <- function(Tev=1.77)
{
	# plot the time shifts expected from an event occuring periodically
	# around Jupiter when it is observed from the Earth
	# the default period Tev = 1.77 days is approximately the period of Io.

	ysyn <- 399 			# days in one synodic Earth-Jupiter year
	evMax <- ysyn/Tev		# so that the plot extends about 1 synodic year
	evSeq <- 1:evMax

	# when the event occurs arounf jupiter is is strictly periodic
	tmOccurs <- evSeq*Tev				# [days]

	d0 <- epochTime(2002, 1, 1)

	d <- d0 + tmOccurs 
	Epos <- earthPosition(d)
	Jpos <- jupiterPosition(d)
	EJdist <- distance(Epos, Jpos)

	vlight <- (3E+5)/(149.6E+6)*60*60*24 			# 173 [AU per day]

	# when the event is observed on Earth, it is delayed
	tmObserved <- tmOccurs - EJdist/vlight	# [days]

	plot(1:evMax, tmOccurs - 1.000*tmObserved,
			xlab="Eclipse Sequence", ylab="Deviation [days]",
			main="Observable Effect", sub=paste("Period ", Tev, " Days"))
}

plot4 <- function(Tev=1.77)
{
	# same as plot3 but with 2 slightly different values of Tev

	ysyn <- 399 			# days in one synodic Earth-Jupiter year
	evMax <- ysyn/Tev		# so that the plot extends about 1 synodic year
	evSeq <- 1:evMax

	# when the event occurs arounf jupiter is is strictly periodic
	tmOccurs <- evSeq*Tev				# [days]

	d0 <- epochTime(2002, 1, 1)

	d <- d0 + tmOccurs 
	Epos <- earthPosition(d)
	Jpos <- jupiterPosition(d)
	EJdist <- distance(Epos, Jpos)

	vlight <- (3E+5)/(149.6E+6)*60*60*24 			# 173 [AU per day]

	# when the event is observed on Earth, it is delayed
	tmObserved <- tmOccurs - EJdist/vlight	# [days]

	m <- matrix(c(tmOccurs - 1.00001*tmObserved, tmOccurs - 0.99999*tmObserved), ncol=2, byrow=F)
	matplot(1:evMax, m,
			xlab="Eclipse Sequence", ylab="Deviation [days]",
			main="Sensitivity of Observable Effect")
}