clear; close all
BOHR2ANG = 0.52917721067; HA2EV = 27.21138602;
defhdf5 = 'results/scf_vext.h5';
defmat = 'results/scf_vext.mat';

mrc = load(defmat);
% bulk
vnab = loadDistArray('results/scf_bulk.h5','/potential/vna'); 
vdhb = loadDistArray('results/scf_bulk.h5','/potential/vdh'); 
velb = vnab.data + vdhb.data;
% defect
vna = loadDistArray(defhdf5,'/potential/vna'); 
vdh = loadDistArray(defhdf5,'/potential/vdh'); 
vel = vna.data + vdh.data;

% parse parameters
avec = mrc.domain.latvec;
fgn = mrc.domain.fgridn;
u = (1:fgn(1))/fgn(1);
v = (1:fgn(2))/fgn(2);
w = (1:fgn(3))/fgn(3);
[u3,v3,w3] = ndgrid(u,v,w);

vext = sawtooth(0.5,u3,v3,w3,3);

tmp = squeeze(mean(mean(reshape(vext, fgn),1),2)); 
ind = 1/6 < w & w < 2/6;
p1 = polyfit(w(ind),tmp(ind)',1);

tmp = squeeze(mean(mean(reshape(vel - velb, fgn),1),2)); 
ind = 1/6 < w & w < 2/6;
p2 = polyfit(w(ind),tmp(ind)',1);

epsi = p1(1)/p2(1);
fprintf('dielectric constant = %f\n',epsi);

% figure(); hold on
% z = w * avec(3,3);
% tmp = squeeze(mean(mean(reshape(vext, fgn),1),2)); 
% plot(z,tmp, '--', 'DisplayName', 'v_{ext}'); hold on
% tmp = squeeze(mean(mean(reshape(vel - velb, fgn),1),2)); 
% plot(z,tmp, '--', 'DisplayName', '\Delta V'); hold on
% vext = sawtooth(0.5/p1(1)*p2(1),u3,v3,w3,3);
% tmp = squeeze(mean(mean(reshape(vext, fgn),1),2)); 
% plot(z,tmp, '--', 'DisplayName', 'v_{fit}'); hold on
% xlabel('z (bohr)')
% ylabel('potential (Ha)')
% legend(); grid on
% saveas(gcf, 'vext_diamond.png')