40 from scipy
import interpolate
42 from typing
import Tuple
50 cobs: (float) : central obstruction
54 d = (1. - cobs * cobs) / nr
55 rad2 = cobs**2 + d / 16. + d * (np.arange(nr, dtype=np.float32))
60 def make_kernels(cobs: float, nr: int, radp: np.ndarray, kl_type: bytes,
61 outscl: float = 3.) -> np.ndarray:
63 This routine generates the kernel used to find the KL modes.
64 The kernel constructed here should be simply a discretization
65 of the continuous kernel. It needs rescaling before it is treated
66 as a matrix for finding the eigen-values. The outer scale
67 should be in units of the diameter of the telescope.
73 cobs : (float): central obstruction
79 kl_type : (bytes) : "kolmo" or "karman"
81 outscl : (float) : outter scale for Von Karman spectrum
88 kers = np.zeros((nth, nr, nr), dtype=np.float32)
89 cth = np.cos((np.arange(nth, dtype=np.float32)) * (2. * np.pi / nth))
90 dth = 2. * np.pi / nth
91 fnorm = -1. / (2 * np.pi * (1. - cobs**2)) * 0.5
94 for j
in range(i + 1):
95 te = 0.5 * np.sqrt(radp[i]**2 + radp[j]**2 - (2 * radp[i] * radp[j]) * cth)
97 if (kl_type == scons.KLType.KOLMO):
99 te = 6.88 * te**(5. / 3.)
101 elif (kl_type == scons.KLType.KARMAN):
103 te = 6.88 * te**(5. / 3.) * (1 - 1.485 * (te / outscl)**
104 (1. / 3.) + 5.383 * (te / outscl)**
105 (2) - 6.281 * (te / outscl)**(7. / 3.))
109 raise TypeError(
"kl type unknown")
111 f = np.fft.fft(te, axis=-1)
112 kelt = fnorm * dth * np.float32(f.real)
113 kers[:, i, j] = kers[:, j, i] = kelt
128 s = np.zeros((nr, nr), dtype=np.float32)
129 for j
in range(nr - 1):
130 rnm = 1. / np.sqrt(np.float32((j + 1) * (j + 2)))
132 s[j + 1, j] = -1 * (j + 1) * rnm
134 rnm = 1. / np.sqrt(nr)
153 azbas = np.zeros((npp, np.int32(1 + nord)), dtype=np.float32)
154 th = np.arange(npp, dtype=np.float32) * (2. * np.pi / npp)
157 for i
in np.arange(1, nord, 2):
158 azbas[:, np.int32(i)] = np.cos((np.int32(i) // 2 + 1) * th)
159 for i
in np.arange(2, nord, 2):
160 azbas[:, np.int32(i)] = np.sin((np.int32(i) // 2) * th)
165 def radii(nr: int, npp: int, cobs: float) -> np.ndarray:
167 This routine generates an nr x npp array with npp copies of the
168 radial coordinate array. Radial coordinate span the range from
169 r=cobs to r=1 with successive annuli having equal areas (ie, the
170 area between cobs and 1 is divided into nr equal rings, and the
171 points are positioned at the half-area mark on each ring). There
172 are no points on the border.
182 cobs: (float) : central obstruction
189 r2 = cobs**2 + (np.arange(nr, dtype=np.float) + 0.) / nr * (1.0 - cobs**2)
191 r = np.transpose(np.tile(rs, (npp, 1)))
200 def polang(r: np.ndarray) -> np.ndarray:
202 This routine generates an array with the same dimensions as r,
203 but containing the azimuthal values for a polar coordinate system.
218 phi1 = np.arange(np1, dtype=np.float) / float(np1) * 2. * np.pi
219 p1, p2 = np.meshgrid(np.ones(nr), phi1)
229 def setpincs(ax: np.ndarray, ay: np.ndarray, px: np.ndarray, py: np.ndarray,
230 cobs: float) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
232 This routine determines a set of squares for interpolating
233 from cartesian to polar coordinates, using only those points
249 cobs: (float) : central obstruction
265 dcar = (ax[nc - 1, 0] - ax[0, 0]) / (nc - 1)
267 rlx = (px - ofcar) / dcar
268 rly = (py - ofcar) / dcar
274 pincx = np.zeros((4, lx.shape[0], lx.shape[1]))
275 pincx[[1, 2], :, :] = lx + 1
276 pincx[[0, 3], :, :] = lx
278 pincy = np.zeros((4, ly.shape[0], ly.shape[1]))
279 pincy[[0, 1], :, :] = ly
280 pincy[[2, 3], :, :] = ly + 1
282 pincw = np.zeros((4, shx.shape[0], shx.shape[1]))
283 pincw[0, :, :] = (1 - shx) * (1 - shy)
284 pincw[1, :, :] = shx * (1 - shy)
285 pincw[2, :, :] = shx * shy
286 pincw[3, :, :] = (1 - shx) * shy
289 axyinap = np.clip(axy, cobs**2. + 1.e-3, 0.999)
294 for z
in range(pincw.shape[0]):
295 for i
in range(pincw.shape[1]):
296 for j
in range(pincw.shape[2]):
297 pincw[z, i, j] = pincw[z, i, j] * axyinap[np.int32(pincx[z, i, j]),
298 np.int32(pincy[z, i, j])]
300 pincw = pincw * np.tile(1.0 / np.sum(pincw, axis=0), (4, 1, 1))
302 return pincx, pincy, pincw
305 def pcgeom(nr, npp, cobs, ncp, ncmar):
307 This routine builds a geom_struct. px and py are the x and y
308 coordinates of points in the polar arrays. cr and cp are the
309 r and phi coordinates of points in the cartesian grids. ncmar
310 allows the possibility that there is a margin of ncmar points
311 in the cartesian arrays outside the region of interest
321 cobs: (float) : central obstruction
349 nused = ncp - 2 * ncmar
351 hw = np.float(ncp - 1) / 2.
353 r =
radii(nr, npp, cobs)
361 np.arange(int(ncp)**2, dtype=np.float) + 1, (int(ncp), int(ncp)), order=
'F')
362 ax = np.float32(ax - 1) % ncp - 0.5 * (ncp - 1)
363 ax = ax / (0.5 * nused)
364 ay = np.transpose(ax)
366 pincx, pincy, pincw =
setpincs(ax, ay, px0, py0, cobs)
369 cr2 = (ax**2 + ay**2)
370 ap = np.clip(cr2, cobs**2 + 1.e-3, 0.999)
372 cr = (cr2 - cobs**2) / (1 - cobs**2) * nr
373 cp = (np.arctan2(ay, ax) + dpi) % dpi
374 cp = (npp / dpi) * cp
376 cr = np.clip(cr, 1.e-3, nr - 1.001)
378 cp = np.clip(cp, 1.e-3, npp - 1.001)
382 return ncp, ncmar, px, py, cr, cp, pincx, pincy, pincw, ap
385 def set_pctr(dim: int, nr, npp, nkl: int, cobs: float, nord, ncmar=
None, ncp=
None):
387 This routine calls pcgeom to build a geom_struct with the
388 right initializations. bas is a gkl_basis_struct built with
437 ncp, ncmar, px, py, cr, cp, pincx, pincy, pincw, ap =
pcgeom(
438 nr, npp, cobs, ncp, ncmar)
439 return ncp, ncmar, px, py, cr, cp, pincx, pincy, pincw, ap
442 def gkl_fcom(kers: np.ndarray, cobs: float, nf: int):
444 This routine does the work : finding the eigenvalues and
445 corresponding eigenvectors. Sort them and select the right
446 one. It returns the KL modes : in polar coordinates : rabas
447 as well as the associated variance : evals. It also returns
448 a bunch of indices used to recover the modes in cartesian
449 coordinates (nord, npo and ordd).
453 kerns : (np.ndarray[ndim= ,dtype=np.float32]) :
455 cobs : (float) : central obstruction
464 fktom = (1. - cobs**2) / nr
466 evs = np.zeros((nr, nt), dtype=np.float32)
482 btemp = (ts.dot(zom).dot(s))[0:nr - 1, 0:nr - 1]
485 v0, newev, vt = np.linalg.svd(
486 fktom * btemp, full_matrices=
True
489 v1 = np.zeros((nr, nr), dtype=np.float32)
490 v1[0:nr - 1, 0:nr - 1] = v0
491 v1[nr - 1, nr - 1] = 1
494 newev = np.concatenate((newev, [0]))
496 evs[:, nxt] = np.float32(newev)
497 kers[nxt, :, :] = np.sqrt(nr) * vs
501 vs, newev, vt = np.linalg.svd(fktom * kers[nxt, :, :], full_matrices=
True)
503 evs[:, nxt] = np.float32(newev)
504 kers[nxt, :, :] = np.sqrt(2. * nr) * vs
505 mxn = max(np.float32(newev))
506 egtmxn = np.floor(evs[:, 0:nxt + 1] > mxn)
508 if ((2 * np.sum(egtmxn) - np.sum(egtmxn[:, 0])) >= nkl):
512 kers = kers[0:nus + 1, :, :]
516 evs = np.reshape(evs[:, 0:nus + 1], nr * (nus + 1), order=
'F')
517 a = np.argsort(-1. * evs)[0:nkl]
527 oind = np.zeros(nkl + 1, dtype=np.int32)
543 tord = (oind) // nr + 1
545 odd = np.arange(nkl, dtype=np.int32) % 2
546 pio = (oind) % nr + 1
549 ordd = 2 * (tord - 1) - np.floor((tord > 1) & (odd)) + 1
553 rabas = np.zeros((nr, nkl), dtype=np.float32)
555 npo = np.zeros(sizenpo, dtype=np.int32)
558 npo[np.int32(ordd[i]) - 1] = npo[np.int32(ordd[i]) - 1] + 1
559 rabas[:, i] = kers[tord[i] - 1, :, pio[i] - 1]
561 return evals, nord, npo, ordd, rabas
580 raise TypeError(
"kl funct order it's so big")
584 ordi = np.int32(ordd[i])
586 azbasi = np.transpose(azbas)
587 azbasi = azbasi[ordi, :]
589 sf1 = np.zeros((nr, npp), dtype=np.float64)
593 sf2 = np.zeros((npp, nr), dtype=np.float64)
597 sf = sf1 * np.transpose(sf2)
616 r = np.arange(nr, dtype=np.float64)
617 phi = np.arange(npp, dtype=np.float64)
618 tab_phi, tab_r = np.meshgrid(phi, r)
619 tab_x = (tab_r / (nr)) * np.cos((tab_phi / (npp)) * 2 * np.pi)
620 tab_y = (tab_r / (nr)) * np.sin((tab_phi / (npp)) * 2 * np.pi)
622 newx = np.linspace(-1, 1, ncp)
623 newy = np.linspace(-1, 1, ncp)
624 tx, ty = np.meshgrid(newx, newy)
626 cd = interpolate.griddata((tab_r.flatten(), tab_phi.flatten()), pol.flatten(),
627 (cr, cp), method=
'cubic')
628 cdf = interpolate.griddata((tab_r.flatten(
"F"), tab_phi.flatten(
"F")),
629 pol.flatten(
"F"), (cr, cp), method=
'cubic')
630 cdxy = interpolate.griddata((tab_y.flatten(), tab_x.flatten()), pol.flatten(),
631 (tx, ty), method=
'cubic')
647 tab_kl = np.zeros((nkl, ncp, ncp), dtype=np.float64)
648 tab_klf = np.zeros((nkl, ncp, ncp), dtype=np.float64)
649 tab_klxy = np.zeros((nkl, ncp, ncp), dtype=np.float64)
653 tab_kl[i, :, :], tab_klf[i, :, :], tab_klxy[i, :, :] =
pol2car(
656 return tab_kl, tab_klf, tab_klxy