Branch function in a sum in Matlab-CodePudding

I want to calculate the amplitude and phase for a single slit diffraction with the initial condition U(x,0)=1, |x|<a and U(x,0)=0, |x|>a, a=2. My code in Matlab is the following

clear all; clc; format long;  22-BPM based on Fresnel integral
Nx=1024; L=100; x=(L/Nx)*(-Nx/2:Nx/2-1); dx=x(2)-x(1); 
tic  

a=2;
syms u(x);
u(x)=piecewise(-a<x<a, 1, x>a, 0, x<-a, 0);
fplot(u);

S=100; zmax=10; z=linspace(0.1,zmax,S);
for m=1:S
    z1=z(m);
    for tr=1:Nx
      field(tr,m)=exp(1i*pi/4)/(2*1i*sqrt(pi*z1))*...
      sum(u.*exp(1i*(x(tr)-x).^2/(4*z1)))*dx;
    end;
end; 

[X,Z]=meshgrid(x,z);
figure(1); surf(X',Z',abs(field).^2); shading interp; view([0 90]); 
axis tight; axis square;
figure(2); surf(X',Z',angle(field)); shading interp; view([0 90]); 
axis tight; axis square;
figure(30); plot(x,abs(field(:,S)),'-b*'); hold on;
figure(31); plot(x,angle(field(:,S)),'-b*'); hold on;

I_onaxis=abs(field(Nx/2 1,:)); % on-axis Intensity of the field
figure(32); plot(z,I_onaxis,'-b*'); hold on;

toc

The error is in line 16, where the sum is.

CodePudding user response：

Don't mix symbolic and numeric variables. Replace the u(x) function with u = abs(x) < a and pre-allocate the matrix field for speed.

clear, clc, format long;  22-BPM based on Fresnel integral
Nx = 1024; 
L  = 100; 
x  = (L/Nx)*(-Nx/2:Nx/2-1); 
dx = x(2)-x(1); 

tic  
a = 2;
u = abs(x) < a;

S = 100; zmax = 10; z = linspace(0.1,zmax,S);
field = zeros(Nx,S);
for m = 1:S
    z1 = z(m);
    for tr = 1:Nx
      field(tr,m) = exp(1i*pi/4)/(2*1i*sqrt(pi*z1))*...
      sum( u .* exp(1i*(x(tr)-x).^2/(4*z1)) ) * dx;
    end
end

[X,Z] = meshgrid(x,z);
figure(1); surf(X',Z',abs(field).^2); shading interp; view([0 90]); 
axis tight; axis square;
figure(2); surf(X',Z',angle(field)); shading interp; view([0 90]); 
axis tight; axis square;
figure(30); plot(x,abs(field(:,S)),'-b*'); hold on;
figure(31); plot(x,angle(field(:,S)),'-b*'); hold on;

I_onaxis=abs(field(Nx/2 1,:)); % on-axis Intensity of the field
figure(32); plot(z,I_onaxis,'-b*'); hold on;

toc