ModuleAnalysis

  1"""
  2Adds various analysis methods.
  3
  4Contains the following functions:
  5    - GetETransmission: Calculates the energy transmission from RayListIn to RayListOut in percent.
  6    - GetResultSummary: Calculate and return FocalSpotSize-standard-deviation and Duration-standard-deviation.
  7    - GetPulseProfile: Retrieve the pulse profile from the given RayList and Detector.
  8    - GetNumericalAperture: Returns the numerical aperture associated with the supplied ray-bundle.
  9    - GetAiryRadius: Returns the radius of the Airy disk.
 10    - GetPlaneWaveFocus: Calculates the approximate polychromatic focal spot of a set of rays.
 11    - GetDiffractionFocus: Calculates the approximate polychromatic focal spot of a set of rays.
 12    - GetClosestSphere: Calculates the closest sphere to the surface of a mirror.
 13    - GetAsphericity: Calculates the maximum distance of the mirror surface to the closest sphere.
 14    - _best_fit_sphere: Calculates the best sphere to fit a set of points.
 15
 16Adds the following methods:
 17    - OpticalChain:
 18        - getETransmission: Calculates the energy transmission from the input to the output of the OpticalChain.
 19        - getResultsSummary: Calculate and return FocalSpotSize-standard-deviation and Duration-standard-deviation.
 20        - getPulseProfile: Retrieve the pulse profile from the output of the OpticalChain.
 21        - getPlaneWaveFocus: Calculates the approximate polychromatic focal spot of the output of the OpticalChain.
 22        - getDiffractionFocus: Calculates the approximate polychromatic focal spot of the output of the OpticalChain.
 23    - Mirror:
 24        - getClosestSphere: Calculates the closest sphere to the surface of a mirror.
 25        - getAsphericity: Calculates the maximum distance of the mirror surface to the closest sphere.
 26
 27Created in July 2024
 28
 29@author: André Kalouguine + Stefan Haessler + Anthony Guillaume
 30"""
 31# %% Module imports
 32import ARTcore.ModuleGeometry as mgeo
 33import ARTcore.ModuleOpticalRay as mray
 34import ARTcore.ModuleProcessing as mp
 35import ARTcore.ModuleOpticalChain as moc
 36import ARTcore.ModuleMirror as mmirror
 37import ART.ModulePlottingMethods as mpm
 38import matplotlib.pyplot as plt
 39from mpl_toolkits.axes_grid1 import make_axes_locatable
 40import numpy as np
 41import math
 42
 43LightSpeed = 299792458000
 44
 45# %% Analysis methods
 46def GetETransmission(RayListIn, RayListOut) -> float:
 47    """
 48    Calculates the energy transmission from RayListIn to RayListOut in percent by summing up the
 49    intensity-property of the individual rays.
 50
 51    Parameters
 52    ----------
 53        RayListIn : list(Ray)
 54            List of incoming rays.
 55        
 56        RayListOut : list(Ray)
 57            List of outgoing rays.
 58
 59    Returns
 60    -------
 61        ETransmission : float
 62    """
 63    ETransmission = 100 * sum(Ray.intensity for Ray in RayListOut) / sum(Ray.intensity for Ray in RayListIn)
 64    return ETransmission
 65
 66def _getETransmission(OpticalChain, IndexIn=0, IndexOut=-1) -> float:
 67    """
 68    Calculates the energy transmission from the input to the output of the OpticalChain in percent.
 69
 70    Parameters
 71    ----------
 72        OpticalChain : OpticalChain
 73            An object of the ModuleOpticalChain.OpticalChain-class.
 74
 75        IndexIn : int, optional
 76            Index of the input RayList in the OpticalChain, defaults to 0.
 77
 78        IndexOut : int, optional
 79            Index of the output RayList in the OpticalChain, defaults to -1.
 80
 81    Returns
 82    -------
 83        ETransmission : float
 84    """
 85    Rays = OpticalChain.get_output_rays()
 86    if IndexIn == 0:
 87        RayListIn = OpticalChain.get_input_rays()
 88    else:
 89        RayListIn = Rays[IndexIn]
 90    RayListOut = Rays[IndexOut]
 91    ETransmission = GetETransmission(RayListIn, RayListOut)
 92    return ETransmission
 93
 94moc.OpticalChain.getETransmission = _getETransmission
 95
 96def GetResultSummary(Detector, RayListAnalysed, verbose=False):
 97    """
 98    Calculate and return FocalSpotSize-standard-deviation and Duration-standard-deviation
 99    for the given Detector and RayList.
100    If verbose, then also print a summary of the results for the given Detector.
101
102    Parameters
103    ----------
104        Detector : Detector
105            An object of the ModuleDetector.Detector-class.
106
107        RayListAnalysed : list(Ray)
108            List of objects of the ModuleOpticalRay.Ray-class.
109
110        verbose : bool
111            Whether to print a result summary.
112
113    Returns
114    -------
115        FocalSpotSizeSD : float
116
117        DurationSD : float
118    """
119    DetectorPointList2DCentre = Detector.get_2D_points(RayListAnalysed)
120    FocalSpotSizeSD = mp.StandardDeviation(DetectorPointList2DCentre)
121    DelayList = Detector.get_Delays(RayListAnalysed)
122    DurationSD = mp.StandardDeviation(DelayList)
123
124    if verbose:
125        FocalSpotSize = mgeo.DiameterPointList(DetectorPointList2DCentre)
126        summarystring = (
127            "At the detector distance of "
128            + "{:.3f}".format(Detector.get_distance())
129            + " mm we get:\n"
130            + "Spatial std : "
131            + "{:.3f}".format(FocalSpotSizeSD * 1e3)
132            + " \u03BCm and min-max: "
133            + "{:.3f}".format(FocalSpotSize * 1e3)
134            + " \u03BCm\n"
135            + "Temporal std : "
136            + "{:.3e}".format(DurationSD)
137            + " fs and min-max : "
138            + "{:.3e}".format(max(DelayList) - min(DelayList))
139            + " fs"
140        )
141
142        print(summarystring)
143
144    return FocalSpotSizeSD, DurationSD
145
146def _getResultsSummary(OpticalChain, Detector = "Focus", verbose=False):
147    """
148    Calculate and return FocalSpotSize-standard-deviation and Duration-standard-deviation
149    for the given Detector and RayList.
150    If verbose, then also print a summary of the results for the given Detector.
151
152    Parameters
153    ----------
154        OpticalChain : OpticalChain
155            An object of the ModuleOpticalChain.OpticalChain-class.
156
157        Detector : Detector or str, optional
158            An object of the ModuleDetector.Detector-class or "Focus" to use the focus detector, defaults to "Focus".
159
160        verbose : bool
161            Whether to print a result summary.
162
163    Returns
164    -------
165        FocalSpotSizeSD : float
166
167        DurationSD : float
168    """
169    if isinstance(Detector, str):
170        Detector = OpticalChain.detectors[Detector]
171    Index = Detector.index
172    RayListAnalysed = OpticalChain.get_output_rays()[Index]
173    FocalSpotSizeSD, DurationSD = GetResultSummary(Detector, RayListAnalysed, verbose)
174    return FocalSpotSizeSD, DurationSD
175
176moc.OpticalChain.getResultsSummary = _getResultsSummary
177
178def GetPulseProfile(Detector, RayList, Nbins=100):
179    """
180    Retrieve the pulse profile from the given RayList and Detector.
181    The pulse profile is calculated by binning the delays of the rays in the RayList
182    at the Detector. The intensity of the rays is also taken into account.
183
184    Parameters
185    ----------
186        Detector : Detector
187            An object of the ModuleDetector.Detector-class.
188
189        RayList : list(Ray)
190            List of objects of the ModuleOpticalRay.Ray-class.
191
192        Nbins : int
193            Number of bins to use for the histogram.
194
195    Returns
196    -------
197        time : numpy.ndarray
198            The bins of the histogram.
199
200        intensity : numpy.ndarray
201            The histogram values.
202    """
203    DelayList = Detector.get_Delays(RayList)
204    IntensityList = [Ray.intensity for Ray in RayList]
205    intensity, time_edges = np.histogram(DelayList, bins=Nbins, weights=IntensityList)
206    time = 0.5 * (time_edges[1:] + time_edges[:-1])
207    return time, intensity
208
209def _getPulseProfile(OpticalChain, Detector = "Focus", Nbins=100):
210    """
211    Retrieve the pulse profile from the output of the OpticalChain.
212    The pulse profile is calculated by binning the delays of the rays in the RayList
213    at the Detector. The intensity of the rays is also taken into account.
214
215    Parameters
216    ----------
217        OpticalChain : OpticalChain
218            An object of the ModuleOpticalChain.OpticalChain-class.
219
220        Detector : Detector or str, optional
221            An object of the ModuleDetector.Detector-class or "Focus" to use the focus detector, defaults to "Focus".
222
223        Nbins : int
224            Number of bins to use for the histogram.
225
226    Returns
227    -------
228        time : numpy.ndarray
229            The bins of the histogram.
230
231        intensity : numpy.ndarray
232            The histogram values.
233    """
234    if isinstance(Detector, str):
235        Detector = OpticalChain.detectors[Detector]
236    Index = Detector.index
237    RayList = OpticalChain.get_output_rays()[Index]
238    time, intensity = GetPulseProfile(Detector, RayList, Nbins)
239    return time, intensity
240
241
242def GetNumericalAperture(RayList: list[mray.Ray], RefractiveIndex: float = 1) -> float:
243    r"""
244    Returns the numerical aperture associated with the supplied ray-bundle 'Raylist'.
245    This is $n\sin\theta$, where $\theta$ is the maximum angle between any of the rays and the central ray,
246    and $n$ is the refractive index of the propagation medium.
247
248    Parameters
249    ----------
250        RayList : list of Ray-object
251            The ray-bundle of which to determine the NA.
252
253        RefractiveIndex : float, optional
254            Refractive index of the propagation medium, defaults to =1.
255
256    Returns
257    -------
258        NA : float
259    """
260    CentralRay = mp.FindCentralRay(RayList)
261    if CentralRay is None:
262        CentralVector = np.array([0, 0, 0])
263        for k in RayList:
264            CentralVector = CentralVector + k.vector
265        CentralVector = CentralVector / len(RayList)
266    else:
267        CentralVector = CentralRay.vector
268    ListAngleAperture = []
269    for k in RayList:
270        ListAngleAperture.append(mgeo.AngleBetweenTwoVectors(CentralVector, k.vector))
271
272    return np.sin(np.amax(ListAngleAperture)) * RefractiveIndex
273
274
275def GetAiryRadius(Wavelength: float, NumericalAperture: float) -> float:
276    r"""
277    Returns the radius of the Airy disk: $r = 1.22 \frac\{\lambda\}\{NA\}$,
278    i.e. the diffraction-limited radius of the focal spot corresponding to a given
279    numerical aperture $NA$ and a light wavelength $\lambda$.
280
281    For very small $NA<10^\{-3\}$, diffraction effects becomes negligible and the Airy Radius becomes meaningless,
282    so in that case, a radius of 0 is returned.
283
284    Parameters
285    ----------
286        Wavelength : float
287            Light wavelength in mm.
288
289        NumericalAperture : float
290
291    Returns
292    -------
293        AiryRadius : float
294    """
295    if NumericalAperture > 1e-3 and Wavelength is not None:
296        return 1.22 * 0.5 * Wavelength / NumericalAperture
297    else:
298        return 0  # for very small numerical apertures, diffraction effects becomes negligible and the Airy Radius becomes meaningless
299
300
301def GetPlaneWaveFocus(OpticalChain, Detector = "Focus", size=None, Nrays=1000, resolution=100):
302    """
303    This function calculates the approximate polychromatic focal spot of a set of rays.
304    To do so, it positions itself in the detector plane and samples a square area.
305    If the size of the area is not given, it will use twice the Airy radius of the system with the largest wavbelength.
306    Otherwise it uses the size.
307    The resolution of the sampling is given by the resolution parameter. So it samples on a grid resolution x resolution.
308
309    To calculate the intensity, it takes Nrays out of the raylist (to subsample if needed).
310    It assimilates every ray to a plane wave, so it calculates a k-vector for each ray (taking into account the wavelength of the ray).
311    It calculates the phase of each ray from the delay and from the intersection position with the detector.
312    The delay from the non-central position is simply sin(alpha)*distance/c, where alpha is the angle between the ray and the normal to the detector
313    and distance is the distance between the intersection point and the central point of the detector. 
314    It then calculates the intensity at each point of the grid by summing the intensity of each plane wave.
315
316    As long as there are not too many rays, this method is faster than doing an FFT.
317    It's also naturally polychromatic.
318    """
319    if isinstance(Detector, str):
320        Detector = OpticalChain.detectors[Detector]
321    Index = Detector.index
322    RayList = OpticalChain.get_output_rays()[Index]
323    if size is None:
324        Wavelengths = [Ray.wavelength for Ray in RayList]
325        Wavelength = max(Wavelengths)
326        NumericalAperture = GetNumericalAperture(RayList)
327        size = 3 * GetAiryRadius(Wavelength, NumericalAperture)
328    X = np.linspace(-size / 2, size / 2, resolution)
329    Y = np.linspace(-size / 2, size / 2, resolution)
330    X, Y = np.meshgrid(X, Y)
331    # We pick Nrays. if there are less than Nrays, we take all of them.
332    PickedRaysGlobal = np.random.choice(RayList, min(Nrays, len(RayList)), replace=False)
333    # We calculate the k-vector of each ray, taking into account the wavelength of the ray.
334    # The units should be in mm^-1
335    PickedRays = [Ray.to_basis(*Detector.basis) for Ray in PickedRaysGlobal]
336    # The rays are now in the reference plane of the detector whose normal is [0,0,1]
337    wavelengths = np.array([Ray.wavelength for Ray in PickedRays])
338    vectors = np.array([Ray.vector for Ray in PickedRays])
339    frequencies = np.array([LightSpeed / Ray.wavelength for Ray in PickedRays]) # mm/s / mm = 1/s
340    k_vectors = 2*np.pi * vectors / wavelengths[:, np.newaxis]
341    angles = np.arccos(np.clip(np.dot(vectors, np.array([0, 0, 1])), -1.0, 1.0))
342    # We calculate the intersection of the rays with the detector plane
343    Intersections = mgeo.IntersectionRayListZPlane(PickedRays)[0]
344    distances = (Intersections - mgeo.Origin[:2]).norm
345    # We calculate the phase of each ray
346    PathDelays = np.array(Detector.get_Delays(PickedRaysGlobal))
347    PositionDelays = np.sum(np.array(Intersections._add_dimension()-mgeo.Origin)*vectors, axis=1)  / LightSpeed * 1e15
348    Delays = (PathDelays+PositionDelays) /  1e15
349    # We have the delays in fs, we can now calculate the phase, taking into account the wavelength of each ray
350    Phases = np.array([2 * np.pi * frequencies[i] * Delays[i] for i in range(len(Delays))])
351    # We also need the intensity of each ray
352    RayIntensities = np.array([Ray.intensity for Ray in PickedRays])
353    # We can now calculate the intensity at each point of the grid
354    Intensity = np.zeros((resolution, resolution))
355    for i in range(len(PickedRays)):
356        Intensity += RayIntensities[i] * np.cos(k_vectors[i][0] * X + k_vectors[i][1] * Y + Phases[i])
357
358    return X, Y, Intensity
359
360moc.getPlaneWaveFocus = GetPlaneWaveFocus
361
362def GetDiffractionFocus(OpticalChain, Detector = "Focus", size=None, Nrays=1000, resolution=100):
363    """
364    This function calculates the approximate polychromatic focal spot of a set of rays.
365    To do so, it positions itself in the detector plane and samples a square area.
366    If the size of the area is not given, it will use twice the Airy radius of the system with the largest wavbelength.
367    Otherwise it uses the size.
368    The resolution of the sampling is given by the resolution parameter. So it samples on a grid resolution x resolution.
369
370    To calculate the intensity, it takes Nrays out of the raylist (to subsample if needed).
371    It assimilates every ray to a plane wave, so it calculates a k-vector for each ray (taking into account the wavelength of the ray).
372    It considers that all the rays are intersecting the detector in the middle and doesn't take into account their phase.
373
374    So it returns the best case scenario for a diffraction limited focus with that numerical aperture
375    """
376    if isinstance(Detector, str):
377        Detector = OpticalChain.detectors[Detector]
378    Index = Detector.index
379    RayList = OpticalChain.get_output_rays()[Index]
380    if size is None:
381        Wavelengths = [Ray.wavelength for Ray in RayList]
382        Wavelength = max(Wavelengths)
383        NumericalAperture = GetNumericalAperture(RayList)
384        size = 3 * GetAiryRadius(Wavelength, NumericalAperture)
385    X = np.linspace(-size / 2, size / 2, resolution)
386    Y = np.linspace(-size / 2, size / 2, resolution)
387    X, Y = np.meshgrid(X, Y)
388    # We pick Nrays. if there are less than Nrays, we take all of them.
389    PickedRaysGlobal = np.random.choice(RayList, min(Nrays, len(RayList)), replace=False)
390    # We calculate the k-vector of each ray, taking into account the wavelength of the ray.
391    # The units should be in mm^-1
392    PickedRays = [Ray.to_basis(*Detector.basis) for Ray in PickedRaysGlobal]
393    # The rays are now in the reference plane of the detector whose normal is [0,0,1]
394    wavelengths = np.array([Ray.wavelength for Ray in PickedRays])
395    vectors = np.array([Ray.vector for Ray in PickedRays])
396    frequencies = np.array([LightSpeed / Ray.wavelength for Ray in PickedRays]) # mm/s / mm = 1/s
397    k_vectors = 2*np.pi * vectors / wavelengths[:, np.newaxis]
398    angles = np.arccos(np.clip(np.dot(vectors, np.array([0, 0, 1])), -1.0, 1.0))
399    # We calculate the intersection of the rays with the detector plane
400    Intersections = mgeo.IntersectionRayListZPlane(PickedRays)[0]
401    distances = (Intersections - mgeo.Origin[:2]).norm
402    # We also need the intensity of each ray
403    RayIntensities = np.array([Ray.intensity for Ray in PickedRays])
404    # We can now calculate the intensity at each point of the grid
405    Intensity = np.zeros((resolution, resolution))
406    for i in range(len(PickedRays)):
407        Intensity += RayIntensities[i] * np.cos(k_vectors[i][0] * X + k_vectors[i][1] * Y)
408
409    return X, Y, Intensity
410
411moc.getDiffractionFocus = GetDiffractionFocus
412
413# %% Asphericity analysis
414# This code is to do asphericity analysis of various surfaces
415# It calculates the closes sphere to the surface of an optical element
416# Then there are two functions. One that simply gives an asphericity value
417# and another one that actually plots the distance the the closest sphere 
418# in much the same way as we plot the MirrorProjection
419
420def BestFitSphere(X,Y,Z):
421    """
422    This function calculates the best sphere to fit a set of points.
423    It uses the least square method to find the center and the radius of the sphere.
424    Cite https://jekel.me/2015/Least-Squares-Sphere-Fit/
425    """
426    A = np.zeros((len(X),4))
427    A[:,0] = X*2
428    A[:,1] = Y*2
429    A[:,2] = Z*2
430    A[:,3] = 1
431
432    #   Assemble the f matrix
433    f = np.zeros((len(X),1))
434    f[:,0] = X**2 + Y**2 + Z**2
435    C, residules, rank, singval = np.linalg.lstsq(A,f)
436
437    #   solve for the radius
438    t = (C[0]*C[0])+(C[1]*C[1])+(C[2]*C[2])+C[3]
439    radius = math.sqrt(t)
440
441    return mgeo.Point(C[:3].flatten()),radius
442
443
444def GetClosestSphere(Mirror, Npoints=1000):
445    """
446    This function calculates the closest sphere to the surface of a mirror.
447    It does so by sampling the surface of the mirror at Npoints points.
448    It then calculates the closest sphere to these points.
449    It returns the radius of the sphere and the center of the sphere.
450    """
451    Points = mpm.sample_support(Mirror.support, Npoints=1000)
452    Points += Mirror.r0[:2]
453    Z = Mirror._zfunc(Points)
454    Points = mgeo.PointArray([Points[:, 0], Points[:, 1], Z]).T
455    spX, spY, spZ = Points[:, 0], Points[:, 1], Points[:, 2]
456    Center, Radius = BestFitSphere(spX, spY, spZ)
457    return Center, Radius
458
459def GetAsphericity(Mirror, Npoints=1000):
460    """
461    This function calculates the maximum distance of the mirror surface to the closest sphere. 
462    """
463    center, radius = GetClosestSphere(Mirror, Npoints)
464    Points = mpm.sample_support(Mirror.support, Npoints=1000)
465    Points += Mirror.r0[:2]
466    Z = Mirror._zfunc(Points)
467    Points = mgeo.PointArray([Points[:, 0], Points[:, 1], Z]).T
468    Points_centered = Points - center
469    Distance = np.linalg.norm(Points_centered, axis=1) - radius
470    Distance*=1e3 # To convert to µm
471    return np.ptp(Distance)
472
473mmirror.Mirror.getClosestSphere = GetClosestSphere
474mmirror.Mirror.getAsphericity = GetAsphericity