knm.py 20.5 KB
 Daniel Brown committed Jul 22, 2014 1 from itertools import combinations_with_replacement as combinations  Daniel Brown committed Jan 06, 2017 2 from pykat.optics.gaussian_beams import BeamParam, HG_mode  Daniel Brown committed Jul 22, 2014 3 from pykat.exceptions import BasePyKatException  Andreas Freise committed Dec 20, 2014 4 from pykat.optics.romhom import u_star_u  Andreas Freise committed Dec 14, 2014 5 from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar  Daniel Brown committed Sep 04, 2014 6 7 from scipy.interpolate import interp2d from scipy.integrate import dblquad  8 from pykat.optics.romhom import ROMWeights  Daniel Brown committed Mar 24, 2015 9 from math import factorial  Andreas Freise committed May 14, 2016 10 from pykat.math.hermite import hermite  Daniel Brown committed Mar 24, 2015 11 12 from scipy.misc import comb from scipy.integrate import newton_cotes  Andreas Freise committed May 14, 2016 13 from pykat.math import newton_weights  Daniel Brown committed Jul 22, 2014 14   Daniel Brown committed Aug 12, 2014 15 import time  Andreas Freise committed Dec 20, 2014 16 import pykat.optics.maps  Daniel Brown committed Jul 22, 2014 17 18 19 20 21 import os.path import numpy as np import pykat import collections import math  Daniel Brown committed Aug 12, 2014 22 import cmath  Daniel Brown committed Aug 01, 2014 23   Daniel Brown committed Sep 04, 2014 24 def makeCouplingMatrix(max_order, Neven=True, Nodd=True, Meven=True, Modd=True):  Daniel Brown committed Aug 01, 2014 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39  max_order = int(max_order) c = [] for n in range(0, max_order+1): for m in range(0, max_order+1): if n+m <= max_order: c.append([n,m]) M = [] for i in c: row = [] for j in c: e = list(i) e.extend(j)  Daniel Brown committed Sep 04, 2014 40 41 42 43 44 45  if not Neven and (e[0]-e[2]) % 2 == 0: continue if not Nodd and (e[0]-e[2]) % 2 == 1: continue if not Meven and (e[1]-e[3]) % 2 == 0: continue if not Modd and (e[1]-e[3]) % 2 == 1: continue  Daniel Brown committed Aug 01, 2014 46 47  row.append(e)  Daniel Brown committed Sep 04, 2014 48 49  M.append(np.array(row).squeeze())  Daniel Brown committed Aug 01, 2014 50   Daniel Brown committed Aug 12, 2014 51  return np.array(M)  Daniel Brown committed Jul 22, 2014 52   Daniel Brown committed Sep 04, 2014 53 54 def adaptive_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, smap=None, delta=(0,0), params={}):  Daniel Brown committed Apr 20, 2015 55  if q1y is None:  Daniel Brown committed Sep 04, 2014 56 57  q1y = q1  Daniel Brown committed Apr 20, 2015 58  if q2y is None:  Daniel Brown committed Sep 04, 2014 59 60 61 62 63 64 65 66 67  q2y = q2 if "epsabs" not in params: params["epsabs"] = 1e-6 if "epsrel" not in params: params["epsrel"] = 1e-6 if "usepolar" not in params: params["usepolar"] = False if len(mode_in) != 2 or len(mode_out) != 2: raise BasePyKatException("Both mode in and out should be a container with modes [n m]")  Andreas Freise committed May 14, 2016 68 69  Hg_in = HG_mode(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1]) Hg_out = HG_mode(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])  Daniel Brown committed Sep 04, 2014 70 71 72 73  Nfuncs = [] Nfuncs.append(0)  Daniel Brown committed Apr 20, 2015 74  if smap is not None:  Daniel Brown committed Sep 04, 2014 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151  if not params["usepolar"]: xlims = (min(smap.x), max(smap.x)) ylims = (min(smap.y), max(smap.y)) def Rfunc(y,x): Nfuncs[-1] += len(x) return (Hg_in.Unm(x+delta[0], y+delta[1]) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).real def Ifunc(y,x): Nfuncs[-1] += len(x) return (Hg_in.Unm(x+delta[0], y+delta[1]) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).imag else: xlims = (0, 2*math.pi) ylims = (0, params["aperture"]) def Rfunc(r, phi): Nfuncs[-1] += len(x) x = r*np.cos(phi) y = r*np.sin(phi) return (r * Hg_in.Unm(x, y) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).real def Ifunc(r, phi): Nfuncs[-1] += len(x) x = r*np.cos(phi) y = r*np.sin(phi) return (r * Hg_in.Unm(x, y) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).imag else: if not params["usepolar"]: _x = 4 * math.sqrt(1+max(mode_in[0],mode_in[1])) * q1.w _y = 4 * math.sqrt(1+max(mode_in[0],mode_in[1])) * q1y.w xlims = (-_x, _x) ylims = (-_y, _y) def Rfunc(y, x): Nfuncs[-1] += len(r) return (Hg_in.Unm(x+delta[0], y+delta[1]) * Hg_out.Unm(x, y).conjugate()).real def Ifunc(y,x): Nfuncs[-1] += len(r) return (Hg_in.Unm(x+delta[0], y+delta[1]) * Hg_out.Unm(x, y).conjugate()).imag else: xlims = (0, 2*math.pi) ylims = (0, params["aperture"]) def Rfunc(r, phi): if hasattr(r, "__len__"): Nfuncs[-1] += len(r) else: Nfuncs[-1] += 1 x = r*np.cos(phi) y = r*np.sin(phi) return (r * Hg_in.Unm(x, y) * Hg_out.Unm(x, y).conjugate()).real def Ifunc(r, phi): if hasattr(r, "__len__"): Nfuncs[-1] += len(r) else: Nfuncs[-1] += 1 x = r*np.cos(phi) y = r*np.sin(phi) return (r * Hg_in.Unm(x, y) * Hg_out.Unm(x, y).conjugate()).imag R, errR = dblquad(Rfunc, xlims[0], xlims[1], lambda y: ylims[0], lambda y: ylims[1], epsabs=params["epsabs"], epsrel=params["epsrel"]) I, errI = dblquad(Ifunc, xlims[0], xlims[1], lambda y: ylims[0], lambda y: ylims[1], epsabs=params["epsabs"], epsrel=params["epsrel"]) params["Nfuncs"] = Nfuncs[0] params["errors"] = (errR, errI) return R + 1j * I  Daniel Brown committed Mar 24, 2015 152 153 def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None, delta=(0,0), params={}, newtonCotesOrder=0):  Daniel Brown committed Jul 22, 2014 154   Daniel Brown committed Mar 24, 2015 155  if Axy is None:  Daniel Brown committed Jul 22, 2014 156 157  Axy == np.ones((len(x), len(y)))  Daniel Brown committed Mar 24, 2015 158  if q1y is None:  Daniel Brown committed Jul 22, 2014 159 160  q1y = q1  Daniel Brown committed Mar 24, 2015 161  if q2y is None:  Daniel Brown committed Jul 22, 2014 162 163 164 165  q2y = q2 if len(mode_in) != 2 or len(mode_out) != 2: raise BasePyKatException("Both mode in and out should be a container with modes [n m]")  Daniel Brown committed Aug 01, 2014 166   Daniel Brown committed Jul 22, 2014 167  dx = abs(x[1] - x[0])  Daniel Brown committed Aug 01, 2014 168 169  dy = abs(y[1] - y[0])  Daniel Brown committed Apr 20, 2015 170  if cache is None:  Andreas Freise committed May 14, 2016 171 172  Hg_in = HG_mode(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1]) Hg_out = HG_mode(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])  Daniel Brown committed Mar 24, 2015 173   Daniel Brown committed Aug 12, 2014 174  U1 = Hg_in.Unm(x+delta[0], y+delta[1])  Daniel Brown committed Aug 01, 2014 175  U2 = Hg_out.Unm(x,y).conjugate()  Daniel Brown committed Mar 24, 2015 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201  if newtonCotesOrder > 0: W = newton_cotes(newtonCotesOrder, 1)[0] if newtonCotesOrder > 1: if (len(x) - len(W)) % newtonCotesOrder != 0: raise ValueError("To use Newton-Cotes order {0} the number of data points in x must ensure: (N_x - ({0}+1)) mod {0} == 0".format(newtonCotesOrder) ) if (len(y) - len(W)) % newtonCotesOrder != 0: raise ValueError("To use Newton-Cotes order {0} the number of data points in y must ensure: (N_y - ({0}+1)) mod {0} == 0".format(newtonCotesOrder) ) wx = np.zeros(x.shape, dtype=np.float64) wy = np.zeros(y.shape, dtype=np.float64) N = len(W) for i in range(0, (len(wx)-1)/newtonCotesOrder): wx[(i*(N-1)):(i*(N-1)+N)] += W for i in range(0, (len(wy)-1)/newtonCotesOrder): wy[(i*(N-1)):(i*(N-1)+N)] += W Wxy = np.outer(wx, wy) if newtonCotesOrder == 0: return dx * dy * np.einsum('ij,ij', Axy, U1*U2) else: return dx * dy * np.einsum('ij,ij', Axy, U1*U2*Wxy)  Daniel Brown committed Aug 01, 2014 202 203 204 205 206 207 208  else: strx = "u1[%i,%i]" % (mode_in[0], mode_out[0]) stry = "u2[%i,%i]" % (mode_in[1], mode_out[1]) return dx * dy * np.einsum('ij,ij', Axy, np.outer(cache[strx], cache[stry]))  Daniel Brown committed Jul 22, 2014 209   Daniel Brown committed Aug 01, 2014 210   Daniel Brown committed Aug 12, 2014 211   Daniel Brown committed Sep 04, 2014 212 def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=(0,0), params={}):  Daniel Brown committed Apr 20, 2015 213  if q1y is None:  Daniel Brown committed Aug 01, 2014 214 215  q1y = q1  Daniel Brown committed Apr 20, 2015 216  if q2y is None:  Daniel Brown committed Aug 01, 2014 217 218 219 220 221 222 223 224 225 226 227 228 229 230  q2y = q2 it = np.nditer(couplings, flags=['refs_ok','f_index']) cache = {} while not it.finished: try: mode_in = [int(it.next()), int(it.next())] mode_out = [int(it.next()), int(it.next())] strx = "u1[%i,%i]" % (mode_in[0], mode_out[0]) stry = "u2[%i,%i]" % (mode_in[1], mode_out[1])  Daniel Brown committed Aug 12, 2014 231 232  #Hg_in = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1]) #Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])  Daniel Brown committed Aug 01, 2014 233 234  if strx not in cache:  Daniel Brown committed Aug 12, 2014 235 236  cache[strx] = u_star_u(q1.z, q2.z, q1.w0, q2.w0, mode_in[0], mode_out[0], x, x+delta[0]) #Hg_in.Un(x) * Hg_out.Un(x).conjugate()  Daniel Brown committed Aug 01, 2014 237 238  if stry not in cache:  Daniel Brown committed Aug 12, 2014 239 240  cache[stry] = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, mode_in[1], mode_out[1], y, y+delta[1]) #Hg_in.Um(y) * Hg_out.Um(y).conjugate()  Daniel Brown committed Aug 01, 2014 241 242 243 244 245 246  except StopIteration: break return cache  Daniel Brown committed Aug 12, 2014 247 248   Daniel Brown committed Jul 25, 2014 249 250 def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None):  Daniel Brown committed Apr 20, 2015 251  if q1y is None:  Daniel Brown committed Jul 25, 2014 252 253  q1y = q1  Daniel Brown committed Apr 20, 2015 254  if q2y is None:  Daniel Brown committed Jul 25, 2014 255 256 257 258 259 260 261 262 263 264 265 266 267  q2y = q2 it = np.nditer(couplings, flags=['refs_ok','f_index']) cache = {} cache["w_ij_Q1Q3"] = weights.w_ij_Q1 + weights.w_ij_Q3 cache["w_ij_Q2Q4"] = weights.w_ij_Q2 + weights.w_ij_Q4 cache["w_ij_Q1Q2"] = weights.w_ij_Q1 + weights.w_ij_Q2 cache["w_ij_Q1Q4"] = weights.w_ij_Q1 + weights.w_ij_Q4 cache["w_ij_Q2Q3"] = weights.w_ij_Q2 + weights.w_ij_Q3 cache["w_ij_Q3Q4"] = weights.w_ij_Q3 + weights.w_ij_Q4 cache["w_ij_Q1Q2Q3Q4"] = weights.w_ij_Q1 + weights.w_ij_Q3 + weights.w_ij_Q2 + weights.w_ij_Q4  Daniel Brown committed Jul 22, 2014 268   Daniel Brown committed Jul 25, 2014 269 270 271 272 273 274 275 276 277  while not it.finished: try: mode_in = [int(it.next()), int(it.next())] mode_out = [int(it.next()), int(it.next())] strx = "x[%i,%i]" % (mode_in[0], mode_out[0]) stry = "y[%i,%i]" % (mode_in[1], mode_out[1]) if strx not in cache:  Daniel Brown committed Aug 01, 2014 278  cache[strx] = u_star_u(q1.z, q2.z, q1.w0, q2.w0, mode_in[0], mode_out[0], weights.EI["xm"].nodes)  Daniel Brown committed Jul 25, 2014 279   Daniel Brown committed Aug 01, 2014 280 281  if stry not in cache: cache[stry] = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, mode_in[1], mode_out[1], weights.EI["ym"].nodes)  Daniel Brown committed Jul 25, 2014 282 283 284 285 286  except StopIteration: break return cache  Daniel Brown committed Aug 12, 2014 287 288 289   Daniel Brown committed Jul 25, 2014 290 def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=None):  Daniel Brown committed Apr 20, 2015 291  if q1y is None:  Daniel Brown committed Jul 25, 2014 292  q1y = q1  Daniel Brown committed Jul 22, 2014 293   Daniel Brown committed Apr 20, 2015 294  if q2y is None:  Daniel Brown committed Jul 25, 2014 295 296 297 298 299 300 301 302 303 304  q2y = q2 # x modes n = mode_in[0] m = mode_out[0] # y modes npr = mode_in[1] mpr = mode_out[1]  305  if isinstance(weights, ROMWeights):  Daniel Brown committed Apr 20, 2015 306  if cache is None:  Daniel Brown committed Apr 10, 2015 307 308  u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EIx.nodes) u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EIy.nodes)  Daniel Brown committed Aug 01, 2014 309   310 311 312 313 314 315 316  w_ij_Q1Q3 = weights.w_ij_Q1 + weights.w_ij_Q3 w_ij_Q2Q4 = weights.w_ij_Q2 + weights.w_ij_Q4 w_ij_Q1Q2 = weights.w_ij_Q1 + weights.w_ij_Q2 w_ij_Q1Q4 = weights.w_ij_Q1 + weights.w_ij_Q4 w_ij_Q2Q3 = weights.w_ij_Q2 + weights.w_ij_Q3 w_ij_Q3Q4 = weights.w_ij_Q3 + weights.w_ij_Q4 w_ij_Q1Q2Q3Q4 = weights.w_ij_Q1 + weights.w_ij_Q2 + weights.w_ij_Q3 + weights.w_ij_Q4  Daniel Brown committed Aug 12, 2014 317   318 319 320  else: strx = "x[%i,%i]" % (mode_in[0], mode_out[0]) stry = "y[%i,%i]" % (mode_in[1], mode_out[1])  Daniel Brown committed Jul 31, 2014 321   322 323  u_x_nodes = cache[strx] u_y_nodes = cache[stry]  Daniel Brown committed Aug 12, 2014 324   325 326 327 328 329 330 331  w_ij_Q1Q3 = cache["w_ij_Q1Q3"] w_ij_Q2Q4 = cache["w_ij_Q2Q4"] w_ij_Q1Q2 = cache["w_ij_Q1Q2"] w_ij_Q1Q4 = cache["w_ij_Q1Q4"] w_ij_Q2Q3 = cache["w_ij_Q2Q3"] w_ij_Q3Q4 = cache["w_ij_Q3Q4"] w_ij_Q1Q2Q3Q4 = cache["w_ij_Q1Q2Q3Q4"]  Daniel Brown committed Aug 12, 2014 332   333  u_xy_nodes = np.outer(u_x_nodes, u_y_nodes)  Daniel Brown committed Aug 12, 2014 334   335 336 337 338  n_mod_2 = n % 2 m_mod_2 = m % 2 npr_mod_2 = npr % 2 mpr_mod_2 = mpr % 2  Daniel Brown committed Aug 12, 2014 339   340 341  if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2: k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)  Daniel Brown committed Aug 12, 2014 342   343 344 345 346 347  elif n_mod_2 != m_mod_2: if npr_mod_2 == mpr_mod_2: k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q3) else: k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q3)  Daniel Brown committed Jul 31, 2014 348   349 350 351 352 353 354 355  elif npr_mod_2 != mpr_mod_2: if n_mod_2 == m_mod_2: k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2) else: k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q3) else:  Daniel Brown committed Apr 20, 2015 356  if cache is None:  Daniel Brown committed Apr 10, 2015 357 358  u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EIx.nodes) u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EIy.nodes)  359 360  w_ij = weights.w_ij  Daniel Brown committed Aug 12, 2014 361  else:  362 363 364 365 366  strx = "x[%i,%i]" % (mode_in[0], mode_out[0]) stry = "y[%i,%i]" % (mode_in[1], mode_out[1]) u_x_nodes = cache[strx] u_y_nodes = cache[stry]  Daniel Brown committed Jul 30, 2014 367   368 369 370 371  u_xy_nodes = np.outer(u_x_nodes, u_y_nodes) k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij)  Daniel Brown committed Aug 12, 2014 372  return k_ROQ  Daniel Brown committed Jul 31, 2014 373   Daniel Brown committed Aug 12, 2014 374 __fac_cache = []  Daniel Brown committed Jul 31, 2014 375   Daniel Brown committed Aug 12, 2014 376 377 def fac(n): global __fac_cache  378 379 380 381  if len(__fac_cache) == 0: return math.factorial(int(n)) else: return __fac_cache[n]  Daniel Brown committed Jul 31, 2014 382   Daniel Brown committed Aug 12, 2014 383 384 385 def m_1_pow(n): if n % 2 == 0: return 1  Daniel Brown committed Jul 31, 2014 386  else:  Daniel Brown committed Aug 12, 2014 387 388 389 390 391 392 393 394  return -1 def __Ss(u, _u, F, _F, d=0): r = 0 for s in range(0, min(u,_u)+1): r += m_1_pow(s) * _F ** (u-s) * _F ** (_u-s) / (fac(2*s+d)*fac(u-s)*fac(_u-s))  Daniel Brown committed Jul 31, 2014 395   Daniel Brown committed Aug 12, 2014 396 397 398 399 400  return r def __S(m, _m, X, _X, F, _F, d=0): if m % 2 == 1:  Daniel Brown committed May 20, 2015 401  lim1 = int((m-1)/2)  Daniel Brown committed Aug 12, 2014 402  else:  Daniel Brown committed May 20, 2015 403  lim1 = int(m/2 )  Daniel Brown committed Aug 12, 2014 404 405  if _m % 2 == 1:  Daniel Brown committed May 20, 2015 406  lim2 = int((_m-1)/2)  Daniel Brown committed Aug 12, 2014 407  else:  Daniel Brown committed May 20, 2015 408  lim2 = int(_m/2)  Daniel Brown committed Aug 12, 2014 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436  r = 0 for u in range(0, lim1+1): for _u in range(0, lim2+1): r += m_1_pow(u) * _X**(m-2*u) * X**(_m-2*_u) / ( fac(m-2*u)*fac(_m-2*_u) ) * __Ss(u, _u, F, _F, d=d) return r def __bayerhelms_kn(n, _n, q1, q2, gamma=0.0): K0 = (q1.zr - q2.zr)/q2.zr K2 = (q1.z - q2.z)/q2.zr K = (K0 + 1j*K2)/2.0 Ktilde = abs(K / (1+K)) if gamma != 0: a = q2.zr * math.sin(gamma) / (cmath.sqrt(1+K.conjugate()) * q2.w0) _X = - a * (q2.z/q2.zr - 1j) X = - a * (q2.z/q2.zr + 1j*(1+2*K.conjugate())) Ex = cmath.exp(-_X*X / 2.0) else: _X = 0.0 X = 0.0 Ex = 1.0  Daniel Brown committed Jul 25, 2014 437   Daniel Brown committed Aug 12, 2014 438 439  _F = K / (2.0 * (1.0+K0)) F = K.conjugate() / 2.0  Daniel Brown committed Jul 25, 2014 440   Daniel Brown committed Aug 12, 2014 441  Sg = __S(n, _n, X, _X, F, _F)  Daniel Brown committed Jul 25, 2014 442   Daniel Brown committed Aug 12, 2014 443 444 445 446 447 448 449 450 451 452 453  if n > 0 and _n > 0: Su = __S(n-1, _n-1, X, _X, F, _F, 1) else: Su = 0 b = m_1_pow(_n) * cmath.sqrt(fac(n) * fac(_n) * (1.0 + K.real)**(n+0.5) * (1.0 + K.conjugate()) ** (-(n+_n+1))) return b * Ex * (Sg - Su) def bayerhelms_HG_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, gamma=(0,0)):  Daniel Brown committed Apr 20, 2015 454  if q1y is None:  Daniel Brown committed Aug 12, 2014 455 456  q1y = q1  Daniel Brown committed Apr 20, 2015 457  if q2y is None:  Daniel Brown committed Aug 12, 2014 458 459 460 461 462 463 464 465 466 467 468 469  q2y = q2 # x modes n = mode_in[0] _n = mode_out[0] # y modes m = mode_in[1] _m = mode_out[1] return __bayerhelms_kn(n,_n, q1, q2, 2*gamma[0]) * __bayerhelms_kn(m, _m, q1y, q2y, 2*gamma[1])  Daniel Brown committed Mar 24, 2015 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 def __sqaure_knm_int(n, _n, R): # This uses the H_n(x) * H_m(x) product identity to reduce the overlap into # a sum of factorial and an integral of a single Hermite with a gaussian function # thus making it easier to solve expR = math.exp(-(R**2)) S = 0 for j in range(0, min(n, _n)+1): _j1 = _n + n - 2*j - 1 if _j1+1 == 0: # for the zeroth order we just have the gaussian integral to solve L = math.sqrt(math.pi) * math.erf(R) elif (_j1+1) % 2 == 1: # if the Hermite is odd then the integral is always 0 as its anti-symmetric L = 0 else: L = 2 * hermite(_j1, 0) - expR * (hermite(_j1, R) - hermite(_j1, -R)) I = 2**j * factorial(j) * comb(n, j) * comb(_n, j) S += I * L return S def square_aperture_HG_knm(mode_in, mode_out, q, R): """ Computes the coupling coefficients for a square aperture. """ # x modes n = mode_in[0] _n = mode_out[0] # y modes m = mode_in[1] _m = mode_out[1]  Andreas Freise committed May 14, 2016 508 509  hg1 = HG_mode(q, n=n, m=m) hg2 = HG_mode(q, n=_n, m=_m)  Daniel Brown committed Mar 24, 2015 510 511 512  kx = hg1.constant_x * hg2.constant_x.conjugate() ky = hg1.constant_y * hg2.constant_y.conjugate()  Daniel Brown committed Apr 10, 2015 513   Daniel Brown committed Mar 24, 2015 514 515 516 517 518 519 520 521 522 523 524  f = q.w / math.sqrt(2) R = R / (q.w / math.sqrt(2)) kx *= f kx *= __sqaure_knm_int(n, _n, R) ky *= f ky *= __sqaure_knm_int(m, _m, R) return kx * ky  Daniel Brown committed Sep 04, 2014 525 526 527  def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=False, profile=False, gamma=(0,0), delta=(0,0), params={}):  Daniel Brown committed Apr 20, 2015 528  if q1y is None:  Daniel Brown committed Jul 25, 2014 529 530  q1y = q1  Daniel Brown committed Apr 20, 2015 531  if q2y is None:  Daniel Brown committed Jul 25, 2014 532 533 534  q2y = q2 assert q1.wavelength == q2.wavelength and q1y.wavelength == q2y.wavelength and q1y.wavelength == q1.wavelength  Daniel Brown committed Jul 22, 2014 535 536 537 538 539 540 541 542  couplings = np.array(couplings) a = couplings.size / 4.0 if int(a) - a != 0: raise BasePyKatException("Iterator should be product of 4, each element of coupling array should be [n,m,n',m']")  Daniel Brown committed Aug 12, 2014 543 544 545 546 547 548 549 550 551 552 553  maxtem = 0 c = couplings.flatten() for i in range(0, c.size/2): maxtem = max(sum(c[i*2:(i*2+2)]), maxtem) global __fac_cache for n in range(0, maxtem+1): __fac_cache.append(math.factorial(n))  Daniel Brown committed Apr 20, 2015 554  if surface_map is not None:  Daniel Brown committed Sep 04, 2014 555  Axy = surface_map.z_xy(wavelength=q1.wavelength)  Daniel Brown committed Aug 12, 2014 556 557 558  x = surface_map.x y = surface_map.y  Daniel Brown committed Jul 22, 2014 559 560 561 562 563 564  K = np.zeros((couplings.size/4,), dtype=np.complex128) it = np.nditer(couplings, flags=['refs_ok','f_index']) i = 0  Daniel Brown committed Aug 12, 2014 565 566 567  if profile: t0 = time.time()  Daniel Brown committed Aug 01, 2014 568   Daniel Brown committed Jul 25, 2014 569  if method == "romhom":  Daniel Brown committed Apr 20, 2015 570  if surface_map is None:  Daniel Brown committed Jul 25, 2014 571 572 573 574  raise BasePyKatException("Using 'romhom' method requires a surface map to be specified") weights = surface_map.ROMWeights  Daniel Brown committed Apr 20, 2015 575  if weights is None:  Daniel Brown committed Jul 25, 2014 576  raise BasePyKatException("The ROM weights need to be generated for this map before use.")  Daniel Brown committed Aug 01, 2014 577 578  cache = __gen_ROM_HG_knm_cache(weights, couplings, q1=q1, q2=q2, q1y=q1y, q2y=q2y)  Daniel Brown committed Jul 25, 2014 579   Daniel Brown committed Aug 01, 2014 580  elif method == "riemann":  Daniel Brown committed Apr 20, 2015 581  if surface_map is None:  Daniel Brown committed Aug 12, 2014 582 583 584  raise BasePyKatException("Using 'riemann' method requires a surface map to be specified") cache = __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=delta)  Daniel Brown committed Jul 25, 2014 585  else:  Daniel Brown committed Aug 01, 2014 586  cache = None  Daniel Brown committed Jul 25, 2014 587  weights = None  Daniel Brown committed Aug 01, 2014 588   Daniel Brown committed Aug 12, 2014 589 590 591 592  if profile: cache_time = time.time() - t0 Ktime = np.zeros((couplings.size/4,), dtype=np.float64)  Daniel Brown committed Aug 01, 2014 593 594 595  if verbose: p = ProgressBar(maxval=couplings.size, widgets=["Knm (%s): " % method, Percentage(), Bar(), ETA()])  Daniel Brown committed Jul 22, 2014 596 597  while not it.finished: try:  Daniel Brown committed Aug 12, 2014 598 599 600  if profile: t0 = time.time()  Daniel Brown committed Jul 22, 2014 601 602  mode_in = [int(it.next()), int(it.next())] mode_out = [int(it.next()), int(it.next())]  Daniel Brown committed Jul 25, 2014 603   Daniel Brown committed Sep 04, 2014 604   Daniel Brown committed Jul 22, 2014 605  if method == "riemann":  Daniel Brown committed Aug 12, 2014 606  K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy, cache=cache, delta=delta)  Daniel Brown committed Jul 25, 2014 607  elif method == "romhom":  Daniel Brown committed Aug 01, 2014 608  K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=cache)  Daniel Brown committed Aug 12, 2014 609 610  elif method == "bayerhelms": K[i] = bayerhelms_HG_knm(mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, gamma=gamma)  Daniel Brown committed Sep 04, 2014 611 612  elif method == "adaptive": K[i] = adaptive_knm(mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, smap=surface_map, delta=delta, params=params)  Daniel Brown committed Jul 22, 2014 613  else:  Daniel Brown committed Jul 25, 2014 614  raise BasePyKatException("method value '%s' not accepted" % method)  Daniel Brown committed Aug 12, 2014 615 616 617 618  if profile: Ktime[i] = time.time() - t0  Daniel Brown committed Jul 22, 2014 619  i +=1  Daniel Brown committed Aug 01, 2014 620 621 622 623  if verbose: p.update(i*4)  Daniel Brown committed Jul 22, 2014 624 625 626 627  except StopIteration: break  Daniel Brown committed Aug 12, 2014 628 629 630 631  if profile: return K.reshape(couplings.shape[:-1]), Ktime.reshape(couplings.shape[:-1]), cache_time else: return K.reshape(couplings.shape[:-1])  Daniel Brown committed Aug 01, 2014 632 633   Daniel Brown committed Mar 24, 2015 634 635 636   Daniel Brown committed Aug 12, 2014 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 def plot_knm_matrix(couplings, knm): import pylab as plt fig = plt.figure() ax = fig.add_subplot(111) cax = ax.imshow(knm, interpolation='nearest') fig.colorbar(cax) numrows, numcols = knm.shape c = couplings[:, 0, :2] c_ = [] for d in c: c_.append("%i,%i"%(d[0], d[1])) ax.set_xticklabels(c_) ax.set_yticklabels(c_) def format_coord(x, y): col = int(x+0.5) row = int(y+0.5) if col>=0 and col=0 and row