% -*- mode: Noweb; noweb-code-mode: python-mode -*- \section{PXD file} \label{sec:pxd-file} <>= from source cimport source from centroiding cimport centroiding cdef extern from "atmosphere.h": cdef cppclass layer: int N_WIDTH_LAYER, N_LENGTH_LAYER float WIDTH_LAYER, LENGTH_LAYER cdef cppclass atmosphere: int N_LAYER float r0 unsigned long N_PHASE_LAYER float *d__phase_screen_LAYER layer *layers void setup(float , float, int, float *, float *, float *, float *) void setup(float , float, int, float *, float *, float *, float *, float, int, float, float) void setup(float , float, int, float *, float *, float *, float *, float, int, float, float, char *, int); void gmt_setup(float , float) void gmt_setup(float , float, float, int, float, float, char *, int) void gmt_setup(float, float, float, int, float, float) void gmt_setup(float , float, int) void gmt_setup(float , float, float, int, float, float, char *, int, int) void gmt_setup(float, float, float, int, float, float, int) void gmt_set_id(int) void gmt_set_eph(float) void cleanup() void info() void save_layer_phasescreens(char *,int) void reset() void get_phase_screen(source *, float , int , float , int , float ) void get_phase_screen(source *, float , int , float , int , float , float) void get_phase_screen(source *, float ) void get_phase_screen(float*, float *, float *, int, source *, float ) void get_phase_screen_gradient(float *, float *, float *, float *, int, source *, float) void get_phase_screen_gradient(centroiding *, int , float , source *, float ) void get_phase_screen_gradient_rolling_shutter( centroiding *, int , float , source *, float, float ) void get_phase_screen_circ_centroids(centroiding *, float, source *, int , float ) void get_phase_screen_circ_uplink_centroids(centroiding *, float, source *, int , float, char ) void rayTracing(float*, float*, float*, int, source*, float) void rayTracing(source*, float, int, float, int, float) <> @ \subsection{Class definitions} \label{sec:class-definitions} <>= from numpy cimport ndarray from utilities cimport cuFloatArray from source cimport Source from centroiding cimport Centroiding from shackHartmann cimport ShackHartmann cdef class Layer: cdef: layer *_c_layer readonly cuFloatArray phase_screen cdef void init(Layer, atmosphere *, int) cdef class AtmosphereAbstract: cdef: atmosphere *_c_atmosphere readonly list layers readonly float L0 readonly int N_LAYER readonly ndarray altitude, xi0, wind_speed, wind_direction @ \section{PYX file} \label{sec:pyx-file} \subsection{Layer} \label{sec:layer} \index{atmosphere!python!Layer} <>= cimport cython cimport numpy as np import numpy as np @cython.boundscheck(False) @cython.wraparound(False) cdef class Layer: def __cinit__(self): self._c_layer = new layer() cdef void init(self, atmosphere *atm, int i_LAYER): self._c_layer = &(atm.layers[i_LAYER]) property N_WIDTH: def __get__(self): return self._c_layer.N_WIDTH_LAYER property N_LENGTH: def __get__(self): return self._c_layer.N_LENGTH_LAYER property WIDTH: def __get__(self): return self._c_layer.WIDTH_LAYER property LENGTH: def __get__(self): return self._c_layer.LENGTH_LAYER @ \subsection{AtmosphereAbstract} \label{sec:atmosphereabstract} \index{atmosphere!python!AtmosphereAbstract} <>= cdef class AtmosphereAbstract: """ Atmosphere base class Attributes ---------- r0 : number The Fried parameter [m]. L0 : number The outer scale of turbulence [m]. N_LAYER : int, optional The number of atmospheric layers, default: 1. altitude : ndarray of floats, optional The atmosphere layer altitudes [m], default: 0.0. xi0 : ndarray of floats, optional The atmosphere layer weights 0<=xi0<=1, the sum must equal 1, default: 1.0. wind_speed : ndarray of floats, optional The atmosphere layer wind speeds [m/s], default: 0.0. wind_direction : ndarray of floats, optional The atmosphere layer wind directions [rd], default: 0.0. """ def __cinit__(self, *args, **kwargs): self._c_atmosphere = new atmosphere() def __init__(self, float r0, float L0, N_LAYER=1, altitude=0.0, xi0=1.0, wind_speed=0.0, wind_direction=0.0): self.L0 = L0 self.N_LAYER = N_LAYER self.altitude = np.array( altitude, dtype=np.float32) self.xi0 = np.array( xi0, dtype=np.float32) self.wind_speed = np.array( wind_speed, dtype=np.float32) self.wind_direction = np.array( wind_direction, dtype=np.float32) <> @ \subsection{Atmosphere} \label{sec:atmosphere} \index{atmosphere!python!Atmosphere} <>= cdef class Atmosphere(AtmosphereAbstract): """ Create an atmosphere Parameters ---------- r0 : number The Fried parameter [m]. L0 : number The outer scale of turbulence [m]. N_LAYER : int, optional The number of atmospheric layers, default: 1. altitude : ndarray of floats, optional The atmosphere layer altitudes [m], default: 0.0. xi0 : ndarray of floats, optional The atmosphere layer weights 0<=xi0<=1, the sum must equal 1, default: 1.0. wind_speed : ndarray of floats, optional The atmosphere layer wind speeds [m/s], default: 0.0. wind_direction : ndarray of floats, optional The atmosphere layer wind directions [rd], default: 0.0. L : float, optional The size of the square phase screen at an altitude of 0km [m], default: 0.0. NXY_PUPIL : int, optional The sampling of the square phase screen at an altitude of 0km [px], default: 0. fov : float, optional The field of view [rd], default: 0.0. duration : float, optional The time length of the phase screens [s], default: 1.0. filename : str, optional The name of the file to save the phase screens of the layers to, default: None. N_DURATION : int, optional The number of phase screens of time length duration, default: 1. Examples -------- >>> import numpy as np >>> import ceo >>> h = np.array([0,1,2.5,5,8,10,12],dtype=np.float32)*1e3 >>> xi0 = np.array([0.3,0.15,0.25,0.1,0.1,0.05,0.05], dtype=np.float32) >>> vs = np.array([10,10,10,10,10,10,10],dtype=np.float32) >>> vo = np.array([0,0,0,0,0,0,0], dtype=np.float32) An atmosphere with r0=15cm and L0=30m is created with >>> atm = Atmosphere(15e-2,30,h,xi0,vs,vo) A 25m square phase screen sampled with 401 pixels across is computed a t=0s with >>> L = 25.0 >>> N = 401 >>> delta = L/(N-1) >>> src =ceo.Source('R',resolution=(N,N)) >>> atm.get_phase_screen(src,delta,N,delta,N,0.0) The phase screen can be computed faster by precomputing the phase screen in the layers of the atmosphere: >>> atm = Atmosphere(15e-2,30,h,xi0,vs,vo,L=L,NXY_PUPIL=N,fov=1.0*ceo.constants.ARCMIN2RAD,duration=1.0) The phase screens are computed for a 1 arc minute field-of-view and for 1s, meaning that a source can be propagated anywhere within this field-of-view and for anytime between 0 and 1s; the fast ray tracing through the precomputed layers is called with >>> atm.ray_tracing(src,delta,N,delta,N,0.0) The precomputed phase screens can be saved to a file by passing the filename to the atmosphere constructor >>> atm = Atmosphere(15e-2,30,h,xi0,vs,vo,L=L,NXY_PUPIL=N,fov=1.0*ceo.constants.ARCMIN2RAD,duration=1.0, filename="myWondrousAtmosphere.bin") If the file already exists the constructor will attempt to load it. The phase screens are loaded into the device memory, so duration must be chosen with care to not overflow the device memory. For longer time sequence of N_DURATIONxduration time length, the full sequence of phase screens is loaded into the device memory. ray_tracing will copy the appropriate time sequence of duration time length from the host to the device based of the time parameters. >>> atm = Atmosphere(15e-2,30,h,xi0,vs,vo,L=L,NXY_PUPIL=N,fov=1.0*ceo.constants.ARCMIN2RAD,duration=10.0, filename="myWondrousAtmosphere.bin",N_DURATION=10) >>> atm.ray_tracing(src,delta,N,delta,N,8.88) See also -------- AtmosphereAbstract : the Atmosphere base class """ def __init__(self, float r0, float L0, N_LAYER=1, altitude=0.0, xi0=1.0, wind_speed=0.0, wind_direction=0.0, float L=0.0, int NXY_PUPIL=0, float fov=0.0, float duration=1.0, str filename=None, int N_DURATION=1): super(Atmosphere,self).__init__(r0,L0, N_LAYER=N_LAYER, altitude=altitude, xi0=xi0, wind_speed=wind_speed, wind_direction=wind_direction) if NXY_PUPIL==0: self._c_atmosphere.setup(r0, self.L0, self.N_LAYER, self.altitude.data, self.xi0.data, self.wind_speed.data, self.wind_direction.data) else: if filename is not None: self._c_atmosphere.setup(r0, self.L0, self.N_LAYER, self.altitude.data, self.xi0.data, self.wind_speed.data, self.wind_direction.data, L, NXY_PUPIL, fov, duration, filename.encode(), N_DURATION) else: self._c_atmosphere.setup(r0, self.L0, self.N_LAYER, self.altitude.data, self.xi0.data, self.wind_speed.data, self.wind_direction.data, L, NXY_PUPIL, fov,duration) self.init() cdef init(self): cdef: int k, N Layer L N = 0 self.layers = [] for k in range(self.N_LAYER): L = Layer() L._c_layer = self._c_atmosphere.layers + k L.phase_screen = cuFloatArray(shape=(L.N_LENGTH,L.N_WIDTH), dev_malloc=False) L.phase_screen._c_gpu.dev_data = self._c_atmosphere.d__phase_screen_LAYER + N self.layers.append( L ) N += L.N_WIDTH*L.N_LENGTH @ \subsection{GmtAtmosphere} \label{sec:gmtatmosphere} \index{atmosphere!python!GmtAtmosphere} <>= cdef class GmtAtmosphere(AtmosphereAbstract): """ Creates an atmosphere with the GMT default Cn2 and wind vector profiles Parameters ---------- r0 : number The Fried parameter [m]. L0 : number The outer scale of turbulence [m]. L : float, optional The size of the square phase screen at an altitude of 0km [m], default: 0.0. NXY_PUPIL : int, optional The sampling of the square phase screen at an altitude of 0km [px], default: 0. fov : float, optional The field of view [rd], default: 0.0. duration : float, optional The time length of the phase screens [s], default: 1.0. filename : str, optional The name of the file to save the phase screens of the layers to, default: None. N_DURATION : int, optional The number of phase screens of time length duration, default: 1. SEED : int, optional The seed of the random generator, default: 2016. ID : int, optional The ID number of the GMT atmosphere model, default: 1. See also -------- AtmosphereAbstract : the Atmosphere base class Atmosphere : the generic atmosphere class """ def __init__(self, float r0, float L0, float L=0.0,int NXY_PUPIL=0,float fov=0.0, float duration=1.0, str filename=None, int N_DURATION=1, int SEED=2016, int ID=1, float EPH=0.0): super(GmtAtmosphere,self).__init__(r0,L0) self._c_atmosphere.gmt_set_id(ID) self._c_atmosphere.gmt_set_eph(EPH) if NXY_PUPIL==0: self._c_atmosphere.gmt_setup(r0, L0, SEED) else: if filename is None: self._c_atmosphere.gmt_setup(r0, L0, L, NXY_PUPIL, fov, duration, SEED) else: self._c_atmosphere.gmt_setup(r0, L0, L, NXY_PUPIL, fov, duration, filename.encode(), N_DURATION, SEED) @ \subsection{JGmtAtmosphere} \label{sec:jgmtatmosphere} \index{atmosphere!python!JGmtAtmosphere} <>= from utilities import JSONAbstract from constants import ARCMIN2RAD class JGmtAtmosphere(JSONAbstract,GmtAtmosphere): """ Creates an atmosphere with the GMT default Cn2 and wind vector profiles Parameters ---------- jsonfile : the JSON file that describes the atmosphere model zipdata : the zip file containing the atmosphere mode cache : the directory where to unzip the atmosphere data Examples -------- >>> import ceo >>> atm = ceo.JGmtAtmosphere(jsonfile='/mnt/bins/gmtAtmosphereL030.json', zipdata='s3://gmto.rconan/gmtAtmosphereL030.zip', cache='/mnt/bins/') >>> import json >>> json.loads(open('/mnt/bins/gmtAtmosphereL030.json').read()) {u'L0': 30, u'duration': 15, u'field-of-view [arcmin]': 20, u'filename': u'/mnt/bins/gmtAtmosphereL030.bin', u'pupil sampling': 346, u'pupil size': 26.0, u'r0': 0.15, u'time sampling': 20} See also -------- JSONAbstract : the JSON file interface class GmtAtmosphere : the atmosphere class with the default GMT atmosphere model """ def __init__(self, *args, **kwargs): print "@(ceo.JGmtAtmosphere)>" JSONAbstract.__init__(self, *args, **kwargs) #self.jprms['filename'] = self.jprms['filename'].encode('utf-8') GmtAtmosphere.__init__(self, **self.jprms) @ <>= def reset(self): self._c_atmosphere.reset() def save_layer_phasescreens(self, bytes filename, int N_DURATION): self._c_atmosphere.save_layer_phasescreens(filename, N_DURATION) def get_phase_screen(self, Source src, float delta_x, int N_x, float delta_y, int N_y, float time, float exponent=0.0): """ Computes the N_x$\times$N_y phase screen of the source wavefront Parameters ---------- src: Source object The light source delta_x: float X-axis spacing N_x: int X-axis sampling delta_y: float Y-axis spacing N_y: int Y-axis sampling time: float The time of the atmosphere snapshot since the beginning of the simulation """ if exponent>0: self._c_atmosphere.get_phase_screen(src._c_source, delta_x, N_x, delta_y, N_y, time, exponent) else: self._c_atmosphere.get_phase_screen(src._c_source, delta_x, N_x, delta_y, N_y, time) def phase_screen(self, Source src, float time): self._c_atmosphere.get_phase_screen(src._c_source, time) def get_phase_values(self, cuFloatArray x, cuFloatArray y, Source src, float time): """ Computes the phase screen of the source wavefront Parameters ---------- x: cuFloatArray X coordinates of the phase screen samples. y: cuFloatArray Y coordinates of the phase screen samples. src: Source object The light source. time: float The time of the atmosphere snapshot since the beginning of the simulation. """ ps = cuFloatArray(shape=x.shape,dev_malloc=True) self._c_atmosphere.get_phase_screen(ps._c_gpu.dev_data, x._c_gpu.dev_data, y._c_gpu.dev_data, x._c_gpu.N, src._c_source, time) return ps def ray_tracing_values(self, cuFloatArray x_PUPIL,cuFloatArray y_PUPIL, int NXY_PUPIL, Source src, float tau): """ Computes the phase screen of the source wavefront Parameters ---------- x: cuFloatArray X coordinates of the phase screen samples. y: cuFloatArray Y coordinates of the phase screen samples. src: Source object The light source. time: float The time of the atmosphere snapshot since the beginning of the simulation. """ ps = cuFloatArray(host_data=np.zeros(x_PUPIL.shape, dtype=np.float32 , order='C')) self._c_atmosphere.rayTracing(x_PUPIL._c_gpu.dev_data, y_PUPIL._c_gpu.dev_data, ps._c_gpu.dev_data, NXY_PUPIL, src._c_source, tau) return ps def ray_tracing(self, Source src, float delta_x, int N_X, float delta_y, int N_Y, float tau): """ Computes the N_x$\times$N_y phase screen of the source wavefront Parameters ---------- src: Source object The light source delta_x: float X-axis spacing N_x: int X-axis sampling delta_y: float Y-axis spacing N_y: int Y-axis sampling time: float The time of the atmosphere snapshot since the beginning of the simulation """ self._c_atmosphere.rayTracing(src._c_source, delta_x, N_X, delta_y, N_Y, tau) def propagate(self, Source src): cdef: float delta int N N = src.rays.N_L delta = src.rays.L/(N-1) self._c_atmosphere.rayTracing(src._c_source, delta, N, delta, N, src.timeStamp) def get_phase_screen_gradient_values(self, double[:] x, double[:] y, Source src, float time): cdef int Nxy Nxy = x.size cdef cuFloatArray sx, sy sx = cuFloatArray(shape=(Nxy,1),dev_malloc=True) sy = cuFloatArray(shape=(Nxy,1),dev_malloc=True) cdef cuFloatArray gx,gy gx = cuFloatArray(host_data=x) gy = cuFloatArray(host_data=y) self._c_atmosphere.get_phase_screen_gradient(sx._c_gpu.dev_data, sy._c_gpu.dev_data, gx._c_gpu.dev_data, gy._c_gpu.dev_data, Nxy, src._c_source, time) return np.asarray(sx.host(),dtype=np.double),np.asarray(sy.host(),dtype=np.double) def get_phase_screen_gradient(self, Centroiding cog, int NL, float d, Source src, float time, float delay=0.0): """ Computes the phase screen gradient of the source wavefront Parameters ---------- cog: Centroiding object The centroids container object. NL: int The size of the square lenslet array. d: float The lenslet array pitch. src: Source object The light source. time: float The time of the atmosphere snapshot since the beginning of the simulation """ if delay>0: self._c_atmosphere.get_phase_screen_gradient_rolling_shutter( cog._c_centroiding, NL, d, src._c_source, time, delay) else: self._c_atmosphere.get_phase_screen_gradient(cog._c_centroiding, NL, d, src._c_source, time) def get_phase_screen_gradientSH(self, ShackHartmann wfs, int NL, float d, Source src, float time): """ Computes the phase screen gradient of the source wavefront Parameters ---------- cog: Centroiding object The centroids container object. NL: int The size of the square lenslet array. d: float The lenslet array pitch. src: Source object The light source. time: float The time of the atmosphere snapshot since the beginning of the simulation """ self._c_atmosphere.get_phase_screen_gradient(&wfs._c_shackHartmann.data_proc, NL, d, src._c_source, time) def get_phase_screen_circ_centroids(self, Centroiding cog, float R, Source src, int N_SRC, float time): self._c_atmosphere.get_phase_screen_circ_centroids(cog._c_centroiding, R, src._c_source, N_SRC, time) def get_phase_screen_circ_uplink_centroids(self, Centroiding cog, float R, Source src, int N_SRC, float time, char focused): self._c_atmosphere.get_phase_screen_circ_uplink_centroids(cog._c_centroiding, R, src._c_source, N_SRC, time,focused) def __dealloc__(self): self._c_atmosphere.cleanup() property r0: def __get__(self): return self._c_atmosphere.r0 def __set__(self,float value): self._c_atmosphere.r0 = value