clear all;
%% load closed-form from prism
xyz    = load('xyz.dat');
x      = xyz(:,1);
y      = xyz(:,2);

%% solution
T      = load('T.dat');   %% our solutions
T2     = load('heath2005_result.dat'); %% reference solutions from prism


lsize=1.5;
fsize=11;
msize= 4;


%% (xx,xy,xz,yx,yy,yz,zx,zy,zz) in T.dat
%% plot tensor (xx,yy,zz)

figure('Position',[0 0 800 600]); 
subplot(3,2,1);
plot(y, T(:,1),'ko', y, T2(:,4), 'k-', 'MarkerSize',  msize, 'LineWidth', lsize);
ylabel('T_{xx}(nT/m)','FontSize', fsize);
set(gca,'fontsize',fsize);

subplot(3,2,2);
plot(y, (T(:,1)-T2(:,4))./T2(:,4)*100, 'ko', 'MarkerSize',  msize, 'LineWidth', lsize);
ylabel('Relative error of T_{xx} (%)','FontSize', fsize);
set(gca,'fontsize',fsize);

subplot(3,2,3);
plot(y, T(:,5), 'ko', y, T2(:,8), 'k-', 'MarkerSize',  msize, 'LineWidth', lsize);
ylabel({'T_{yy}(nT/m)'},'FontSize', fsize);
set(gca,'fontsize',fsize);

subplot(3,2,4);
plot(y, (T(:,5)-T2(:,8))./T2(:,8)*100, 'ko', 'MarkerSize',  msize, 'LineWidth', lsize);
ylabel('Relative error of T_{yy} (%)','FontSize', fsize);
set(gca,'fontsize',fsize);

subplot(3,2,5);
plot(y, T(:,9), 'ko',  y, T2(:,12), 'k-', 'MarkerSize',  msize, 'LineWidth', lsize);
xlabel('y(m)','FontSize', fsize);
ylabel({'T_{zz}(nT/m)'},'FontSize', fsize);
ylim([-5000,1000]);
hl=legend('Our closed-form solutions', 'Closed-form solutions for prism');
set(hl, 'Box', 'off', 'location', 'Best','FontSize', fsize) ;
set(gca,'fontsize',fsize);

subplot(3,2,6);
plot(y, (T(:,9)-T2(:,12))./T2(:,12)*100, 'ko', 'MarkerSize',  msize, 'LineWidth', lsize);
xlabel('y(m)','FontSize', fsize);
ylabel('Relative error of T_{zz} (%)','FontSize', fsize);
set(gca,'fontsize',fsize);


%% plot tensor (yz)
figure('Position',[0 0 800 250]); 
subplot(1,2,1);
plot(y, T(:,6),'ko', y, T2(:,9), 'k-', 'MarkerSize',  msize, 'LineWidth', lsize);
ylabel({'T_{yz}(nT/m)'},'FontSize', fsize);
ylim([-1500,2000]);
hl=legend('Our closed-form solutions', 'Closed-form solutions for prism');
set(hl, 'Box', 'off', 'location', 'Best','FontSize', fsize) ;
xlabel('y(m)','FontSize', fsize);
set(gca,'fontsize',fsize);

subplot(1,2,2);
plot(y, (T(:,6)-T2(:,9))./T2(:,9)*100, 'ko', 'MarkerSize',  msize, 'LineWidth', lsize);
ylabel('Relative error of T_{yz} (%)','FontSize', fsize);
xlabel('y(m)','FontSize', fsize);
set(gca,'fontsize',fsize);