ModuleDefects

  1"""
  2Provides classes for different deformations of optical surfaces.
  3
  4
  5Created on Fri Mar 10 2023
  6
  7@author: Semptum + Stefan Haessler 
  8"""
  9# %% Modules
 10
 11import numpy as np
 12#import cupy as cp
 13from abc import ABC, abstractmethod
 14import scipy
 15import math
 16from ARTcore.ModuleZernike import zernike_gradient
 17import logging
 18
 19logger = logging.getLogger(__name__)
 20
 21# %% Abstract class
 22
 23
 24class Defect(ABC):
 25    """
 26    Abstract class representing a defect on the optical surface.
 27    """
 28
 29    @abstractmethod
 30    def RMS(self):
 31        pass
 32
 33    @abstractmethod
 34    def PV(self):
 35        pass
 36
 37class MeasuredMap(Defect):
 38    """
 39    Class describing a defect map measured experimentally (for instance using a Hartmann wavefront sensor). Note that the provided image should cover the entire support. 
 40    """
 41    def __init__(self, Support, Map):
 42        self.deformation = Map
 43        self. Support = Support
 44        rect = Support._CircumRect() # TODO fix this so it just uses a rectangle with known size
 45        X = np.linspace(-rect[0], rect[0], num=self.deformation.shape[0])
 46        Y = np.linspace(-rect[1], rect[1], num=self.deformation.shape[1])
 47        self.DerivX, self.DerivY = np.gradient(self.deformation, rect/self.deformation.shape)
 48        DerivInterpX = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(self.DerivX), method="linear")
 49        DerivInterpY = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(self.DerivY), method="linear")
 50        SurfInterp = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(self.deformation), method="linear")
 51
 52        self.rms = np.std(self.deformation)
 53        self.DerivInterp = lambda x: np.array([DerivInterpX(x[:2]), DerivInterpY(x[:2])])
 54        self.SurfInterp = lambda x: SurfInterp(x[:2])
 55    def get_normal(self, Point):
 56        Grad = self.DerivInterp(Point)
 57        dX, dY =  Grad.flatten()
 58        norm = np.linalg.norm([dX, dY, 1])
 59        dX /= norm
 60        dY /= norm
 61        return np.array([dX, dY, np.sqrt(1 - dX**2 - dY**2)])
 62
 63    def get_offset(self, Point):
 64        return self.SurfInterp(Point)
 65
 66    def RMS(self):
 67        return self.rms
 68
 69    def PV(self):
 70        pass
 71
 72class Fourrier(Defect):
 73    """
 74    
 75    """
 76    def __init__(self, Support, RMS, slope=-2, smallest=0.1, biggest=None):
 77        # The sizes are the wavelength in mm
 78        rect = Support._CircumRect()  # TODO fix this so it just uses a rectangle with known size
 79        if biggest is None:
 80            biggest = np.max(rect)
 81            
 82        k_max = 2 / smallest
 83        k_min = 2 / biggest
 84        ResX = int(round(k_max * rect[0] / 2))+1
 85        ResY = int(round(k_max * rect[1]))
 86
 87        kXX, kYY = np.meshgrid(
 88            np.linspace(0, k_max, num=ResX, dtype='float32', endpoint=False),
 89            np.linspace(-k_max, k_max, num=ResY, dtype='float32', endpoint=False),
 90            sparse=True)
 91
 92        maskedFFT = np.ma.masked_outside(np.sqrt(kXX**2 + kYY**2), k_min, k_max)
 93
 94        FFT = maskedFFT**slope * np.exp(
 95            1j *np.random.uniform(0, 2 * np.pi, size=maskedFFT.shape).astype('float32')
 96            )
 97        FFT = FFT.data*(1-FFT.mask)
 98
 99        deformation = np.fft.irfft2(np.fft.ifftshift(FFT, axes=0))
100        RMS_factor = RMS/np.std(deformation)
101        deformation *= RMS_factor
102
103        DerivX = np.fft.irfft2(np.fft.ifftshift(FFT * 1j * kXX * RMS_factor, axes=0))*np.pi/2
104        kY = np.concatenate((kYY[kYY.shape[0]//2:],kYY[:kYY.shape[0]//2]))
105        DerivY = np.fft.irfft2(np.fft.ifftshift(FFT * 1j * RMS_factor, axes=0)*kY)*np.pi/2
106        #del FFT
107
108        X = np.linspace(-rect[0]/2, rect[0]/2, num=(ResX-1)*2) # Because iRfft
109        Y = np.linspace(-rect[1]/2, rect[1]/2, num=ResY)
110        
111        DerivInterpX = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(DerivX), method="linear")
112        DerivInterpY = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(DerivY), method="linear")
113        SurfInterp = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(deformation), method="linear")
114
115        self.DerivX = DerivX
116        self.DerivY = DerivY
117        self.rms = np.std(deformation)
118        self.deformation = deformation
119        self.DerivInterp = lambda x: np.array([DerivInterpX(x[:2]), DerivInterpY(x[:2])])
120        self.SurfInterp = lambda x: SurfInterp(x[:2])
121
122    def get_normal(self, Point):
123        """
124        Calculates and returns the surface normal at the given Point.
125        It uses the derivative interpolators to compute the partial derivatives of the surface
126        deformation and returns a normalized vector representing the surface normal.
127        """
128        dX, dY = self.DerivInterp(Point)
129        norm = np.linalg.norm([dX, dY, 1])
130        dX /= norm
131        dY /= norm
132        return np.array([dX, dY, np.sqrt(1 - dX**2 - dY**2)])
133
134    def get_offset(self, Point):
135        """
136        Calculates and returns the surface offset at the given Point.
137        It uses the surface deformation interpolator to determine the displacement
138        at the specified point on the surface.
139        """
140        return self.SurfInterp(Point)
141
142    def RMS(self):
143        """
144        Returns the root mean square (RMS) value of the surface deformation (self.rms).
145        """
146        return self.rms
147
148    def PV(self):
149        pass
150
151
152class Zernike(Defect):
153    def __init__(self, Support, coefficients):
154        self.coefficients = coefficients
155        self.max_order = np.max([k[0] for k in coefficients])
156        self.support = Support
157        self.R = Support._CircumCirc()
158
159    def get_normal(self, Point):
160        x, y, z = Point / self.R
161        val, derivX, derivY = zernike_gradient([x], [y], self.max_order)
162        dX = 0
163        dY = 0
164        for k in self.coefficients:
165            dX += self.coefficients[k] * derivX[k][0][1][0]
166            dY += self.coefficients[k] * derivY[k][0][1][0]
167        dX /= self.R
168        dY /= self.R
169        return np.array([-dX, -dY, 1])
170
171    def get_offset(self, Point):
172        x, y, z = Point / self.R
173        val, derivX, derivY = zernike_gradient([x], [y], self.max_order)
174        Z = 0
175        for k in self.coefficients:
176            Z += self.coefficients[k] * val[k][0][1][0]
177        return Z
178
179    def RMS(self):
180        return np.sqrt(np.sum([i**2 for i in self.coefficients.values()]))
181
182    def PV(self):
183        pass