% -*- mode: Noweb; noweb-code-mode: python-mode -*- \section{PXD file} \label{sec:pxd-file} <>= from numpy cimport ndarray from utilities cimport mask, cuFloatArray, cuFloatComplexArray, cuDoubleArray, cuIntArray, MaskAbstract, gpu_t, float2 from rayTracing cimport Coordinates cdef extern from "utilities.h": ctypedef double rtd ctypedef struct vector: rtd x rtd y rtd z cdef extern from "source.h": ctypedef struct ray: vector coordinates vector directions vector surface_normal rtd optical_path_length <> <> @ \section{Ray bundle} \label{sec:ray-bundle-1} \index{source!python!Bundle} \subsection{PXD file} \label{sec:pxd-file-1} <>= cdef cppclass bundle: int N_RAY, N_BUNDLE, N_L int *d__piston_mask rtd refractive_index, L rtd *d__sphere_radius rtd *d__sphere_distance double rot_angle vector *d__origin vector *d__sphere_origin ray *d__ray ray *d__chief_ray void setup(rtd, int, int , vector, int) void setup(rtd, int , vector, int) void setup_free(int, double *, double *, vector) void setup_free(double, double, int, double *, double *, vector) void cleanup() void to_z_plane(rtd) void to_sphere(vector) void to_sphere(rtd , rtd) void to_focal_plane(rtd , rtd) void get_coordinates(double *) void get_chief_coordinates(double *) void get_sphere_origins(double *) void get_directions(double *) void get_chief_directions(double *) void get_optical_path_length(double *) void get_chief_optical_path_length(double *) void get_optical_path_difference(double *) void get_optical_path_difference(double *, float , int , float , int ) void get_vignetting(double *) void gmt_truss_onaxis() void gmt_m2_baffle() # static routines void get_n_iteration(int *) @ \subsubsection{Class definition} \label{sec:class-definition} <>= # Bundle cdef class Bundle: cdef: bundle *_c_bundle cuDoubleArray _coordinates_ cuDoubleArray _chief_coordinates_ cuDoubleArray _directions_ cuDoubleArray _chief_directions_ cuDoubleArray _optical_path_length_ cuDoubleArray _chief_optical_path_length_ cuDoubleArray _optical_path_difference_ cuDoubleArray _vignetting_ cuDoubleArray _sphere_origins_ cuIntArray _n_iteration_ public cuDoubleArray sphere_radius public cuDoubleArray sphere_distance public cuIntArray _piston_mask_ @ \subsection{PYX file} \label{sec:pyx-file} <>= # polar bundle from scipy.optimize import brentq from ceo.constants import ARCSEC2RAD cdef class Bundle: """ A class to represent a bundle of rays This class is usually instanciated by the Source class Parameters ---------- src : Source The Source object that contains the pointer to CEO bundle structure Attributes ---------- N_BUNDLE : int, read only The number of ray bundle coordinates : cuDoubleArray, read only The coordinate vectors of the rays as a (N_RAYxN_BUNDLE)x3 array chief_coordinates : cuDoubleArray, read only The coordinate vector of the chief ray directions : cuDoubleArray, read only The direction vectors of the rays as a (N_RAYxN_BUNDLE)x3 array chief_directions : cuDoubleArray, read only The direction vector of the chief ray optical_path_length : cuDoubleArray, read only The optical path lengths of the rays as a N_BUNDLExN_RAY array chief_optical_path_length : cuDoubleArray, read only The optical path length of the chief ray optical_path_difference : cuDoubleArray, read only The optical path differences of the rays as a N_BUNDLExN_RAY array vignetting : cuDoubleArray, read only The vignetting flag of the rays as a N_BUNDLExN_RAY array sphere_radius : cuDoubleArray, read only The radius of the reference sphere used to compute the optical path difference sphere_distance : cuDoubleArray, read only The distance from the last surface to the reference sphere used to compute the optical path difference piston_mask : numpy ndarray A 7 columns array where each column is a mask corresponding to 1 segment L : rtd, read only The size in meter of the box encompassing the ray bundle N_L : int, read only The sampling of the box encompassing the ray bundle (the number of rays is N_L$\times$N_L """ def __cinit__(self, *args, **kwargs): self._c_bundle = new bundle() def __init__(self): self._coordinates_ = cuDoubleArray( shape=(self._c_bundle.N_RAY*self._c_bundle.N_BUNDLE,3),dev_malloc=True) self._chief_coordinates_ = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,3),dev_malloc=True) self._directions_ = cuDoubleArray( shape=(self._c_bundle.N_RAY*self._c_bundle.N_BUNDLE,3),dev_malloc=True) self._chief_directions_ = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,3),dev_malloc=True) self._optical_path_length_ = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,self._c_bundle.N_RAY),dev_malloc=True) self._chief_optical_path_length_ = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,1),dev_malloc=True) self._optical_path_difference_ = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,self._c_bundle.N_RAY),dev_malloc=True) self._vignetting_ = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,self._c_bundle.N_RAY),dev_malloc=True) self._sphere_origins_ = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,3),dev_malloc=True) self._n_iteration_ = cuIntArray( shape=(self._c_bundle.N_BUNDLE,self._c_bundle.N_RAY),dev_malloc=True) self.sphere_radius = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,1),dev_malloc=False) self.sphere_radius._c_gpu.dev_data = self._c_bundle.d__sphere_radius self.sphere_distance = cuDoubleArray( shape=(self._c_bundle.N_BUNDLE,1),dev_malloc=False) self.sphere_distance._c_gpu.dev_data = self._c_bundle.d__sphere_distance self._piston_mask_ = cuIntArray(shape=(self._c_bundle.N_BUNDLE,self._c_bundle.N_RAY)) self._piston_mask_._c_gpu.dev_data = self._c_bundle.d__piston_mask; <> @ This a special [[Bundle]] class that is solely used by the [[Source]] class \index{source!python!SourceBundle} <>= cdef class SourceBundle(Bundle): """ A class to represent a bundle of rays This class is instanciated by the Source class Parameters ---------- src : Source The Source object that contains the pointer to CEO bundle structure See also -------- Bundle: the parent class for all ray bundle classes """ def __init__(self,Source src): self._c_bundle = &(src._c_source.rays) super(SourceBundle,self).__init__() @ This ray bundle used user--defined ray coordinates: \index{source!python!FreeBundle} <>= cdef class FreeBundle(Bundle): """ A class to represent a bundle of rays Parameters ---------- x : ndarray the x coordinates of the rays in the entrance pupil y : ndarray the y coordinates of the rays in the entrance pupil origin : list the point of origin of the chief ray See also -------- Bundle: the parent class for all ray bundle classes """ def __init__(self, ndarray x, ndarray y, list origin, double zenith=0.0, double azimuth=0.0): cdef vector __origin __origin.x = origin[0] __origin.y = origin[1] __origin.z = origin[2] self._c_bundle.setup_free(zenith, azimuth, x.size, x.data, y.data, __origin) super(FreeBundle,self).__init__() def __dealloc__(self): self._c_bundle.cleanup() @ <>= """ def __dealloc__(self): self._c_bundle.cleanup() """ def to_z_plane(self,rtd z_plane): """ Propagates the rays to a given plane Parameters ---------- z_plane : double The z coordinates of the plane """ self._c_bundle.to_z_plane(z_plane) def to_focal_plane(self, rtd focal_plane_distance, rtd focal_plane_radius): """ Propagates the rays to the focal plane Parameters ---------- focal_plane_distance : double, optional The location of the focal plane on the optical axis focal_plane_radius : double, optional The radius of curvature of the focal plane """ self._c_bundle.to_focal_plane(focal_plane_distance, focal_plane_radius) def to_sphere(self, sphere_origin=None, rtd focal_plane_distance=0.0, rtd focal_plane_radius=0.0): """ Computes the optical path difference with respect to the reference sphere Parameters ---------- sphere_origin : list of double, optional A 3 elements list containing the x,y and z coordinates of the sphere center; default: None focal_plane_distance : double, optional The location of the focal plane on the optical axis; default: 0.0 focal_plane_radius : double, optional The radius of curvature of the focal plane; default: 0.0 """ cdef vector origin if sphere_origin is not None: origin.x = sphere_origin[0] origin.y = sphere_origin[1] origin.z = sphere_origin[2] self._c_bundle.to_sphere(origin) if focal_plane_distance!=0.0: self._c_bundle.to_sphere(focal_plane_distance, focal_plane_radius) def to_focal_plane(self, rtd focal_plane_distance): self._c_bundle.to_z_plane(focal_plane_distance) def get_optical_path_difference(self, float delta_x, int N_x, float delta_y, int N_y): opd = cuDoubleArray(shape=(N_x,N_y),dev_malloc=True) self._c_bundle.get_optical_path_difference(opd._c_gpu.dev_data, delta_x, N_x, delta_y, N_y) return opd def ee80(self,spaxel="circle"): """ Computes the size of the geometric EE80 square Parameters ---------- spaxel : char, optional The spaxel geometry, square or circle; default: circle Returns ------- float The size of the geometric EE80 square in meter """ assert(spaxel=="circle" or spaxel=="square","spaxel is either a circle or a square!") V = self.vignetting.host().flatten()==1 xyz = self.coordinates.host() x = np.abs( xyz[V,0] ) y = np.abs( xyz[V,1] ) r = np.hypot(x,y) def ensquared(hd): idx = np.logical_and(x<=hd , y<=hd) return float(idx.sum())/V.sum() def encircled(hd): idx = r<=hd return float(idx.sum())/V.sum() if spaxel=="circle": f = lambda x: encircled(x)-0.8 b = r.max() else: f = lambda x: ensquared(x)-0.8 b = np.maximum(x.max(),y.max()) try: return 2*brentq(f,0,b) except ValueError: return np.float('inf') def gmt_truss_onaxis(self): self._c_bundle.gmt_truss_onaxis() def gmt_m2_baffle(self): self._c_bundle.gmt_m2_baffle() property N_BUNDLE: def __get__(self): return self._c_bundle.N_BUNDLE property origin: def __get__(self): cdef gpu_t[vector] *_c_gpu cdef vector origin _c_gpu = new gpu_t[vector]() _c_gpu.setup(1) _c_gpu.dev_data = self._c_bundle.d__origin _c_gpu.host_data = &origin _c_gpu.dev2host() return origin.x,origin.y,origin.z def __set__(self,val): cdef gpu_t[vector] *_c_gpu cdef vector origin origin.x = val[0] origin.y = val[1] origin.z = val[2] _c_gpu = new gpu_t[vector]() _c_gpu.setup(1) _c_gpu.dev_data = self._c_bundle.d__origin _c_gpu.host_data = &origin _c_gpu.host2dev() property rot_angle: def __get__(self): return self._c_bundle.rot_angle def __set__(self,val): self._c_bundle.rot_angle = val property N_L: def __get__(self): return self._c_bundle.N_L property coordinates: def __get__(self): self._c_bundle.get_coordinates(self._coordinates_._c_gpu.dev_data) return self._coordinates_ property chief_coordinates: def __get__(self): self._c_bundle.get_chief_coordinates(self._chief_coordinates_._c_gpu.dev_data) return self._chief_coordinates_ property directions: def __get__(self): self._c_bundle.get_directions(self._directions_._c_gpu.dev_data) return self._directions_ property chief_directions: def __get__(self): self._c_bundle.get_chief_directions(self._chief_directions_._c_gpu.dev_data) return self._chief_directions_ property optical_path_length: def __get__(self): self._c_bundle.get_optical_path_length(self._optical_path_length_._c_gpu.dev_data) return self._optical_path_length_ property chief_optical_path_length: def __get__(self): self._c_bundle.get_chief_optical_path_length(self._chief_optical_path_length_._c_gpu.dev_data) return self._chief_optical_path_length_ property optical_path_difference: def __get__(self): self._c_bundle.get_optical_path_difference(self._optical_path_difference_._c_gpu.dev_data) return self._optical_path_difference_ property vignetting: def __get__(self): self._c_bundle.get_vignetting(self._vignetting_._c_gpu.dev_data) return self._vignetting_ property sphere_origins: def __get__(self): self._c_bundle.get_sphere_origins(self._sphere_origins_._c_gpu.dev_data) return self._sphere_origins_ property n_iteration: def __get__(self): self._c_bundle.get_n_iteration(self._n_iteration_._c_gpu.dev_data) return self._n_iteration_ property piston_mask: def __get__(self): P = self._piston_mask_.host() return [ np.array( [P[kk,:].flatten()==k for k in range(1,8)] ) for kk in range(self._piston_mask_.shape[0])] property refractive_index: def __get__(self): return self._c_bundle.refractive_index def __set__(self,rtd value): self._c_bundle.refractive_index = value property L: def __get__(self): return self._c_bundle.L @ \section{Complex amplitude} \label{sec:complex-amplitude-1} \index{source!python!Complex\_amplitude} \subsection{PXD file} \label{sec:pxd-file-2} <>= cdef cppclass complex_amplitude: int N int N_PX float *phase float *amplitude mask *M void reset() void reset(complex_amplitude *) void reset_phase() void reset_phase(complex_amplitude *) void masked(mask *) void masked() void add_phase(float , float *) void rms(float *) void finite_difference(float *, float *, int, float); void finite_difference(float *, float *, int, float, mask *); void gradient_average(float *, float *, int , float) void gradient_average(float *, float *, float) void segments_gradient_average(float *, float * , float, int *) void segments_gradient_averageFast(float *, float * , float, int *) @ \subsubsection{Class definition} \label{sec:class-definition-1} <>= cdef class Complex_amplitude: cdef: complex_amplitude *_c_complex_amplitude readonly cuFloatArray phase readonly cuFloatArray amplitude @ \subsection{PYX file} \label{sec:pyx-file-1} <>= import numpy as np cimport numpy as np cimport cython @cython.boundscheck(False) @cython.wraparound(False) cdef class Complex_amplitude: """ A class to represent a wavefront complex amplitude Parameters ---------- Source : Source a CEO Source object Attributes ----------- phase : float The complex amplitude angle amplitude : float, read only the complex amplitude magnitude See also -------- Source: a class for astronomical sources """ def __cinit__(self,Source src=None): self._c_complex_amplitude = new complex_amplitude() if src is not None: self._c_complex_amplitude = &(src._c_source.wavefront) self.__alloc__() def __alloc__(self, tuple shape=None): if shape is None: shape = (self._c_complex_amplitude.N, self._c_complex_amplitude.N_PX/\ self._c_complex_amplitude.N) self.phase = cuFloatArray(shape=shape) self.phase._c_gpu.dev_data = self._c_complex_amplitude.phase self.amplitude = cuFloatArray(shape=shape) self.amplitude._c_gpu.dev_data = self._c_complex_amplitude.amplitude def rms(self,int units_exponent=0): """ Computes the rms of the wavefront phase Parameters ---------- units_exponent: int, optional 3 for km, 0 for meter, -2 for cm, -3 for mm, -6 for micron, -9 for nm, etc; default: 0 Returns ------- float the wavefront phase rms """ cdef float[::1] rms = np.zeros(self._c_complex_amplitude.N, dtype=np.single) self._c_complex_amplitude.rms(&rms[0]) return np.array(rms)*10**-units_exponent def reset(self, cuFloatArray phase = None, Complex_amplitude wavefront=None): """ Reset the wavefront amplitude to 1 and phase to 0 or to the given phase Parameters ---------- phase : cuFloatArray, optional The GPU array the wavefront phase is reset to; default to None See also -------- cuFloatArray : a class for GPU host and device float data """ if wavefront is None: self._c_complex_amplitude.reset() if phase is not None: self._c_complex_amplitude.add_phase(1,phase._c_gpu.dev_data) else: self._c_complex_amplitude.reset(wavefront._c_complex_amplitude) def reset(self, Complex_amplitude wavefront=None): if wavefront is None: self._c_complex_amplitude.reset_phase() else: self._c_complex_amplitude.reset_phase(wavefront._c_complex_amplitude) def addPhase(self, cuFloatArray phase): self._c_complex_amplitude.add_phase(1,phase._c_gpu.dev_data) def axpy(self, float alpha, cuFloatArray phase): self._c_complex_amplitude.add_phase(alpha,phase._c_gpu.dev_data) def __iadd__(self, x): cdef cuFloatArray data if isinstance(x,(np.ndarray,list,tuple,float)): data = cuFloatArray(host_data=x) if isinstance(x,cuDoubleArray): data = x.__single__() if isinstance(x,cuFloatArray): data = x self._c_complex_amplitude.add_phase(1,data._c_gpu.dev_data) return self def __isub__(self, x): cdef cuFloatArray data if isinstance(x,(np.ndarray,list,tuple,float)): data = cuFloatArray(host_data=x) if isinstance(x,cuDoubleArray): data = x.__single__() if isinstance(x,cuFloatArray): data = x self._c_complex_amplitude.add_phase(-1,data._c_gpu.dev_data) return self def finiteDifference(self, int NL, float d, MaskAbstract M=None): """ Computes the average finite difference of the wavefront Parameters ---------- NL : int The linear size of the lenslet array d : float The lenslet pitch M : MaskAbstract, optional The valid lenslet mask Returns ------- cuFloatArray The wavefront finite differences in an array of size [N_SRCxNLx2,NL] """ cdef cuFloatArray sxy sxy = cuFloatArray(shape=(self._c_complex_amplitude.N*NL*2,NL), dev_malloc=True) if M is None: self._c_complex_amplitude.finite_difference(sxy._c_gpu.dev_data, sxy._c_gpu.dev_data + NL*NL, NL, d) else: self._c_complex_amplitude.finite_difference(sxy._c_gpu.dev_data, sxy._c_gpu.dev_data + NL*NL, NL, d, M._c_mask) return sxy def gradientAverage(self, int NL=1, float d=1.0): """ Computes the average finite difference of the wavefront Parameters ---------- NL : int The linear size of the lenslet array d : float The lenslet pitch Returns ------- cuFloatArray The wavefront gradient average in an array of size [N_SRCxNLx2,NL] """ cdef cuFloatArray gsxy cdef float[:] sxy if NL>1: gsxy = cuFloatArray(shape=(self._c_complex_amplitude.N*NL*2,NL), dev_malloc=True) self._c_complex_amplitude.gradient_average(gsxy._c_gpu.dev_data, gsxy._c_gpu.dev_data + NL*NL, NL, d) return gsxy else: sxy = np.zeros(2*self._c_complex_amplitude.N,dtype=np.single) self._c_complex_amplitude.gradient_average(&sxy[0], &sxy[self._c_complex_amplitude.N], d) return np.asarray(sxy) def gradientAverageFast(self, float d): """ Computes the average finite difference of the wavefront Parameters ---------- d : float The lenslet pitch Returns ------- cuFloatArray The wavefront gradient average in an array of size [N_SRCx2,1] """ cdef float[:] sxy sxy = np.zeros(2*self._c_complex_amplitude.N,dtype=np.single) self._c_complex_amplitude.gradient_average(&sxy[0], &sxy[self._c_complex_amplitude.N], d) return np.asarray(sxy) @ \section{Source} \label{sec:source} \index{source!python!Source} \subsection{PXD file} \label{sec:pxd-file-3} <>= cdef cppclass source: int N_SRC const char *photometric_band float fwhm, magnitude float2 *d__otf complex_amplitude wavefront bundle rays void setup(const char *,float , float , float ) void setup(const char *,float , float , float , int ) void setup(const char *,float *, float *, float , int , int ) void setup(const char *,float *, float *, float , int , rtd, int, vector) void setup(const char *,float *,double *, double *, float , int , rtd, int, vector) void setup_chief(const char *,float *,double *, double *, float , int , rtd, int, vector, vector) void setup(const char *,float *,double *, double *, float , int , int, double *, double *, vector) void cleanup() void reset_rays() void reset_rays(int) void opd2phase() void opd2phase(int) void info() float n_photon() float n_photon(float) float n_background_photon(float) float wavelength() float spectral_bandwidth() void update_directions(double *, double *, int) void update_magnitude(float *, int) void optical_transfer_function(float2 *) @ \subsubsection{Class definition} \label{sec:class-definition-2} <>= from cpython.object cimport Py_GT cdef class Source: cdef: source *_c_source readonly int N_SRC readonly int n, m readonly int size double[:] _zenith_, _azimuth_ float[:] _magnitude_ readonly float height readonly Bundle rays readonly Complex_amplitude wavefront readonly sphere_distance readonly tuple OPTICALPATH public double samplingTime, timeStamp public dict parameters public bytes band @ \subsection{PYX file} \label{sec:pyx-file-2} <>= cdef class Source: """ A class to represent an astronomical source Parameters ---------- photometric_band : char The sources photometric band zenith : list, tuple or numpy array, optional The sources zenith angles [rd], defaults to 0. azimuth : list, tuple or numpy array, optional The sources azimuth angles [rd], defaults to 0. height : float, optional The sources altitude [m], defaults to infinity. resolution : tuple A 2 element tuple with the sampling [n,m] in pixel of the source wavefront complex amplitude. fwhm : float, optional The fwhm of the source intensity distribution in detector pixel unit (before binning), defaults to None. rays_box_size : float, optional The size of the ray bundle [m], defaults to 25.5. rays_box_sampling : int, optional The linear sampling of the ray bundle rays_origin : list, optional A 3 element list, with the first 2 being the (x,y) origin of the ray bundle and the 3rd the ray bundle starting z coordinates, defaults to [0,0,25] samplingTime : double the source propagation sampling time; default: 1s Attributes ---------- n : int wavefront resolution m : int wavefront resolution size : int number of sources zenith : ndarray zenith angle azimuth : ndarray azimuth angle thetax : ndarray tan(zenith)cos(azimuth) thetay : ndarray tan(zenith)sin(azimuth) height : float source height ray : Bundle geometric ray bundle phase : cuFloatArray wavefront phase wavelength : float, read only the wavelength of the source spectral_bandwidth : float, read only the spectral bandwidth associated to the source nPhoton : float, read only number of photon [m^-2.s^-1] magnitude : float, write only star magnitude wavefront : Complex_amplitude the wavefront complex amplitude sphere_distance : float the distance to the reference sphere used to compute the OPD OPTICALPATH : tuple (readonly) the source optical path samplingTime : double the source propagation sampling time; default: 1s timeStamp : double the time associated with the current propagation, the "samplingTime" is added to "timeStamp" after each call to the unary + operator parameters : dict the object parameters in a dictionary Examples -------- >>> import math >>> import numpy as np >>> import ceo An on-axis source in K band is simply defined with: >>> src = Source("K") The resolution of the source wavefront is specifed with >>> n = 256 >>> src = Source("K",resolution=(n,n)) For a 20" off-axis and 10th magnitude source, zenith and azimuth are set with >>> src = Source("K",zenith=20*math.pi/180/3600,azimuth=math.pi/4,resolution=(n,n),magnitude=10) A Laser guide star constellation of 6 sources is defined with >>> zen = np.ones(6)*30*math.pi/180/3600, >>> azi = np.linspace(0,5,6)*2*math.pi/6 >>> lgs = ceo.Source("R",zenith=zen,azimuth=azi,height=90e3,resolution=(n,n)) If the source is to be used for ray tracing through an optical model like the GMT, then >>> lgs = ceo.Source("R",zenith=zen,azimuth=azi,height=90e3, rays_box_size=25.5, rays_box_sampling=n, rays_origin=[0,0,25]) """ def __cinit__(self, *args, **kwargs): self._c_source = new source() self.parameters = {'args':args,'kwargs':kwargs} def __init__(self,photometric_band, zenith=0, azimuth=0, magnitude=None, height=float("inf"), resolution = (0,0), fwhm = None, rays_box_size=None, rays_box_sampling=None, rays_origin=[0.0,0.0,25.0], chief_ray_origin=None, samplingTime = 1.0, int N_RAY=0, double[:] rays_x=None, double[:] rays_y=None): cdef float[:] z,a if rays_box_sampling is not None: resolution = [rays_box_sampling, rays_box_sampling] self.n = resolution[0] self.m = resolution[1] self._zenith_ = np.array( zenith, dtype=np.float64, ndmin=1) self._azimuth_ = np.array( azimuth, dtype=np.float64, ndmin=1) if magnitude is None: self._magnitude_ = np.zeros_like( self._zenith_ , dtype=np.float32) else: self._magnitude_ = np.array( magnitude, dtype=np.float32, ndmin=1) self.height = float(height); self.size = self._zenith_.size self.band = photometric_band.encode() cdef vector _origin_, _chief_origin_ if N_RAY>0: _origin_.x = rays_origin[0] _origin_.y = rays_origin[1] _origin_.z = rays_origin[2] self._c_source.setup(self.band, &self._magnitude_[0], &self._zenith_[0], &self._azimuth_[0], height, self.size, N_RAY, &rays_x[0], &rays_y[0], _origin_) else: if rays_box_size is None: z = np.array( zenith, dtype=np.float32, ndmin=1) a = np.array( azimuth, dtype=np.float32, ndmin=1) self._c_source.setup(self.band, &z[0],&a[0], height, self.size, np.prod(resolution)) else: _origin_.x = rays_origin[0] _origin_.y = rays_origin[1] _origin_.z = rays_origin[2] if chief_ray_origin is None: self._c_source.setup(self.band, &self._magnitude_[0], &self._zenith_[0], &self._azimuth_[0], height, self.size, rays_box_size, rays_box_sampling, _origin_) else: _chief_origin_.x = chief_ray_origin[0] _chief_origin_.y = chief_ray_origin[1] _chief_origin_.z = chief_ray_origin[2] self._c_source.setup_chief(self.band, &self._magnitude_[0], &self._zenith_[0], &self._azimuth_[0], height, self.size, rays_box_size, rays_box_sampling, _origin_, _chief_origin_) #if magnitude is not None: # self._c_source.magnitude = magnitude if fwhm is not None: self._c_source.fwhm = fwhm self.wavefront = Complex_amplitude(self) self.rays = SourceBundle(self) self.N_SRC = self._c_source.N_SRC self.sphere_distance = None self.samplingTime = samplingTime self.timeStamp = 0.0 def __dealloc__(self): self._c_source.cleanup() def __rshift__(Source x, tuple y): """ Set the optical path Parameters ---------- x : Source The Source object y : tuple A list of CEO object the source is propagating through """ x.OPTICALPATH = y #print "OPTICAL PATH: "+str(map(lambda z: z.__class__.__name__,x.OPTICALPATH)) def __richcmp__(Source x, tuple y, int op): """ Append to the optical path Parameters ---------- x : Source The Source object y : tuple A list of CEO object the source is propagating through """ if op==Py_GT: x.OPTICALPATH = x.OPTICALPATH + y #print "OPTICAL PATH: "+str(map(lambda z: z.__class__.__name__,x.OPTICALPATH)) else: assert False def __invert__(self): """ Reset the source """ self.reset() def __pos__(self): """ Reset the source and propagate the source through the optical path """ assert self.OPTICALPATH is not None, "OPTICALPATH is not set!" self.reset() #map(lambda x: x.propagate(self), self.OPTICALPATH) l = [x.propagate(self) for x in self.OPTICALPATH] self.timeStamp += self.samplingTime def __iadd__(self, int x): cdef int k for k in range(x): self.__pos__() return self def opd2phase(self, int FLAG=1): """ Transfer the OPD from ray tracing to the wavefront phase Parameters ---------- FLAG : int (optional) If FLAG is 0, the rays vignetting mask is not reset, this is used for sequential ray tracing (default: 1) """ self._c_source.opd2phase(FLAG) def updateDirections(self, zenith, azimuth): """ Update the zenith and azimuth direction angles Parameters ---------- zenith : float, tuple, list or numpy array The zenith angle array in radian azimuth : float, tuple, list or numpy array The azimuth angle array in radian """ self._zenith_ = np.array(zenith ,dtype=np.float64,ndmin=1) self._azimuth_ = np.array(azimuth,dtype=np.float64,ndmin=1) self._c_source.update_directions(&(self._zenith_[0]),&(self._azimuth_[0]),self.N_SRC) def nBackgroundPhoton(self, float backgroundMagnitude): return self._c_source.n_background_photon(backgroundMagnitude) property phase: def __get__(self): x = cuFloatArray(shape=(self.n*self.size,self.m)) x._c_gpu.dev_data = self._c_source.wavefront.phase return x def __set__(self,cuFloatArray val): self._c_source.wavefront.add_phase(1,val._c_gpu.dev_data) property amplitude: def __get__(self): x = cuFloatArray(shape=(self.n*self.size,self.m)) x._c_gpu.dev_data = self._c_source.wavefront.amplitude return x property wavelength: def __get__(self): return self._c_source.wavelength() property spectral_bandwidth: def __get__(self): return self._c_source.spectral_bandwidth() property nPhoton: def __get__(self): cdef int k NPH = np.zeros(self._magnitude_.size) for k in range(self._magnitude_.size): NPH[k] = self._c_source.n_photon(self._magnitude_[k]) return NPH property zenith: def __get__(self): return np.asarray(self._zenith_) property azimuth: def __get__(self): return np.asarray(self._azimuth_) property magnitude: def __get__(self): return np.asarray(self._magnitude_) def __set__(self, list magnitude): self._magnitude_ = np.array(magnitude,dtype=np.float32,ndmin=1) self._c_source.update_magnitude(&(self._magnitude_[0]), self._magnitude_.size) property theta_x: def __get__(self): return np.tan(self._zenith_)*np.cos(self._azimuth_) property theta_y: def __get__(self): return np.tan(self._zenith_)*np.sin(self._azimuth_) property fwhm: def __get__(self): return self._c_source.fwhm def __set__(self,float value): self._c_source.fwhm = value def masked(self, MaskAbstract tel): """ Apply the binary mask of the telescope pupil to the source wavefront Parameters ---------- mask : Mask The binary mask structure """ self._c_source.wavefront.masked(tel._c_mask) def masked(self): """ Apply the binary mask of the telescope pupil to the source wavefront """ self._c_source.wavefront.masked() def reset(self): """ Reset the wavefront amplitude to 1 and phase to 0 and re--initialize the ray bundle """ self._c_source.wavefront.reset() self._c_source.reset_rays() def resetRays(self, int FLAG=1): """ Re--initialize the ray bundle Parameters ---------- FLAG : int (optional) If FLAG is 0, the rays vignetting mask is not reset, this is used for sequential ray tracing (default: 1) """ self._c_source.reset_rays(FLAG) def piston(self, where='pupil',int units_exponent=0): """ Get the piston corresponding to either the pupil or the segments Parameters ---------- where : string, optional Either "pupil" for a piston on the full pupil or "segments" for pistons of each segment; default: pupil units_exponent : int, optional Multiply the piston values by 10^-units_exponent; default: 0 Returns ------- ndarray of float The vector of piston values in meters as either an N_SRC vector for "pupil" or a N_SRCx7 array for "segments" Examples -------- >>> import ceo >>> gs = ceo.Source("V",rays_box_size=25.5, rays_box_sampling=256, rays_origin=[0,0,25]) >>> gmt = ceo.GMT_MX(25.5,256) >>> gmt.propagate(gs) The piston over the entire pupil is obtained with: >>> gs.piston(where="pupil") The 7 segment pistons are retrieved with >>> gs.piston(where="segments") See also -------- GMT_MX : a class embedding GMT M1 and M2 classes GMT_M1 : a class for GMT M1 model """ assert where=="pupil" or where=="segments", "where parameter is either ""pupil"" or ""segments""" cdef: int k tuple shape float[:,::1] amplitude, phase shape = (self.size,self.n*self.m) amplitude = self.wavefront.amplitude.host(shape=shape) phase = self.wavefront.phase.host(shape=shape) if where=='pupil': return np.sum(phase,axis=1)/np.sum(amplitude,axis=1)*10**-units_exponent if where=='segments': # Q = mask*amplitude.T # return np.dot(Q,np.reshape( self.wavefront.phase.host() , (-1,) ) ).flatten()/np.sum(Q,axis=1) ps = np.zeros((self.size,7),dtype=np.float32) for k in range(self.size): Q = self.rays.piston_mask[k]*amplitude[k,:] ps[k,:] = np.dot(Q,phase[k,:])/np.sum(Q,axis=1) return np.array(ps)*10**-units_exponent def phaseRms(self,where='pupil',int units_exponent=0): """ Get the phase rms corresponding to either the pupil or the segments Parameters ---------- where : string, optional Either "pupil" for a piston on the full pupil or "segments" for pistons of each segment; default: pupil units_exponent : int, optional Multiply the rms phase values by 10^-units_exponent; default: 0 Returns ------- ndarray of float The vector of phase rms values in meters as either an N_SRC vector for "pupil" or a N_SRCx7 array for "segments" """ assert where=="pupil" or where=="segments", "where parameter is either ""pupil"" or ""segments""" cdef: int k, s tuple shape shape = (self.size,self.n*self.m) phase = self.wavefront.phase.host(shape=shape) if where=='pupil': return self.wavefront.rms(units_exponent) if where=='segments': rms = np.zeros((self.size,7),dtype=np.float32) for k in range(self.size): Q = self.rays.piston_mask[k] for s in range(7): rms[k,s] = np.std(phase[k,Q[s,:]]) return np.array(rms)*10**-units_exponent def segmentsWavefrontGradientSlow(self): """ Computes the average angle of arrival (tip and tilt) of each segment Returns cuFloatArray The average angle of arrival of each segment in an array of size [N_SRCx7x2,1] as [cx,cy] """ cdef cuFloatArray sxy sxy = cuFloatArray( host_data=np.zeros((7*self.N_SRC*2,1) ) ) self._c_source.wavefront.segments_gradient_average(sxy._c_gpu.dev_data, sxy._c_gpu.dev_data+7, self._c_source.rays.L, self._c_source.rays.d__piston_mask) return sxy def segmentsWavefrontGradient(self): """ Computes the average angle of arrival (tip and tilt) of each segment Returns ------- cuFloatArray The average angle of arrival of each segment in an array of size [N_SRCx7x2,1] as [cx,cy] """ cdef float[:,:] sxy sxy = np.zeros((2*self.N_SRC*7,1),dtype=np.single) self._c_source.wavefront.segments_gradient_averageFast(&sxy[0,0], &sxy[7*self.N_SRC,0], self._c_source.rays.L, self._c_source.rays.d__piston_mask) return np.asarray(sxy) def opticalTransferFunction(self): cdef cuFloatComplexArray otf cdef int N N = self.rays.N_L*2-1 otf = cuFloatComplexArray(shape=(N,N),dev_malloc=True) self._c_source.optical_transfer_function(otf._c_gpu.dev_data) return otf @ \subsubsection{JSource} \label{sec:jsource} \index{source!python!JSource} <>= from utilities import JSONAbstract class JSource(JSONAbstract,Source): """ docstring for JSource """ def __init__(self, *args, **kwargs): print "@(ceo.JSource)>" JSONAbstract.__init__(self, *args, **kwargs) Source.__init__(self,self.jprms["band"], zenith = self.jprms["zenith"], azimuth = self.jprms["azimuth"], height = np.float(self.jprms["height"]), fwhm = self.jprms["fwhm"], magnitude = self.jprms["magnitude"], rays_box_size = self.jprms["pupil size"], rays_box_sampling = self.jprms["pupil sampling"], rays_origin = self.jprms["rays origin"]) @ \section{PSSn} \label{sec:source-pssn} \index{source!python!PSSn} <>= cdef cppclass pssn: int N int N_OTF int N_O float2 *d__W float2 *d__O float2 *d__O0 float2 *d__C float2 *buffer void setup(source *, float, float) void cleanup() void otf(source *) float eval() float eval(float *) float oeval() float oeval(float *) @ \subsubsection{Class definition} \label{sec:class-definition-pssn} <>= # PSSn cdef class PSSn: cdef: pssn *_c_pssn public cuFloatComplexArray W public cuFloatComplexArray O0 public cuFloatComplexArray O public cuFloatComplexArray C public cuFloatComplexArray buffer @ \subsection{PYX file} \label{sec:pyx-file-pssn} <>= #PSSn cdef class PSSn: def __cinit__(self, Source src, float r0, float L0): self._c_pssn = new pssn() self._c_pssn.setup(src._c_source,r0,L0) self.W = cuFloatComplexArray( shape=(self._c_pssn.N_OTF*self._c_pssn.N,self._c_pssn.N_OTF), dev_malloc=False) self.W._c_gpu.dev_data = self._c_pssn.d__W self.O0 = cuFloatComplexArray( shape=(self._c_pssn.N_OTF*self._c_pssn.N,self._c_pssn.N_OTF), dev_malloc=False) self.O0._c_gpu.dev_data = self._c_pssn.d__O0 self.O = cuFloatComplexArray( shape=(self._c_pssn.N_OTF*self._c_pssn.N,self._c_pssn.N_OTF), dev_malloc=False) self.O._c_gpu.dev_data = self._c_pssn.d__O self.C = cuFloatComplexArray( shape=(self._c_pssn.N_OTF*self._c_pssn.N,self._c_pssn.N_OTF), dev_malloc=False) self.C._c_gpu.dev_data = self._c_pssn.d__C self.buffer = cuFloatComplexArray( shape=(self._c_pssn.N_OTF*self._c_pssn.N,self._c_pssn.N_OTF), dev_malloc=False) self.buffer._c_gpu.dev_data = self._c_pssn.buffer def __dealloc__(self): self._c_pssn.cleanup() def otf(self,Source src): self._c_pssn.otf(src._c_source) def eval(self): return self._c_pssn.eval() def evals(self): cdef float[:] results results = np.zeros(self._c_pssn.N,dtype=np.single) self._c_pssn.eval(&results[0]) return np.asarray(results) def oeval(self): return self._c_pssn.oeval() def oevals(self): cdef float[:] results results = np.zeros(self._c_pssn.N,dtype=np.single) self._c_pssn.oeval(&results[0]) return np.asarray(results) @property def N_O(self): return self._c_pssn.N_O @N_O.setter def N_O(self,int value): self._c_pssn.N_O = value @ \subsection{PXD file} \label{sec:pxd-file-10} \section{Ray tracing routines} \label{sec:ray-tracing-routines-1} <>= cdef void transform_to_S(source *, conic *) cdef void transform_to_S(source *, aperture *) cdef void transform_to_R(source *, conic *) cdef void transform_to_R(source *, aperture *) cdef void intersect(source *, conic *) cdef void reflect(source *) cdef void thin_lens(source *) @ <>= # ray tracing def Transform_to_S(Source src, Conic F): transform_to_S(src._c_source, F._c_conic) def Transform_to_S_from_A(Source src, Aperture A): transform_to_S(src._c_source, A._c_aperture) def Transform_to_R(Source src, Conic F): transform_to_R(src._c_source, F._c_conic) def Transform_to_R_from_A(Source src, Aperture A): transform_to_R(src._c_source, A._c_aperture) def Intersect(Source src, Conic F): intersect(src._c_source, F._c_conic) def Reflect(Source src): reflect(src._c_source)