I am currently participating in research on transient absorption spectroscopy and four wave mixing. In the experimental design, an extreme ultraviolet (XUV) laser pulse and infrared (IR) laser pulse are sent through a finite gas sample with a certain time delay between them. My current task is to develop a function that, given data on the system and input pulses, provides the transmission spectrum of the XUV pulse electric field after exiting the sample.

The Procedure

The gist of the procedure is to take two matrices describing the sample susceptibility, combine them into a single matrix, diagonalize said matrix, normalize the left and right eigenvector matrices according to each other, apply an exponential function to the eigenvalue matrix, and finally compute the exiting spectrum given an input spectrum.

The Code

function E_out = FiniteSampleTransmissionXUV(w, X_0, X_nd, E_in, a, N, tau)
% Function takes as input:
% > w - A vector that describes the frequency spectrum range
% > X_0 - A vector that describes the susceptibility 
%   without the IR pulse
% > X_nd - A matrix that describes the susceptibility with the
%   IR pulse
% > E_in - A vector that describes the frequency spectrum of the 
%   incoming XUV pulse electric field
% > a - A constant that describes the intensity of the IR pulse
% > N - A constant that describes the optical depth of the
%   sample
% > tau - A constant that describes the time delay between the 
%   IR pulse and the XUV pulse
% Function provides output:
% > A vector that describes the frequency spectrum of the
%   XUV pulse electric field at the end of the sample

    % determines number of frequency steps from input frequency range
    n = size(w);    % constant
    % determines frequency step size
    delta_w = 1 / (n(2) - 1);   % constant
    % create matrix sqrt(w_i*w_j)
    sqrt_w = w.^(1/2).' * w.^(1/2);     % nxn matrix

    % combine X_0 and X_nd into total suscptibility matrix
    X_ij = (a * delta_w * sqrt_w .* X_nd) + ...
            diag(diag(sqrt_w).' .* (X_0));      % nxn matrix
    % diagonalize X_ij where sum(U_R_i^2) = 1
    [U_R, LAMBDA, U_L] = eig(X_ij);     % nxn matrices
    % attain the function that scales U_L'*U_R
    F = sum(U_L' * U_R, 1);     % row vector
    sqrt_F = F.^2;
    % scale U_R and U_L so that U_L'*U_R = 1
    U_R_bar = U_R ./ sqrt_F;    % nxn matrix
    U_L_bar = U_L ./ sqrt_F;    % nxn matrix

    % apply exponential transform to eigenvalues      % diagonal nxn matrix
    exp_LAMBDA = diag(exp(1i * (2*pi*N / (3*10^8)) * diag(LAMBDA)));
    % create phase shift due to the time delay
    tau_phase = exp(1i * w * tau);  % row vector
    % recombine transformed susceptibility matrix
    ULAMBDAU = U_R_bar * exp_LAMBDA * U_L_bar';     % nxn matrix
    % apply effects of finite propagation and pulse
    % time delay to input electric field spectrum
    E_out = ULAMBDAU * (tau_phase .* E_in).';   % vector


The following is a script to test and demonstrate the function.

n = 100;    % number of frequency steps
w = linspace(0,1,n);    % linearly spaced 1D array

X_0 = (1 + 1i) * (sin(5*pi*w)).^2;

X_nd = (1 + 1i) * ((sin(1*pi*w')) * sin(1*pi*w)).^2;

E_in = GaussianSpectrum(w, 1, 0.5, 0.1);

a = 1;
N = 10^8;
tau = 0;
E_out = FiniteSampleTransmissionXUV(w, X_0, X_nd, E_in, a, N, tau);

plot(w, E_in,'b')
hold on
plot(w, abs(E_out),'r')
hold off

function E_w_x0 = GaussianSpectrum(w, E_0, w_0, sigma)
% creates a quasi-gaussian spectrum defined by the input
% frequency range w, amplitude E_0, center frequency w_0,
% and spectral width sigma

    E_w_x0 = E_0 * sin(pi*w) .* ...
             (exp(-((w - w_0) / (2^(1/2) * sigma)).^2));


In running the test, we see the frequency spectrum of the XUV pulse before entering the sample in blue with the exiting spectrum in red. One sees that frequencies of the XUV pulse have been absorbed, while new frequencies have been generated.

enter image description here

Please let me know what violates best practices and how the function could be refactored to be more robust and clean. Thank you!


1 Answer 1


Very nice code, it’s hard to find anything to comment on it.

3*10^8 is better written as 3e8 (IEEE float format), then MATLAB doesn’t need to do any computations, though I doubt this would ever be a significant portion of the computation time. It is up to you what you find most readable.

In your test code, 2^(1/2) is more commonly written as sqrt(2). Again, what is more readable is up to you.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.