% Experiment 4 of ECE 474
% This prog is based on my derivation, which is a modified version of Eq. (1.3.17) of Salehi.
close all;clear;clc;clf reset
n1 = 1.50;n2 = 1.47;delta = (n1 - n2)/n1;a = 25e-6;
z = 0:0.1:50;z = z*a;x0 = 0.2*a*1;y0 = 0.2*a*1;tetx0 = 0.08736;tety0 = 0.0873;
zf = 0:0.02:50;zf = zf*a*100;
x_z = x0*cos(sqrt(2*delta)*z/a) + a*tetx0*sin(sqrt(2*delta)*z/a)/sqrt(2*delta);
y_z = y0*cos(sqrt(2*delta)*z/a) + a*tety0*sin(sqrt(2*delta)*z/a)/sqrt(2*delta);

set(gcf,'Renderer','Zbuffer');set(gcf,'Color',[1 1 1]);
plot3(x_z,y_z,z,'-r','LineWidth',4);set(gca,'FontSize',13);%colormap colorcube;%view([-21.5 60.0])
%title('3D Plot of Characteristic Equation for HE and EH Modes','FontSize',12,'FontWeight','bold');
grid on
axis([-a a -a a min(z) max(z)]);hold on
plot3(a*cos(sqrt(2*delta)*zf/a),a*sin(sqrt(2*delta)*zf/a),zf/100,'--g','LineWidth',0.1);
axis square
xlabel('\itx\rm\bf axis','FontSize',12,'FontWeight','bold');
ylabel('\ity\rm\bf axis','FontSize',12,'FontWeight','bold');
zlabel('\itz\rm\bf axis','FontSize',12,'FontWeight','bold');
n1str = ['\itn\rm\bf_1 = ' num2str(n1,'%3.2f')];n2str = ['\itn\rm\bf_2 = ' num2str(n2,'%3.2f')];
astr = ['\ita\rm\bf = ' num2str(a/1e-6,'%3.0f') ' \mum'];
x0str = ['\itx\rm\bf_0\it_n = \rm\bf ' num2str(x0/a)];y0str = ['\ity\rm\bf_0\it_n = \rm\bf ' num2str(y0/a)];
tetx0str = ['\it\theta\rm\bf_0\it_x = \rm\bf ' num2str(tetx0*180/pi,'%3.0f') '^0'];
tety0str = ['\it\theta\rm\bf_0\it_y = \rm\bf ' num2str(tety0*180/pi,'%3.0f') '^0'];
legc = {n1str;' ';n2str;' ';astr;' ';x0str;' ';y0str;' ';tetx0str;' ';tety0str};
text(max(x_z)*0.8,max(y_z),max(z)*0.8,legc,'FontSize',14,'FontWeight','bold','BackgroundColor','white');
arctan = atan(2*(x0*tetx0 +y0*tety0)/ ...
   (sqrt(2*delta)/a*(x0^2 + y0^2) - a/sqrt(2*delta)*(tetx0^2 + tety0^2)));
zp = 0.5*a/sqrt(2*delta)*arctan;
x_zp = x0*cos(sqrt(2*delta)*zp/a) + a*tetx0*sin(sqrt(2*delta)*zp/a)/sqrt(2*delta);
y_zp = y0*cos(sqrt(2*delta)*zp/a) + a*tety0*sin(sqrt(2*delta)*zp/a)/sqrt(2*delta);
r_min = sqrt(x_zp^2 + y_zp^2)
arctan = arctan + pi;zp = 0.5*a/sqrt(2*delta)*arctan;
x_zp = x0*cos(sqrt(2*delta)*zp/a) + a*tetx0*sin(sqrt(2*delta)*zp/a)/sqrt(2*delta);
y_zp = y0*cos(sqrt(2*delta)*zp/a) + a*tety0*sin(sqrt(2*delta)*zp/a)/sqrt(2*delta);
r_max = sqrt(x_zp^2 + y_zp^2)
