% Matlab script file to model apparent resistivity
% measured from a Schlumberger array over an earth model 
% that consists of two layers over a half space
%
%	Written by:
%		T. Boyd (4/96)


% First setup default values for all of the input parameters
% If the script has already been run do not overwrite the
% current values

if exist('savee1') == 0
	savee1 = 2.0;
end
if exist('savee2') == 0
	savee2 = 5.0;
end
if exist('saver1') == 0
	saver1 = 500.0;
end
if exist('saver2') == 0
	saver2 = 10.0;
end
if exist('saver3') == 0
	saver3 = 2500.0;
end

% Prompt user for input data file name, open data file, and
% get observed apparent resistivity observations. Input file is assumed to
% contain a series of electrode spacings and apparent resistivities,
% one per line, in the form 'spacing avalue'.

if exist('fid') == 0 % Data file has not been opened or read
	[file_name,path_name] = uigetfile('*.dat','Apparent Resistivity Data File',300,300);
	ctype = computer;	% Get computer type
		% First test for flavors of computers and construct path
	if ctype(1:2) == 'PC'
		full_name = sprintf('%s\%s',path_name,file_name);	% PC
	elseif ctype(1:3) == 'MAC'
		full_name = sprintf('%s:%s',path_name,file_name);	% MAC
	elseif ctype(1:3) == 'VAX'
		full_name = sprintf('%s:[%s]',path_name,file_name);	% VMS
	else
		full_name = sprintf('%s/%s',path_name,file_name);	% UNIX
	end
	fid = fopen(full_name,'r');
	if fid ~= -1
		[A, ndata] = fscanf(fid,'%f %f');
		dist = A(1:2:ndata);	% Array of electrode spacings
		obsar = A(2:2:ndata);	% Array of apparent resistivities
		ndata = round(ndata/2);	% Number of spacings at which resistivities were collected
		clear A; % release original data array
		fclose(fid);
		fid = 1;
	else
		error('No data file selected');
		clear;
	end
end

% Now prompt user for input

pstring = sprintf('Thickness of first layer (given to nearest half meter) [%.1f]: ',savee1);
e1 = input(pstring);
if isempty(e1)
	e1 = savee1;
end
savee1 = e1;

pstring = sprintf('Thickness of second layer (given to nearest half meter) [%.1f]: ',savee2);
e2 = input(pstring);
if isempty(e2)
	e2 = savee2;
end
savee2 = e2;

pstring = sprintf('Resistivity of first layer (ohm/m) [%.1f]: ',saver1);
r1 = input(pstring);
if isempty(r1)
	r1 = saver1;
end
saver1 = r1;

pstring = sprintf('Resistivity of second layer (ohm/m) [%.1f]: ',saver2);
r2 = input(pstring);
if isempty(r2)
	r2 = saver2;
end
saver2 = r2;

pstring = sprintf('Resistivity of underlying halfspace (ohm/m) [%.1f]: ',saver3);
r3 = input(pstring);
if isempty(r3)
	r3 = saver3;
end
saver3 = r3;


% Set depth scale length and convert depths to integer increments of this length

eo = 0.5;
if e1 < 1.0
	e1 = 1.0;
end
e1 = round(e1 / eo) * eo; % assure that e1 and e2 are multiples of eo
e2 = round(e2 / eo) * eo;

% compute reflection coefficients and constants used to compute Q

k1 = (r2 - r1) / (r2 + r1);
k2 = (r3 - r2) / (r3 + r2);
n1 = 0;
n2 = floor(e1 / eo);
n3 = floor(e2 / eo);
n4 = floor( (e1 + e2) / eo);

% set convergence criteria and compute maximum number of terms necessary for convergence

conver = 5.0; % Acceptable error in Ohm-meters
nterms = floor(sqrt(dist(ndata)^3 / (8.0 * conver)));
if nterms < 2*n4
	nterms = 2*n4
end
if nterms > 1000.0
	fprintf(1,'Solution truncated at 1000 terms\n');
	fprintf(1,'It may be inaccurate at large electrode spacings\n');
	nterms = 1000.0;
end

% compute p and h

p = (1:1:nterms) * 0.0;
h = (1:1:nterms) * 0.0;
p(n2) = k1;
p(n4) = k2;
h(n3) = k1 * k2;

% now compute q

q = (1:1:nterms)* 0.0 + p;
for m = 2:nterms
	upper = min(n4,m-1);
	for i = 1:upper
		q(m) = q(m) + (p(i) - h(i)) * q(m-i);
	end
end

% now compute apparent resistivity

ar = (1:1:ndata) * 0.0;
for i = 1:ndata
	l = dist(i);
	voltage = 0.0;
	for m = 1:nterms
		voltage = voltage + (q(m) * (1.0 + 4.0 * m^2 / (l/eo)^2)^(-1.5));
	end
	ar(i) = r1 * (1.0 + 2.0 * voltage);
end
% now compute RMS error for the model
	smr = 0;
for i = 1:ndata
	error = (ar(i) - obsar(i))^2;
	smr = smr + error;
	end
	rms = sqrt(smr)/ndata;
pstring = sprintf('r1 = %.1f r2 = %.1f r3 = %.1f z1 = %.1f z2 = %.1f RMS =%.1f',r1,r2,r3,e1,e2,rms);


% Plot the result
loglog(dist,obsar,'ro', dist,ar)

title(pstring)
xlabel('Electrode Spacing (AB/2) (m)')
ylabel('Apparent Resistivity (Ohm - m)')