close all

%  Parameters

L = 75;        % Depth of soil column [mm]


% Discretization in space

M = 30+1; % Number of points for space variable z
param.M=M;

dz = L/(M-1); % Step size in space

z = (0:M-1)'*dz; % Discretized points along the depth
param.z=z;

% Data points at 5cm, 10cm, 15cm, pellet position, position of uniform
% initial concentration

param.l10 = find(abs(z-10)==min(abs(z-10)),1);
param.l35 = find(abs(z-35)==min(abs(z-35)),1);
param.l50 = find(abs(z-50)==min(abs(z-50)),1);
param.lpel = find(abs(z-17.5)==min(abs(z-17.5)),1);
param.lbrk2 = find(abs(z-27.5)==min(abs(z-27.5)),1);

l10 =param.l10;
l35 =param.l35;
l50 =param.l50;
lpel=param.lpel;
lbrk2=param.lbrk2;

% Set up a function to solve an ODE obtained by discretizing space

e = ones(M,1);
A = spdiags([e -2*e e], -1:1, M, M);
A(1,1) = -1; A(end,end) = -1; % No flux at the boundaries
A = A/(dz^2); % The second derivative operator


A = [A,sparse(M,M);sparse(M,2*M)]; % Expanding to accomodate the two parts yl(concentration in liquid) and ys(concentration in soil)
param.A=A;

% Other matrices needed in the equation

J1 = [-speye(M),sparse(M,M);speye(M),sparse(M,M)]; % source/sink corresponding to beta1
param.J1=J1;        
J2 = [sparse(M,M),speye(M);sparse(M,M),-speye(M)]; % source/sink corresponding to beta2
param.J2=J2;

T=importdata('Sampling_times.txt'); % Data points in time, discard first 4 points
T=T(5:end);
param.T=T;

% Observed data and standard deviation for three depths, discard first 4
% readings
C11=importdata('CorC11.txt');
C11=C11(5:end);
param.C11=C11;
C12=importdata('CorC12.txt');
C12=C12(5:end);
param.C12=C12;
C13=importdata('CorC13.txt');
C13=C13(5:end);
param.C13=C13;
SD11=importdata('SDCorC11.txt');
SD11=SD11(5:end);
param.SD11=SD11;
SD12=importdata('SDCorC12.txt');
SD12=SD12(5:end);
param.SD12=SD12;
SD13=importdata('SDCorC13.txt');
SD13=SD13(5:end);
param.SD13=SD13;

% Optimization for finding best fit for parameters
options = optimset('TolX', 1e-20, 'MaxFunEvals', 10000, 'Display', 'Iter', 'FinDiffType', 'central','UseParallel',false,'TolFun', 1.0000e-020);%,'PlotFcns',@optimplotfval);
        lb=[0,0.47,0,0,1,1];
        AA= [0,0,-1.213e-3,1,0,0];
        b= [0];
        Aeq = [];
        beq = [];
        ub= [1,10000,1,1,900,900];
        nonlcon = [];
        [x_est,fval] = fmincon(@(x)sol6(x,param),[0.07,8.5,0.011,1.3e-05,9,144], AA,b,Aeq,beq,lb,ub,nonlcon,options);
        
% Estimated parameters to solve y=(D*A+beta1*J1+beta2*J2)x using
% initial value y0.
        
E=x_est(1)*A+x_est(3)*J1+x_est(4)*J2;
y0=x_est(2)*100*[(exp(-((z(1:lpel-1)-17.5).^2)/x_est(5)));(exp(-((z(lpel:lbrk2-1)-17.5).^2)/x_est(6)));exp(-((z(lbrk2)-17.5).^2)/x_est(6))];
y0=[y0;.015*(x_est(3)/x_est(4))+0*y0];

[t,y] = ode45(@(t,y)E*y,[0,335],y0); 

% Time points vector for plots
        
kdp=zeros(size(T)); % Indices corresponding to the data points
for j=1:size(kdp,1)
    kdp(j)=find(abs(t-T(j))==min(abs(t-T(j))),1);
end


% Relative error of the fit
pc_err=((norm(y(kdp,l10)-C11))^2+(norm(y(kdp,l35)-C12))^2+(norm(y(kdp,l50)-C13))^2)/((norm(C11))^2+(norm(C12))^2+(norm(C13))^2);
pc_err=100*pc_err^(1/2)


% Data fit comparision for C11
figure;
plot(t,y(:,l10));
hold on
errorbar(t(kdp),C11,SD11,'o')
ylim([0,max([y(kdp,l10);C11])])
hold off
xlabel('Time [h]','FontSize',20)
ylabel('Concentration [ppm]','FontSize',20)
legend({'model','data C11'},'FontSize',20)

% Data fit comparision for C12
figure;
plot(t,y(:,l35));
hold on
errorbar(t(kdp),C12,SD12,'o')
ylim([0,max([y(kdp,l35);C12])])
hold off
xlabel('Time [h]','FontSize',20)
ylabel('Concentration [ppm]','FontSize',20)
legend({'model','data C12'},'FontSize',20)

% Data fit comparision for C13
figure;
plot(t,y(:,l50));
hold on
errorbar(t(kdp),C13,SD13,'o')
ylim([0,max([y(kdp,l50);C13])])
hold off
xlabel('Time [h]','FontSize',20)
ylabel('Concentration [ppm]','FontSize',20)
legend({'model','data C13'},'FontSize',20)

% Plot of initial concentration or movie of concentration change through
% time
for i=1:1
figure(4)
pause(0.1)
figure(4)
plot(z,y(i,1:M))
% ylim([-1 max(max(y(1,:)))]);
end
hold off
xlabel('Depth [mm]','FontSize',20)
ylabel('Initial concentration [ppm]','FontSize',20)

% Soil digestion data
D=-10*importdata('Depth.txt');
S=importdata('SoilData.txt');

% Comparision of final adsorbed P with soil digestion data
figure;
plot(z,1e-2*y(end,M+1:end))
% ylim([-1 max(max(1e-3*y(end,:)))]);
hold on
plot(D,S)
hold off
xlabel('Depth [mm]','FontSize',20)
ylabel('Final concentration [ppm]','FontSize',20)
legend({'model','data Bound P'},'FontSize',20,'Location','east')