% Solution to the Diffusion equation
% for a given fixed boundary allowed to
% change through time 

clear

% parameters
D = 61000;

parta = [7.8 -3.41 -0.01 0.749 1.449 5.019 2.207 1.489 2.129];
partb = [-0.532 -0.397 -0.299 -0.226];
ho = [parta partb]; % change in delta head elevation through time (parsed)

% variables
t = 0.1;
x = 1:10000; % all possible integer distances downdelta
y = 25:25:10000; % distances where values are selected from

%matrices
ha = zeros(1,10000);
hb = zeros(15,400);
hc = zeros(16,400);
mdlcore = zeros(13,1);


for i=1:2
	ha(1,x) = ho(1,1)*erfc(x/sqrt(4*D*(t)))-6;
	hc(i,:) = ha(1,y);
	if t == 0.1
	     t = 20;
	else t = t + 20; % time step is 20 years
    end
end


for i=3:14
	for tau = 1:i-1
	   ha(1,x) = ho(1,tau)*erfc(x/sqrt(4*D*(t-tau))); % calc each pulse at fixed t
	   hb(tau,:) = ha(1,y);
	end
	hc(i,:) = sum(hb)-6; % add all pulses together
	t = t + 20;
end

hd = hc';

save test hd /ascii

% location of Auger 2 is at x=3175 --> y(127,1)

auger = [0.5 0.5 0.6 0.8 1.2 3.46 2.1 1.7 2.2 0.2 0.2 0.2 0.2];
auger=auger';


mdlelev = hd(127,2:14);

mdlcore(1,1) = (30.48/20)*(mdlelev(1,1)+6);

for j=1:12
   mdlcore(j+1,1) = (30.48/20).*(mdlelev(1,j+1)-mdlelev(1,j));
end

check = [auger mdlcore]