e0 = 8.85e-12;
De = 17;
K2 = 7.4e-12;
d = 20;
E= [0.037 0.077 0.117 0.177 0.257];
%E = 0.037;
%E = (0:0.001:0.01);
z = (0:0.001:20);
zc = 2.9;
value = sqrt(e0*De/(K2));

nume =  -pi.*exp(-value.*E.*zc);
denom = 4.*(exp(-value.*E.*zc) + exp(value.*E.*zc));

nume2 = -pi.*exp(-value.*E.*zc);
denom2 = 8.*(cosh(value.*E.*zc));


g = figure();
hold on
for i = 1:length(E)
    k2 = (-pi.*exp(-value.*E(i).*zc))/(8.*(cosh(value.*E(i).*zc)));
    k1 = -(pi./4 + k2);
    f = (pi/4) + k1.*exp(-value.*E(i).*z) + k2.*exp(value.*E(i).*z);
    plot(z,f,'Color','black','LineWidth',1.5)
end

ylabel('Twist angle, \phi [rad]','FontSize', 18);
xlabel('Height from MM surface, z [µm]','FontSize', 18);

yticks([0 pi/8 pi/4]);
yticklabels({'0','\pi/8','\pi/4'});
xlim([0,2.9]);
ylim([0,pi/4]);

ax = gca;
ax.FontSize = 15;
ax.YLabel.FontSize = 20;
ax.XLabel.FontSize = 20;
ax.LineWidth = 1.2;
ax.FontName = 'Arial';

box on


ax = gca; exportgraphics(g,'DistanceMM_Efield.png','Resolution',600)
hold off