File:Gunners’ equations.svg: Difference between revisions

#! /usr/bin/R -f
q4 <- -10/51; q3 <- -40/357; q2 <- 8045/2499; q1 <- 220/119; q0 <- 220/833
Q <- expression(q4*x^4 + q3*x^3 + q2*x^2 + q1*x + q0)
Y <- function(x) {eval(parse(text=paste("sqrt(",Q,")")),list(x=x))}
Y_ <- function(x) {eval(D(parse(text=paste("sqrt(",Q,")")),"x"),list(x=x))}

x1 <- -pi; y1 <- Y(x1); y1_ <- Y_(x1)

a3 <- -8*y1
a2 <- 4*y1_^2 + 12*q3*x1 + 24*q4*x1^2 - 4*y1_^2 + 4*q2
a1 <- -4*q3*y1_ - 16*q4*x1*y1_ + 8*q4*y1
a0 <- q3^2 - 4*q4*q3*x1 - 8*q4^2*x1^2 + 4*q4*y1_^2 - 4*q4*q2

a_roots <- polyroot(c(a0,a1,a2,a3))

a <- a_roots[1]
b <- y1_ - 2*a*x1
c <- y1 + a*x1^2 - x1*y1_
x <- -(2*a*b - q3 +2*a^2*x1 - 2*q4*x1) / (2*(a^2 - q4))
y <- Y(x)
y_ <- Y_(x)

x_data <- seq(0,sqrt(16.5),0.001)

svg(); plot(main=c("“gunners’ equations” solved to fit a parabolic trajectory",
                    "over a quartic curve given one of two points of tangency"),
    x=c(rev(x_data),-x_data,-rev(x_data),x_data), xlab="x",
    y=c(Y(rev(x_data)),Y(-x_data),-Y(-rev(x_data)),-Y(x_data)), ylab="y",
    type="l", lwd=3, col="brown", asp=1, ylim=c(-4.25,5.75))
abline(0,0)
abline(v=c(x1, x, x1 + (x-x1)/4, x1 + (x-x1)/3, x1 + 2*(x-x1)/3, x1 + 3*(x-x1)/4), col="yellow")
abline(y1 - x1*(y - y1)/(x - x1) - (x - x1)^2*a/4, (y - y1)/(x - x1), col="yellow")
abline(y1 - x1*(y - y1)/(x - x1) - (x - x1)^2*a/3, (y - y1)/(x - x1), col="yellow")
abline(y1 - x1*(y - y1)/(x - x1) - (x - x1)^2*a/2, (y - y1)/(x - x1), col="yellow")
abline(y1 - x1*(y - y1)/(x - x1) - (x - x1)^2*a, (y - y1)/(x - x1), col="yellow")
abline(y1 - x1*(y1_+(y - y1)/(x - x1))/2, (y1_+(y - y1)/(x - x1))/2, col="orange")
abline(y - x*(y_+(y - y1)/(x - x1))/2, (y_+(y - y1)/(x - x1))/2, col="orange")
abline(v=(x + x1)/2, col="purple")
abline(y1 - x1*(y - y1)/(x - x1), (y - y1)/(x - x1), col="purple")

x_vals <- seq(-7,7,0.001)
lines(x_vals, a*x_vals^2 + b*x_vals +c, lwd=2, col="green")
abline(y1 - x1*y1_, y1_, col="red")
abline(y - x*y_, y_, col="red")

text(0,c(5.9, 5.5, 5.1, 4.7),labels=c("y = a*x^2 + b*x + c",
    paste("a = ", a),paste("b = ", b), paste("c = ", c)), col="green")
text(0,-4.3,labels="y^2 = -10/51*x^4 - 40/357*x^3 + 8045/2499*x^2 + 220/119*x + 220/833", col="brown")
text(x1,c(y1-0.2,y1-0.5),pos=1,labels=c(paste("x1 = ", x1),paste("y1 = ", y1)))
text(x-0.9,c(y-0.3,y-0.6),pos=1,labels=c(paste("x = ", x),paste("y = ", y)))

points(c(x1,x),c(y1,y))
points(c((x + x1)/2, (x + x1)/2, (x + x1)/2),
    c((y + y1)/2, (y + y1)/2 - (x - x1)^2*a/4 , (y + y1)/2 - (x - x1)^2*a/2))

c(a*x1^2 + b*x1 + c - y1, 2*a*x1 + b - y1_, a*x^2 + b*x + c - y, 2*a*x + b - y_)
c(a3*a^3 + a2*a^2 + a1*a + a0, a, b, c, x)