MatLab: Plotting Output Electric Field Spectrum

Goal

This code combines two matrices, diagonalizes that matrix, and finds the right and left eigenvectors, right and left eigenvectors are normalized with respect to each other, then it is multiplied by an exponential function that is a function of the Eigen values and plots the output spectrum with respect to omega(w: frequency)

Code

 %%Main V2

%has an input electric field and outputs the transmitted

%define some global variables
n = 100; %step count
a = linspace(1,100,n); %our dummy parameter
p = 1; %p is the index value used of a
N = 10^8; %N is the density of particles needs to be changed
c = 3*10^8; %m/s
k = (2*pi)/c;
%tau = linspace(0,1,n);
tau = 0; % time delay

%the Frequency grid
w = linspace(0,1,n);
wp = linspace(0,1,n);

%%define a diagonal X0 s.t it is zero on the bounds
chi_0 = (1+1i)*(sin(5*pi*w).^2);
%chi_0 get diagonalized in EigVect function

%%define a 3D function same as above
%X_nd non diagonal
for i=1:n
for j=1:n
X_nd(i,j) = (1+1i)*sin(5*pi.*w(i).*wp(j)).^2;
end
end

%%Diagonalize Xij and Normalize EigVectors

[U_r, L, U_L] = EigVect(chi_0,X_nd,a(p),w,n);
%get the right/left eigent vectors and eigen values L
[U_rbar,U_Lbar] = Norm(U_r,U_L);
%normalize the right/left eigen vectors

%%Define the incoming electric field

E_0 = 1; %amplitude
w_0 = 0.5; %choose value from 0.3-0.7 doesn't matter for now
sigma = 0.1;

%This for loop creates a vector of E
for i = 1:n
E(i) = E_0.*sin(pi*w(i))*exp(-((w(i)-w_0)/(sigma*sqrt(2))).^2);
end

%%Compute the transmitted spectrum E(omega,tau;a,N) at the exit of the sample

%exponential function of lambda
Exp_lambda = exp(N*1i*k*diag(L));
%time delay
tau_phase = exp(1i*tau.*w);

UXU = U_rbar*diag(Exp_lambda)*(U_Lbar');

%Output Electric field
E_out = UXU*(tau_phase.*E).';

hold on
plot(w,E, 'b');
plot(w,abs(E_out), 'r'); %for output plot absolute value
legend('Input','Output')
title(['a = ',num2str(a(p)),', N = ',num2str(N),', tau = ',num2str(tau)])

%%Functions
%Normalization of the eigenvectors w.r.t one another
function [A_bar,B_bar] = Norm(A,B)
f = sum(B'*A,1);
sqrt_f = sqrt(f);

A_bar = A./sqrt_f;
B_bar = B./sqrt_f;
end

%input a chi_0 and a non diagonal portion return right and left eigvectors
%and the diagonal eigen values

function [A_bar, L, B_bar] = EigVect(x,y,a,w,n)
X_0 = diag(x);

dw = 1/(n-1);

for i = 1:n
for j = 1:n
X_ij(i,j) = (a*dw*y(i,j)+X_0(i,j))*sqrt(w(i)*w(j));
end
end
%y is X_nd

[A_bar, L, B_bar] = eig(X_ij);
%returns the non-normalized right and left eigenvectors of X_ij
%and the diagonal eigenvalue matrix

end


I am looking for a general code review as well as any suggestions to bring the computation time down. Thank you!

• Do you have any specific concerns, or are you looking for a generic review? May 31, 2021 at 15:23
• Are you in the same class as this person? codereview.stackexchange.com/questions/261316/… — You could learn a lot from reviewing that code! May 31, 2021 at 17:26
• I am looking for a generic code review as well as some shortcuts to speed up the computation time.
– Mate
May 31, 2021 at 18:08