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);


