% -*- mode: Noweb; noweb-code-mode: python-mode -*- \section{PXD file} \label{sec:pxd-file} <>= cimport numpy as np from utilities cimport rtd, vector, cuFloatArray, cuDoubleArray from source cimport source, bundle, Bundle, Source cdef extern from "rayTracing.h": <> <> @ \section{Coordinate system} \label{sec:coordinate-system-1} \index{rayTracing!python!Coordinate\_system} \subsection{PXD file} \label{sec:pxd-file-1} <>= cdef extern from "utilities.h": cdef cppclass gpu_t[T]: T *dev_data T *host_data int N void setup() void dev_malloc() void free_dev() void dev2host() void host2dev() void reset() cdef cppclass coordinate_system: vector *origin vector *euler_angles rtd *R char *tag void update() @ \subsubsection{Class definitions} \label{sec:class-definitions} <>= cimport numpy as np cdef class Quaternion: cdef: public np.ndarray qw, qx, qy, qz cdef class Coordinates: cdef: gpu_t[vector] *_c_gpu double[:,:] data vector *v cdef void init(Coordinates,vector *) cdef class Coordinate_system: cdef: coordinate_system *_c_coordinate_system double[:,:,::1] _R_ int N public Coordinates origin, euler_angles public bytes tag <> cdef void init(Coordinate_system, coordinate_system *) @ with <>= double[:,:,:] __data_log__ int __log__, __log_idx__ @ \subsection{PYX file} \label{sec:pyx-file} <>= import numpy as np # coordinate_system cdef class Coordinate_system: """ A class for the coordinate system origins and euler angles Attributes ---------- N : int The number of coordinate systems origin : Coordinates The origin of the coordinate system(s) euler_angles : Coordinates The euler angles of the coordinates system(s) R : ndarray The rotation natrix of the coordinate systems as a [N,3,3] array """ def __cinit__(self, int N=1): self.N = N self._c_coordinate_system = new coordinate_system() self.origin = Coordinates((self.N,3)) self.euler_angles = Coordinates((self.N,3)) <> cdef void init(self, coordinate_system *cs): self._c_coordinate_system = cs self.origin.v = cs.origin self.euler_angles.v = cs.euler_angles self.tag = cs.tag def update(self, double[:,:] Txyz=None, double[:,:] Rxyz=None, **kwargs): """ Updates the coordinate systems on the device Parameters ---------- Txyz : ndarray The coordinates translation updates in the $7\times3$ matrix form Rxyz : ndarray The coordinates rotation updates in the $7\times3$ matrix form """ if Txyz is not None: self.origin[:] = Txyz if Rxyz is not None: self.euler_angles[:] = Rxyz self._c_coordinate_system.update() if self.__log__>0: k = self.__log_idx__%self.__log__ self.__log_idx__ += 1 self.__data_log__[:,:3,k] = self.origin.data self.__data_log__[:,3:,k] = self.euler_angles.data def reset(self): """ Resets the coordinate systems to its nominal position """ self.origin[:] = 0.0 self.euler_angles[:] = 0.0 self._c_coordinate_system.update() def __getitem__(self,key): if key=="origin" or key=="Txyz": return self.origin[:] elif key=="euler_angles" or key=="Rxyz": return self.euler_angles[:] else: raise KeyError("Available keys are: ""origin"",""Txyz"" or ""euler_angles"",""Rxyz""") def __setitem__(self,key,value): if key=="origin" or key=="Txyz": self.origin[:] = value elif key=="euler_angles" or key=="Rxyz": self.euler_angles[:] = value else: raise KeyError("Available keys are: ""origin"",""Txyz"" or ""euler_angles"",""Rxyz""") property R: def __get__(self): self._R_ = self._c_coordinate_system.R return np.array(self._R_) <> @ with <>= self.__log__ = 0 self.__log_idx__ = 0 self.__data_log__ = np.empty((0,0,0)) @ and <>= <> self.__data_log__ = np.zeros((self.N,6,self.__log__)) <> @ <>= property log: def __get__(self): return self.__log__ def __set__(self,value): self.__log__ = value self.__log_idx__ = 0 @ <>= property data_log: def __get__(self): return np.asarray(self.__data_log__) @ \subsubsection{Coordinates} \label{sec:coordinates} \index{rayTracing!python!Coordinates} <>= # Coordinates cdef class Coordinates: """ An interface between an array of device vectors and an array of host coordinates Examples -------- >>> import ceo >>> gmt = ceo.GMT_MX(25.5,256, M1_radial_order=8, M2_radial_order=14) The origin of M1 segment in the motion coordinate system is a Coordinates object >>> print gmt.M1.motion_CS.origin to retrieve the origin values: >>> print gmt.M1.motion_CS.origin[:] Segment 1 is moved along the y-axis of 1mm: >>> gmt.M1.motion_CS.origin[0,1] = 1e-3 All the segments are resetted to their nominal position with gmt.M1.motion_CS.origin[:] = 0 """ def __cinit__(self,tuple shape): self._c_gpu = new gpu_t[vector]() self.data = np.zeros(shape, dtype=np.float64) cdef void init(self,vector *v): self._c_gpu.host_data = self.v self._c_gpu.dev_data = v def __getitem__(self,index): cdef int k self._c_gpu.dev2host() for k in range(self.data.shape[0]): self.data[k,0] = self.v[k].x self.data[k,1] = self.v[k].y self.data[k,2] = self.v[k].z return np.array( self.data.__getitem__(index) , dtype=np.float64) def __setitem__(self,index,value): self.data.__setitem__(index,np.asarray(value, dtype=np.float64) ) cdef int k for k in range(self.data.shape[0]): self.v[k].x = self.data[k,0] self.v[k].y = self.data[k,1] self.v[k].z = self.data[k,2] #self._c_gpu.host2dev() def __repr__(self): data = np.array(self[:]) return data.__repr__() def __str__(self): max_abs_data = np.max(np.abs(self[:])) if max_abs_data!=0: exp = np.floor(np.log10(max_abs_data)) else: exp = 0 s = np.array_str(np.array(self.data)*(10**-exp),precision=3,suppress_small=True) return "10^%dX\n"%exp + s #return "@(CEO)>Coordinates:\nX,Y,Z coordinates in meter: 10^%dX\n"%exp + s property shape: def __get__(self): return self.data.shape property size: def __get__(self): return self.data.size @ \subsubsection{Quaternions} \label{sec:quaternions} \index{rayTracing!python!Quaternions} <>= # Quaternion cdef class Quaternion: def __cinit__(self, qw=np.array(0.0,ndmin=1), qx=np.array(0.0,ndmin=1), qy=np.array(0.0,ndmin=1), qz=np.array(0.0,ndmin=1), np.ndarray origins = None, np.ndarray euler_angles = None): self.qw = np.array(qw) self.qx = np.array(qx) self.qy = np.array(qy) self.qz = np.array(qz) if origins is not None: self.from_origins(origins) if euler_angles is not None: self.from_euler_angles(euler_angles) def from_euler_angles(self, np.ndarray euler_angles): if euler_angles.ndim==1: euler_angles = euler_angles[None,:] ca = np.cos(euler_angles[:,0]*0.5) sa = np.sin(euler_angles[:,0]*0.5) cb = np.cos(euler_angles[:,1]*0.5) sb = np.sin(euler_angles[:,1]*0.5) cc = np.cos(euler_angles[:,2]*0.5) sc = np.sin(euler_angles[:,2]*0.5) self.qw = ca*cb*cc + sa*sb*sc self.qx = sa*cb*cc - ca*sb*sc self.qy = ca*sb*cc + sa*cb*sc self.qz = ca*cb*sc - sa*sb*cc def from_origins(self, np.ndarray origins): if origins.ndim==1: origins = origins[None,:] self.qx = origins[:,0] self.qy = origins[:,1] self.qz = origins[:,2] self.qw = np.zeros_like(self.qx) def norm(self): red = self.qw**2 + self.qx**2 + self.qy**2 + self.qz**2 return np.sqrt(red) def conj(self): Q = Quaternion(self.qw,-self.qx,-self.qy,-self.qz) return Q def __add__(x,y): if isinstance(x,Quaternion) and isinstance(y,Quaternion): Q1 = x Q2 = y return Quaternion(Q1.qw+Q2.qw,Q1.qx+Q2.qx,Q1.qy+Q2.qy,Q1.qz+Q2.qz) else: return NotImplemented def __sub__(x,y): if isinstance(x,Quaternion) and isinstance(y,Quaternion): Q1 = x Q2 = y return Quaternion(Q1.qw-Q2.qw,Q1.qx-Q2.qx,Q1.qy-Q2.qy,Q1.qz-Q2.qz) else: return NotImplemented def __mul__(x,y): if isinstance(x,Quaternion) and isinstance(y,Quaternion): Q1 = x Q2 = y a1 = Q1.qw b1 = Q1.qx c1 = Q1.qy d1 = Q1.qz a2 = Q2.qw b2 = Q2.qx c2 = Q2.qy d2 = Q2.qz return Quaternion( a1*a2-b1*b2-c1*c2-d1*d2, a1*b2+b1*a2+c1*d2-d1*c2, a1*c2-b1*d2+c1*a2+d1*b2, a1*d2+b1*c2-c1*b2+d1*a2) else: return NotImplemented def to_euler_angles(self): q0 = self.qw q1 = self.qx q2 = self.qy q3 = self.qz a = np.arctan2(2.0*(q0*q1+q2*q3),1-2.0*(q1*q1+q2*q2)) b = np.arcsin(2.0*(q0*q2-q3*q1)) c = np.arctan2(2.0*(q0*q3+q1*q2),1-2.0*(q2*q2+q3*q3)) return np.concatenate((a[:,None],b[:,None],c[:,None]),axis=1) def __getitem__(self,index): wxyz = np.zeros((self.qw.size,4)) wxyz[:,0] = self.qw wxyz[:,1] = self.qx wxyz[:,2] = self.qy wxyz[:,3] = self.qz return wxyz.__getitem__(index) def __setitem__(self,index,value): wxyz = np.zeros((self.qw.size,4)) wxyz[:,0] = self.qw wxyz[:,1] = self.qx wxyz[:,2] = self.qy wxyz[:,3] = self.qz wxyz[index] = value self.qw = wxyz[:,0] self.qx = wxyz[:,1] self.qy = wxyz[:,2] self.qz = wxyz[:,3] def __repr__(self): return [self.qw,self.qx,self.qy,self.qz].__repr__() def __str__(self): return "@(CEO)>Quaternion:" + \ "\nqw:" + np.array_str(self.qw,precision=3,suppress_small=True) + \ "\nqx:" + np.array_str(self.qx,precision=3,suppress_small=True) + \ "\nqy:" + np.array_str(self.qy,precision=3,suppress_small=True) + \ "\nqz:" + np.array_str(self.qz,precision=3,suppress_small=True) @ \section{Conic surface} \label{sec:surface-1} \index{rayTracing!python!Conic} \subsection{PXD file} \label{sec:pxd-file-2} <>= cdef cppclass conic: rtd refractive_index, c, k vector origin vector *d__origin coordinate_system ref_frame void setup(rtd , rtd , vector , vector , vector, rtd) void setup(rtd , rtd , vector , vector , vector, rtd, int, rtd *) void cleanup() void trace(bundle *) @ \subsubsection{Class definition} \label{sec:class-definition} <>= # conic cdef class Conic: cdef: conic *_c_conic readonly Coordinate_system ref_frame readonly Coordinates origin public bint coord_break public dict material @ \subsection{PYX file} \label{sec:pyx-file-3} <>= from cython cimport boundscheck, wraparound from refractors import glass_index # conic cdef class Conic: """ A class defining a conic surface Parameters ---------- c : rtd The vertex curvature k : rtd The conic constant origin : list, optional The location of the conic surface [m]; default: [0.0,0.0,0.0] euler_angles : list, optional The rotation along the x,y and z axis [rd]; default: [0.0,0.0,0.0] conic_origin : list, optional The coordinates of the center of the conic surface; default: [0.0,0.0,0.0] material : str, optional The material type of the surface e.g. "mirror", "BK7", etc.. ; default: '' refractive_index : rtd, optional The refractive index of the surface, either 0 for a dummy surface, -1 for a reflective surface, >0 for a refractive surface; default: 0.0 coord_break : bint, optional The coordinate break surface type flag, True of False; default: False asphere_a : rtd[:], optional The coefficient of an even asphere surface; default: numpy.empty(0) Attributes ---------- ref_frame : Coordinate_system, read only The coordinate system of the surface origin : Coordinates, read only The location of the origin of the conic surface within the surface coordinate system coord_break : bint The coordinate break surface type flag, True of False; default: False material : dict The material type of the surface e.g. "mirror", "BK7", etc.. """ def __cinit__(self,rtd c, rtd k, list origin = [0.0,0.0,0.0], list euler_angles = [0.0,0.0,0.0], list conic_origin = [0.0,0.0,0.0], str material = '', rtd refractive_index = 0.0, bint coord_break = False, rtd[:] asphere_a = np.empty(0)): self._c_conic = new conic() cdef vector __origin __origin.x = origin[0] __origin.y = origin[1] __origin.z = origin[2] cdef vector __euler_angles __euler_angles.x = euler_angles[0] __euler_angles.y = euler_angles[1] __euler_angles.z = euler_angles[2] cdef vector __conic_origin __conic_origin.x = conic_origin[0] __conic_origin.y = conic_origin[1] __conic_origin.z = conic_origin[2] with boundscheck(False), wraparound(False): self._c_conic.setup(c, k, __origin, __euler_angles, __conic_origin, refractive_index, asphere_a.size, &asphere_a[0]) self.ref_frame = Coordinate_system(1) self.ref_frame.init( &(self._c_conic.ref_frame) ) self.origin = Coordinates((1,3)) self.origin.v = &(self._c_conic.origin) self.origin.init(self._c_conic.d__origin) self.coord_break = coord_break self.material = {'name': material, 'formula': 15, 'c': []} def __dealloc__(self): self._c_conic.cleanup() def __str__(self): origin = "\n . Origin: "+self.ref_frame.origin.__str__() euler_angles = "\n . Euler angles: "+self.ref_frame.euler_angles.__str__() if self.coord_break: return ". coordinate break"+origin+euler_angles else: return ". material: "+self.material["name"]+\ " [c=%g,k=%g]"%(self._c_conic.c,self._c_conic.k)+\ origin+euler_angles def trace(self, Bundle rays): """ Traces the rays of the source to the conic optical surface Parameters ---------- src : Source A Source object See also -------- Source : the celestial source class """ self._c_conic.trace(rays._c_bundle) def refractive_index(self, float wavelength=-1.0): if self.material['name'] == '': return self._c_conic.refractive_index else: if wavelength>0: return glass_index(self.material['formula'], wavelength*1e6, None,None, self.material['c']) @ \section{Aperture} \label{sec:aperture-1} \index{rayTracing!python!Aperture} \subsection{PXD file} \label{sec:pxd-file-3} <>= cdef cppclass aperture: coordinate_system ref_frame void setup(rtd, rtd, int, vector , vector) void setup_GMT_M1(rtd, int) void cleanup() void vignetting(bundle *) @ \subsection{PYX file} \label{sec:pyx-file-4} <>= # aperture cdef class Aperture: """ A class defining a circular aperture for ray tracing Parameters ---------- D : rtd The aperture diameter [m] ri : rtd The central obscuration ratio D_px : int The sampling of the diameter [pixel] origin : list, optional The location of the aperture [m]; default: [0.0,0.0,0.0] euler_angles : list, optional The rotation along the x,y and z axis [rd]; default: [0.0,0.0,0.0] Attributes ---------- ref_frame : Coordinate_system, read only The coordinate system of the aperture """ cdef: aperture *_c_aperture public Coordinate_system ref_frame def __cinit__(self, rtd D, rtd ri, int D_px, origin=[0.0,0.0,0.0], euler_angles=[0.0,0.0,0.0]): self._c_aperture = new aperture() cdef vector __origin __origin.x = origin[0] __origin.y = origin[1] __origin.z = origin[2] cdef vector __euler_angles __euler_angles.x = euler_angles[0] __euler_angles.y = euler_angles[1] __euler_angles.z = euler_angles[2] self._c_aperture.setup(D, ri, D_px, __origin, __euler_angles) self.ref_frame = Coordinate_system(1) self.ref_frame.init( &(self._c_aperture.ref_frame) ) def __dealloc__(self): self._c_aperture.cleanup() def vignetting(self, Bundle rays): """ Vignets the rays of the source according to the aperture geometry Parameters ---------- src : Source A Source object See also -------- Source : the celestial source class """ self._c_aperture.vignetting(rays._c_bundle) @ \section{Zernike} \label{sec:zernike} \index{rayTracing!python!ZernikeS} \subsection{PXD file} \label{sec:pxd-file-4} <>= cdef cppclass zernikeS: rtd *a void setup(int , rtd *, vector , vector , int) void cleanup() void surface(rtd *, rtd *, rtd *, int , int) void update(rtd *) void surface_derivatives(rtd *, rtd *, rtd *, rtd *, int ) void surface_and_derivatives(rtd *, rtd *, rtd *, rtd *, rtd *, int); void projection(float *, rtd *, rtd *, int) rtd zernike_surface(rtd , rtd , unsigned int , rtd *) @ \subsubsection{Class definition} \label{sec:class-definition-1} <>= from numpy cimport ndarray from utilities cimport cuDoubleArray cdef class ZernikeS: cdef: zernikeS *_c_zernikeS public int max_n, n_mode, N_SURF public ndarray a public ndarray Zmat, invZmat, projZmat double [:,::1] _a_ <> @ \subsection{PYX file} \label{sec:pyx-file-5} <>= # Zernike surface cdef class ZernikeS: """ A class to represent a surface as the weighted sum of Zernike polynomials Parameters ---------- max_n : int The largest of the radial orders a : ndarray, optional The Zernike polynomial coefficients as a N_SURFxn_modes array; default to all 0 N_SURF : int, optional The number of Zernike surfaces; default: 1 origin : list, optional Origin of the coordinate system of the Zernike polynomials; default: [0,0,0] euler_angles : list, optional Euler angles of the coordinate system of the Zernike polynomials; default: [0,0,0] Attributes ---------- max_n : int The largest of the radial orders n_mode : int The total number of Zernike modes N_SURF : int The number of Zernike surfaces a : ndarray The Zernike polynomial coefficients as a N_SURFxn_modes array Examples -------- >>> import ceo >>> import numpy as np >>> import matplotlib as plt Lets define a Zernike surface made of the first 15 Zernike polynomials (4th radial order) >>> Z = ceo.ZernikeS(4) and setting it to be a pure focus mode: >>> Z.a[3] = 1.0 >>> Z.update() The surface is computed with >>> npx = 512 >>> u = np.linspace(-1,1,npx) >>> x,y = np.meshgrid(u,u, indexing='xy') >>> r = np.hypot(x,y) >>> o = np.arctan2(y,x) >>> cuo = ceo.cuDoubleArray(host_data=o) >>> cur = ceo.cuDoubleArray(host_data=r) >>> S = Z.surface(cur,cuo) and its derivative with >>> (dSdx,dSdy) = Z.surface_derivatives(cur,cuo) >>> fig, (ax1,ax2,ax3) = plt.subplots(ncols=3,sharey=True) >>> fig.set_size_inches(20,4.5) >>> h1 = ax1.imshow(S.host(),interpolation=None,origin='lower') >>> plt.colorbar(h1,ax=ax1) >>> h2 = ax2.imshow(dSdx.host(),interpolation=None,origin='lower') >>> plt.colorbar(h2,ax=ax2) >>> h3 = ax3.imshow(dSdy.host(),interpolation=None,origin='lower') >>> plt.colorbar(h3,ax=ax3) """ def __cinit__(self, int _max_n_, ndarray _a_=None, origin=[0.0,0.0,0.0], euler_angles=[0.0,0.0,0.0], int N_SURF=1): self._c_zernikeS = new zernikeS() self.max_n = _max_n_ self.n_mode = (_max_n_+1)*(_max_n_+2)/2 self.N_SURF = N_SURF if _a_ is None: self.a = np.zeros((self.N_SURF,self.n_mode), dtype=np.float64) else: self.a = _a_ cdef vector __origin __origin.x = origin[0] __origin.y = origin[1] __origin.z = origin[2] cdef vector __euler_angles __euler_angles.x = euler_angles[0] __euler_angles.y = euler_angles[1] __euler_angles.z = euler_angles[2] self._c_zernikeS.setup( self.max_n, self.a.data, __origin, __euler_angles, self.N_SURF) <> def __dealloc__(self): self._c_zernikeS.cleanup(); def update(self, rtd[:,:] modes=None, **kwargs): """ Updates the Zernike surface based on the Zernike coefficients """ cdef: double[:,:] __a__ if modes is not None: self.a[:] = modes self._c_zernikeS.update(self.a.data) if self.__log__>0: k = self.__log_idx__%self.__log__ self.__log_idx__ += 1 __a__ = self.a self.__data_log__[:,:,k] = __a__ def reset(self): """ Resets the Zernike coefficients to zero and update the Zernike surface """ self.a = np.zeros((self.N_SURF,self.n_mode)) self._c_zernikeS.update(self.a.data) def __invert__(self): self.reset() def __ixor__(self, tuple t): self.a[t[0],t[1]] = t[2] self.update() return self def surface(self, cuDoubleArray r, cuDoubleArray o, int surf_id=0, int[:] surf_ids=None, cuDoubleArray out=None): """ Computes the Zernike surface in polar coordinates Parameters ---------- r : cuDoubleArray The normalized radius o : cuDoubleArray The azimuth surf_id : int, optional The surface index; default: 0 surf_ids : int[:], optional The surface index vector; default: None out : cuDoubleArray, optional The Zernike surface boutput array; default: None Returns ------- S : cuDoubleArray The Zernike surface """ cdef int k cdef cuDoubleArray S if out is None: S = cuDoubleArray(shape=r.shape, dev_malloc = True) self._c_zernikeS.surface(S._c_gpu.dev_data, r._c_gpu.dev_data, o._c_gpu.dev_data, r._c_gpu.N, surf_id) return S else: if surf_ids is None: surf_ids = np.arange(self.N_SURF,dtype=np.int32) assert out.size>=surf_ids.size*r.size, "out array is too small!" for k in range(surf_ids.size): self._c_zernikeS.surface(out._c_gpu.dev_data + r.size*k, r._c_gpu.dev_data, o._c_gpu.dev_data, r._c_gpu.N, surf_ids[k]) def surface_derivatives(self, cuDoubleArray r, cuDoubleArray o): """ Computes the Zernike surface x and y derivatives in polar coordinates Parameters ---------- r : cuDoubleArray The normalized radius o : cuDoubleArray The azimuth Returns ------- dSdx : cuDoubleArray The Zernike surface x derivative dSdy : cuDoubleArray The Zernike surface y derivative """ cdef cuDoubleArray dSdx, dSdy dSdx = cuDoubleArray(shape=r.shape, dev_malloc = True) dSdy = cuDoubleArray(shape=r.shape, dev_malloc = True) self._c_zernikeS.surface_derivatives(dSdx._c_gpu.dev_data, dSdy._c_gpu.dev_data, r._c_gpu.dev_data, o._c_gpu.dev_data, r._c_gpu.N) return dSdx,dSdy def surface_and_derivatives(self, cuDoubleArray r, cuDoubleArray o): """ Computes the Zernike surface x and y derivatives in polar coordinates Parameters ---------- r : cuDoubleArray The normalized radius o : cuDoubleArray The azimuth Returns ------- S : cuDoubleArray The Zernike surface dSdx : cuDoubleArray The Zernike surface x derivative dSdy : cuDoubleArray The Zernike surface y derivative """ cdef cuDoubleArray S, dSdx, dSdy S = cuDoubleArray(shape=r.shape, dev_malloc = True) dSdx = cuDoubleArray(shape=r.shape, dev_malloc = True) dSdy = cuDoubleArray(shape=r.shape, dev_malloc = True) self._c_zernikeS.surface_and_derivatives(S._c_gpu.dev_data, dSdx._c_gpu.dev_data, dSdy._c_gpu.dev_data, r._c_gpu.dev_data, o._c_gpu.dev_data, r._c_gpu.N) return S,dSdx,dSdy def projection(self, cuFloatArray phase, cuDoubleArray r, cuDoubleArray o): """ Projects the phase onto a set on Zernike polynomials Parameters ---------- phase : cuFloatArray The wavefront phase r : cuDoubleArray The radius of the coordinates o : cuDoubleArray The azimuth of the coordinates """ self._c_zernikeS.projection(phase._c_gpu.dev_data, r._c_gpu.dev_data, o._c_gpu.dev_data, r.size) self._a_ = self._c_zernikeS.a self.a = np.asarray(self._a_) def fitting_init(self, Source gs, int alphaId=0): P = np.rollaxis( np.array(gs.rays.piston_mask ),0,3) ## Find center coordinates (in pixels) of each segment mask u = np.arange(gs.n) v = np.arange(gs.m) x,y = np.meshgrid(u,v) x = x.reshape(1,-1,1) y = y.reshape(1,-1,1) xc = np.sum(x*P,axis=1)/P.sum(axis=1) yc = np.sum(y*P,axis=1)/P.sum(axis=1) ## Preliminary estimation of radius (in pixels) of each segment mask (assuming that there is no central obscuration) Rs = np.sqrt(P.sum(axis=1)/np.pi) ## Polar coordinates rho = np.hypot( x - xc[:,np.newaxis,:], y - yc[:,np.newaxis,:]) #temporal rho vector theta = np.arctan2( y - yc[:,np.newaxis,:], x - xc[:,np.newaxis,:]) * P ## Estimate central obscuration area of each segment mask ObsArea = np.sum(rho < 0.9*Rs[:,np.newaxis,:] * ~P.astype('bool'), axis=1) ## Improve estimation of radius of each segment mask Rs = np.sqrt( (P.sum(axis=1)+ObsArea) / np.pi) ## Normalize rho vector (unitary radius) rho = rho / Rs[:,np.newaxis,:] * P #final rho vector # Build a Zernike Influence-function Matrix for all segments #alphaId = 0 # only on-axis direction supported... self.Zmat = np.zeros((gs.n*gs.m,self.n_mode,7)) for segId in range(1,8): self.reset() cutheta = cuDoubleArray(host_data=theta[segId-1,:,alphaId].reshape(gs.m,gs.n)) curho = cuDoubleArray(host_data= rho[segId-1,:,alphaId].reshape(gs.m,gs.n)) for k in range(self.n_mode): self.a[0,k] = 1 self.update() S = self.surface(curho,cutheta).host(shape=(gs.m*gs.n,1))*P[segId-1,:,alphaId].reshape(-1,1) self.Zmat[:,k,segId-1] = S.flatten() self.a[0,k] = 0 print 'Zernike Inflence Function Matrix:' print np.shape(self.Zmat) #Pseudo-inverse of Zmat self.invZmat = np.zeros((self.n_mode,gs.m*gs.n,7)) for segId in range(1,8): self.invZmat[:,:,segId-1] = np.linalg.pinv(self.Zmat[:,:,segId-1]) #print 'inverse of Zernike Influence Function Matrix:' #print self.invZmat.shape #Projection matrix self.projZmat = np.zeros((self.n_mode,gs.m*gs.n,7)) for segId in range(1,8): self.projZmat[:,:,segId-1] = self.Zmat[:,:,segId-1].T / np.sum(P[segId-1,:,alphaId]) def projectionS(self, phase): cdef: int segId a = np.zeros((self.n_mode,7)) for segId in range(1,8): a[:,segId-1] = np.dot(self.projZmat[:,:,segId-1], phase.reshape(-1)) return a def fitting(self, phase): #Project final residual WF onto segment Zernike modes cdef: int segId arec = np.zeros((self.n_mode, 7)) for segId in range(1,8): arec[:,segId-1] = np.dot(self.invZmat[:,:,segId-1], phase.reshape(-1)) return arec <> def Zernike_Surface(rtd r, rtd o, int max_n, np.ndarray a): return zernike_surface(r, o, max_n, a.data) @ with @ and <>= <> self.__data_log__ = np.zeros((self.N_SURF,self.n_mode,self.__log__), dtype=np.float64) <> @ \section{Ray tracing routines} \label{sec:ray-tracing-routines} \subsection{PXD file} \label{sec:pxd-file-5} <>= cdef void transform_to_S(bundle *, conic *) cdef void transform_to_S(bundle *, aperture *) cdef void transform_to_R(bundle *, conic *) cdef void transform_to_R(bundle *, aperture *) cdef void intersect(bundle *, conic *) cdef void reflect(bundle *) cdef void refract(bundle *, rtd) cdef void thin_lens(bundle *) @ \subsection{PYX file} \label{sec:pyx-file-6} \subsubsection{Transform\_to\_S} \label{sec:transform_to_s} \index{rayTracing!python!Transform\_to\_S} <>= # ray tracing def Transform_to_S(Bundle rays, Conic F): """ Transforms the rays of the source into the coordinate system of the conic object Parameters ---------- src : Source A Source object F : Conic A Conic object See also -------- Source : the celestial source class Conic : the optics conic surface class """ transform_to_S(rays._c_bundle, F._c_conic) #def Transform_to_S_from_A(Source src, Aperture A): # transform_to_S(src.rays._c_bundle, A._c_aperture) @ \subsubsection{Transform\_to\_R} \label{sec:transform_to_r} \index{rayTracing!python!Transform\_to\_R} <>= def Transform_to_R(Bundle rays, Conic F): """ Transforms the rays of the source into the coordinate system of the rays This is the inverse of Transform_to_S Parameters ---------- src : Source A Source object F : Conic A Conic object See also -------- Source : the celestial source class Conic : the optics conic surface class """ transform_to_R(rays._c_bundle, F._c_conic) #def Transform_to_R_from_A(Source src, Aperture A): # transform_to_R(src.rays._c_bundle, A._c_aperture) @ \subsubsection{Intersect} \label{sec:intersect} \index{rayTracing!python!Intersect} <>= def Intersect(Bundle rays, Conic F): """ Intersection of the rays of the source with the surface represented by the conic object Parameters ---------- src : Source A Source object F : Conic A Conic object See also -------- Source : the celestial source class Conic : the optics conic surface class """ intersect(rays._c_bundle, F._c_conic) @ \subsubsection{Reflect} \label{sec:reflect} \index{rayTracing!python!Reflect} <>= def Reflect(Bundle rays): """ Reflection of the rays of the source from the surface represented by the conic object Parameters ---------- src : Source A Source object See also -------- Source : the celestial source class """ reflect(rays._c_bundle) @ \subsubsection{Refract} \label{sec:refract} \index{rayTracing!python!Refract} <>= def Refract(Bundle rays, rtd mu): """ Refraction of the rays of the source through the surface represented by the conic object Parameters ---------- src : Source A Source object mu : rtd The ratio of the index of refraction associated to the incident and refracted rays See also -------- Source : the celestial source class """ refract(rays._c_bundle, mu) @ \subsubsection{Trace} \label{sec:trace} \index{rayTracing!python!Trace} <>= from ceo cimport GMT_M1, GMT_M2, GmtMirrors from GMTLIB import GMT_MX def Trace(Bundle rays, list S, bint global_CS=True, float wavelength=550e-9): cdef: int k, n rtd n_S, mu list xyz Conic _S_ n = len(S) #src.reset() xyz = [ rays.coordinates.host() ] for k in range(n): #print 'Material refractive index: %f'%src.rays.refractive_index if isinstance(S[k],Aperture): S[k].vignetting(rays) elif isinstance(S[k],(GMT_M1,GMT_M2)): S[k].trace(rays) elif isinstance(S[k],(GmtMirrors)): S[k].M2.blocking(rays) S[k].M1.trace(rays) S[k].M2.trace(rays) else: _S_ = S[k] Transform_to_S(rays,_S_) if not _S_.coord_break: Intersect(rays,_S_) n_S = _S_.refractive_index(wavelength) if n_S!=0: if n_S==-1: Reflect(rays) else: mu = rays.refractive_index/n_S if mu!=1.0: Refract(rays,mu) rays.refractive_index = n_S if global_CS: Transform_to_R(rays,_S_) xyz.append( rays.coordinates.host() ) return xyz #chief() @ <>= _S_ = S[idx-1] ceo.Transform_to_S(src,_S_) if not _S_.coord_break: ceo.Intersect(src,_S_) n_S = _S_.refractive_index if n_S!=0: if n_S==-1: ceo.Reflect(src) else: _mu_ = src.rays.refractive_index/n_S ceo.Refract(src,_mu_) src.rays.refractive_index = n_S