44 Zx = np.array(np.tile(np.arange(n), (n, 1)), dtype=np.float64) + 1
47 Xx = np.array(np.zeros(n) + n + 1, dtype=np.float64)
48 Xy = np.array(np.arange(n) + 1, dtype=np.float64)
50 ns = int(np.log2(n + 1) + 1)
51 stencil = np.zeros((n, n))
54 for i
in range(2, ns):
55 stencil[::2**(i - 1), 2**(i - 2)] = 1
56 stencil.itemset(2**(i - 1) - 1, 1)
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)
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)
68 stencil = np.roll(stencil, n // 2, axis=0)
69 stencil = np.fliplr(stencil)
71 istencil = np.where(stencil.flatten() != 0)[0]
73 return Zx, Zy, Xx, Xy, istencil
79 stencil = np.zeros((n, n))
80 ns = int(np.log2(n + 1) + 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)
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)
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)
96 return np.sum(stencil)
100 """ Compute_size2(np.ndarray[ndim=1, dtype=np.int64_t] size)
102 Compute the size of a stencil, given the screen size
106 size: (np.ndarray[ndim=1,dtype=np.int64_t]) :screen size
108 stsize = np.zeros(len(size), dtype=np.int64)
110 for i
in range(len(size)):
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.
119 The known values z are the values of the phase screen that are pointed by
120 the stencil indexes (istencil)
124 size2 = istencil.shape[0]
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))
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)
136 xz = -
phase_struct((Xx_r.T - Zx_r)**2 + (Xy_r.T - Zy_r)**2, L0)
138 xz += np.resize(tmp, (size2, size)).T + np.resize(tmp2, (size, size2))
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).
150 Xx_r = np.resize(Xx, (size, size))
151 Xy_r = np.resize(Xy, (size, size))
153 tmp = np.resize(
phase_struct((Zxn - Xx)**2 + (Zyn - Xy)**2, L0), (size, size))
155 xx = -
phase_struct((Xx_r - Xx_r.T)**2 + (Xy_r - Xy_r.T)**2, L0)
165 def Czz(n, Zx, Zy, ist, L0):
166 """ Czz computes the covariance matrix of the already known phase values z.
168 The known values z are the values of the phase screen that are pointed by
169 the stencil indexes (istencil)
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))
179 phase_struct((Zx[0, n - 1] - Zx_s)**2 + (Zy[0, n - 1] - Zy_s)**2, L0),
182 zz = -
phase_struct((Zx_r - Zx_r.T) ** 2 + (Zy_r - Zy_r.T) ** 2, L0) + \
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.
195 The method used is described by Fried & Clark in JOSA A, vol 25, no 2, p463, Feb 2008.
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.
201 SEE ALSO: extrude createStencil Cxx Cxz Czz
209 zz =
Czz(n, Zx, Zy, istencil, L0)
212 xz =
Cxz(n, Zx, Zy, Xx, Xy, istencil, L0)
215 xx =
Cxx(n, Zx[0, n - 1], Zy[0, n - 1], Xx, Xy, L0)
217 U, s, V = np.linalg.svd(zz)
222 zz1 = np.dot(np.dot(U, np.diag(s1)), V)
233 bbt = xx - np.dot(A, xz.T)
236 U1, l, V1 = np.linalg.svd(bbt)
239 B = np.dot(U1, np.sqrt(np.diag(l)))
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]:]
247 isty = (n * n - 1) - isty
249 istencil = (n * n - 1) - istencil
251 return np.asfortranarray(A.astype(np.float32)), np.asfortranarray(
252 B.astype(np.float32)), istencil.astype(np.uint32), isty.astype(np.uint32)
255 def extrude(p, r0, A, B, istencil):
256 """ DOCUMENT p1 = extrude(p,r0,A,B,istencil)
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
262 The method used is described by Fried & Clark in JOSA A, vol 25, no 2, p463, Feb 2008.
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.
270 AB, n, A, B, istencil;
272 p1 = extrude(p,r0,A,B,istencil);
275 SEE ALSO: AB() createStencil() Cxx() Cxz() Czz()
278 amplitude = r0**(-5. / 6)
280 z = p.flatten()[istencil]
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
295 return 6.88 * r**(5. / 6.)
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))
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.
312 k1 = 0.1716613621245709486
313 dprf0 = (2 * np.pi / L0) * r
315 Xlim = 0.75 * 2 * np.pi
316 ilarge = np.where(dprf0 > Xlim)
317 ismall = np.where(dprf0 <= Xlim)
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])
332 if (ilarge[0].size > 0):
334 if (ismall[0].size > 0):
337 return k1 * L0**(5. / 3.) * res
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.
349 - This function makes a floating point interrupt for x=0
350 and should not be used in this case.
352 - Works only for x>0.
356 k2 = 1.00563491799858928388289314170833
357 k3 = 1.25331413731550012081
358 a1 = 0.22222222222222222222
359 a2 = -0.08641975308641974829
360 a3 = 0.08001828989483310284
362 res = k2 - k3 * np.exp(-x) * x**(1 / 3.) * (1.0 + x_1 * (a1 + x_1 * (a2 + x_1 * a3)))
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) ,
374 Setting x22 = (x/2)^2, setting uda = (1/2)^a, and multiplying by x^a,
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
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.
385 The first term of the series has also been skipped, as it
386 vanishes with another term in the expression of Dphi.
394 Ga = 2.01126983599717856777
395 Gma = -3.74878707653729348337
396 s = np.zeros(x.shape)
397 for n
in range(k + 1):
420 screen_size : screen size (in pixels)
422 r0 : total r0 @ 0.5 microns
424 A, B, istx, isty =
AB(screen_size, L0)
425 phase = np.zeros((screen_size, screen_size))
433 for i
in range(2 * screen_size):
434 phase =
extrude(phase, r0, A, B, istx)
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)
446 creates a phase screen and fill it with turbulence
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
456 screen = np.zeros((screen_size, screen_size),
458 for i
in range(2 * screen_size):
459 screen =
extrude(screen, r0 / pupixsize, A, B, ist)