% starVelocities.m
% Solution to problem 1 of control 2, Computacion I, 2016/02/08. 
% J.M.Soler, Feb.2016

% Set velocity distribution parameters (teacher only)
clear all
% vg = [330,490];     % mean velocity of each galaxy
% dg = [65,40];       % velocity dispersion of each galaxy
% ns = [3500,1500];   % number of stars in each galaxy

% Find star velocities (teacher only)
% nt = sum(ns);
% ig = 1+(rand(1,nt)<ns(2)/nt);         % galaxy of each star
% vs = vg(ig) + dg(ig).*randn(1,nt);    % velocity of each star
% save('starVelocities.dat','vs','-ascii');

% Get star velocities
vs = load('starVelocities.dat');
nt = numel(vs);

% Find histogram of star velocities
nx = 100;
[y,x] = hist(vs,nx);
dx = x(2)-x(1);       % bin interval

% Set initial parameters for fit (estimated from histogram plot)
vg = [300,500];       % mean velocity of each galaxy
dg = [50,50];         % velocity dispersion of each galaxy
hg = [120,80];        % max number of stars per bin

% Fit distribution with two gaussians
initParams = [vg,dg,hg];
[params,~,~,cov] = nlinfit(x,y,@twoGaussians,initParams);
yfit = twoGaussians(params,x);
vg = params(1:2);
dg = params(3:4);
hg = params(5:6);
ns = hg.*sqrt(2*pi*dg.^2)/dx;
dvg = sqrt(diag(cov(1:2,1:2)))';
ddg = sqrt(diag(cov(3:4,3:4)))';
dhg = sqrt(diag(cov(5:6,5:6)))';
dns = ns.*(dhg./hg+ddg./dg);   % error propagation

% Print results of fit
fprintf('galaxy 1 mean velocity (km/s)       = %6.1f +/- %4.1f\n',vg(1),dvg(1));
fprintf('galaxy 2 mean velocity (km/s)       = %6.1f +/- %4.1f\n',vg(2),dvg(2));
fprintf('galaxy 1 velocity dispersion (km/s) = %6.1f +/- %4.1f\n',dg(1),ddg(1));
fprintf('galaxy 2 velocity dispersion (km/s) = %6.1f +/- %4.1f\n',dg(2),ddg(2));
fprintf('galaxy 1 number of stars            = %6.1f +/- %4.1f\n',ns(1),dns(1));
fprintf('galaxy 2 number of stars            = %6.1f +/- %4.1f\n',ns(2),dns(2));

% Plot histogram of star velocities and fit
bar(x,y), hold on
plot(x,yfit,'r')
xlabel('star velocity (km/s)')
ylabel('number of stars')
