
function toterr = borestore1(x)
D = x(1);
S = x(2);
load('amps_phase_filtdata.mat');
row = 9; 
col = 1;
A1 = table2array(amplitudes(row,col));
Ph1 = table2array(cor_relphases(row,col));

amp1 = A1; %amplitude from data
phase1 = Ph1; %phase from data

%user inputs%
rw = .075;  %well radius (.076m) minus equipment if applicable (.0167 outside radius pvc)
r = 39.89;  % distance from pumping well (m)
Q0 = .0057; %formation injection amplitude 
P = 2.04;  % period (minutes)

w = 2*pi/P;     % omega

%dimensionless parameters
Cs = .0001*pi*rw^2; 
Co = pi*rw^2;
rD = r/rw;
Cds = Cs/(2*pi*rw^2*S);
Cdo = Co/(2*pi*rw^2*S);
freq = w*rw^2/D; %dimensionless freq
iw = 1i*freq;
        

z1 = (Cdo*sqrt(iw)*besselk(0,sqrt(iw)) + besselk(1,sqrt(iw)))/(Cdo*sqrt(iw));
z2 = (Cds*sqrt(iw)*besselk(0,sqrt(iw)) + besselk(1,sqrt(iw)))/(Cds*sqrt(iw));
g = (besselk(0,rD*sqrt(iw)))^2 - (z1*z2);
hDO = -(besselk(0,rD*sqrt(iw))*besselk(1,sqrt(iw)))/(iw^(3/2)*Cdo*g);

amp = (Q0/(2*pi*D*S))*abs(hDO/Cds); 
phase = angle(hDO/Cds);

errphase = ((phase1 - phase)/phase1)^2; %error in model phase vs data phase
erramp = ((amp1 - amp)/amp1)^2;  %error in model amp vs data amp
toterr = errphase + erramp;      % total error

end
