COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
shesha.util.iterkolmo Namespace Reference

Stencil and matrices computation for the creation of a turbulent screen. More...

Functions

def create_stencil (n)
 TODO docstring. More...
 
def stencil_size (n)
 TODO docstring. More...
 
def stencil_size_array (size)
 Compute_size2(np.ndarray[ndim=1, dtype=np.int64_t] size) More...
 
def Cxz (n, Zx, Zy, Xx, Xy, istencil, L0)
 Cxz computes the covariance matrix between the new phase vector x (new column for the phase screen), and the already known phase values z. More...
 
def Cxx (n, Zxn, Zyn, Xx, Xy, L0)
 Cxx computes the covariance matrix of the new phase vector x (new column for the phase screen). More...
 
def Czz (n, Zx, Zy, ist, L0)
 Czz computes the covariance matrix of the already known phase values z. More...
 
def AB (n, L0, deltax, deltay, rank=0)
 DOCUMENT AB, n, A, B, istencil This function initializes some matrices A, B and a list of stencil indexes istencil for iterative extrusion of a phase screen. More...
 
def extrude (p, r0, A, B, istencil)
 DOCUMENT p1 = extrude(p,r0,A,B,istencil) More...
 
def phase_struct (r, L0=None)
 TODO docstring. More...
 
def rodconan (r, L0)
 The phase structure function is computed from the expression Dphi(r) = k1 * L0^(5. More...
 
def asymp_macdo (x)
 Computes a term involved in the computation of the phase struct function with a finite outer scale according to the Von-Karman model. More...
 
def macdo_x56 (x, k=10)
 Computation of the function f(x) = x^(5/6)*K_{5/6}(x) using a series for the esimation of K_{5/6}, taken from Rod Conan thesis : K_a(x)=1/2 \sum_{n=0}^\infty \frac{(-1)^n}{n!} \left(\Gamma(-n-a) (x/2)^{2n+a} + \Gamma(-n+a) (x/2)^{2n-a} \right) , with a = 5/6. More...
 
def create_screen_assist (screen_size, L0, r0)
 
def create_screen (r0, pupixsize, screen_size, L0, A, B, ist)
 DOCUMENT create_screen screen = create_screen(r0,pupixsize,screen_size,&A,&B,&ist) More...
 

Detailed Description

Stencil and matrices computation for the creation of a turbulent screen.

Author
COMPASS Team https://github.com/ANR-COMPASS
Version
5.0.0
Date
2020/05/18

This file is part of COMPASS https://anr-compass.github.io/compass/

Copyright (C) 2011-2019 COMPASS Team https://github.com/ANR-COMPASS All rights reserved. Distributed under GNU - LGPL

COMPASS is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or any later version.

COMPASS: End-to-end AO simulation tool using GPU acceleration The COMPASS platform was designed to meet the need of high-performance for the simulation of AO systems.

The final product includes a software package for simulating all the critical subcomponents of AO, particularly in the context of the ELT and a real-time core based on several control approaches, with performances consistent with its integration into an instrument. Taking advantage of the specific hardware architecture of the GPU, the COMPASS tool allows to achieve adequate execution speeds to conduct large simulation campaigns called to the ELT.

The COMPASS platform can be used to carry a wide variety of simulations to both testspecific components of AO of the E-ELT (such as wavefront analysis device with a pyramid or elongated Laser star), and various systems configurations such as multi-conjugate AO.

COMPASS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with COMPASS. If not, see https://www.gnu.org/licenses/lgpl-3.0.txt.

Function Documentation

◆ AB()

def shesha.util.iterkolmo.AB (   n,
  L0,
  deltax,
  deltay,
  rank = 0 
)

DOCUMENT AB, n, A, B, istencil This function initializes some matrices A, B and a list of stencil indexes istencil for iterative extrusion of a phase screen.

The method used is described by Fried & Clark in JOSA A, vol 25, no 2, p463, Feb 2008. The iteration is : x = A(z-zRef) + B.noise + zRef with z a vector containing "old" phase values from the initial screen, that are listed thanks to the indexes in istencil.

SEE extrude createStencil Cxx Cxz Czz

Definition at line 202 of file iterkolmo.py.

202  """
203 
204  # if (rank == 0):
205  # print("create stencil and Z,X matrices")
206  Zx, Zy, Xx, Xy, istencil = create_stencil(n)
207  # if (rank == 0):
208  # print("create zz")
209  zz = Czz(n, Zx, Zy, istencil, L0)
210  # if (rank == 0):
211  # print("create xz")
212  xz = Cxz(n, Zx, Zy, Xx, Xy, istencil, L0)
213  # if (rank == 0):
214  # print("create xx")
215  xx = Cxx(n, Zx[0, n - 1], Zy[0, n - 1], Xx, Xy, L0)
216 
217  U, s, V = np.linalg.svd(zz)
218  s1 = s
219  s1[s.size - 1] = 1
220  s1 = 1. / s1
221  s1[s.size - 1] = 0
222  zz1 = np.dot(np.dot(U, np.diag(s1)), V)
223  # if (rank == 0):
224  # print("compute zz pseudo_inverse")
225  # zz1=np.linalg.pinv(zz)
226 
227  # if (rank == 0):
228  # print("compute A")
229  A = np.dot(xz, zz1)
230 
231  # if (rank == 0):
232  # print("compute bbt")
233  bbt = xx - np.dot(A, xz.T)
234  # if (rank == 0):
235  # print("svd of bbt")
236  U1, l, V1 = np.linalg.svd(bbt)
237  # if (rank == 0):
238  # print("compute B")
239  B = np.dot(U1, np.sqrt(np.diag(l)))
240 
241  test = np.zeros((n * n), np.float32)
242  test[istencil] = np.arange(A.shape[1]) + 1
243  test = np.reshape(test, (n, n), "C")
244  isty = np.argsort(test.T.flatten("C")).astype(np.uint32)[n * n - A.shape[1]:]
245 
246  if (deltay < 0):
247  isty = (n * n - 1) - isty
248  if (deltax < 0):
249  istencil = (n * n - 1) - istencil
250 
251  return np.asfortranarray(A.astype(np.float32)), np.asfortranarray(
252  B.astype(np.float32)), istencil.astype(np.uint32), isty.astype(np.uint32)
253 
254 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ asymp_macdo()

def shesha.util.iterkolmo.asymp_macdo (   x)

Computes a term involved in the computation of the phase struct function with a finite outer scale according to the Von-Karman model.

The term involves the MacDonald function (modified bessel function of second kind) K_{5/6}(x), and the algorithm uses the asymptotic form for x ~ infinity.

Warnings
  • This function makes a floating point interrupt for x=0 and should not be used in this case.
  • Works only for x>0.

Definition at line 354 of file iterkolmo.py.

354  # k2 is the value for
355  # gamma_R(5./6)*2^(-1./6)
356  k2 = 1.00563491799858928388289314170833
357  k3 = 1.25331413731550012081 # sqrt(pi/2)
358  a1 = 0.22222222222222222222 # 2/9
359  a2 = -0.08641975308641974829 # -7/89
360  a3 = 0.08001828989483310284 # 175/2187
361  x_1 = 1. / x
362  res = k2 - k3 * np.exp(-x) * x**(1 / 3.) * (1.0 + x_1 * (a1 + x_1 * (a2 + x_1 * a3)))
363  return res
364 
365 
366 def macdo_x56(x, k=10):
Here is the caller graph for this function:

◆ create_screen()

def shesha.util.iterkolmo.create_screen (   r0,
  pupixsize,
  screen_size,
  L0,
  A,
  B,
  ist 
)

DOCUMENT create_screen screen = create_screen(r0,pupixsize,screen_size,&A,&B,&ist)

creates a phase screen and fill it with turbulence r0 total r0 @ 0.5m pupixsize pupil pixel size (in meters) screen_size screen size (in pixels) A A array for future extrude B B array for future extrude ist istencil array for future extrude

Definition at line 454 of file iterkolmo.py.

454 
455  # AB, screen_size, A, B, ist,L0 # initialisation for A and B matrices for phase extrusion
456  screen = np.zeros((screen_size, screen_size),
457  dtype=np.float32) # init of first phase screen
458  for i in range(2 * screen_size):
459  screen = extrude(screen, r0 / pupixsize, A, B, ist)
460 
461  return screen
Here is the call graph for this function:

◆ create_screen_assist()

def shesha.util.iterkolmo.create_screen_assist (   screen_size,
  L0,
  r0 
)

Definition at line 424 of file iterkolmo.py.

424  A, B, istx, isty = AB(screen_size, L0)
425  phase = np.zeros((screen_size, screen_size))
426 
427  print(stencil_size(screen_size))
428 
429  # pl.ion()
430  # pl.imshow(phase, animated=True)
431  # pl.show()
432 
433  for i in range(2 * screen_size):
434  phase = extrude(phase, r0, A, B, istx)
435  # pl.clf()
436  # pl.imshow(phase, cmap='Blues')
437  # pl.draw()
438 
439  return phase
440 
441 
442 def create_screen(r0, pupixsize, screen_size, L0, A, B, ist):
Here is the call graph for this function:

◆ create_stencil()

def shesha.util.iterkolmo.create_stencil (   n)

TODO docstring.

Definition at line 43 of file iterkolmo.py.

43  """
44  Zx = np.array(np.tile(np.arange(n), (n, 1)), dtype=np.float64) + 1
45  Zy = Zx.T
46 
47  Xx = np.array(np.zeros(n) + n + 1, dtype=np.float64)
48  Xy = np.array(np.arange(n) + 1, dtype=np.float64)
49 
50  ns = int(np.log2(n + 1) + 1)
51  stencil = np.zeros((n, n))
52 
53  stencil[:, 0] = 1
54  for i in range(2, ns):
55  stencil[::2**(i - 1), 2**(i - 2)] = 1
56  stencil.itemset(2**(i - 1) - 1, 1)
57  # stencil[0,2**(i-1)-1]=1
58 
59  # i=ns
60  stencil.itemset(2**(ns - 1) - 1, 1)
61  for i in range(0, n, 2**(ns - 1)):
62  stencil.itemset(2**(ns - 2) + i * n, 1)
63  # i=ns+1
64  stencil.itemset(2**(ns) - 1, 1)
65  for i in range(0, n, 2**(ns)):
66  stencil.itemset(2**(ns - 1) + i * n, 1)
67 
68  stencil = np.roll(stencil, n // 2, axis=0)
69  stencil = np.fliplr(stencil)
70 
71  istencil = np.where(stencil.flatten() != 0)[0]
72 
73  return Zx, Zy, Xx, Xy, istencil
74 
75 
Here is the caller graph for this function:

◆ Cxx()

def shesha.util.iterkolmo.Cxx (   n,
  Zxn,
  Zyn,
  Xx,
  Xy,
  L0 
)

Cxx computes the covariance matrix of the new phase vector x (new column for the phase screen).

Definition at line 148 of file iterkolmo.py.

148  """
149  size = Xx.shape[0]
150  Xx_r = np.resize(Xx, (size, size))
151  Xy_r = np.resize(Xy, (size, size))
152 
153  tmp = np.resize(phase_struct((Zxn - Xx)**2 + (Zyn - Xy)**2, L0), (size, size))
154 
155  xx = -phase_struct((Xx_r - Xx_r.T)**2 + (Xy_r - Xy_r.T)**2, L0) # + \
156  # tmp+tmp.T
157 
158  xx += tmp + tmp.T
159 
160  xx *= 0.5
161 
162  return xx
163 
164 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Cxz()

def shesha.util.iterkolmo.Cxz (   n,
  Zx,
  Zy,
  Xx,
  Xy,
  istencil,
  L0 
)

Cxz computes the covariance matrix between the new phase vector x (new column for the phase screen), and the already known phase values z.

The known values z are the values of the phase screen that are pointed by the stencil indexes (istencil)

Definition at line 121 of file iterkolmo.py.

121  """
122 
123  size = Xx.shape[0]
124  size2 = istencil.shape[0]
125 
126  Xx_r = np.resize(Xx, (size2, size))
127  Xy_r = np.resize(Xy, (size2, size))
128  Zx_s = Zx.flatten()[istencil]
129  Zy_s = Zy.flatten()[istencil]
130  Zx_r = np.resize(Zx_s, (size, size2))
131  Zy_r = np.resize(Zy_s, (size, size2))
132 
133  tmp = phase_struct((Zx[0, size - 1] - Xx)**2 + (Zy[0, size - 1] - Xy)**2, L0)
134  tmp2 = phase_struct((Zx[0, size - 1] - Zx_s)**2 + (Zy[0, size - 1] - Zy_s)**2, L0)
135 
136  xz = -phase_struct((Xx_r.T - Zx_r)**2 + (Xy_r.T - Zy_r)**2, L0)
137 
138  xz += np.resize(tmp, (size2, size)).T + np.resize(tmp2, (size, size2))
139 
140  xz *= 0.5
141 
142  return xz
143 
144 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Czz()

def shesha.util.iterkolmo.Czz (   n,
  Zx,
  Zy,
  ist,
  L0 
)

Czz computes the covariance matrix of the already known phase values z.

The known values z are the values of the phase screen that are pointed by the stencil indexes (istencil)

Definition at line 170 of file iterkolmo.py.

170  """
171 
172  size = ist.shape[0]
173  Zx_s = Zx.flatten()[ist]
174  Zx_r = np.resize(Zx_s, (size, size))
175  Zy_s = Zy.flatten()[ist]
176  Zy_r = np.resize(Zy_s, (size, size))
177 
178  tmp = np.resize(
179  phase_struct((Zx[0, n - 1] - Zx_s)**2 + (Zy[0, n - 1] - Zy_s)**2, L0),
180  (size, size))
181 
182  zz = -phase_struct((Zx_r - Zx_r.T) ** 2 + (Zy_r - Zy_r.T) ** 2, L0) + \
183  tmp + tmp.T
184 
185  zz *= 0.5
186 
187  return zz
188 
189 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ extrude()

def shesha.util.iterkolmo.extrude (   p,
  r0,
  A,
  B,
  istencil 
)

DOCUMENT p1 = extrude(p,r0,A,B,istencil)

Extrudes a phase screen p1 from initial phase screen p. p1 prolongates p by 1 column on the right end. r0 is expressed in pixels

The method used is described by Fried & Clark in JOSA A, vol 25, no 2, p463, Feb 2008. The iteration is : x = A(z-zRef) + B.noise + zRef with z a vector containing "old" phase values from the initial screen, that are listed thanks to the indexes in istencil.

Examples n = 32; AB, n, A, B, istencil; p = array(0.0,n,n); p1 = extrude(p,r0,A,B,istencil); pli p1

SEE AB() createStencil() Cxx() Cxz() Czz()

Definition at line 277 of file iterkolmo.py.

277 
278  amplitude = r0**(-5. / 6)
279  n = p.shape[0]
280  z = p.flatten()[istencil]
281  zref = p[0, n - 1]
282  z -= zref
283  newColumn = np.dot(A, z) + np.dot(B, np.random.normal(0, 1, n) * amplitude) + zref
284  p1 = np.zeros((n, n), dtype=np.float32)
285  p1[:, 0:n - 1] = p[:, 1:]
286  p1[:, n - 1] = newColumn
287 
288  return p1
289 
290 
291 def phase_struct(r, L0=None):
Here is the caller graph for this function:

◆ macdo_x56()

def shesha.util.iterkolmo.macdo_x56 (   x,
  k = 10 
)

Computation of the function f(x) = x^(5/6)*K_{5/6}(x) using a series for the esimation of K_{5/6}, taken from Rod Conan thesis : K_a(x)=1/2 \sum_{n=0}^\infty \frac{(-1)^n}{n!} \left(\Gamma(-n-a) (x/2)^{2n+a} + \Gamma(-n+a) (x/2)^{2n-a} \right) , with a = 5/6.

Setting x22 = (x/2)^2, setting uda = (1/2)^a, and multiplying by x^a, this becomes : x^a * Ka(x) = 0.5 $ -1^n / n! [ G(-n-a).uda x22^(n+a) + G(-n+a)/uda x22^n ] Then we use the following recurrence formulae on the following quantities : G(-(n+1)-a) = G(-n-a) / -a-n-1 G(-(n+1)+a) = G(-n+a) / a-n-1 (n+1)! = n! * (n+1) x22^(n+1) = x22^n * x22 and at each iteration on n, one will use the values already computed at step (n-1). The values of G(a) and G(-a) are hardcoded instead of being computed.

The first term of the series has also been skipped, as it vanishes with another term in the expression of Dphi.

Definition at line 388 of file iterkolmo.py.

388 
389  a = 5. / 6.
390  fn = 1. # initialisation factorielle 0!=1
391  x2a = x**(2. * a)
392  x22 = x * x / 4. # (x/2)^2
393  x2n = 0.5 # init (1/2) * x^0
394  Ga = 2.01126983599717856777 # Gamma(a) / (1/2)^a
395  Gma = -3.74878707653729348337 # Gamma(-a) * (1/2.)^a
396  s = np.zeros(x.shape) # array(0.0, dimsof(x));
397  for n in range(k + 1): # (n=0; n<=k; n++) {
398  dd = Gma * x2a
399  if n:
400  dd += Ga
401  dd *= x2n
402  dd /= fn
403  # addition to s, with multiplication by (-1)^n
404  if (n % 2):
405  s -= dd
406  else:
407  s += dd
408  # prepare recurrence iteration for next step
409  if (n < k):
410  fn *= n + 1 # factorial
411  Gma /= -a - n - 1 # gamma function
412  Ga /= a - n - 1 # idem
413  x2n *= x22 # x^n
414 
415  return s
416 
417 
418 def create_screen_assist(screen_size, L0, r0):
419  """
420  screen_size : screen size (in pixels)
421  L0 : L0 in pixel
422  r0 : total r0 @ 0.5 microns
423  """
Here is the caller graph for this function:

◆ phase_struct()

def shesha.util.iterkolmo.phase_struct (   r,
  L0 = None 
)

TODO docstring.

Definition at line 294 of file iterkolmo.py.

294  if L0 is None:
295  return 6.88 * r**(5. / 6.)
296  else:
297  return rodconan(np.sqrt(r), L0)
298 
299 
300 def rodconan(r, L0):
Here is the call graph for this function:
Here is the caller graph for this function:

◆ rodconan()

def shesha.util.iterkolmo.rodconan (   r,
  L0 
)

The phase structure function is computed from the expression Dphi(r) = k1 * L0^(5.

/3) * (k2 - (2.pi.r/L0)^5/6 K_{5/6}(2.pi.r/L0))

For small r, the expression is computed from a development of K_5/6 near 0. The value of k2 is not used, as this same value appears in the series and cancels with k2. For large r, the expression is taken from an asymptotic form.

Definition at line 309 of file iterkolmo.py.

309 
310  # k1 is the value of :
311  # 2*gamma_R(11./6)*2^(-5./6)*pi^(-8./3)*(24*gamma_R(6./5)/5.)^(5./6);
312  k1 = 0.1716613621245709486
313  dprf0 = (2 * np.pi / L0) * r
314 
315  Xlim = 0.75 * 2 * np.pi
316  ilarge = np.where(dprf0 > Xlim)
317  ismall = np.where(dprf0 <= Xlim)
318 
319  res = r * 0.
320  """
321  # TODO those lines have been changed (cf trunk/yoga_ao/yorick/yoga_turbu.i l 258->264)
322  if((ilarge[0].size > 0)and(ismall[0].size == 0)):
323  res[ilarge] = asymp_macdo(dprf0[ilarge])
324  # print("simulation atmos with asymptotic MacDonald function")
325  elif((ismall[0].size > 0)and(ilarge[0].size == 0)):
326  res[ismall] = -macdo_x56(dprf0[ismall])
327  # print("simulation atmos with x56 MacDonald function")
328  elif((ismall[0].size > 0)and(ilarge[0].size > 0)):
329  res[ismall] = -macdo_x56(dprf0[ismall])
330  res[ilarge] = asymp_macdo(dprf0[ilarge])
331  """
332  if (ilarge[0].size > 0):
333  res[ilarge] = asymp_macdo(dprf0[ilarge])
334  if (ismall[0].size > 0):
335  res[ismall] = -macdo_x56(dprf0[ismall])
336 
337  return k1 * L0**(5. / 3.) * res
338 
339 
340 def asymp_macdo(x):
Here is the call graph for this function:
Here is the caller graph for this function:

◆ stencil_size()

def shesha.util.iterkolmo.stencil_size (   n)

TODO docstring.

Definition at line 78 of file iterkolmo.py.

78  """
79  stencil = np.zeros((n, n))
80  ns = int(np.log2(n + 1) + 1)
81 
82  stencil[:, 0] = 1
83  for i in range(2, ns):
84  stencil[::2**(i - 1), 2**(i - 2)] = 1
85  stencil.itemset(2**(i - 1) - 1, 1)
86 
87  # i=ns
88  stencil.itemset(2**(ns - 1) - 1, 1)
89  for i in range(0, n, 2**(ns - 1)):
90  stencil.itemset(2**(ns - 2) + i * n, 1)
91  # i=ns+1
92  stencil.itemset(2**(ns) - 1, 1)
93  for i in range(0, n, 2**(ns)):
94  stencil.itemset(2**(ns - 1) + i * n, 1)
95 
96  return np.sum(stencil)
97 
98 
Here is the caller graph for this function:

◆ stencil_size_array()

def shesha.util.iterkolmo.stencil_size_array (   size)

Compute_size2(np.ndarray[ndim=1, dtype=np.int64_t] size)

Compute the size of a stencil, given the screen size

:parameters:

size (np.ndarray[ndim=1,dtype=np.int64_t]) :screen size

Definition at line 107 of file iterkolmo.py.

107  """
108  stsize = np.zeros(len(size), dtype=np.int64)
109 
110  for i in range(len(size)):
111  stsize[i] = stencil_size(size[i])
112  return stsize
113 
114 
Here is the call graph for this function: