%hubble
% Loads data downloaded from the NASA Extragalactic Database (NED).
% (2008) https://ned.ipac.caltech.edu/level5/NED05D/ned05D_6.html	
% Note top level website is: https://ned.ipac.caltech.edu/	
% Galaxy distance in Mega-parsecs is plotted against recessional velocitiy
% from Earth in km/s. From this we can deduce the Hubble law:
%
% v = H0 x d
%
% v is the recessional velocity
% H0 is the Hubble constant (at the present Epoch! It will vary with the
% age of the Universe, depending on the parameters of the Friedmann
% equations).
% d is distance.
%
% This function creates a plot of v vs d and overlays the 2020 best
% estimate of the Hubble law (H0 = 72.1 km/s /Mpc
% https://en.wikipedia.org/wiki/Hubble%27s_law
%
% Note since the Universe is expanding, this law would be the same
% regardless of the position of observation.
%
% LAST UPDATED by Andy French April 2017

function hubble

%Hubble constant (2016) km/s /Mpc
H0 = 72.1;

%Conversion from Mega-parsecs to metres
Mpc = 3.086e22;

%Fontsize
fsize = 18;

%dmax /MPc
dmax = 200;

%vmax /km/s
vmax = 14000;

%

%Load NED data
[num,txt,raw] = xlsread('Hubble law raw data.xls');

%Extract distance /Mpc and radial velocitity columns
s=size(raw);
d = [];
v = [];
for n=4:s(1)
    d = [d,raw{n,5}];
    v = [v,raw{n,16}];
end
save hubble d v

%Plot v vs d
figure('color',[1 1 1],'name','Hubble')
plot(d,v,'r+','markersize',3); hold on;
plot( [0,dmax],[0, H0*dmax],'k-','linewidth',2 );
set(gca,'fontsize',fsize);
xlabel('Distance from Earth /Mpc')
ylabel('Radial velocity /kms^{-1}')
title({[num2str(length(d)),' galaxies from NASA Extragalactic Database (2008)',]...
    ['Hubble law v = H_0d where: H_0 = ',num2str(H0),' kms^{-1}/Mpc']})
xlim([0,dmax]); ylim([0,vmax]); grid on;

%Print PNG and close
print(gcf,'Hubble.png','-dpng','-r300');
close(gcf)

%End of code