COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
iterkolmo.py
1 
37 
38 import numpy as np
39 
40 
41 def create_stencil(n):
42  """ TODO: docstring
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 
76 def stencil_size(n):
77  """ TODO: docstring
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 
99 def stencil_size_array(size):
100  """ Compute_size2(np.ndarray[ndim=1, dtype=np.int64_t] size)
101 
102  Compute the size of a stencil, given the screen size
103 
104  :parameters:
105 
106  size: (np.ndarray[ndim=1,dtype=np.int64_t]) :screen size
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 
115 def Cxz(n, Zx, Zy, Xx, Xy, istencil, L0):
116  """ Cxz computes the covariance matrix between the new phase vector x (new
117  column for the phase screen), and the already known phase values z.
118 
119  The known values z are the values of the phase screen that are pointed by
120  the stencil indexes (istencil)
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 
145 def Cxx(n, Zxn, Zyn, Xx, Xy, L0):
146  """ Cxx computes the covariance matrix of the new phase vector x (new
147  column for the phase screen).
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 
165 def Czz(n, Zx, Zy, ist, L0):
166  """ Czz computes the covariance matrix of the already known phase values z.
167 
168  The known values z are the values of the phase screen that are pointed by
169  the stencil indexes (istencil)
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 
190 def AB(n, L0, deltax, deltay, rank=0):
191  """ DOCUMENT AB, n, A, B, istencil
192  This function initializes some matrices A, B and a list of stencil indexes
193  istencil for iterative extrusion of a phase screen.
194 
195  The method used is described by Fried & Clark in JOSA A, vol 25, no 2, p463, Feb 2008.
196  The iteration is :
197  x = A(z-zRef) + B.noise + zRef
198  with z a vector containing "old" phase values from the initial screen, that are listed
199  thanks to the indexes in istencil.
200 
201  SEE ALSO: extrude createStencil Cxx Cxz Czz
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 
255 def extrude(p, r0, A, B, istencil):
256  """ DOCUMENT p1 = extrude(p,r0,A,B,istencil)
257 
258  Extrudes a phase screen p1 from initial phase screen p.
259  p1 prolongates p by 1 column on the right end.
260  r0 is expressed in pixels
261 
262  The method used is described by Fried & Clark in JOSA A, vol 25, no 2, p463, Feb 2008.
263  The iteration is :
264  x = A(z-zRef) + B.noise + zRef
265  with z a vector containing "old" phase values from the initial screen, that are listed
266  thanks to the indexes in istencil.
267 
268  Examples
269  n = 32;
270  AB, n, A, B, istencil;
271  p = array(0.0,n,n);
272  p1 = extrude(p,r0,A,B,istencil);
273  pli, p1
274 
275  SEE ALSO: AB() createStencil() Cxx() Cxz() Czz()
276  """
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):
292  """ TODO: docstring
293  """
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):
301  """ The phase structure function is computed from the expression
302  Dphi(r) = k1 * L0^(5./3) * (k2 - (2.pi.r/L0)^5/6 K_{5/6}(2.pi.r/L0))
303 
304  For small r, the expression is computed from a development of
305  K_5/6 near 0. The value of k2 is not used, as this same value
306  appears in the series and cancels with k2.
307  For large r, the expression is taken from an asymptotic form.
308  """
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):
341  """ Computes a term involved in the computation of the phase struct
342  function with a finite outer scale according to the Von-Karman
343  model. The term involves the MacDonald function (modified bessel
344  function of second kind) K_{5/6}(x), and the algorithm uses the
345  asymptotic form for x ~ infinity.
346 
347  Warnings :
348 
349  - This function makes a floating point interrupt for x=0
350  and should not be used in this case.
351 
352  - Works only for x>0.
353  """
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):
367  """ Computation of the function
368  f(x) = x^(5/6)*K_{5/6}(x)
369  using a series for the esimation of K_{5/6}, taken from Rod Conan thesis :
370  K_a(x)=1/2 \sum_{n=0}^\infty \frac{(-1)^n}{n!}
371  \left(\Gamma(-n-a) (x/2)^{2n+a} + \Gamma(-n+a) (x/2)^{2n-a} \right) ,
372  with a = 5/6.
373 
374  Setting x22 = (x/2)^2, setting uda = (1/2)^a, and multiplying by x^a,
375  this becomes :
376  x^a * Ka(x) = 0.5 $ -1^n / n! [ G(-n-a).uda x22^(n+a) + G(-n+a)/uda x22^n ]
377  Then we use the following recurrence formulae on the following quantities :
378  G(-(n+1)-a) = G(-n-a) / -a-n-1
379  G(-(n+1)+a) = G(-n+a) / a-n-1
380  (n+1)! = n! * (n+1)
381  x22^(n+1) = x22^n * x22
382  and at each iteration on n, one will use the values already computed at step (n-1).
383  The values of G(a) and G(-a) are hardcoded instead of being computed.
384 
385  The first term of the series has also been skipped, as it
386  vanishes with another term in the expression of Dphi.
387  """
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  """
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):
443  """ DOCUMENT create_screen
444  screen = create_screen(r0,pupixsize,screen_size,&A,&B,&ist)
445 
446  creates a phase screen and fill it with turbulence
447  r0 : total r0 @ 0.5m
448  pupixsize : pupil pixel size (in meters)
449  screen_size : screen size (in pixels)
450  A : A array for future extrude
451  B : B array for future extrude
452  ist : istencil array for future extrude
453  """
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
shesha.util.iterkolmo.create_screen
def create_screen(r0, pupixsize, screen_size, L0, A, B, ist)
DOCUMENT create_screen screen = create_screen(r0,pupixsize,screen_size,&A,&B,&ist)
Definition: iterkolmo.py:454
shesha.util.iterkolmo.stencil_size
def stencil_size(n)
TODO docstring.
Definition: iterkolmo.py:78
shesha.util.iterkolmo.create_stencil
def create_stencil(n)
TODO docstring.
Definition: iterkolmo.py:43
shesha.util.iterkolmo.Cxx
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).
Definition: iterkolmo.py:148
shesha.util.iterkolmo.AB
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 ind...
Definition: iterkolmo.py:202
shesha.util.iterkolmo.stencil_size_array
def stencil_size_array(size)
Compute_size2(np.ndarray[ndim=1, dtype=np.int64_t] size)
Definition: iterkolmo.py:107
shesha.util.iterkolmo.rodconan
def rodconan(r, L0)
The phase structure function is computed from the expression Dphi(r) = k1 * L0^(5.
Definition: iterkolmo.py:309
shesha.util.iterkolmo.asymp_macdo
def asymp_macdo(x)
Computes a term involved in the computation of the phase struct function with a finite outer scale ac...
Definition: iterkolmo.py:354
shesha.util.iterkolmo.phase_struct
def phase_struct(r, L0=None)
TODO docstring.
Definition: iterkolmo.py:294
shesha.util.iterkolmo.macdo_x56
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},...
Definition: iterkolmo.py:388
shesha.util.iterkolmo.create_screen_assist
def create_screen_assist(screen_size, L0, r0)
Definition: iterkolmo.py:424
shesha.util.iterkolmo.Czz
def Czz(n, Zx, Zy, ist, L0)
Czz computes the covariance matrix of the already known phase values z.
Definition: iterkolmo.py:170
shesha.util.iterkolmo.extrude
def extrude(p, r0, A, B, istencil)
DOCUMENT p1 = extrude(p,r0,A,B,istencil)
Definition: iterkolmo.py:277
shesha.util.iterkolmo.Cxz
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),...
Definition: iterkolmo.py:121