#https://en.wikipedia.org/wiki/Hubble%27s_law
#Objects observed in deep space (extragalactic space, 10 megaparsecs (Mpc) or more) are found to have a Doppler shift interpretable as relative velocity away from Earth;
#This Doppler-shift-measured velocity, of various galaxies receding from the Earth, is approximately proportional to their distance from the Earth for galaxies up to a few hundred megaparsecs away.[1][2]
dist = c(0.032,0.034,0.214,0.263,0.275,0.275,0.45,0.5,0.5,
         0.63,0.8,0.9, 0.9,0.9,0.9,1,1.1,1.1,1.4,1.7,2,2,2,2)
vel = c(170, 290, -130, -70,-185,-220, 200, 290, 270, 200, 
        300, -30, 650, 150, 500, 920, 450, 500, 500, 960, 500, 
        850, 800, 1090)
hubble = data.frame(dist, vel)
hubble
#
#make a plot of velocity (Y) versus distance (x)
# like Excel the x's come first
#
plot(hubble$dist, hubble$vel,
     main = "Hubble's Data",
     xlab = "Distance (Megaparsec)",
     ylab = "Velocity (km/s)")
#
# if you are going to work a lot with a particular set of data 
# you can "attach" it 
# then instead of saying hubble$dist you can just say dist 
# because of the way we constructed our data frame here 
# it is unnecessary 
#
# we can get a "linear model" of the data 
# also known as fitting the data to a line 
# also known as linear regression 
# Excel calls it a TrendLine
#
# the method is lm() and SORRY this time it wants y first 
# and x second with a TILDE in between
#
hub_fit = lm(hubble$vel ~ hubble$dist)
#
# summary got us a lot of info about a set of data 
# summary lso gets us a lot of info about a model 
summary(hub_fit)
# 
# we can add the fit-line (Trendline) to the plot 
# using abline()
#
abline(hub_fit)
#
# You can get help on the lm (linear model) 
#
help(lm)
#
# to display equation extract slope 
# round it or it displays too many decimal places 
# paste together the equation 
# use text to display it 
# note we give the equation an x and y coordinate  
#
intercept=round(hub_fit$coefficients["(Intercept)"],2)
slope=round(hub_fit$coefficients["hubble$dist"],2)
fit_eq = paste("y = ", slope, "x + ", intercept)
text(0.25,800,fit_eq )