%% Analyzing model output
load('name'); load(basename);
load underneath  % extra parameters from mf_adapt

%% Grid
[xGr,yGr,zGr,xm,ym,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);

%% Get unformatted data
C=readMT3D('MT3D001.UCN');

NT=length(C);

%% Plot
figure;  hold on; grid on; leg{NT+1}='dummy'; % allocate

clr=repmat('rbgk',1,100);

set(gca,'xlim',[10 130]);
plot(STCONC(:,1,1),ym-ym(1),'k');
leg{1}='initial';
for it=1:NT
    fprintf(' %g',C(it).values(2000,1,1));
    plot(C(it).values(:,1,1),ym-ym(1),clr(it));
    leg{it+1}=sprintf('%10g y',C(it).time);  
end
fprintf('\n');
legend(leg);

xlabel('temp [C]');
ylabel('z [m]');
title('Temperature distribution in geothermal system');

set(gca,'xlim',[0 150],'ylim',[-4000 0]);

%% Compare with numerical and analytical solutions (see text)
% In fact we store "temperature", i.e. heat/rhoc

sigma=sqrt( 2*D/R*[C.time] );      % sigma of continuous spreading
sigma0=sqrt(2*D/R*tp);             % sigma at the end of the system's life

% Analytical temperature due to a pulse according to the total heat
% extracted at t=tp (see theory in the text)
TC=(T0+G*L)+DTEMP*(sigma0*sqrt(8/pi)+W)./(sigma*sqrt(2*pi));

% plot the analytical temperatures on the curves of the numerical simulatin
% at the center of the aquifer
for it=1:NT;
    plot(TC(it),-L,['p',clr(it)]); 
end

%% Compare numerical and analytical heat stored
figure; hold on; grid on;

H0=DTEMP*(sigma0*sqrt(8/pi)+W);  %initial heat/rhoc stored at tp

for it=1:NT
    %compares to see wheater the point match the 45degree line
    plot(TC(it),C(it).values(2000,1,1),[clr(it),'o']);
    HT=sum((C(it).values(:,1,1)-(T0+G*(ym(1)-ym))).*Dy);
    % HT and H0 should be the same for any point in time except at the 
    % end due to loss of heat through the fixe temperature boundaries
    fprintf('%15g%15g%15g%15g\n',TC(it),C(it).values(2000,1,1),HT,H0);
end
xlabel('TC(it)'); ylabel('C(it).values(2000,1,1)');
xlim=get(gca,'xlim');
set(gca,'xlim',[0 xlim(end)], 'ylim',[0 xlim(end)]);
line([0 xlim(end)],[0 xlim(end)]);
title('Compare numerical and analytical solution for temperature in center of aquifer');
