# What is wrong in my script?

2 views (last 30 days)
Athira T Das on 14 Jul 2022
Answered: Walter Roberson on 25 Jul 2022
clc; clear all; close all;
syms m1 m2 s1 s2 j1 j2
I = 0;
lambda = 1060*10^-9;
M=0
z=linspace(0.00001,8000)
wo = 0.02;
C = 10^(-7);
k=2*pi/lambda;
po=(0.545*C^2*k^2*z).^(-(3/5));
b=0.1;
T1=0;
T2=0;
T3=0;
T4=0;
T5=0;
T6=0;
delta= ((1i*k)./(2*z))+(1./(wo.^2))+(1./(po.^2));
delta_star= subs(delta, 1i, -1i);
etta = delta_star - (1./(delta.*po.^4));
X=((k./(2.*z)).^2).*((1./(2.*1i.*sqrt(delta))).^M).*(1./(16.*delta.*etta));
alpha = (1i.*k./(2.*z)).*(1./(delta.*po.^2)-1);
beta_p = (b./(2.*wo)).*(1./(delta.*po.^2)+1);
beta_n = (b./(2.*wo)).*(1./(delta.*po.^2)-1);
A = (alpha.^2./etta)-((k.^2)-(4.*z.^2.*delta));
B1_p = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta);
B1_n = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta);
B2_p = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta));
B2_n = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta));
C1_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
C1_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
e1=(exp(C1_p).*hermiteH(m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(m2+j2,((-1i.*beta_p)./sqrt(etta)))));
e2=(exp(C1_p).*hermiteH(M-m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M-m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(M-m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M-m2+j2,((-1i.*beta_p)./sqrt(etta)))));
O = zeros(size(z));
for m1 =0:M
T2=0;
for m2=0:M
T3=0;
for s1=0:m1/2
T4=0;
for j1=0:m1-(2*s1)
T5=0;
for s2=0:(M-m1)/2
T6=0;
for j2=0:(M-m1-(2*s2))
T6=T6+(factorial(M-m1-(2.*s2))/(factorial(j2)*factorial(M-m1-(2.*s2)-j2))).*((1./(2.*1i.*sqrt(etta))).^(M-m2+j2)).*((1./(po.^2)).^j2).*((b./2.*wo).^(M-m1-2*s2-j2)).*e2;
end
T5=T5+(((factorial(M-m1).*(-1).^s2)./(factorial(s2).*factorial(M-m1-2.*s2))).*(((2.*1i)./(delta.^(0.5))).^(M-m1-2.*s2))).*T6;
end
if (m1-(2*s2))<0
break
end
T4=T4+(((factorial(m1-(2.*s1)))/(factorial(j1)*factorial(m1-(2.*s1)-j1))).*((1./(2.*1i.*sqrt(etta))).^(m2+j1)).*((1./(po.^2)).^j1).*((b./2.*wo).^(m1-2.*s1-j1))).*e1.*T5;
end
T3=T3+(((factorial(m1).*(-1).^s1)./(factorial(s1).*factorial(m1-(2.*s1)))).*(((2.*1i)./(sqrt(delta))).^(m1-2.*s1))).*T4;
end
T2=T2+nchoosek(M,m2)*((-1i)^(M-m2))*T3;
end
T1=T1+nchoosek(M,m1)*((1i)^(M-m1))*T2;
end
I0 = -real(X.*T1)
I_numerical = double(vpa(I0))
I_o=abs(I_numerical)
writematrix(I_o,"Normalized_intensity_M=0_b0.1_if.xlsx")
writematrix(z,"Propagation_distance.xlsx")
Athira T Das on 15 Jul 2022
@TorstenYes it is from reliable source

Walter Roberson on 25 Jul 2022
Your code defines e1 and e2 in terms of symbolic variables m2 and j2 . Later it does for loops that assign values to m2 and j2. However, assigning a value to m2 and j2 at the MATLAB level does not affect the symbolic variables by the same name. You need to subs() the numeric values in.
format long g
syms m2 s1 s2 j1 j2
I = 0;
lambda = 1060*10^-9;
M=0;
z=linspace(0.00001,8000);
wo = 0.02;
C = 10^(-7);
k=2*pi/lambda;
po=(0.545*C^2*k^2*z).^(-(3/5));
b=0.1;
T1=0;
T2=0;
T3=0;
T4=0;
T5=0;
T6=0;
delta= ((1i*k)./(2*z))+(1./(wo.^2))+(1./(po.^2));
delta_star= subs(delta, 1i, -1i);
etta = delta_star - (1./(delta.*po.^4));
X=((k./(2.*z)).^2).*((1./(2.*1i.*sqrt(delta))).^M).*(1./(16.*delta.*etta));
alpha = (1i.*k./(2.*z)).*(1./(delta.*po.^2)-1);
beta_p = (b./(2.*wo)).*(1./(delta.*po.^2)+1);
beta_n = (b./(2.*wo)).*(1./(delta.*po.^2)-1);
A = (alpha.^2./etta)-((k.^2)-(4.*z.^2.*delta));
B1_p = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta);
B1_n = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta);
B2_p = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta));
B2_n = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta));
C1_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
C1_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
e1=(exp(C1_p).*hermiteH(m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(m2+j2,((-1i.*beta_p)./sqrt(etta)))));
e2=(exp(C1_p).*hermiteH(M-m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M-m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(M-m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M-m2+j2,((-1i.*beta_p)./sqrt(etta)))));
O = zeros(size(z));
for m1 =0:M
T2=0;
for m2=0:M
T3=0;
for s1=0:m1/2
T4=0;
for j1=0:m1-(2*s1)
T5=0;
for s2=0:(M-m1)/2
T6=0;
for j2=0:(M-m1-(2*s2))
T6=T6+(factorial(M-m1-(2.*s2))/(factorial(j2)*factorial(M-m1-(2.*s2)-j2))).*((1./(2.*1i.*sqrt(etta))).^(M-m2+j2)).*((1./(po.^2)).^j2).*((b./2.*wo).^(M-m1-2*s2-j2)).*subs(e2);
end
T5=T5+(((factorial(M-m1).*(-1).^s2)./(factorial(s2).*factorial(M-m1-2.*s2))).*(((2.*1i)./(delta.^(0.5))).^(M-m1-2.*s2))).*T6;
end
if (m1-(2*s2))<0
break
end
T4=T4+(((factorial(m1-(2.*s1)))/(factorial(j1)*factorial(m1-(2.*s1)-j1))).*((1./(2.*1i.*sqrt(etta))).^(m2+j1)).*((1./(po.^2)).^j1).*((b./2.*wo).^(m1-2.*s1-j1))).*subs(e1).*T5;
end
T3=T3+(((factorial(m1).*(-1).^s1)./(factorial(s1).*factorial(m1-(2.*s1)))).*(((2.*1i)./(sqrt(delta))).^(m1-2.*s1))).*T4;
end
T2=T2+nchoosek(M,m2)*((-1i)^(M-m2))*T3;
end
T1=T1+nchoosek(M,m1)*((1i)^(M-m1))*T2;
end
I0 = -real(X.*T1);
I_numerical = double(vpa(I0));
I_o=abs(I_numerical);
writematrix(I_o,"Normalized_intensity_M=0_b0.1_if.xlsx")
writematrix(z,"Propagation_distance.xlsx")

### Categories

Find more on Dialog Boxes in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by