% series.m
% anything after a "%" is a comment

clear

load epsilon.txt; % Text file with m strain data, 1 value each line
load a.txt; % Text file with m slit depths
load C.txt; % Text file with m x n C matrix
load P.txt; % Text file with m x n P matrix

m=length(epsilon); 
N=size(C,2); % maximum value of n

for n=1:N; % loop over possible expansion orders
	Cn=C(:,1:n);  % Take submatrix of C for fit order < N
      Pn=P(:,1:n);  % Take submatrix of P for fit order < N    
	A=pinv(Cn)*epsilon; % least squares fit using MATLAB functionality
           % could use A=(((Cn'*Cn)^-1)*Cn')*epsilon  to look like Eq. 5
	sigma(:,n)=Pn*A;  % stresses per Eq. 3
	epsilon_fit=Cn*A; % fit strains per Eq. 4
	u_epsilon_bar(n)=norm(epsilon_fit-epsilon)/sqrt(m-n); 
% avg strain uncert. & misfit per Eq. 13 (just for plotting)
      u_epsilon=sqrt(m/(m-n))*(epsilon-epsilon_fit); 
            % Individual strain uncertainties, Eq. 14
      B=pinv(Cn);  % Eq. 5. Again, could replace with B=(((Cn'*Cn)^-1)*Cn')
      V=B*diag(u_epsilon.^2)*B'; % Eq. 11
	s(:,n)=sqrt(diag(Pn*V*Pn')); % Eq. 8
end

for n=2:N-1; 
     s_model(:,n)=std(sigma(:,n-1:n+1)')'; % Model error per Eq. 16
     s_total(:,n)=sqrt(s(:,n).^2+s_model(:,n).^2); 
           % Total error per Eq. 17
     s_total_bar(n)=norm(s_total(:,n))/sqrt(m); % Avg total error
end

% select minimum total error :
[minimum_s_total_bar,ntemp]=min(s_total_bar(2:N-1)); 
n_star=ntemp+1 % Because “min” starts at 2
        
% echo results to screen and make plots:
sprintf('%s','    x    stress    +-')
fprintf('%8.3f %8.2f %8.3f\n',[a sigma(:,n_star) s_total(:,n_star)]')

figure(1);
plot([2:N-1],s_total_bar(2:N-1),'-o')
title('uncertainty'),xlabel('n');

figure(2);
plot([1:N],u_epsilon_bar,'-o')
title('average strain misfit'),xlabel('n');

figure(3);
errorbar(a,sigma(:,n_star),s_total(:,n_star))
title(sprintf('n = %i',n_star)),xlabel('x'),ylabel('stress');