Added plots, fixed python files

Python/trap_forces_axial.py

@@ -14,13 +14,14 @@ sns.set()
 
 # Intensity profile graphics
 rho = np.linspace(-r_max, r_max, 500)
-fig1 = plt.figure(1, figsize=(10, 6))
+plt.figure(1, figsize=(10, 6))
 I = gauss(rho)
 I0 = np.max(I)
 plt.plot(rho, I / I0)
 plt.fill_between(rho, I / I0, 0, alpha=0.3)
-plt.xlabel('r, m', fontsize=18)
+plt.title("Forces along Y axis")
-plt.ylabel('I(r)', fontsize=18)
+plt.xlabel('$r, m$', fontsize=18)
+plt.ylabel('$I(r)$', fontsize=18)
 
 # Integration
 def q_res_g(z, func):
@@ -47,11 +48,12 @@ axial = axial_g + axial_s
 z = -z[::-1]
 
 # Graphics
-fig2 = plt.figure(2, figsize=(10, 6))
+plt.figure(2, figsize=(10, 6))
 plt.plot(z, axial_g, '-.', lw=1, label='$F_{g}$')
 plt.plot(z, axial_s, '--', lw=1, label='$F_{s}$')
 plt.plot(z, axial, lw=1, label='$F_{t}$')
-plt.xlabel('r, m', fontsize=18)
+plt.xlabel('$r, m$', fontsize=18)
-plt.ylabel('F, N', fontsize=18)
+plt.ylabel('$F, N$', fontsize=18)
 plt.legend(fontsize=18)
+
 plt.show()

Python/trap_forces_efficiencies.py

@@ -14,12 +14,11 @@ t_deg = t * 180 / np.pi
 pol = np.pi / 4
 
 plt.figure(figsize=(13, 8))
-plt.plot(t_deg, q_s(t, pol), '--', label='Q_s')
+plt.plot(t_deg, q_s(t, pol), '--', label='$Q_s$')
-plt.plot(t_deg, -q_g(t, pol), '-.', label='Q_g')
+plt.plot(t_deg, -q_g(t, pol), '-.', label='$Q_g$')
-plt.plot(t_deg, q_mag(t, pol), label='Q_t')
+plt.plot(t_deg, q_mag(t, pol), label='$Q_t$')
-plt.grid()
 plt.xlabel(r'$\theta$, deg', fontsize=18)
-plt.ylabel('Q', fontsize=18)
+plt.ylabel('$Q$', fontsize=18)
 plt.legend(fontsize=18)
-plt.title('Beam efficiencies', fontsize=20)
+plt.title('Ray efficiencies', fontsize=20)
 plt.show()

Python/trap_forces_transverse.py

@@ -13,13 +13,14 @@ sns.set()
 
 # Intensity profile graphics
 rho = np.linspace(-r_max, r_max, 500)
-fig1 = plt.figure(1, figsize=(10, 6))
+plt.figure(1, figsize=(10, 6))
 I = gauss(rho)
 I0 = np.max(I)
 plt.plot(rho, I / I0)
 plt.fill_between(rho, I / I0, 0, alpha=0.3)
-plt.xlabel('r, m', fontsize=18)
+plt.title("Forces along Z axis")
-plt.ylabel('I(r)', fontsize=18)
+plt.xlabel('$r, m$', fontsize=18)
+plt.ylabel('$I(r)$', fontsize=18)
 
 
 # Integration
@@ -46,11 +47,12 @@ transverse = transverse_g + transverse_s
 
 
 # Graphics
-fig2 = plt.figure(2, figsize=(10, 6))
+plt.figure(2, figsize=(10, 6))
 plt.plot(y, transverse_g, '-.', lw=1, label='$F_{g}$')
 plt.plot(y, transverse_s, '--', lw=1, label='$F_{s}$')
 plt.plot(y, transverse, lw=1, label='$F_{t}$')
-plt.xlabel('r, m', fontsize=18)
+plt.xlabel('$r, m$', fontsize=18)
-plt.ylabel('F, m', fontsize=18)
+plt.ylabel('$F, N$', fontsize=18)
 plt.legend()
+
 plt.show()

README.md

@@ -8,7 +8,9 @@ For dielectric spheres with a diameter much larger than the wavelength, geometri
 
 A simple geometrical-optical model of a single-beam gradient trap is shown in Fig. 1.
 
-![](images/fig_1.png)
+|![](images/fig_1.png)|
+|---|
+|(A) A single-beam gradient optical trap in the consideration of geometric optics with the focus $f$ located on the axis $Z$. (B) The geometry of the incident beam contributing to the forces $F_g$ and $F_s$|
 
 The trap consists of an incident parallel beam of a better arbitrary mode composition and polarization, which falls into a high-aperture micro lens and focuses in focus $f$. The maximum angle of deflection of the rays at the output of a micro lens is determined by its numerical aperture
 
@@ -28,7 +30,6 @@ $$F_s = \frac{n_1 P}{c} \left[ 1 + R \cos{2\theta} - \frac{T^2[\cos(2\theta - 2r
 
 $$F_g = \frac{n_1 P}{c} \left[ R \sin{2\theta} - \frac{T^2[\sin(2\theta - 2r) + R\sin{2\theta}]}{1 + R^2 + 2R\cos{2r}} \right] = \frac{n_1 P}{c} Q_g$$
 
-
 where $F_s$ – scattering force [N];
 
 $F_g$ – gradient force (return force) [N];
@@ -59,8 +60,20 @@ It can be seen that the forces acting on the particle depend on the polarization
 
 Let's plot the efficiencies $Q$ of the forces obtained. We introduce the relative refractive index $n = \frac{n_2}{n_1}$. For the calculation, we will use the values $NA = 1.25$, from where $\phi_{max} = 70^\circ$; $n_1 = 1.33$; $𝑛 = 1.2$, from where $n_2 = 1.6$ (polystyrene). The value of $f$ and the radius of the particle do not need to be introduced, because in geometrical-optical consideration, the result does not depend on their values.
 
+|![](images/fig_2.png)|
+|---|
+|Fig. 2. Efficiency of forces acting on a particle at $n = 1.2$|
+
 In Fig. 2 $Q = \sqrt{Q_s^2 +Q_g^2}$ is the efficiency of the resulting force acting on the particle. This graph shows the contribution of each ray incident on a particle at a certain angle. It can be seen that the maximum value of $Q_g$ is approximately at the value of the angle of incidence of $70^\circ$, which shows the need to use high values of $NA$. For comparison, we present graphs for other $n$ in Fig. 3 and 4.
 
+|![](images/fig_3.png)|
+|---|
+|Fig. 3. Efficiency of forces acting on a particle at $n = 1.1$|
+
+|![](images/fig_4.png)|
+|---|
+|Fig. 4. Efficiency of forces acting on a particle at $n = 1.8$|
+
 It can be seen that at $n\rightarrow 1$ the required maximum angle of incidence becomes too high, and at $n\rightarrow 1.4$ $Q_s$ is equal to or exceeds $Q_g$ in most of the angle range, which indicates the difficulty of obtaining a stable capture.
 
 ## Force along Z-axis
@@ -92,6 +105,10 @@ The efficiencies of forces acting along the $Z$ $Q_{tz}$, $Q_{gz}$, $Q_{sz}$ axi
 
 Let's build graphs based on the data given above (Fig. 5).
 
+|![](images/fig_5.png)|
+|---|
+|Fig. 5. Forces arising from the displacement of a particle along the axis $Z$|
+
 It can be seen from the graph that when the particle is displaced from the equilibrium position, the forces that return the particle back grow. The force acquires maximum values in the coordinates $z = -1.06\ \mu m$ and $z = 1.02\ \mu m$. The stability of the capture and the rigidity of the trap depend on the specific values of these forces. It can also be noticed that the equilibrium position of the particle is slightly shifted relative to zero and that due to symmetry in this case, the force does not depend on the polarization of light.
 
 ## Force along Y-axis
@@ -136,7 +153,11 @@ $$F_{gY}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{gY
 
 $$F_{sY}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{sY} r d\beta dr$$
 
-Using the same data, we will plot the efficiencies of forces along the $Y$ axis.
+Using the same data, we will plot the efficiencies of forces along the $Y$ axis (Fig. 7).
+
+|![](images/fig_7.png)|
+|---|
+|Fig. 7. Forces arising from the displacement of a particle along the axis $Y$|
 
 It can be seen that the returning force acts symmetrically when the particle is displaced, and the equilibrium position is already at zero.
 
@@ -146,7 +167,7 @@ In the case of a Gaussian beam (the main $TEM_{00}$ mode), the intensity distrib
 
 $$I(r) = \exp\left(-\frac{2r^2}{w_0}\right)$$
 
-Graphs of forces acting along the $Z$ and $Y$ axes under the condition $\varepsilon = 1.0\ (w_0 = r_{max})$ are shown in Fig. 7 and 8.
+Graphs of forces acting along the $Z$ and $Y$ axes under the condition $\varepsilon = 1.0\ (w_0 = r_{max})$ are shown in Fig. 5 and 7.
 
 Calculate the distributions for other values of $\varepsilon$ and summarize the data in Table 1. In this table, the value $A = 1 - \exp\left(-\frac{2r_{max}^2}{w_0}\right)$ is the fraction of power that falls on the input pupil of the micro lens. 
 
@@ -170,7 +191,7 @@ Everything described above is designed in the form of Python and Matlab files:
 
 - Forces along $Z$ axis: [trap_forces_axial.py](Python/trap_forces_axial.py) and [trap_forces_axial.m](Matlab/trap_forces_axial.m)
 
-- Forces along $Y$ axis: [trap_forces_transverse.py](Python/trap_forces_transverse.py) and [trap_forces_transverse.m](Matlab/trap_forces_transverse.m)
+- Forces along $Y$ axis: [trap_forces_transverse.py](Python/trap_forces_transverse.py) and [trap_forces_transverse.m](Matlab/trap_forces_transverse.m) (**takes more time to calculate!**)
 
 Launch one of these files, and the graphs will be plotted. You can change initial constants and do your experiments.
 

README_RU.md

@@ -8,7 +8,7 @@
 
 |![](images/fig_1.png)|
 |---|
-|Fig. 6. (A) Однолучевая градиентная оптическая ловушка в рассмотрении геометрической оптики с фокусом $f$, расположенным на оси $Z$. (B) Геометрия падающего луча, дающего вклад в силы $F_g$ и $F_S$|
+|Рис. 1. (A) Однолучевая градиентная оптическая ловушка в рассмотрении геометрической оптики с фокусом $f$, расположенным на оси $Z$. (B) Геометрия падающего луча, дающего вклад в силы $F_g$ и $F_S$|
 
 Ловушка состоит из падающего параллельного пучка лучше произвольного модового состава и поляризации, который попадает в высокоапертурный микрообъектив и фокусируется в фокусе $f$. Максимальный угол отклонения лучей на выходе микрообъектива определяется его числовой апертурой
 
@@ -60,8 +60,20 @@ $$n_1 \sin{\theta} = n_2 \sin{r}$$
 Построим графики для эффективностей $Q$ полученных сил. Введем относительный показатель преломления $n = \frac{n_2}{n_1}$. Для расчета будем использовать значения $NA = 1.25$, откуда $\phi_{max} = 70^\circ$;
 $n_1 = 1.33$; $𝑛 = 1.2$, откуда $n_2 = 1.6$ (полистирол). Величину $f$ и радиус частицы вводить не нужно, т.к. в геометрооптическом рассмотрении от их значений результат не зависит.
 
+|![](images/fig_2.png)|
+|---|
+|Рис. 2. Эффективности сил, действующих на частицу при $n = 1.2$|
+
 На рис. 2 $Q = \sqrt{Q_s^2 + Q_g^2}$ – эффективность результирующей силы, действующей на частицу. Этот график показывает вклад каждого луча, падающего на частицу под определенным углом. Видно, что максимальное значение $Q_g$ находится примерно при значении угла падения $70^\circ$, что показывает необходимость использования высоких значений $NA$. Для сравнения приведем графики для других $n$ на рис. 3 и 4.
 
+|![](images/fig_3.png)|
+|---|
+|Рис. 3. Эффективности сил, действующих на частицу при $n = 1.1$|
+
+|![](images/fig_4.png)|
+|---|
+|Рис. 4. Эффективности сил, действующих на частицу при $n = 1.8$|
+
 Видно, что при $n \rightarrow 1$ необходимый максимальный угол падения становится слишком высок, а при $n \rightarrow 1.4\$ $Q_s$ равна или превышает $Q_g$ в большей части диапазона углов, что говорит о сложности получения стабильного захвата.
 
 ## Сила вдоль оси Z
@@ -93,6 +105,10 @@ $$F_{sZ}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{sZ
 
 Построим графики на основе данных, приведенных выше (рис. 5).
 
+|![](images/fig_5.png)|
+|---|
+|Рис. 5. Силы, возникающие при смещении частицы вдоль оси $Z$|
+
 Из графика видно, что при смещении частицы из положения равновесия растут силы, которые возвращают частицу обратно. Максимальные значения сила приобретает в координатах $z = −1.06\ \mu m$ и $z = 1.02\ \mu m$. От конкретных значений этих сил зависит стабильность захвата и жесткость ловушки. Также можно заметить, что положение равновесия частицы немного смещено относительно нуля и то, что благодаря симметрии в данном случае сила не зависит от поляризации света.
 
 ## Сила вдоль оси Y
@@ -137,7 +153,11 @@ $$F_{gY}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{gY
 
 $$F_{sY}^\Sigma = \frac{1}{\pi r_{max}} \int_{0}^{2\pi} \int_{0}^{r_{max}} F_{sY} r d\beta dr$$
 
-По тем же данным построим графики эффективностей сил вдоль оси $Y$.
+По тем же данным построим графики эффективностей сил вдоль оси $Y$ (рис. 7).
+
+|![](images/fig_7.png)|
+|---|
+|Рис. 7. Силы, возникающие при смещении частицы вдоль оси $Y$|
 
 Видно, что возвращающая сила симметрично действует при смещении частицы, и положение равновесия уже находится в нуле.
 
@@ -171,7 +191,7 @@ $$I(r) = \exp\left(-\frac{2r^2}{w_0}\right)$$
 
 - Силы вдоль оси $Z$: [trap_forces_axial.py](Python/trap_forces_axial.py) и [trap_forces_axial.m](Matlab/trap_forces_axial.m)
 
-- Силы вдоль оси $Y$: [trap_forces_transverse.py](Python/trap_forces_transverse.py) и [trap_forces_transverse.m](Matlab/trap_forces_transverse.m)
+- Силы вдоль оси $Y$: [trap_forces_transverse.py](Python/trap_forces_transverse.py) и [trap_forces_transverse.m](Matlab/trap_forces_transverse.m) (**занимает больше времени на вычисление!**)
 
 Запустите один из этих файлов, и будут построены графики. Вы можете изменить начальные константы и провести свои эксперименты.