% eda08_08
% factor analysis on Mid Atlantic Ridge Rocks dataset
% Q-mode factors with geographical plots

% load data
D = load('rocks_with_lat_lon.txt');
sio2 = D(:,1);   % SiO2
tio2 = D(:,2);   % TiO2
als03 = D(:,3);  % Al203
feot = D(:,4);   % FeO-total
mgo = D(:,5);    % MgO
cao = D(:,6);    % CaO
na20 = D(:,7);   % Na2O
k20 = D(:,8);    % K2O
lat = D(:,9);    % lat
lon = D(:,10);   % lon
Ns = size(D);
N = Ns(1);
M = Ns(2)-2;

% compute factors and factor loadings using singular value decompostion
S=D(:,1:M);
[U, SIGMA, V] = svd(S',0);
sigma = diag(SIGMA);
Ns = length(sigma);
F = V';
C = U*SIGMA;

P = 4;
CP = C(:,1:P);
FP = F(1:P,:);

% map bounds
left = -60;
right = 10;
bottom = -70;
top = 70;

% plot singular values
figure(1);
clf;
subplot(1,3,1);
set(gca,'LineWidth',2);
hold on;
axis( [left, right, bottom, top] );
nf = 2;
fmin = mean(FP(nf,:))-2*std(FP(nf,:));
fmax = mean(FP(nf,:))+2*std(FP(nf,:));
fmin = -2*std(FP(nf,:));
fmax = 2*std(FP(nf,:));
for i=[1:N]
    h = 3;
    c = (FP(nf,i)-fmin)/(fmax-fmin);
    if( c<0 )
        c=0;
    elseif ( c>1)
        c=1;
    end
    rectangle('Position',[lon(i), lat(i), h, h], 'FaceColor', [c,c,c]);
end
xlabel('lon, deg');
ylabel('lat, deg');
title('Factor 2');

subplot(1,3,2);
set(gca,'LineWidth',2);
hold on;
axis( [left, right, bottom, top] );
nf = 3;
fmin = -2*std(FP(nf,:));
fmax = 2*std(FP(nf,:));
for i=[1:N]
    h = 3;
    c = (FP(nf,i)-fmin)/(fmax-fmin);
    if( c<0 )
        c=0;
    elseif ( c>1)
        c=1;
    end
    rectangle('Position',[lon(i), lat(i), h, h], 'FaceColor', [c,c,c]);
end
xlabel('lon, deg');
ylabel('lat, deg');
title('Factor 3');

subplot(1,3,3);
set(gca,'LineWidth',2);
hold on;
axis( [left, right, bottom, top] );
nf = 4;
fmin = -2*std(FP(nf,:));
fmax = 2*std(FP(nf,:));
for i=[1:N]
    h = 3;
    c = (FP(nf,i)-fmin)/(fmax-fmin);
    if( c<0 )
        c=0;
    elseif ( c>1)
        c=1;
    end
    rectangle('Position',[lon(i), lat(i), h, h], 'FaceColor', [c,c,c]);
end
xlabel('lon, deg');
ylabel('lat, deg');
title('Factor 4');

% plot coastlines on maps
for myfig = [1:3]
subplot(1,3,myfig);
coastfile = 'coastdata.txt';
fd = fopen(coastfile);
N = 0;
for i=[1:100000]
    line = fgetl(fd);
    if ( ~ischar(line) )
        break;
    end
    if( line(1) == '>' )
        if( N == 0 )
            continue;
        end
        plot( mylon, mylat, 'k.', 'LineWidth', 2 );
        N=0;
    else
        c = sscanf(line,'%f%f');
        N=N+1;
        mylat(N) = c(2);
        mylon(N) = c(1);
    end
end
fclose(fd);
end

    
