diff --git a/FinalCode/BootStrapFunction.m b/FinalCode/BootStrapFunction.m
index 7a57550..121ddd5 100644
--- a/FinalCode/BootStrapFunction.m
+++ b/FinalCode/BootStrapFunction.m
@@ -1,5 +1,7 @@
 function BootStrap = BootStrapFunction(data,Nboot,FieldNames)
-PointEstResults = WIPsubfnFitModel(data);
+
+PointEstResults  = FitProcessModel(data);
+
 % initialize bootstrap values based on the size of the results structure.
 % FieldNames to bootstrap
 
@@ -38,7 +40,7 @@
         Samp = [Samp1; Samp2];
     end
     temp.data = temp.data(Samp,:);
-    Results = WIPsubfnFitModel(temp);
+    Results = FitProcessModel(temp);
     BootStrap.beta(:,:,i) = Results.beta;
     BootStrap.B(:,:,i) = Results.B;
     BootStrap.Paths{i} = Results.Paths;
diff --git a/UnusedCode/WIPsubfnCalculateBCaCI.m b/FinalCode/CalculateBCaCI.m
similarity index 100%
rename from UnusedCode/WIPsubfnCalculateBCaCI.m
rename to FinalCode/CalculateBCaCI.m
diff --git a/UnusedCode/WIPsubfnCalculateBCaLimits.m b/FinalCode/CalculateBCaLimits.m
similarity index 100%
rename from UnusedCode/WIPsubfnCalculateBCaLimits.m
rename to FinalCode/CalculateBCaLimits.m
diff --git a/UnusedCode/WIPsubfnConvertPathsToMatrix.m b/FinalCode/ConvertPathsToMatrix.m
similarity index 85%
rename from UnusedCode/WIPsubfnConvertPathsToMatrix.m
rename to FinalCode/ConvertPathsToMatrix.m
index bbeef31..d58a348 100644
--- a/UnusedCode/WIPsubfnConvertPathsToMatrix.m
+++ b/FinalCode/ConvertPathsToMatrix.m
@@ -1,4 +1,4 @@
-function data = WIPsubfnConvertPathsToMatrix(Paths,PathNumber)
+function data = ConvertPathsToMatrix(Paths,PathNumber)
 % number of bootstraps
 Nrepeats = size(Paths,2);
 
diff --git a/FinalCode/CreateBCaCI.m b/FinalCode/CreateBCaCI.m
new file mode 100644
index 0000000..e2e6805
--- /dev/null
+++ b/FinalCode/CreateBCaCI.m
@@ -0,0 +1,66 @@
+function BCaCI = CreateBCaCI(Results,BootStrap,JackKnife,Thresholds)
+% from the point estimate, the bootstrap and the jacknife results return
+% all of the bias-corrected, accelerated confidence intervals.
+
+% Therefore, the whatever bootstrap estimates are calculated this function
+% creates the BCa confidence intervals for it.
+
+% Find the parameters in the bootstrap
+FieldNames = fieldnames(BootStrap);
+
+Nthresh = length(Thresholds);
+NPaths = size(Results.Paths,1);
+% prepare the structure to hold all of the BCa confidence intervals
+BCaCI = {};
+for i = 1:length(FieldNames)
+    Value = getfield(Results,FieldNames{i});
+    if iscell(Value)
+        BlankValue = zeros(numel(Results.Paths{1}),1,2,NPaths,Nthresh);
+    else
+        BlankValue = zeros([size(Value) 2 Nthresh]);
+    end
+    BCaCI = setfield(BCaCI,FieldNames{i},BlankValue);
+end
+%% put the confidence intervals into this structure
+% cycle over the input fieldnames
+for i = 1:length(FieldNames)
+    if ~isempty(strmatch(FieldNames{i},'Paths'))
+        CurrentBCaCI = zeros(size(getfield(BCaCI,FieldNames{i})));
+        % cycle over all the input thresholds
+        for t = 1:Nthresh
+            CurrentAlpha = Thresholds(t);
+            for j = 1:NPaths
+                % This takes that date from each path and flattens it
+                BootStrapData = ConvertPathsToMatrix(BootStrap.Paths,j);
+                JackKnifeData = ConvertPathsToMatrix(JackKnife.Paths,j);
+                
+                PointEstimate = reshape(Results.Paths{j},numel(Results.Paths{j}),1);
+                
+                % The flattended data allows it to be used by this
+                % program which calculates the adjusted alpha limits used for
+                % determining the confidence intervals
+                [Alpha1 Alpha2 Z p] = CalculateBCaLimits(JackKnifeData,PointEstimate, BootStrapData,CurrentAlpha);
+                
+                % find the confidence intervals for these adjusted
+                % alpha limits
+                CurrentBCaCI(:,:,:,j,t) = CalculateBCaCI(BootStrapData,Alpha1,Alpha2,PointEstimate);
+            end
+        end
+        % put the confidence interval data back into the structure
+        % TO DO: somehow UNFLATTEN the data
+        BCaCI = setfield(BCaCI,FieldNames{i},CurrentBCaCI);
+    else
+        CurrentBCaCI = zeros(size(getfield(BCaCI,FieldNames{i})));
+        for t = 1:Nthresh
+            CurrentAlpha = Thresholds(t);
+            BootStrapData = getfield(BootStrap,FieldNames{i});
+            JackKnifeData = getfield(JackKnife,FieldNames{i});
+            PointEstimate = getfield(Results,FieldNames{i});
+            [Alpha1 Alpha2 Z p] = CalculateBCaLimits(JackKnifeData,PointEstimate, BootStrapData,CurrentAlpha);
+          
+            % find the confidence intervals
+            CurrentBCaCI(:,:,:,t) = CalculateBCaCI(BootStrapData,Alpha1,Alpha2,PointEstimate);
+        end
+        BCaCI = setfield(BCaCI,FieldNames{i},CurrentBCaCI);
+    end
+end
\ No newline at end of file
diff --git a/FinalCode/CycleOverVoxelsProcessBootstrap.m b/FinalCode/CycleOverVoxelsProcessBootstrap.m
index 97a6b97..bcba68f 100644
--- a/FinalCode/CycleOverVoxelsProcessBootstrap.m
+++ b/FinalCode/CycleOverVoxelsProcessBootstrap.m
@@ -8,6 +8,7 @@
 % Prepare the output structure
 Results = cell(Nvoxels,1);
 for i = 1:Nvoxels
+    fprintf(1,'%d of %d voxels\n',i,Nvoxels);
     % Extract the data for this voxel
     OneVoxelModel = ExtractDataFromVoxel(ModelInfo,i);
     % Perform the analysis for this voxel
diff --git a/FinalCode/ExampleSetup.m b/FinalCode/ExampleSetup.m
index 0056c58..a4bbb6a 100644
--- a/FinalCode/ExampleSetup.m
+++ b/FinalCode/ExampleSetup.m
@@ -51,6 +51,8 @@
 BEHAVIOR = randn(Nsub,1);
 AgeGroup = round(rand(Nsub,1));
 
+COVARIATES = randn(Nsub,2);
+
 fprintf(1,'Done preparing data.\n');
 
 %% General Setup
@@ -85,7 +87,7 @@
 % correctsed accelerated confidence intervals determined from bootstrap
 % resampling.
 % Are there any voxel-wise bootstrap resamplings?
-Nboot = 5;
+Nboot = 100;
 
 % An alternative to voxel-wise statistics is a map-wise statistics based
 % off of the maximum statistic approach from a series of permutation
@@ -144,11 +146,20 @@
 DataHeader.dt = [16 0];
 ModelInfo.DataHeader = DataHeader;
 
-%% Model Setup
+%% Model Specific Setup Number One
+% THe first model is to perform a simple mediation analysis with the
+% voxel-wise brain data as the mediation between a dichotomous variable
+% (Age group) and a continuous variable (Behavior). 
+
+%     M
+%    / \
+%   /   \
+%  X     Y
+
 Model1 = ModelInfo;
 
 % Name of the output folder
-Tag = 'ExampleModel1_perm';
+Tag = 'ExampleModel1_boot';
 Model1.Tag = Tag;
 
 % The first model is a basic mediation model testing whether the effect of age
@@ -167,7 +178,7 @@
 % VAR_2 = VAR_1
 % VAR_3 = VAR_1 + VAR_2
 
-Direct = zeros(Nvar);
+Direct = zeros(Model1.Nvar);
 Direct(1,[2 3]) = 1;
 Direct(2,3) = 1;
 
@@ -181,13 +192,13 @@
 
 % For this model there are no interactions so this is kept as a zero
 % matrix.
-Inter = zeros(Nvar);
+Inter = zeros(Model1.Nvar);
 
 % The paths to be tested are included in another matrix the same size of
 % the DIRECT matrix. One difference here is that the PATHS matrix may have
 % a third dimension to test multiple paths within a single model. The steps
 % along a path are specified with integers, i.e. step 1, step 2 ...
-Paths = zeros(Nvar);
+Paths = zeros(Model1.Nvar);
 Paths(1,2) = 1;
 Paths(2,3) = 2;
 
@@ -195,7 +206,8 @@
 Model1.Direct = Direct;
 Model1.Inter = Inter;
 Model1.Paths = Paths;
-%%
+
+%% Run the analysis
 ResultsFolder = PrepareDataForProcess(Model1);
 % The writing of the images should ideally be perfomed by the cluster once
 % the analyses have completed. Therefore, a check is needed or better yet a
@@ -203,6 +215,114 @@
 % e.g. qsub -hold_jid job1,job2 -N job3 ./c.sh
 WriteOutResults(ResultsFolder)
 
+
+%% Model Specific Setup Number Two
+% The second model is to perform a simple mediation analysis with the
+% voxel-wise brain data as the mediation between a dichotomous variable
+% (Age group) and a continuous variable (Behavior). This analysis now
+% includes covariates also.
+
+
+Model2 = ModelInfo;
+
+% Add the covariate names
+Model2.Names = [Model2.Names 'Cov1' 'Cov2'];
+
+% Add the covariate data
+Model2.data = [Model2.data COVARIATES(:,1) COVARIATES(:,2)];
+
+% Update the number of variables value 
+Model2.Nvar = length(Model2.Names);
+
+% Name of the output folder
+Tag = 'ExampleModel2_boot';
+Model2.Tag = Tag;
+
+% Specify the simple mediation model as above
+
+Direct = zeros(Model2.Nvar);
+Direct(1,[2 3]) = 1;
+Direct(2,3) = 1;
+
+% Now specify the covariates in the models
+Direct([4,5],2) = 1;
+Direct([4,5],3) = 1;
+
+% For this model there are no interactions so this is kept as a zero
+% matrix.
+Inter = zeros(Model2.Nvar);
+
+% The paths do not change
+Paths = zeros(Model2.Nvar);
+Paths(1,2) = 1;
+Paths(2,3) = 2;
+
+
+Model2.Direct = Direct;
+Model2.Inter = Inter;
+Model2.Paths = Paths;
+
+ResultsFolder = PrepareDataForProcess(Model2);
+
+% The writing of the images should ideally be perfomed by the cluster once
+% the analyses have completed. Therefore, a check is needed or better yet a
+% wait command for the cluster job: 
+% e.g. qsub -hold_jid job1,job2 -N job3 ./c.sh
+WriteOutResults(ResultsFolder)
+
+
+%% Model Specific Setup Number Three
+% The third model is to perform a moderated-mediation analysis with the
+% voxel-wise brain data as the mediation between a dichotomous variable
+% (Age group) and a continuous variable (Behavior). 
+% Now there is an interaction between the brain measure and age group in
+% predicting the behavior variable
+
+%     M
+%    / \
+%   /  /\
+%  /  /  \
+% X --    Y
+
+Model3 = ModelInfo;
+
+% Update the number of variables value 
+Model3.Nvar = length(Model3.Names);
+
+% Name of the output folder
+Tag = 'ExampleModel3_boot';
+Model3.Tag = Tag;
+
+% Specify the simple mediation model as above
+
+Direct = zeros(Model3.Nvar);
+Direct(1,[2 3]) = 1;
+Direct(2,3) = 1;
+
+
+% For this model there is interaction
+Inter = zeros(Model3.Nvar);
+Inter([1 2],3) = 1
+
+
+% The paths do not change
+Paths = zeros(Model3.Nvar);
+Paths(1,2) = 1;
+Paths(2,3) = 2;
+
+
+Model3.Direct = Direct;
+Model3.Inter = Inter;
+Model3.Paths = Paths;
+
+ResultsFolder = PrepareDataForProcess(Model3);
+
+% The writing of the images should ideally be perfomed by the cluster once
+% the analyses have completed. Therefore, a check is needed or better yet a
+% wait command for the cluster job: 
+% e.g. qsub -hold_jid job1,job2 -N job3 ./c.sh
+WriteOutResults(ResultsFolder)
+
 %%
 
 % TODO 
@@ -224,22 +344,3 @@
 % transient times of large amounts of data in memory. I am afraid to
 % overload a single computer even transiently.
 % 
-% 
-% 
-% WIPsubfnPrepareDataForProcessPERMUTE
-% % THis performs the permutation testing approach
-% WIPsubfnVoxelWiseProcessPERMUTE(InDataPath,count,Nperm)
-% % This applies the BCaCI to each voxel
-% Results = WIPsubfnProcessData(Model)
-%     
-%     BootStrap = WIPsubfnBootStrap(Model,Model.Nboot,FieldNames);
-% 
-%     % Perform the jack-knife step
-%     JackKnife = WIPJackKnife(Model,Results,FieldNames);
-%     
-%     % Calculate the BCaci values for each parameter
-%     Results.BCaCI = WIPsubfnCreateBCaCI(Results,BootStrap,JackKnife,Model.Thresh);
-% 
-% % THis does the actual regression model fitting.
-% Results = WIPsubfnFitModel(data)
-% 
diff --git a/FinalCode/JackKnifeFunction.m b/FinalCode/JackKnifeFunction.m
index 2c43130..1e201df 100644
--- a/FinalCode/JackKnifeFunction.m
+++ b/FinalCode/JackKnifeFunction.m
@@ -22,7 +22,8 @@
     Include = 1:Model.Nsub;
     Include(i) = 0;
     tempModel.data = Model.data(find(Include),:);
-    tempResults = WIPsubfnFitModel(tempModel);
+    tempResults = FitProcessModel(tempModel);
+    
     % store all the parameter estimates in the JackKnife structure
     JackKnife.beta(:,:,i) = tempResults.beta;
     JackKnife.B(:,:,i) = tempResults.B;
diff --git a/FinalCode/OneVoxelProcessBootstrap.m b/FinalCode/OneVoxelProcessBootstrap.m
index d9272a4..2b85451 100644
--- a/FinalCode/OneVoxelProcessBootstrap.m
+++ b/FinalCode/OneVoxelProcessBootstrap.m
@@ -22,7 +22,7 @@
     JackKnife = JackKnifeFunction(Model,Results,FieldNames);
     
     % Calculate the BCaci values for each parameter
-    Results.BCaCI = WIPsubfnCreateBCaCI(Results,BootStrap,JackKnife,Model.Thresholds);
+    Results.BCaCI = CreateBCaCI(Results,BootStrap,JackKnife,Model.Thresholds);
   
 end
 
diff --git a/FinalCode/PrepareDataForProcess.m b/FinalCode/PrepareDataForProcess.m
index 99d8745..a37e13c 100644
--- a/FinalCode/PrepareDataForProcess.m
+++ b/FinalCode/PrepareDataForProcess.m
@@ -71,8 +71,19 @@
     % This analysis is run on the host computer
     switch ModelType
         case 'bootstrap'
-            % Need to split the data and cycle over voxels
-            Results = CycleOverVoxelsProcessBootstrap(ModelInfo);
+            % Save the data
+            InDataPath = fullfile(DataFolder,'ModelInfo');
+            Str = ['save ' InDataPath ' ModelInfo '];
+            eval(Str);
+            % Process the data
+            Parameters = CycleOverVoxelsProcessBootstrap(ModelInfo);
+            
+            % Save the results
+            ResultsPath = fullfile(OutFolder,'Results');
+            mkdir(ResultsPath)
+            Str = sprintf('save %s %s',fullfile(ResultsPath,sprintf('Bootstrap_count%04d_%dSamp',1,ModelInfo.Nboot)),'Parameters');
+            eval(Str)
+            
         case 'permutation'
             % The permutation analysis requires calculation of the point
             % estimate and then calculation of all of the permutations. The
diff --git a/FinalCode/WriteOutParameterMaps.m b/FinalCode/WriteOutParameterMaps.m
index e7ac980..26f26a1 100644
--- a/FinalCode/WriteOutParameterMaps.m
+++ b/FinalCode/WriteOutParameterMaps.m
@@ -34,8 +34,8 @@ function WriteOutParameterMaps(Tag,Parameters,ModelInfo)
                 FileName = sprintf('%s%sX',FileName,ModelInfo.Names{InterVar(j)});
             end
             FileName = sprintf('%s_%s.nii',FileName(1:end-1),Tag);
-            I = zeros(ModelInfor.DataHeader.dim);
-            I(ModelInfo.Indices) = squeeze(temp(Nvar+2,i,:));
+            I = zeros(ModelInfo.DataHeader.dim);
+            I(ModelInfo.Indices) = squeeze(temp(ModelInfo.Nvar+2,i,:));
             % Create the header for this image
             Vo = ModelInfo.DataHeader;
             Vo.fname = fullfile(ModelInfo.ResultsPath,FileName);
diff --git a/FinalCode/WriteOutResults.m b/FinalCode/WriteOutResults.m
index bfd0eb1..b225e8f 100644
--- a/FinalCode/WriteOutResults.m
+++ b/FinalCode/WriteOutResults.m
@@ -67,17 +67,89 @@ function WriteOutResults(ResultsFolder)
                 PointEstimate(:,:,kk,i) = Parameters{i}.Paths{kk};
             end
         end
+        
+        WriteOutPermutationPaths(ModelInfo,MaxPermPaths,MinPermPaths,PointEstimate)
     case 'bootstrap'
+        % locate the results files
+        F = dir(fullfile(ResultsFolder,'Results','Bootstrap*.mat'));
+        NFiles = length(F);
+        
+        % Check to see if the analyses are completed
+        if ~(NFiles == ModelInfo.NJobSplit)
+            errordlg('This analysis is not complete.');
+        end
+        
+        % load a single results fle to determine the size of the paths
+        load(fullfile(ResultsFolder,'Results',F(1).name))
+        [n m] = size(Parameters{1}.Paths{1});
+        PointEstimate = zeros(m,n,length(Parameters));
+        % cycle over number of voxels
+        for i = 1:length(Parameters)
+            PointEstimate(:,:,i) = Parameters{i}.Paths{:};
+            % cycle over thresholds
+            for j = 1:length(ModelInfo.Thresholds)
+                BCaCIUpper(:,:,j,i) = Parameters{i}.BCaCI.Paths(:,:,j,1,1);
+                BCaCILower(:,:,j,i) = Parameters{i}.BCaCI.Paths(:,:,j,1,2);
+            end
+        end
+        
+        WriteOutBootstrapPaths(ModelInfo,PointEstimate,BCaCIUpper,BCaCILower,m,n)
 end
 
-%% WRITE OUT ALL IMAGES
+
+
+
+
+
+%% WRITE OUT ALL IMAGES from the regression models
 WriteOutParameterMaps('beta',Parameters,ModelInfo)
 WriteOutParameterMaps('B',Parameters,ModelInfo)
 WriteOutParameterMaps('t',Parameters,ModelInfo)
 
+%%
+function WriteOutBootstrapPaths(ModelInfo,PointEstimate,BCaCIUpper,BCaCILower,m,n)
+%% WRITE OUT THE PATH IMAGES FOR THE BOOTSTRAP TEST
 
-%% WRITE OUT THE PATH IMAGES FOR THE PERMUTATION TEST
+for i = 1:m % cycle over path number
+    for j = 1:n % cycle over probed level in the path
+        
+        % write unthresholed path estimate
+        Vo = ModelInfo.DataHeader;
+        Vo.fname = (fullfile(ModelInfo.ResultsPath,sprintf('Path%d_level%d.nii',i,j)));
+        % Prepare the data matrix
+        I = zeros(ModelInfo.DataHeader.dim);
+        % Extract the path data values
+        temp = squeeze(PointEstimate(i,j,:));
+        I(ModelInfo.Indices) = temp;
+        % Write the images
+        spm_write_vol(Vo,I);
+    end
+end
+% Write out the thresholded path images
+for k = 1:length(ModelInfo.Thresholds)
+    for i = 1:m % cycle over path number
+        for j = 1:n % cycle over probed level in the path
+            
+            % Find those locations where the product ofthe confidence
+            % intervals is greater then zero. These are locatons where the
+            % CI do not include zero.
+            temp = squeeze(BCaCILower(j,i,k,:).*BCaCIUpper(j,i,k,:)) > 0;
+            Vo = ModelInfo.DataHeader;
+            Vo.fname = (fullfile(ModelInfo.ResultsPath,sprintf('Path%d_level%d_%0.3f.nii',i,j,ModelInfo.Thresholds(k))));
+            % Prepare the data matrix
+            I = zeros(ModelInfo.DataHeader.dim);
 
+            I(ModelInfo.Indices) = temp;
+            % Write the images
+            spm_write_vol(Vo,I);
+        end
+    end
+end
+
+
+
+function WriteOutPermutationPaths(ModelInfo,MaxPermPaths,MinPermPaths,PointEstimate)
+%% WRITE OUT THE PATH IMAGES FOR THE PERMUTATION TEST
 % cycle over the thresholds requested
 for j = 1:length(ModelInfo.Thresholds)
     % Find the number in a sorted list of permutations that corresponds to
@@ -94,9 +166,9 @@ function WriteOutResults(ResultsFolder)
     end
     
     for kk = 1:o % cycle over the number of paths
-        % cycle over ...
-        for i = 1:m 
-            % sort the max and min permutaion results
+        % cycle over the probed value??
+        for i = 1:m
+            % sort the max and min permutation results
             sMax(:,i) = sort(squeeze(MaxPermPaths(i,:,kk,:)),'descend');
             sMin(:,i) = sort(squeeze(MinPermPaths(i,:,kk,:)));
             % find the permutation value based on the sorted values
@@ -115,6 +187,7 @@ function WriteOutResults(ResultsFolder)
             spm_write_vol(Vo,I);
             
             % Find the locations that exceed the two-tailed threshold
+% >>>> I THINK THERE IS A BUG HERE WITH THE SECOND INDEX BEING SET TO 1 <<<
             temp(find((PointEstimate(i,1,kk,:) < Mx(i))&(PointEstimate(i,1,kk,:) >0))) = 0;
             temp(find((PointEstimate(i,1,kk,:) > Mn(i))&(PointEstimate(i,1,kk,:) <0))) = 0;
             
diff --git a/UnusedCode/WIPsubfnCreateBCaCI.m b/UnusedCode/WIPsubfnCreateBCaCI.m
deleted file mode 100644
index 712ce10..0000000
--- a/UnusedCode/WIPsubfnCreateBCaCI.m
+++ /dev/null
@@ -1,55 +0,0 @@
- function BCaCI = WIPsubfnCreateBCaCI(Results,BootStrap,JackKnife,Thresholds)
- % from the point estimate, the bootstrap and the jacknife results return
- % all of the bias-corrected, accelerated confidence intervals
-    FieldNames = fieldnames(BootStrap);
-    Nthresh = length(Thresholds);
-    NPaths = size(Results.Paths,1);
-    BCaCI = {};
-    for i = 1:length(FieldNames)
-        Value = getfield(Results,FieldNames{i});
-        if iscell(Value)
-            BlankValue = zeros(numel(Results.Paths{1}),1,2,NPaths,Nthresh);
-        else
-            BlankValue = zeros([size(Value) 2 Nthresh]);
-        end
-        BCaCI = setfield(BCaCI,FieldNames{i},BlankValue);
-    end
-    %% put the confidence intervals into this structure
-    % cycle over the input fieldnames
-    for i = 1:length(FieldNames)
-        if ~isempty(strmatch(FieldNames{i},'Paths'))
-            CurrentBCaCI = zeros(size(getfield(BCaCI,FieldNames{i})));
-            % cycle over all the input thresholds
-            for t = 1:Nthresh
-                CurrentAlpha = Thresholds(t);
-                for j = 1:NPaths
-                    % This takes that date from each path and flattens it
-                    BootStrapData = WIPsubfnConvertPathsToMatrix(BootStrap.Paths,j);
-                    JackKnifeData = WIPsubfnConvertPathsToMatrix(JackKnife.Paths,j);
-                    PointEstimate = reshape(Results.Paths{j},numel(Results.Paths{j}),1);
-                    % The flattended data allows it to be used by this
-                    % program which calculates the adjusted alpha limits used for
-                    % determining the confidence intervals
-                    [Alpha1 Alpha2 Z p] = WIPsubfnCalculateBCaLimits(JackKnifeData,PointEstimate, BootStrapData,CurrentAlpha);
-                    % find the confidence intervals for these adjusted
-                    % alpha limits
-                    CurrentBCaCI(:,:,:,j,t) = WIPsubfnCalculateBCaci(BootStrapData,Alpha1,Alpha2,PointEstimate);
-                end
-            end
-            % put the confidence interval data back into the structure
-            % TO DO: somehow UNFLATTEN the data
-            BCaCI = setfield(BCaCI,FieldNames{i},CurrentBCaCI);
-        else
-            CurrentBCaCI = zeros(size(getfield(BCaCI,FieldNames{i})));
-            for t = 1:Nthresh
-                CurrentAlpha = Thresholds(t);
-                BootStrapData = getfield(BootStrap,FieldNames{i});
-                JackKnifeData = getfield(JackKnife,FieldNames{i});
-                PointEstimate = getfield(Results,FieldNames{i});
-                [Alpha1 Alpha2 Z p] = WIPsubfnCalculateBCaLimits(JackKnifeData,PointEstimate, BootStrapData,CurrentAlpha);
-                % find the confidence intervals
-                CurrentBCaCI(:,:,:,t) = WIPsubfnCalculateBCaci(BootStrapData,Alpha1,Alpha2,PointEstimate);
-            end
-            BCaCI = setfield(BCaCI,FieldNames{i},CurrentBCaCI);
-        end
-    end
\ No newline at end of file