% Unit 4 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%=================================================================
% plot the positions of galaxies from 'GalaxiesInTheUniverse.dat'
%=================================================================
clear all
close all
clc
format long

%% OPTION 1: read in the data using load()
load('Galaxies.dat');
data = Galaxies;
clear Galaxies;

%% OPTION 2: read in the data using fscanf() (Note: there are 4 columns)
fid      = fopen('Galaxies.dat','r');
Ncolumns = 4;
data     = fscanf(fid,'%f %f %f %f');
Nrows    = fix(length(data)/Ncolumns);
data     = reshape(data,Ncolumns,Nrows)';
fclose(fid);

%% plot positions
figure(1)
S=(data(:,4))/min(data(:,4));
C=log10(data(:,4));
scatter3(data(:,1),data(:,2),data(:,3),S)
xlabel('x [Mpc]')
ylabel('y [Mpc]')
zlabel('z [Mpc]')

% generate the histogram and plot on a log-log scale
figure(2)
[Nobj,Mass] = hist(data(:,4));

loglog((Mass),(Nobj),'o')
xlabel('galaxy mass [solar masses]')
ylabel('number of galaxies')
legend('data from cosmological simulation')
ylim([0.1 1e5])

%% in case cftools are installed
% fit power-law model to data
foptions = fitoptions('Method','NonlinearLeastSquares','StartPoint',[-4 60]);
fobject  = fittype('m.*x+c', 'options',foptions);
[pfit,pof] = fit(log10(Mass'),log10(Nobj'),fobject);

m  = pfit.m
c  = pfit.c

hold on
loglog((Mass),10.^(m.*log10(Mass)+c),'r')
legend('cosmological simulation','power-law model')
