function []=ForcePole_BM()

% Calculates the body mass of an animal from the force pole output from the
% primate leaping dataset

%#############################
%## -- CODE STARTS BELOW -- ##
%#############################
warning('off','MATLAB:nargchk:deprecated')
close all

% Get path to local toolbox for saving user options'

if ismac
    lTbxPath=toolboxdir('local');
else
    lTbxPath=userpath;
end

% Gravitational acceleration (at Rootstown, according to Wolfram)
AccG=9.80169;

%##########################################
%## ASK USER TO SELECT VIDEO DATA FOLDER ##
%##########################################

if exist(fullfile(lTbxPath,'LeapForcePole_DataDir.mat'),'file')==2
    load(fullfile(lTbxPath,'LeapForcePole_DataDir.mat'));
else
    DataDir=uigetdir('',...
        'Choose the Experiment Data folder for the CMZ Leaping Project...');
    if DataDir==0 return; end
    saveDataDir=questdlg(...
        'Would you like to save this location for future use?',...
        'Save Processed Data Location','Yes','No','Yes');
    if strcmp(saveDataDir,'Yes')
        save(fullfile(lTbxPath,'LeapForcePole_DataDir.mat'),'DataDir');
    end
end

%################################
%## SELECT THE DATE TO ANALYZE ##
%################################

ExptFldr=uigetdir(DataDir,'Choose the experiment to analyze...');
if ExptFldr==0 return; end

%############################################
%## LIST THE AVAILABLE TRIALS FOR ANALYSIS ##
%############################################
TrialList=dir(fullfile(ExptFldr,'*.txt')); 
TrialList=char(TrialList.name);

if isempty(TrialList)
    uiwait(msgbox('The selected experiment contains no usable data.',...
        'Aborting','modal'));
    return
end

TrialList=cellstr(TrialList);
TrialList=strrep(TrialList,'.txt','');

%###########################
%## EXPERIMENT PARAMETERS ##
%###########################
ForceRateList={'50','100','200','500','600'};
FRanswer=listdlg('PromptString','Force plate sampling rate:',...
    'ListString',ForceRateList,...
    'ListSize',[150 75],...
    'Name','Force Rate Selection');
if isempty(FRanswer) return; end
ForceRate=str2double(ForceRateList{FRanswer});


[bmt,ok]=listdlg('PromptString',...
                'Select Body Mass Trial',...
                'Name','Select Trials','ListString',TrialList);
if ok==0 return; end

%#########################
%## CALCULATE BODY MASS ##
%#########################
BMtrial=TrialList{bmt};
BMdata=importdata(fullfile(ExptFldr,[BMtrial,'.txt']));
if size(BMdata,2)>6
    BMdata=BMdata(:,[3,9]); % Just vertical force
else
    BMdata=BMdata(:,3); % Just vertical force
end

% Baseline the data
BMtime=(1:size(BMdata,1))/ForceRate;
BMdata=baseline(BMtime,BMdata);

% Choose loading plateau
BMdata=sum(BMdata,2);
fg=figure;
scsz=get(0,'screensize');
set(fg,'position',[50 50 scsz(3:4)-100],'nextplot','replacechildren',...
    'numbertitle','off','toolbar','none','menu','none')
set(fg,'name','Total Vertical Force')
plot(BMtime,BMdata)           

yline(0,'--k','LineWidth',2);

% Zoom in on jump
title('Zoom in on the body mass event',...
    'color',[.9 0 0],'fontsize',14,'fontweight','bold')
[ZoomTime,~]=ginput(2);
xlim(ZoomTime);

title('Select the period of stable loading to calculate body weight.',...
    'color',[.9 0 0],'fontsize',14,'fontweight','bold')
[x1,~]=ginput(2);
s=find(BMtime>=x1(1),1,'first');
e=find(BMtime<=x1(2),1,'last');
delete(fg)

% Calculate body weight as mean of vertical force and body mass
% as body weight divided by AccG
BodyWght=mean(BMdata(s:e));
BodyMass=BodyWght/AccG;
fprintf('Body mass for this trial is %3.1f grams\n',BodyMass*1000);

%#############################
%## BEGIN SUBFUNCTIONS HERE ##
%#############################.

function [FrcData_BL,BLidx]=baseline(Time,ForceData,kr,kc,BLidx)
if nargin==2
    kr=size(ForceData,1);
    kc=size(ForceData,2);
    BLidx=[];
elseif nargin==3
    kc=size(ForceData,2);
    BLidx=[];
elseif nargin==4
    BLidx=[];
end

if isempty(BLidx)
    scsz=get(0,'screensize');
    fg=figure;
    set(fg,'position',[50 50 scsz(3:4)-100],'nextplot',...
        'replacechildren','numbertitle','off',...
        'toolbar','none','menu','none')
    BL='No';
    while strcmp(BL,'No')
        set(fg,'name','Raw Channels')
        plot(Time,ForceData);
        xlim([Time(1) Time(kr)])
        yline(0,'--k','LineWidth',2);
        title('Select The Baseline Time Period PRIOR to The Active Period.',...
            'color',[.9 0 0],'fontsize',14,'fontweight','bold')
        [x1,~]=ginput(2);
        s1=find(Time>=x1(1),1,'first');
        e1=find(Time<=x1(2),1,'last');

        title('Select The Baseline Time Period AFTER The Active Period.',...
            'color',[0 0 .9])
        [x2,~]=ginput(2);
        s2=find(Time>=x2(1),1,'first');
        e2=find(Time<=x2(2),1,'last');

        set(fg,'name','Baselined Channels')
        bslns=mean(ForceData([s1:e1,s2:e2],:),1);
        FrcData_BL=ForceData-repmat(bslns,kr,1); % All Baselined Data Channels
        plot(Time,FrcData_BL);
        xlim([Time(1) Time(kr)])
        yline(0,'--k','LineWidth',2);

        BL=questdlg('Are you satisfied with this baseline setting?',...
            'Confirm Baseline','Yes','No','Yes');

        BLidx=[s1 e1 s2 e2];
    end
    delete(fg)
else
    bslns=mean(ForceData([BLidx(1):BLidx(2),BLidx(3):BLidx(4)],:),1);
    FrcData_BL=ForceData-repmat(bslns,kr,1); % All Baselined Data Channels
end