ModuleSource

55def PointSource(S: np.ndarray, Axis: np.ndarray, Divergence: float, NbRays: int, Wavelength=None):
56    """
57    Return a list of rays simulating rays from a point source.
58
59    ![Illustration of a point source.](pointsource.svg)
60
61    Parameters
62    ----------
63        S : np.ndarray
64            Coordinate vector of the source point (in mm).
65
66        Axis : np.ndarray
67            Axis vector of the cone (main direction of the rays)
68
69        Divergence : float
70            Cone apex *half*-angle in *radians* (similar to gaussian beam).
71
72        NbRays : int
73            Number of rays in the list.
74
75        Wavelength : float, optional
76            Wavelength of the light, set as attribute to all rays.
77
78    """
79    RayList = _Cone(Divergence, NbRays, Wavelength=Wavelength)
80    RayList = mgeo.RotationRayList(RayList, np.array([0, 0, 1]), Axis)
81    RayList = mgeo.TranslationRayList(RayList, S)
82    return RayList

 86def ExtendedSource(S: np.ndarray, Axis: np.ndarray, Diameter: float, Divergence: float, NbRays: int, Wavelength=None):
 87    """
 88    Return a list of rays simulating rays from an extended array of point sources, distributed over a disk with Diameter.
 89
 90    Parameters
 91    ----------
 92        S : np.ndarray
 93            Coordinate vector of the source point (in mm).
 94
 95        Axis : np.ndarray
 96            Axis vector of the cone (main direction of the rays)
 97
 98        Diameter : float
 99            Diameter of the extended source
100
101        Divergence : float
102            Cone apex *half*-angle in *radians* (similar to gaussian beam).
103
104        NbRays : int
105            Number of rays in the list.
106
107        Wavelength : float, optional
108            Wavelength of the light, set as attribute to all rays.
109
110    """
111    NbPointSources = max(5, int(250 * Diameter))
112    NbPointSources = min(NbPointSources, 100)
113
114    MatrixXY = mgeo.SpiralVogel(NbPointSources, Diameter / 2)  # the positions of the point sources
115
116    NbRaysPerPointSource = max(100, int(NbRays / NbPointSources))
117    RayList = []
118    PointSourceRayList = _Cone(Divergence, NbRaysPerPointSource, Wavelength=Wavelength)
119    for k in range(NbPointSources):
120        ShiftedPointSourceRayList = mgeo.TranslationRayList(PointSourceRayList, [MatrixXY[k, 0], MatrixXY[k, 1], 0])
121        for l in range(NbRaysPerPointSource):
122            ShiftedPointSourceRayList[l]._number += (
123                k * NbRaysPerPointSource
124            )  # we allow ourselves exceptionally to modify the "internal" _number attribute
125        RayList.extend(ShiftedPointSourceRayList)
126
127    RayList = mgeo.RotationRayList(RayList, np.array([0, 0, 1]), Axis)
128    RayList = mgeo.TranslationRayList(RayList, S)
129    return RayList

133def PlaneWaveDisk(Centre: np.ndarray, Axis: np.ndarray, Radius: float, NbRays: int, Wavelength=None):
134    """
135    Return a list of rays, all parallel and distributed over a round cross section according to a Vogel-spiral (see illustration). Simulating a 'round' collimated beam.
136
137    ![Illustration of a Vogel spiral.](doc_VogelSpiral.png)
138
139    Parameters
140    ----------
141        Centre : np.ndarray
142            Coordinate vector of the center of the source disk (in mm).
143
144        Axis : np.ndarray
145            Axis vector of the the plane wave (vector of all rays)
146
147        Radius : float
148            Radius of the source ray-bundle.
149
150        NbRays : int
151            Number of rays in the list.
152
153        Wavelength : float, optional
154            Wavelength of the light, set as attribute to all rays.
155
156    """
157    MatrixXY = mgeo.SpiralVogel(NbRays, Radius)
158    RayList = []
159    num = 0
160    for k in range(NbRays - 1):
161        x = MatrixXY[k, 0]
162        y = MatrixXY[k, 1]
163        RayList.append(mray.Ray(np.array([x, y, 0]), np.array([0, 0, 1]), Number=num, Wavelength=Wavelength))
164        num = num + 1
165    RayList = mgeo.RotationRayList(RayList, np.array([0, 0, 1]), Axis)
166    RayList = mgeo.TranslationRayList(RayList, Centre)
167    return RayList

171def PlaneWaveSquare(Centre: np.ndarray, Axis: np.ndarray, SideLength: float, NbRays: int, Wavelength=None):
172    """
173    Return a list of rays, all parallel and distributed over a square cross section awith SideLength. Simulating a 'square' collimated beam.
174
175    Parameters
176    ----------
177        Centre : np.ndarray
178            Coordinate vector of the center of the source disk (in mm).
179
180        Axis : np.ndarray
181            Axis vector of the the plane wave (vector of all rays)
182
183        SideLength : float
184            Side length of the square source ray-bundle.
185
186        NbRays : int
187            Number of rays in the list.
188
189        Wavelength : float, optional
190            Wavelength of the light, set as attribute to all rays.
191
192    """
193    RayList = []
194    x = np.linspace(-SideLength / 2, SideLength / 2, int(np.sqrt(NbRays)))
195    y = np.linspace(-SideLength / 2, SideLength / 2, int(np.sqrt(NbRays)))
196    RayList.append(mray.Ray(np.array([0, 0, 0]), np.array([0, 0, 1]), Number=0))
197    num = 1
198    for i in x:
199        for j in y:
200            if abs(x) > 1e-4 and abs(y) > 1e-4:
201                RayList.append(mray.Ray(np.array([i, j, 0]), np.array([0, 0, 1]), Number=num, Wavelength=Wavelength))
202                num += 1
203    RayList = mgeo.RotationRayList(RayList, np.array([0, 0, 1]), Axis)
204    RayList = mgeo.TranslationRayList(RayList, Centre)
205    return RayList

217def ApplyGaussianIntensityToRayList(RayList, IntensityFraction=1 / np.e**2):
218    """
219    Apply a Gaussain intensity profile to a ray list, with intensity =1 at the center and falling
220    to 'IntensityFraction' at the edge (arbitrary units).
221
222    Parameters
223    ----------
224        RayList : list of Ray-objects
225            Coordinate vector of the center of the source disk (in mm).
226
227        IntensityFraction : float, optional
228            Relative intensity in (0,1), reached at the edge of the ray bundle. Default is 1/e^2.
229
230    """
231    if IntensityFraction >= 1 or IntensityFraction <= 0:
232        # raise ValueError('IntensityFraction should be between 0 and 1!')
233        print(
234            "When applying a Gaussian intensity profile to a ray list, the IntensityFraction should be between 0 and 1! I'm setting it to 1/e^2."
235        )
236        IntensityFraction = 1 / np.e**2
237
238    Axis = mp.FindCentralRay(RayList).vector
239    Divergence = 0
240    for l in RayList:
241        # find the largest angle with the central ray out of all the rays in the list
242        Divergence = max(Divergence, mgeo.AngleBetweenTwoVectors(Axis, l.vector))
243
244    if (
245        Divergence > 1e-12
246    ):  # we're dealing with a point source and not a plane wave, so we can apply Gaussian intensity profile as function of angle
247        for k in RayList:
248            Angle = mgeo.AngleBetweenTwoVectors(Axis, k.vector)
249            Intensity = np.exp(-2 * (np.tan(Angle) / Divergence) ** 2 * -0.5 * np.log(IntensityFraction))
250            k.intensity = Intensity
251    else:  # otherwise we apply it as function of ray distance from the Axis
252        MaxDist = 0
253        for l in RayList:
254            MaxDist = max(MaxDist, np.linalg.norm(l.point))
255        for k in RayList:
256            Intensity = np.exp(-2 * (np.linalg.norm(k.point) / MaxDist) ** 2 * -0.5 * np.log(IntensityFraction))
257            k.intensity = Intensity
258
259    return RayList