clc;
clear;

load('matlab.mat');         % Contains all of the gridboxes, containing temperature anomalies
month = 'July';             % These readings are all from this month
[numOfGridBoxes, numOfDates] = size(DATA);

% This magical one line does it all!!!
[U,S,V] = svd(DATA,0); % U gridbox (temperatures), V time ('), S covariance

numOfCoef = min(numOfDates, numOfGridBoxes);        % The number of eigenvalues

% Variance explained by each Eigenvalue
varExplained = zeros(numOfCoef,3);                 
varExplained(:,2) = diag(S.^2) / sum(diag(S.^2));

% Let's see if it really is our data.
DATA(1:20,1:12)             % small segment of data,
reDATA = U*S*V';            % reBuild the Data
reDATA(1:20,1:12)           % notice how its the same as tho original DATA
pause;
%%
% Cummulative variance explained
for counter = 1:numOfCoef
    varExplained(counter,1) = counter;
    varExplained(counter,3) = sum(varExplained(1:counter,2));
end
varExplained
pause;

% Let's get ready to graph it all up. This is just some preprocessing
[numOfLat, dummy] = size(LAT);
[numOfLon, dummy] = size(LON);
meshX = zeros(numOfLat, numOfLon);
meshY = zeros(numOfLat, numOfLon);
GRID  = zeros(numOfLat, numOfLon);
gridCounter = 1;
for counterLat = 1:numOfLat
    for counterLon = 1:numOfLon
        meshX(counterLat,counterLon) = LON(counterLon);
        meshY(counterLat,counterLon) = LAT(counterLat);
        GRID(counterLat,counterLon) = gridCounter;
        gridCounter = gridCounter + 1;
    end
end

% First Graph
figDATA = figure;
figure(figDATA);
% cummulative percent variance explained
subplot(2,1,1);
plot(varExplained(1:50,1), varExplained(1:50,3));
title('Cumulative Sum of Variance Explained');

% percent explained by each EOF
subplot(2,1,2);
plot(varExplained(1:50,1), varExplained(1:50,2));
title('Normalized Eigenvalues of Covariance Matrix');
pause;
close(figDATA);
% Graphs of the first 15 EOFs (looking at the time components)
for eofCounter = 0:2
    % folder and filename prep
    figDATA = figure;
    figure(figDATA);
    %plotting
    for counter = 1:5
        subplot(5,1,counter);
        plot(1:numOfDates,V(:,counter+eofCounter*5));
        title(['EOF ', num2str(counter+eofCounter*5)]);
        xlabel('Time');
        ylabel('Temp Anom');
    end
    pause;
    close(figDATA);
end

% Graphs of the first 5 EOFs (Spacial components)
for eofCounter = 1:5

    EOF = zeros(numOfLat, numOfLon);
    loopCounter = 1;
    for latCounter = 1:1:numOfLat
        for lonCounter = 1:numOfLon
            EOF(latCounter,lonCounter) = U(loopCounter,eofCounter);
            loopCounter = loopCounter + 1;
        end   
    end

    % Contour Plot
    figDATA = figure;
    figure(figDATA);
    [c,h] = contourm(meshY,meshX,EOF);
    clabel(c,h)
    title(['Contour plot of ', month, ' EOF ', num2str(eofCounter)])
    xlabel('Longitude');
    ylabel('Latitude');
    pause;
    close(figDATA);
    % surface contour plot (the cool graph)
    figDATA = figure;
    figure(figDATA);
    surfc(meshX, meshY, EOF);
    title(['Surface plot of ', month, ' EOF ', num2str(eofCounter)]);
    xlabel('Longitude');
    ylabel('Latitude');
    zlabel('temperature');
    % view(-15,25);
    colormap hsv;
    colorbar;
    legend;
    pause;
    close(figDATA);
end