function [Aw_r2 Af_r2 Ah_r2 Aw Af Ah A_all2]=kh_get_pure_terms(tz,tn,sz,sn,Ntz,Nsz,rp)
%
% Script to calculate the proportion of variance explained by each of the
% three pure subduction processes (Aw, Af, Ah).
%
% All input terms are density levels*profiles (ie 40*m)
%
% Written by Kieran Helm (28-11-2007)
%
% Includes SVD check (as on p 121 7.31 of menke1984).  Use with
% subtraction method of calc ntz and nts (21-12-2007)
% A_all2=proportion of variance explained by all 3 terms

warning off
[rows cols]=size(tz);
Aw_r2=NaN(rows,cols);
Af_r2=NaN(rows,cols);
Ah_r2=NaN(rows,cols);
A_all2=NaN(rows,cols);

for c=1:cols % 
        for r=1:rows % 
            x=1/(rp(r,c)-1);
%             x=[sz(r,c)-tz(r,c)]/(rp(r,c)-1); % added as test (18-12-2007)
            M=[-(rp(r,c)-1) 0 -rp(r,c);...
                1 rp(r,c) 0;...
                rp(r,c) rp(r,c) rp(r,c);...
                0 rp(r,c)-1 -1;...
                1 rp(r,c) 0;...
                1 1 1];
            M=x*M; % added coefficient x out the front (24-9-2007)
            d=[tz(r,c); tn(r,c); Ntz(r,c); sz(r,c); sn(r,c); Nsz(r,c)];
            
%             Aw(r,c)=M(:,1)\d;    % (Af=Ah=0)
%             Af(r,c)=M(:,2)\d;    % (Aw=Ah=0)
%             Ah(r,c)=M(:,3)\d;    % (Aw=Af=0)
            Aw(r,c)=inv(M(:,1)'*M(:,1))*M(:,1)'*d; % (Af=Ah=0) full method
            Af(r,c)=inv(M(:,2)'*M(:,2))*M(:,2)'*d; % (Aw=Ah=0)
            Ah(r,c)=inv(M(:,3)'*M(:,3))*M(:,3)'*d; % (Aw=Af=0)
            
            % SVD to check closure:  p==2 
            if isnan(rp(r,c))==0
                [U,S,V]=svd(M);
                A_all=V(:,1:2)*inv(S(1:2,1:2))*U(:,1:2)'*d; 
                A_allest=M*A_all;
                A_all2(r,c)=(A_allest'*A_allest)/(d'*d);
            end
            
            % Calculate R^2
            Aw_dest=M*[Aw(r,c); 0; 0];
            Af_dest=M*[0; Af(r,c); 0];
            Ah_dest=M*[0; 0; Ah(r,c)];
            Aw_r2(r,c)=(Aw_dest'*Aw_dest)/(d'*d);
            Af_r2(r,c)=(Af_dest'*Af_dest)/(d'*d);
            Ah_r2(r,c)=(Ah_dest'*Ah_dest)/(d'*d);
        end
end

