Repository: sympy/sympy
Pull request: physics : Added refractive index in Gaussian Beam Parameters and refactored code
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #19169 
Closes #19169 

#### Brief description of what is fixed or changed
According to earlier implementation ,when transmission to different medium would take place Rayleigh length `z_r` would change, while wavelength stays the same. Hence beam waist, `w_0` changed which is not physical.
Now as can be verified from the test cases when refractive index is doubled Rayleigh length proportionately gets doubled but beam waist remains the same.
Also code /docstring has been refactored and made more understandable for users/contributors.

1. waist is replaced by beam waist which is the fundamental parameter used .

All the information I have added/changed is validated by these online resources -
https://en.wikipedia.org/wiki/Gaussian_beam
https://phyweb.physics.nus.edu.sg/~l2000/pc2193/Experiments/lab.pdf

#### Other comments
Note : User must be aware of when the medium is changed . Users reading code must recognize what `w` is referring to at different parts of the code .
A test has been modified for better understanding of users.

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* physics.optics
  * Added refractive index to beam parameters.
<!-- END RELEASE NOTES -->


Problem statement:
Gaussian Optics / Refractive index not considered
Dear sympy maintainers,

please correct me, if I am wrong, but I fear that the beam waist in 

class sympy.physics.optics.gaussopt.BeamParameter

is not correctly computed.

From the source:

def w_0(self):
        """
        The beam waist (minimal radius).

        See Also
        ========

        w : the beam radius at `1/e^2` intensity

        Examples
        ========

        >>> from sympy.physics.optics import BeamParameter
        >>> p = BeamParameter(530e-9, 1, w=1e-3)
        >>> p.w_0
        0.00100000000000000
        """
        return sqrt(self.z_r/pi*self.wavelen)

After  transmission through a surface with the refractive index changing, the Rayleigh length z_r would change, while wavelength stays the same. According to this implementation, w_0 changes, which is not physical.

If I might help to solve this, I would be happy to contribute. 
However, I have not a very good understanding of sympy, and this code is interfaced of course with sympy.

Best regards,
Lukas

Hints from issue discussion:
I don't know the physics so well but you are of course welcome to correct it. Many contributors can help with the sympy aspects but I'm not sure how many would know the physics for this case.
Here the formula implemented is incorrect .This will fix itself once the correct formula is used .The relation between Rayleigh length `z_r` and beam waist `w_0` involves the refractive index too .
Check out beam parameters from wikipedia-
https://en.wikipedia.org/wiki/Gaussian_beam

I'll open up a pull request on this issue soon .

Commit message: Positional_Arguments_Fixed

Patch to review:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index e027d9ad0f46..c74e4620ec8f 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -4912,7 +4912,7 @@ def test_sympy__physics__optics__waves__TWave():
 
 def test_sympy__physics__optics__gaussopt__BeamParameter():
     from sympy.physics.optics import BeamParameter
-    assert _test_args(BeamParameter(530e-9, 1, w=1e-3))
+    assert _test_args(BeamParameter(530e-9, 1, w=1e-3, n=1))
 
 
 def test_sympy__physics__optics__medium__Medium():
diff --git a/sympy/physics/optics/gaussopt.py b/sympy/physics/optics/gaussopt.py
index 48d519e92fe0..56037ff8162f 100644
--- a/sympy/physics/optics/gaussopt.py
+++ b/sympy/physics/optics/gaussopt.py
@@ -487,6 +487,7 @@ class BeamParameter(Expr):
     z : the distance to waist, and
     w : the waist, or
     z_r : the rayleigh range.
+    n : the refractive index of medium.
 
     Examples
     ========
@@ -506,9 +507,9 @@ class BeamParameter(Expr):
     >>> from sympy.physics.optics import FreeSpace
     >>> fs = FreeSpace(10)
     >>> p1 = fs*p
-    >>> p.w.n()
+    >>> p.w_z.n()
     0.00101413072159615
-    >>> p1.w.n()
+    >>> p1.w_z.n()
     0.00210803120913829
 
     See Also
@@ -526,18 +527,19 @@ class BeamParameter(Expr):
     # subclass it. See:
     # https://groups.google.com/d/topic/sympy/7XkU07NRBEs/discussion
 
-    def __new__(cls, wavelen, z, z_r=None, w=None):
+    def __new__(cls, wavelen, z, z_r=None, w=None, n=1):
         wavelen = sympify(wavelen)
         z = sympify(z)
+        n = sympify(n)
 
         if z_r is not None and w is None:
             z_r = sympify(z_r)
         elif w is not None and z_r is None:
-            z_r = waist2rayleigh(sympify(w), wavelen)
-        else:
-            raise ValueError('Constructor expects exactly one named argument.')
+            z_r = waist2rayleigh(sympify(w), wavelen, n)
+        elif z_r is None and w is None:
+            raise ValueError('Must specify one of w and z_r.')
 
-        return Expr.__new__(cls, wavelen, z, z_r)
+        return Expr.__new__(cls, wavelen, z, z_r, n)
 
     @property
     def wavelen(self):
@@ -551,6 +553,10 @@ def z(self):
     def z_r(self):
         return self.args[2]
 
+    @property
+    def n(self):
+        return self.args[3]
+
     @property
     def q(self):
         """
@@ -582,9 +588,10 @@ def radius(self):
         return self.z*(1 + (self.z_r/self.z)**2)
 
     @property
-    def w(self):
+    def w_z(self):
         """
-        The beam radius at `1/e^2` intensity.
+        The radius of the beam w(z), at any position z along the beam.
+        The beam radius at `1/e^2` intensity (axial value).
 
         See Also
         ========
@@ -597,7 +604,7 @@ def w(self):
 
         >>> from sympy.physics.optics import BeamParameter
         >>> p = BeamParameter(530e-9, 1, w=1e-3)
-        >>> p.w
+        >>> p.w_z
         0.001*sqrt(0.2809/pi**2 + 1)
         """
         return self.w_0*sqrt(1 + (self.z/self.z_r)**2)
@@ -605,12 +612,12 @@ def w(self):
     @property
     def w_0(self):
         """
-        The beam waist (minimal radius).
+         The minimal radius of beam at `1/e^2` intensity (peak value).
 
         See Also
         ========
 
-        w : the beam radius at `1/e^2` intensity
+        w_z : the beam radius at `1/e^2` intensity (axial value).
 
         Examples
         ========
@@ -620,7 +627,7 @@ def w_0(self):
         >>> p.w_0
         0.00100000000000000
         """
-        return sqrt(self.z_r/pi*self.wavelen)
+        return sqrt(self.z_r/(pi*self.n)*self.wavelen)
 
     @property
     def divergence(self):
@@ -678,7 +685,7 @@ def waist_approximation_limit(self):
 # Utilities
 ###
 
-def waist2rayleigh(w, wavelen):
+def waist2rayleigh(w, wavelen, n=1):
     """
     Calculate the rayleigh range from the waist of a gaussian beam.
 
@@ -697,7 +704,7 @@ def waist2rayleigh(w, wavelen):
     pi*w**2/wavelen
     """
     w, wavelen = map(sympify, (w, wavelen))
-    return w**2*pi/wavelen
+    return w**2*n*pi/wavelen
 
 
 def rayleigh2waist(z_r, wavelen):
diff --git a/sympy/physics/optics/tests/test_gaussopt.py b/sympy/physics/optics/tests/test_gaussopt.py
index ed099d254433..523f7e564208 100644
--- a/sympy/physics/optics/tests/test_gaussopt.py
+++ b/sympy/physics/optics/tests/test_gaussopt.py
@@ -57,8 +57,8 @@ def test_gauss_opt():
     assert streq(N(p.z_r), Float(5.92753330865999))
     fs = FreeSpace(10)
     p1 = fs*p
-    assert streq(N(p.w), Float(0.00101413072159615))
-    assert streq(N(p1.w), Float(0.00210803120913829))
+    assert streq(N(p.w_z), Float(0.00101413072159615))
+    assert streq(N(p1.w_z), Float(0.00210803120913829))
 
     w, wavelen = symbols('w wavelen')
     assert waist2rayleigh(w, wavelen) == pi*w**2/wavelen
@@ -90,8 +90,13 @@ def test_gauss_opt():
     z, l, w = symbols('z l r', positive=True)
     p = BeamParameter(l, z, w=w)
     assert p.radius == z*(pi**2*w**4/(l**2*z**2) + 1)
-    assert p.w == w*sqrt(l**2*z**2/(pi**2*w**4) + 1)
+    assert p.w_z == w*sqrt(l**2*z**2/(pi**2*w**4) + 1)
     assert p.w_0 == w
     assert p.divergence == l/(pi*w)
     assert p.gouy == atan2(z, pi*w**2/l)
     assert p.waist_approximation_limit == 2*l/pi
+
+    p = BeamParameter(530e-9, 1, w=1e-3, n=2)
+    assert streq(p.q, 1 + 3.77358490566038*I*pi)
+    assert streq(N(p.z_r), Float(11.8550666173200))
+    assert streq(N(p.w_0), Float(0.00100000000000000))


Code review comments:
diff_hunk: @@ -582,9 +588,10 @@ def radius(self):
         return self.z*(1 + (self.z_r/self.z)**2)
 
     @property
-    def w(self):
+    def w_z(self):
path: sympy/physics/optics/gaussopt.py
text: @user1:
PR is not letting me unresolve conversation so I am adding comment here:

I don't find the word "spot" in SymPy -- where is the spot size talked about? It is not yet clear why this should be changed from `w` to `w_z`.

@author:
I'll unresolve that conversation ,thanks !

@author:
Unresolved and addressed that change !