Integration of a multivariate symbolic expression

I am trying to evaluate the imanigary part of a correlation function.

What I've got so far is:

% Defining the constants
cm2au = 1/219474.6305;
fs2au = 1/0.0242;
kB = 3.16681e-6; % au/K

% Temperature
T = 300 ; 
% T = 0.001
beta = 1./(kB*T) ;

%Coupling constant 
lambda = 2 * cm2au ;
% lambda = 600 * cm2au ;
gamma = 53.08 *cm2au ;

% Spectral density
syms omega
J = 2*lambda .* (omega*gamma/(omega.^2 .* gamma^2)) ;
fplot(J)
title('Spectral density')
ylabel('J(\omega)')
xlabel('\omega')

% Correlation function
ReC_pos = J * (1+ 1/(exp(beta*omega)-1)) ;  % ReC(omega) for omega > 0
ReC_neg = ReC_pos * exp(-beta*omega) ;       % ReC(-omega) = ReC(omega) * exp(-beta*omega)

% Imaginary part
syms omega_prime 
term1 = subs(ReC_pos, omega, omega_prime) / (omega - omega_prime);
term2 = subs(ReC_neg, omega, omega_prime) / (omega + omega_prime);
ImC_1 = int(term1,omega_prime,[0 Inf],PrincipalValue=true) ;
ImC_2 = int(term2,omega_prime,[0 Inf],PrincipalValue=true) ;
ImC = ImC_1 +ImC_2 ;

% Plot both positive and negative parts
figure
hold on
fplot(ReC_pos, [0 500*cm2au], 'b', 'DisplayName', 'ReC(\omega), \omega > 0')
fplot(subs(ReC_neg, omega, -omega), [-500*cm2au 0], 'r', 'DisplayName', 'ReC(-\omega), \omega < 0')
hold off
title('Correlation function')
ylabel('C(\omega)')
xlabel('\omega')
legend('Location', 'best')

The real part looks right but I can't seem to evaluate the integral of the imaginary part. It should look like a complex Lorentzian.

[Torsten]

There are several parameters missing to run your code.

[Sy Dat]

I'm terribly sorry. I added the parameters.

[Torsten]

There are still parameters missing (see above). Matlab interprets "gamma" as a reference to the gamma function and thus misses the argument.

Better test the code before posting.

[Sy Dat]

Sorry, now, it should be complete.