diff --git a/EstimatesForAccordionLattice.m b/EstimatesForAccordionLattice.m
index 34f681c..19cac75 100644
--- a/EstimatesForAccordionLattice.m
+++ b/EstimatesForAccordionLattice.m
@@ -43,6 +43,7 @@ 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));
@@ -74,22 +75,36 @@ title(['\bf Upper bound = ' num2str(round(max(nu_z * 1E-3),2)) ' kHz ; \bf Lower
 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);
-TrapDepth       = ((8 * a .* Powers) ./ (pi * waist_y * waist_z)) / (BoltzmannConstant * 1E-6); % in µK
+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, TrapDepth, LineWidth=2.0)
-xlim([0.0 5.5]);
+plot(Powers, TrapDepthsToPlot, LineWidth=2.0)
+xlim([0.0 5.25]);
 xlabel('Powers (W)', FontSize=16)
-ylabel('Trap depth (µK)', FontSize=16)
-title(['\bf Upper bound = ' num2str(round(max(TrapDepth),2)) ' µK ; \bf Lower bound = ' num2str(round(min(TrapDepth),2)) ' µK'], 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 Recoil Energy - All energy scales in an optical lattice are naturally parametrized by the lattice recoil energy
+%% 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);
@@ -124,5 +139,25 @@ 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
\ No newline at end of file