%% Physical constants
PlanckConstant               = 6.62607015E-34;
PlanckConstantReduced        = 6.62607015E-34/(2*pi);
FineStructureConstant        = 7.2973525698E-3;
ElectronMass                 = 9.10938291E-31;
GravitationalConstant        = 6.67384E-11;
ProtonMass                   = 1.672621777E-27;
AtomicMassUnit               = 1.660539066E-27; 
BohrRadius                   = 5.2917721067E-11; 
BohrMagneton                 = 9.274009994E-24; 
BoltzmannConstant            = 1.38064852E-23;
StandardGravityAcceleration  = 9.80665;
SpeedOfLight                 = 299792458;
StefanBoltzmannConstant      = 5.670373E-8;
ElectronCharge               = 1.602176634E-19;
VacuumPermeability           = 1.25663706212E-6; 
DielectricConstant           = 8.8541878128E-12;
ElectronGyromagneticFactor   = -2.00231930436153;
AvogadroConstant             = 6.02214076E23;
ZeroKelvin                   = 273.15;
GravitationalAcceleration    = 9.80553;
VacuumPermittivity           = 1 / (SpeedOfLight^2 * VacuumPermeability);
HartreeEnergy                = ElectronCharge^2 / (4 * pi * VacuumPermittivity * BohrRadius);
AtomicUnitOfPolarizability   = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3

% Dy specific constants
Dy164Mass                    = 163.929174751*1.660539066E-27;
Dy164IsotopicAbundance       = 0.2826;
DyMagneticMoment             = 9.93*9.274009994E-24;

%% Lattice spacing 

Wavelength     = 532e-9;
theta          = linspace(1.5, 18.0, 100);
LatticeSpacing = Wavelength ./ (2.*sin((theta*pi/180)/2));

figure(1);
set(gcf,'Position',[100 100 950 750])
plot(theta, LatticeSpacing * 1E6, LineWidth=2.0)
xlim([0 19]);
ylim([0.5 21]);
xlabel('Angle (deg)', FontSize=16)
ylabel('Lattice spacing (µm)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(LatticeSpacing * 1E6),1)) ' µm ; \bf Lower bound = ' num2str(round(min(LatticeSpacing * 1E6),1)) ' µm'], FontSize=16)
grid on

%% Scaling of vertical trap frequency with lattice spacing
Wavelength      = 532e-9;
a               = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
Power           = 5;
waist_y         = 250E-6;
waist_z         = 50E-6;
thetas          = linspace(1.5, 18.0, 100);
LatticeSpacings = zeros(1, length(thetas));
Omega_z         = zeros(1, length(thetas));

for idx = 1:length(thetas)
    theta                = 0.5 * thetas(idx) .* pi/180;
    LatticeSpacings(idx) = Wavelength ./ (2.*sin(theta));
    Omega_z(idx)         = sqrt(((16 * a * Power) / (pi * Dy164Mass * waist_y * waist_z)) * ...
                           ((2 * (cos(theta)/waist_z)^2) + ((Wavelength * sin(theta)/pi)^2 * ...
                           ((1/waist_y^4) + (1/waist_z^4))) + (pi / LatticeSpacings(idx))^2));
    
end

nu_z         = Omega_z ./ (2*pi);

figure(2);
set(gcf,'Position',[100 100 950 750])
plot(LatticeSpacings * 1E6, nu_z * 1E-3, LineWidth=2.0)
xlim([0.5 21]);
xlabel('Lattice spacing (µm)', FontSize=16)
ylabel('Trap frequency (kHz)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(nu_z * 1E-3),2)) ' kHz ; \bf Lower bound = ' num2str(round(min(nu_z * 1E-3),2)) ' kHz'], FontSize=16)
grid on

%% Scaling of trap depth with power
a                               = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
waist_y                         = 250E-6;
waist_z                         = 50E-6;
Powers                          = linspace(0.1, 5, 100);
TrapDepths                      = ((8 * a .* Powers) ./ (pi * waist_y * waist_z));
TrapDepthsInHz                  = TrapDepths ./ PlanckConstant; % in Hz
TrapDepthsInmicroK              = TrapDepths ./ BoltzmannConstant; % in µK
TwoPhotonRecoilEnergy           = (2*PlanckConstantReduced*2*pi/Wavelength)^2 / (2 * Dy164Mass);
TrapDepthsInUnitsOfRecoilEnergy =  TrapDepths ./ TwoPhotonRecoilEnergy;

% TrapDepthsToPlot = TrapDepthsInHz * 1E-3;
% units            = ' kHz';

% TrapDepthsToPlot = TrapDepthsInmicroK * 1E6;
% units            = ' µK';

TrapDepthsToPlot = TrapDepthsInUnitsOfRecoilEnergy;
units            = ' E_r';


figure(3);
set(gcf,'Position',[100 100 950 750])
plot(Powers, TrapDepthsToPlot, LineWidth=2.0)
xlim([0.0 5.25]);
xlabel('Powers (W)', FontSize=16)
ylabel(['Trap depth (' units ' )'], FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(TrapDepthsToPlot),2)) units '; \bf Lower bound = ' num2str(round(min(TrapDepthsToPlot),2)) units], FontSize=16)
grid on

%% Scaling of the lattice recoil Energy - All energy scales in an optical lattice are naturally parametrized by the lattice recoil energy

LatticeSpacing = linspace(2E-6, 20E-6, 100);
RecoilEnergy   = PlanckConstant^2 ./ (8 .* Dy164Mass .* LatticeSpacing.^2);

figure(4);
set(gcf,'Position',[100 100 950 750])
semilogy(LatticeSpacing * 1E6, RecoilEnergy/PlanckConstant, LineWidth=2.0)
xlim([0.5 21]);
xlabel('Lattice spacing (µm)', FontSize=16)
ylabel('Recoil Energy (Hz)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(RecoilEnergy / PlanckConstant),1)) ' Hz; Lower bound = ' num2str(round(min(RecoilEnergy / PlanckConstant),1)) ' Hz'], FontSize=16)
grid on

%% Interference pattern spacing in ToF - de Broglie wavelength associated with the relative motion of atoms

ExpansionTime             = linspace(1E-3, 20.0E-3, 100);

figure(5);
set(gcf,'Position',[100 100 950 750])
labels = [];

for ls = [2E-6:2E-6:5E-6 6E-6:6E-6:20E-6]
    InteferencePatternSpacing = (PlanckConstant .* ExpansionTime) ./ (Dy164Mass * ls);
    plot(ExpansionTime*1E3, InteferencePatternSpacing* 1E6, LineWidth=2.0, DisplayName=['\bf Lattice spacing = ' num2str(round(max(ls * 1E6),1)) ' µm'])
    hold on
end
xlim([0 22]);
xlabel('Free expansion time (milliseconds)', FontSize=16)
ylabel('Interference pattern period (µm)', FontSize=16)
title('\bf Interference of condensates - Fringe period is the de Broglie wavelength associated with the relative motion of atoms', FontSize=12)
legend(labels, 'Location','NorthWest', FontSize=12);
grid on
legend show

%% Scaling of frequency of oscillation in the first order in Kapitza-Dirac scattering 

a                         = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
waist_y                   = 250E-6;
waist_z                   = 50E-6;
Powers                    = linspace(0.001, 0.4, 100);
TrapDepths                = ((8 * a .* Powers) ./ (pi * waist_y * waist_z));
TwoPhotonRecoilEnergy     = (2*PlanckConstantReduced*2*pi/Wavelength)^2 / (2 * Dy164Mass);
RabiOscillationFrequency  = (1/PlanckConstantReduced) .* (sqrt(TrapDepths.^2/2 + TwoPhotonRecoilEnergy^2));

TrapDepthsInHz                  = TrapDepths./ PlanckConstant;
TwoPhotonRecoilEnergyInHz       = TwoPhotonRecoilEnergy / PlanckConstant;
TrapDepthsInUnitsOfRecoilEnergy = TrapDepthsInHz ./ TwoPhotonRecoilEnergyInHz;

figure(6);
set(gcf,'Position',[100 100 950 750])
plot(TrapDepthsInUnitsOfRecoilEnergy, RabiOscillationFrequency .* 1E-3, LineWidth=2.0)
xlim([0 4]);
xlabel('Trap depths (E_r)', FontSize=16)
ylabel('Rabi oscillation frequency (kHz)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(RabiOscillationFrequency .* 1E-3),1)) ' kHz; Lower bound = ' num2str(round(min(RabiOscillationFrequency .* 1E-3),1)) ' kHz'], FontSize=16)
grid on

%% Rabi oscillations in the first order in Kapitza-Dirac scattering 

a                               = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
waist_y                         = 250E-6;
waist_z                         = 50E-6;
Power                           = 0.2;
TrapDepth                       = ((8 * a .* Power) ./ (pi * waist_y * waist_z));
TwoPhotonRecoilEnergy           = (2*PlanckConstantReduced*2*pi/Wavelength)^2 / (2 * Dy164Mass);
TrapDepthsInUnitsOfRecoilEnergy = TrapDepth ./ TwoPhotonRecoilEnergy;
RabiOscillationFrequency        = (1/PlanckConstantReduced) .* (sqrt(TrapDepth.^2/2 + TwoPhotonRecoilEnergy^2));

alpha = TwoPhotonRecoilEnergy / PlanckConstantReduced;
beta  = TrapDepth / PlanckConstantReduced;
C     = beta^2 / ((2*beta^2)  + (4*alpha^2));

PulseDurations                  = linspace(1E-6, 150E-6, 1000);
PopulationInFirstOrders         = C .* (sin(0.5 .* PulseDurations .* (RabiOscillationFrequency)));

figure(6);
set(gcf,'Position',[100 100 950 750])
plot(PulseDurations .* 1E6, PopulationInFirstOrders, LineWidth=2.0, DisplayName=['\bf Power = ' num2str(Power) ' W / Trap depth = ' num2str(round(TrapDepthsInUnitsOfRecoilEnergy, 1)) ' E_r'])
xlabel('Pulse duration (µs)', FontSize=16)
ylabel('Fraction of atoms in first order', FontSize=16)
title('\bf Expected Rabi oscillation', FontSize=16)
grid on
legend(FontSize=16)