%constants
e0 = 8.85e-12;
De = 17;
K2 = 7.4e-12;
d = 20;
zc = 2.9;


%simplified
value = sqrt(e0*De/(K2));
E = (0:0.00001:0.5);

%%functions
%k2 = k2_out(E);
%k1 = -(pi./4 + k2);
%f = (pi/4) + k1.*exp(-value.*E.*zc) + k2.*exp(value.*E.*zc);
k2 = arrayfun(@k2_out,E);
k1 = -(pi./4 + k2);
f = arrayfun(@(a,b,c) (pi/4) + a.*exp(-value.*b.*zc) + c.*exp(value.*b.*zc),k1,E,k2);

%Plotting
g = figure();
hold on
plot(E,f,'Color','black','LineWidth',1.5)
ylabel('Twist angle, \phi [rad]','FontSize', 18);
xlabel('Electric field, E [V/µm]','FontSize', 18);
%ylabel('Transmittance','FontSize', 18);
%xlabel('Frequency [THz]','FontSize', 18);
yticks([0 pi/8 pi/4]);
yticklabels({'0','\pi/8','\pi/4'});
%xlim([0,0.5]);
ylim([0,0.8]);
yline(pi/6,'--','Color','black','LineWidth',1.2);
xline(0.1348,'--','Color','black','LineWidth',1.2);


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,'Efield_Twist_Solution.png','Resolution',600)
hold off

%%Functions

function k = k2_out(ElectricField)
    e0 = 8.85e-12;
    De = 17;
    K2 = 7.4e-12;
    value = sqrt(e0*De/(K2));
    zc = 2.9;
    k = (-pi.*exp(-value.*ElectricField.*zc))/(8.*(cosh(value.*ElectricField.*zc)));
end