function [x,data] = HowISolvedIt()
%% Load data
[East,Central,Mountain,Pacific] = GetTheData();

%% Combine data to huge array
data = [MatrixForm(East, -5); 
    MatrixForm(Central, -6); 
    MatrixForm(Mountain, -7); 
    MatrixForm(Pacific, -8)];

%% Scatter plots for Sunrise time (converted to UTC)
figure;
subplot(2,2,1);
scatter(data(:,1), data(:,6) - data(:,4), '.');
xlabel('Latitude'); ylabel('Sunrise (UTC)');
subplot(2,2,2);
scatter(data(:,2), data(:,6) - data(:,4), '.');
xlabel('Longitude'); ylabel('Sunrise (UTC)');
subplot(2,2,3);
scatter(data(:,3), data(:,6) - data(:,4), '.');
xlabel('Day'); ylabel('Sunrise (UTC)');
subplot(2,2,4);
scatter(data(:,4), data(:,6) - data(:,4), '.');
xlabel('Timezone'); ylabel('Sunrise (UTC)');

%% Scatter plots for Length of day (in hours)
figure;
subplot(2,2,1);
scatter(data(:,1), data(:,5),'.'); 
xlabel('Latitude'); ylabel('Length of day (hrs)');
subplot(2,2,2);
scatter(data(:,2), data(:,5),'.');
xlabel('Longitude'); ylabel('Length of day (hrs)');
subplot(2,2,3);
scatter(data(:,3), data(:,5),'.');
xlabel('Day'); ylabel('Length of day (hrs)');
subplot(2,2,4);
scatter(data(:,4), data(:,5),'.');
xlabel('Timezone'); ylabel('Length of day (hrs)');

%% Call LSQNONLIN
[x1,resnorm1,res1,exitflag1,output1,lambda1] = ...
    lsqnonlin(@(x) MySunrise(x,data), [6, 3/50, 1.5, 2*pi/365, 0.01, 0.2, 2*pi/200, 0.2, 2*pi/200]);
[x2,resnorm2,res2,exitflag2,output2,lambda2] = ...
    lsqnonlin(@(x) MyLengthOfDay(x,data), [12, 3.5/44, 1, 2*pi/365, 160*pi/365]);
x = [x1,x2];

fprintf('Max norm error in sunrise time = %g mins\n', 60*max(abs(res1)));
fprintf('Max norm error in sunset time = %g mins\n', 60*max(abs(res1+res2)));
fprintf('LSQNONLIN took %d and %d iterations for sunrise and length of day respectively\n',...
    output1.iterations, output2.iterations);
fprintf('Parameters found:\n');
disp(x);

%% Plot residuals for sunrise with factors
figure;
subplot(2,2,1);
scatter(data(:,1), res1, '.');
xlabel('Latitude'); ylabel('Sunrise residual');
subplot(2,2,2);
scatter(data(:,2), res1, '.');
xlabel('Longitude'); ylabel('Sunrise residual');
subplot(2,2,3);
scatter(data(:,3), res1, '.');
xlabel('Day'); ylabel('Sunrise residual');
subplot(2,2,4);
scatter(data(:,4), res1, '.');
xlabel('Timezone'); ylabel('Sunrise residual');

% Plot residuals for length of day with factors
figure;
subplot(2,2,1);
scatter(data(:,1), res2, '.');
xlabel('Latitude'); ylabel('Length of Day residual');
subplot(2,2,2);
scatter(data(:,2), res2, '.');
xlabel('Longitude'); ylabel('Length of Day residual');
subplot(2,2,3);
scatter(data(:,3), res2, '.');
xlabel('Day'); ylabel('Length of Day residual');
subplot(2,2,4);
scatter(data(:,4), res2, '.');
xlabel('Timezone'); ylabel('Length of Day residual');
end

% MySunrise used for LSQNONLIN 
function y = MySunrise( x, data )
% y = Rise - (Timezone + A + B*Long + (C - cos(J*Lat))*cos(D*Day) + K*cos(L*Day) + M*sin(N*Day))
y = data(:,6) - (data(:,4) + x(1) + x(2) .* data(:,2) + (x(3) ...
    - cos(x(5)*data(:,1)) ...
    ) .* cos (x(4) .* data(:,3)) ...
    + x(6)*cos(x(7) .* data(:,3)) ...
    + x(8)*sin(x(9) .* data(:,3)));
end

% MyLengthOfDay used for LSQNONLIN
function y = MyLengthOfDay( x, data )
% y = Hours - (E + (F*Lat - G)*sin(H*Day - I))
y = data(:,5) - (x(1) + (x(2) .* data(:,1) - x(3)) .* sin (x(4) .* data(:,3) - x(5)));
end

%% Internal helper functions for converting data to a large matrix
function data = MatrixForm(Dataset, timezone)
data = [];
for k=1:length(Dataset)
%     Name = Dataset(k).Name; 
    Lat = Dataset(k).Lat; 
    Long = Dataset(k).Long;
    Hours = Dataset(k).Set - Dataset(k).Rise;
    nDays = length(Hours);
    one = ones(nDays,1);
    % Data for sunrise is time-zone corrected (so it is in UTC)
    % NOTE: Assumes that we don't use Daylights Saving Time!
    data = [data;Lat*one, Long*one, (1:nDays)', ...
        timezone*one, Hours, ...
        Dataset(k).Rise, Dataset(k).Set];
end
end