function Control1_2_JMS(model)
% Solution to problem 2 of Control 1
% JMS, Oct.2022

if model=='A'
  d0 = 0.2;
  V0 = 1;
elseif model=='B'
  d0 = 0.2;
  V0 = 0.5;
elseif model=='C'
  d0 = 0.2;
  V0 = 1.5;
elseif model=='D'
  d0 = 0.2;
  V0 = 2;
else
  error('unknown model')
end
p0 = 1/(10*sqrt(2));
pt = p0*[ [1;1;1],[-1;-1;1],[1;-1;-1],[-1;1;-1] ];
dmin = 0;
dmax = 3*d0;
dd = 0.01*d0;

figure(11)
d = dmin:dd:dmax;
plot(d,Potencial_C1_JMS(d,d0,V0),'LineWidth',1.5)
xlabel ('d (nm)')
ylabel('Energy (eV)')
grid on
set(gca,'FontSize',15)

Dij = Distancia_JMS(pt);
Eij = Potencial_C1_JMS(Dij,d0,V0);
Etot = 0.5*sum(sum(Eij));
fprintf('Etot = %f eV\n',Etot)

x = linspace(0,10*p0,100);
for ix = 1:numel(x)
  p = [x(ix);0;0];
  Dij = Distancia_JMS([pt,p]);
  Eij = Potencial_C1_JMS(Dij,d0,V0);
  E(ix) = 0.5*sum(sum(Eij));
end

figure(12)
plot(x,E,'LineWidth',1.5)
xlabel ('d (nm)')
ylabel('Energy (eV)')
grid on
set(gca,'FontSize',15)

end
