%% 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;
TrapDepth       = ((8 * a * Power) / (pi * waist_y * waist_z)) / (BoltzmannConstant * 1E-6); % in µK
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 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(3);
set(gcf,'Position',[100 100 950 750])
semilogy(LatticeSpacing * 1E6, RecoilEnergy/PlanckConstant, LineWidth=2.0, DisplayName=['\bf Max = ' num2str(round(max(RecoilEnergy / PlanckConstant),1)) ' Hz; Min = ' num2str(round(min(RecoilEnergy / PlanckConstant),1)) ' Hz'])
xlim([0.5 21]);
xlabel('Lattice spacing (µm)', FontSize=16)
ylabel('Recoil Energy (Hz)', FontSize=16)
title('\bf Scaling of Recoil Energy - All energy scales in an optical lattice are naturally parametrized by the lattice recoil energy', FontSize=12)
grid on
legend(FontSize=12)

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

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

figure(4);
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