close all
clear all
clc

%%

g = 9.81; % [m/s^2]
M = 40;   % [kg]

F_N = M*g; % [N]

IPS.ct_per_m = 7.87 * 1000 * 2 * (2^16/40);
%%

chamberName = {'itmx';'itmy';'ham4';'ham5'};

% Displacement [ct]
                %  Unloaded  Loaded
pos_ct(1).itmx = [ 5981.94  4448.44;  % V1
                   3871.45  2349.48;  % V2
                   1922.08  746.751;  % V3
                   2975.34  1770.58]; % V4
            
pos_ct(2).itmx = [ 5969.13  4457.84;  % V1
                   3856.28  2350.21;  % V2
                   1911.90   757.02;  % V3
                   2954.00  1795.10]; % V4
    
pos_ct(1).itmy = [-7699.07 -9096.09;  % V1
                  -5174.46 -6507.43;  % V2
                  -1100.12 -2402.80;  % V3
                  -5964.78 -7267.80]; % V4
              
pos_ct(2).itmy = [-7802.56 -9051.71;  % V1
                  -5259.05 -6552.60;  % V2
                  -1078.30 -2439.72;  % V3
                  -5959.87 -7222.42]; % V4
    
pos_ct(1).ham4 = [-1505.14 -3079.57;  % V1
                   -481.11 -1895.79;  % V2
                  -1620.80 -2833.71;  % V3
                  -1215.04 -2530.35]; % V4
              
pos_ct(2).ham4 = [-1642.29 -3060.41;  % V1
                   -524.20 -1918.00;  % V2
                  -1576.87 -2776.95;  % V3
                  -1245.68 -2458.46]; % V4      

pos_ct(1).ham5 = [ 1697.28   280.52;  % V1
                   1731.27   432.54;  % V2
                  -3512.46 -4763.14;  % V3
                  -1586.59 -2966.36]; % V4
              
pos_ct(2).ham5 = [ 1352.71   262.20;  % V1
                   1450.48   391.35;  % V2
                  -3502.74 -4773.14;  % V3
                  -1605.18 -2958.40]; % V4
    
springRate = [];              
for iChamber = 1:4
%     disp(['Chamber: ' chamberName{iChamber}])
    for iTrial = 1:3
        if iTrial < 3
            pos_m(iTrial).(chamberName{iChamber}) = pos_ct(iTrial).(chamberName{iChamber}) / IPS.ct_per_m;
            
            % mean displacement of each vertical sensor is equivalent to Z
            % displacement
            disp_m(iTrial).(chamberName{iChamber}) = abs(mean(pos_m(iTrial).(chamberName{iChamber})(:,1)) ...
                                                         - mean(pos_m(iTrial).(chamberName{iChamber})(:,2)));
            
            springRate_Npm(iTrial).(chamberName{iChamber}) = F_N / disp_m(iTrial).(chamberName{iChamber});
        else
            disp_m(iTrial).(chamberName{iChamber}) = abs(mean(pos_m(1).(chamberName{iChamber})(:,2)) ...
                                                         - mean(pos_m(2).(chamberName{iChamber})(:,1)));
            springRate_Npm(iTrial).(chamberName{iChamber}) = F_N / disp_m(iTrial).(chamberName{iChamber});
        end
%         disp(['Trial: ' num2str(iTrial)]);
%         disp(['Disp: ' num2str(disp_m(iTrial).(chamberName{iChamber}))])
%         disp(['Rate: ' num2str(springRate_Npm(iTrial).(chamberName{iChamber}))])
    end
    
    springRate_mean.(chamberName{iChamber}) = mean([springRate_Npm(:).(chamberName{iChamber})]);
    stringRate_unc.(chamberName{iChamber}) = std([springRate_Npm(:).(chamberName{iChamber})]) / sqrt(3);
    disp([upper(chamberName{iChamber}) ' ' num2str(springRate_mean.(chamberName{iChamber}),3) ' (' num2str(stringRate_unc.(chamberName{iChamber}),3) ')'])
    springRate = [springRate;springRate_mean.(chamberName{iChamber}) stringRate_unc.(chamberName{iChamber})];
end

%%
figure(1)
errorbar(1:4,springRate(:,1),springRate(:,2),'*')
set(gca,'XTick',1:4,'XTickLabel',chamberName)
xlabel('Chamber')
ylabel('Spring Rate [N/m]')
title('2015-02-03 H1 HEPI Spring Rates')

FillPage('w')
IDfig(', J. Kissel, H. Radkins, J. Warner')
saveas(gcf,'2015-02-03_HEPI_Spring_Rates.pdf')
    
    