clear; clc; %% Select the output flag for the CCD camera= 1: NIRCCD; 2: VISCCD CCD_flag=1; if CCD_flag==1 Nw=256; % Number of pixels of the NIR CCD elseif CCD_flag==2 Nw=1340; % Number of pixels of the Visible light CCD end %% Output file name for simulation result SI_data_output='Dimer_2DCS.txt'; LO_spectrum_output='LOfield.txt' ; %% Parameter assignments of laser source w_0=229. ; % [THz] Central frequency of the exciting pulse t_d=70. ; % e-1 pulse width [fs] delw=440./t_d ; % e-1 spectral width [THz] of the laser pulse %% Model parameter definition and assignment for the Heterodimer % For more information, please see arXiv:1307.5855v1. w_a=226. ; % [THz] Energy level at site a w_b=232. ; % [THz] Energy level at site b mu_a=0.8 ; % [D] Transition moment at site a mu_b=1.6 ; % [D] Transition moment at site b delE=0.5*(w_b-w_a) ; J=3. ; % J/h [THz]: Dimer's coupling strength avew=0.5*(w_a+w_b) ; theta2=atan(J/delE) ; theta=0.5*theta2 ; w_1=avew-delE/cos(theta2); % Energy of the excited level 1 of the dimer w_2=avew+delE/cos(theta2); % Energy of the excited level 2 of the dimer w_f1=w_2 ;% Energy from the excited level 1 to the final level of the dimer w_f2=w_1 ;% Energy from the excited level 2 to the final level of the dimer omega_1 = w_1*(pi/500.) ; omega_2 = w_2*(pi/500.) ; omega_f1=omega_2 ; omega_f2=omega_1 ; mu_1g=mu_a*cos(theta)+mu_b*sin(theta) ; % Transition moment from g to 1 mu_2g=-mu_a*sin(theta)+mu_b*cos(theta) ; % Transition moment from g to 2 mu_f1=mu_a*sin(theta)+mu_b*cos(theta) ; % Transition moment from 1 to f mu_f2=mu_a*cos(theta)-mu_b*sin(theta) ; % Transition moment from 2 to f gam=3.0 ; % Damping constant [THz] of the excited levels decay=gam/1000. ; % [fs-1] % Terms of the echo field term_d1=2.*mu_1g^4 ; term_d2=2.*mu_1g^2*mu_2g^2 ; term_c1=2.*mu_2g^4 ; term_c2=term_d2 ; term_2q1=-mu_f1*mu_1g*(mu_f1*mu_1g+mu_f2*mu_2g) ; term_2q2=-mu_f2*mu_2g*(mu_f2*mu_2g+mu_f1*mu_1g) ; %% Some scalings % Lambda: wavelength [in nm] and the corresponding freq [in THz] % c=299.792458 [nm/fs], freq[THz]=1000*c[nm/fs]/lambda[nm]. % omega[rad/fs] = 2000*pi*freq[THz] c = 299.792458; % [nm/fs] chng_f2nm = @(x) 1000*c./x; % Change from [THz] to [nm] %% Prepare the pulse spectrum of the exciting pulse % freq =frequency [THz], lambda=wavelength [nm] delf=4.0*delw/Nw ; freq=(w_0+2.*delw)-(0:Nw-1)*delf ; % frequency mesh [THz] for SI % omega = freq.*(pi/500.) ; % Angular frequency mesh [rad/fs] for SI lambda = chng_f2nm(freq); % wavelength mesh [nm] %% Simulation of Echo field E(w_ab, T=0, w_em) E_echo=zeros(Nw,Nw); for i=1:Nw fq=-freq(i) ; arg1=term_d1*1./(fq+w_1+sqrt(-1)*gam); arg2=term_d2*1./(fq+w_2+sqrt(-1)*gam); arg3=term_c1*1./(fq+w_2+sqrt(-1)*gam); arg4=term_c2*1./(fq+w_1+sqrt(-1)*gam); arg5=term_2q1*1./(fq+w_1+sqrt(-1)*gam); arg6=term_2q2*1./(fq+w_2+sqrt(-1)*gam); for j=1:Nw G_d(j)=1./(freq(j)-w_1+sqrt(-1)*gam) ; G_c(j)=1./(freq(j)-w_2+sqrt(-1)*gam) ; G_2q1(j)=1./(freq(j)-w_f1+sqrt(-1)*gam) ; G_2q2(j)=1./(freq(j)-w_f2+sqrt(-1)*gam) ; E_echo(j,i)=(arg1+arg2)*G_d(j)+(arg3+arg4)*G_c(j)+arg5*G_2q1(j)+arg6*G_2q2(j) ; end end Echo=zeros(Nw,Nw); for i=1:Nw Echo(:,i)= abs(E_echo(:,i))./max(abs(E_echo(:,1))) ; end %% Plotting the 2D data f_min=min(freq); f_max=max(freq); figure(1) contourf(freq,-freq,Echo,20); set(gca,'YDir','normal'); axis([f_min f_max -f_max -f_min],'square'); set(gca,'FontSize',16); title('\fontsize{18}\bfDirect Calc. Abs(Echo)'); xlabel('\fontsize{18}Emission Frequency (THz)'); ylabel('\fontsize{18}Excitation Frequency (THz)'); colormap(jet) ; cb=colorbar; set(cb,'FontSize',16);