qhat <- 8.6733
rho <- c(0.97949, 0.96973, 0.95996)
qi <- c(18.629, 12.790, 9.9702)

chi <- 3.841

r2 <- rho^2;
r1 <- rho;
r0 <- c(1,1,1) 

A <- cbind(r2,r1,r0)

print(qi - qhat)

y <- qi - qhat - chi

#plot(rho,y)

x <- solve(A,y)

a <- x[1]
b <- x[2]
c <- x[3]

print( (-b - sqrt(b^2-4*a*c))/(2*a))
print( (-b + sqrt(b^2-4*a*c))/(2*a))

#rho <- c(0.87988, 0.86035, 0.79980)
#qi <- c(10.194, 10.634, 11.459)
#rho <- c(0.87988, 0.87012, 0.86035)
#qi <- c(10.194, 10.425, 10.634)
#rho <- c(0.78027, 0.79004, 0.79980)
#qi <- c(11.609, 11.537, 11.459)
#rho <- c(0.88965, 0.87988, 0.87012)
#qi <- c(9.9399, 10.194, 10.425)
rho <- c(0.71973, 0.64941, 0.60059)
qi <- c(11.939, 12.289, 12.544)

r2 <- rho^2;
r1 <- rho;
r0 <- c(1,1,1)

A <- cbind(r2,r1,r0)

print(qi - qhat)

y <- qi - qhat - chi

plot(rho,y)

x <- solve(A,y)

a <- x[1]
b <- x[2]
c <- x[3]

print( (-b - sqrt(b^2-4*a*c))/(2*a))
print( (-b + sqrt(b^2-4*a*c))/(2*a))