function [sds,sdskew] = qb(tnspec,btrain,newspec,cnter,radfrac,sensitiv)
% function definition of qb with parameters
% tnspec= which is training spectra
% btrain =bootstrap replicates calculated using routine replica
% newspec= sample spectrum
% cnter= center of calibration set calculated using routine replica
% radfrac= fraction of points in hypercylinder
% sensitiv= sensitivity
%
% Copyright 2003 Robert A. Lodder
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
% 
% Contact Information located at www.pharm.uky.edu or search the World Wide Web
% for Analytical Spectroscopy Research Group.
 

b = size(btrain,1);   %finds the number of rows in btrain
qdist = zeros(1,b);   %creates a null row matrix
s02 = sqrt(sum((newspec-cnter).^2));  %computing the squareroot of sum of the suares
s0r = sqrt(sum((btrain-repmat(cnter,b,1)).^2,2)); %repmat creates a bX1 tilings of cnter
s2r = sqrt(sum((btrain-repmat(newspec,b,1)).^2,2));
sub = (s02+s0r+s2r)/2;   
area = sqrt(sub.*(sub-s02).*(sub-s0r).*(sub-s2r));
radial = (2*area)/s02;
project = sqrt(s0r.^2-radial.^2);
locs = find((s02.^2+s0r.^2) < s2r.^2);  %finds the indices where s02*s02 + s0r*s0r < s2r*s2r(getting the locus)
project(locs) = project(locs)*-1; 
qdist = project;
qrr = sort(radial);  %sorts the elements in ascending order
radii = qrr(radfrac*b);
qdist(find(radial > radii)) = 0; % setting the elments of qdist to zero where radial > radii
qdist = sort(qdist(find(qdist))); % sorts all non zero elememts in qdist
lindex = floor(0.16*length(qdist)); % lower limit found by rouding to nearest integer
uindex = floor(0.84*length(qdist)); % upper limit found by rounding to nearest integer
if(length(qdist) < 50)
    '** Need more replicates in hypercylinder **'
end
sd = std(qdist)*sqrt(size(tnspec,1)); %std returns the standard deviation 
sds = sqrt(sum((cnter-newspec).^2))/sd;  %calculation of the standard deviation distances
% BIAS ADJUSTMENT
alpha = normcdf(-1,0,1); % computes the cumulative distribution function with mean 0 and standard deviation 1
za = norminv(alpha,0,1); % inverse of the cumulative ditribution function with mean o and standard deviation 1
tcenter = median(tnspec); % finding the median 
cs02 = s02;
cs0r = sqrt(sum((tcenter-cnter).^2));
cs2r = sqrt(sum((tcenter-newspec).^2));
csub = (cs02+cs0r+cs2r)/2;
carea = sqrt(csub*(csub-cs02)*(csub-cs0r)*(csub-cs2r));
cradial = (2*carea)/cs02;
cproject = sqrt(cs0r^2-cradial^2);
if((s02^2+cs0r^2) > cs2r^2)
    cproject = -cproject;
end
n = length(qdist); % finds the length of the vector
if(floor(n/2) == n/2)
    md = (qdist(n/2)+qdist(n/2+1))/2;
else
    md = qdist(floor(n/2+0.5));
end
cproject = cproject*sensitiv + md;
fdist = qdist-cproject;
index = 1:length(fdist); 
if(cproject > max(qdist))
    zelement = length(qdist)-1;
elseif(cproject < min(qdist))
    zelement = 1;
else
    rootloc = find(abs(fdist)==min(abs(fdist)));
    zelement = rootloc(1);
end
z0 = norminv(zelement/length(qdist),0,1);
if(abs(2*z0) > abs(za))
    error(' -- Decrease skew sensitivity. --');
end
sensitiv = abs(sensitiv);
lowind = floor(normcdf(2*z0+za,0,1)*length(qdist));
upind = floor(normcdf(2*z0-za,0,1)*length(qdist));
if(lowind < 2)
    '** Warning ** Too few replicates'
end
if(upind > length(qdist)-2)
    '** Warning ** Too few replicates'
end
if(lowind < 1)
    lowind = 1;
end
if(upind > length(qdist))
    upind = length(qdist);
end
lowlim = qdist(lowind);
uplim = qdist(upind);
euc = sqrt(sum((cnter-newspec).^2));
fac = abs(norminv(alpha));
erd = sqrt(size(tnspec,1));
if(abs(2*z0)>abs(fac))
    '** Warning ** SKEW CORRECTION exceeds replicates'
end
sdskew = euc/((uplim/fac)*erd);