model = mphload('std_phc_pardiso2.mph')

ModelUtil.showProgress(true)

% resolution in k of band diagram, ie each "direction" will consist of N points
N = 90

% redefining a here, since transporting it to the model is too painful for now.
a = 4.4e-7

% reciprocal grid vectors
b1 = 4*pi/(sqrt(3)*a)*[0 -1]
b2 = 2*pi/a * [1 -1/sqrt(3)]

M_point = b1/2
K_point = dot(b2, b1)/norm(b1) * (b1 + b2) / norm(b1 + b2)

bands = []

k_sweep = zeros(3 * N + 1, 2)

for i = 1:1:N,
 i
 k_sweep(i,:) = M_point / N * (i-1)
end
for i = 1:1:N,
 i
 k_sweep(N+i,:) = M_point + (K_point-M_point) / N * (i-1)
end
for i = 1:1:N,
 i
 k_sweep(2*N+i,:) = K_point * (1-(i-1)/N)
end
k_sweep(3*N+1, :) = [0,0]

for i = 1:1:(3*N+1),
 i
 %this shortcut can be taken due to symmetry
 model.param.set('k_x',k_sweep(i,1))
 model.param.set('k_y',k_sweep(i,2))
 model.study('std1').run()
 
 %the real meaning of (normally small) imaginary frequencies for the exact
 %(!) solution (well, FEM, but still) is yet to be known and ignored here.
 bands = cat(2, bands, mphglobal(model,'real(solid.freq)'))
 plot(bands');
 drawnow;
end

 
 
% model.save('C:\Users\noichl\mppcache.mph')




% mphplot(model, 'ph2')
% mphgeom(model, 'geom1')
% mphmesh(model)
% mphnavigator