knm.py 20.6 KB
 Daniel Brown committed Jul 22, 2014 1 from itertools import combinations_with_replacement as combinations  Andreas Freise committed Dec 20, 2014 2 from pykat.optics.gaussian_beams import beam_param, HG_beam  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 10 11 12 from math import factorial from pykat.maths.hermite import hermite from scipy.misc import comb from scipy.integrate import newton_cotes  Daniel Brown committed Jul 22, 2014 13   Daniel Brown committed Aug 12, 2014 14 import time  Andreas Freise committed Dec 20, 2014 15 import pykat.optics.maps  Daniel Brown committed Jul 22, 2014 16 17 18 19 20 import os.path import numpy as np import pykat import collections import math  Daniel Brown committed Aug 12, 2014 21 import cmath  Daniel Brown committed Aug 01, 2014 22   Daniel Brown committed Sep 04, 2014 23 def makeCouplingMatrix(max_order, Neven=True, Nodd=True, Meven=True, Modd=True):  Daniel Brown committed Aug 01, 2014 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38  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 39 40 41 42 43 44  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 45 46  row.append(e)  Daniel Brown committed Sep 04, 2014 47 48  M.append(np.array(row).squeeze())  Daniel Brown committed Aug 01, 2014 49   Daniel Brown committed Aug 12, 2014 50  return np.array(M)  Daniel Brown committed Jul 22, 2014 51   Daniel Brown committed Sep 04, 2014 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 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 def adaptive_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, smap=None, delta=(0,0), params={}): if q1y == None: q1y = q1 if q2y == None: 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]") 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]) Nfuncs = [] Nfuncs.append(0) if smap != None: 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 151 152 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 153   Daniel Brown committed Mar 24, 2015 154  if Axy is None:  Daniel Brown committed Jul 22, 2014 155 156  Axy == np.ones((len(x), len(y)))  Daniel Brown committed Mar 24, 2015 157  if q1y is None:  Daniel Brown committed Jul 22, 2014 158 159  q1y = q1  Daniel Brown committed Mar 24, 2015 160  if q2y is None:  Daniel Brown committed Jul 22, 2014 161 162 163 164  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 165   Daniel Brown committed Jul 22, 2014 166  dx = abs(x[1] - x[0])  Daniel Brown committed Aug 01, 2014 167 168 169 170 171  dy = abs(y[1] - y[0]) if cache == None: 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 Mar 24, 2015 172   Daniel Brown committed Aug 12, 2014 173  U1 = Hg_in.Unm(x+delta[0], y+delta[1])  Daniel Brown committed Aug 01, 2014 174  U2 = Hg_out.Unm(x,y).conjugate()  Daniel Brown committed Mar 24, 2015 175 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  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 201 202 203 204 205 206 207  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 208   Daniel Brown committed Aug 01, 2014 209   Daniel Brown committed Aug 12, 2014 210   Daniel Brown committed Sep 04, 2014 211 def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=(0,0), params={}):  Daniel Brown committed Aug 01, 2014 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229  if q1y == None: q1y = q1 if q2y == None: 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 230 231  #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 232 233  if strx not in cache:  Daniel Brown committed Aug 12, 2014 234 235  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 236 237  if stry not in cache:  Daniel Brown committed Aug 12, 2014 238 239  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 240 241 242 243 244 245  except StopIteration: break return cache  Daniel Brown committed Aug 12, 2014 246 247   Daniel Brown committed Jul 25, 2014 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None): if q1y == None: q1y = q1 if q2y == None: 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 267   Daniel Brown committed Jul 25, 2014 268 269 270 271 272 273 274 275 276  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 277  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 278   Daniel Brown committed Aug 01, 2014 279 280  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 281 282 283 284 285  except StopIteration: break return cache  Daniel Brown committed Aug 12, 2014 286 287 288   Daniel Brown committed Jul 25, 2014 289 290 291 def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=None): if q1y == None: q1y = q1  Daniel Brown committed Jul 22, 2014 292   Daniel Brown committed Jul 25, 2014 293 294 295 296 297 298 299 300 301 302 303  if q2y == None: q2y = q2 # x modes n = mode_in[0] m = mode_out[0] # y modes npr = mode_in[1] mpr = mode_out[1]  304 305 306 307  if isinstance(weights, ROMWeights): if cache == None: u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["x"].nodes) u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["y"].nodes)  Daniel Brown committed Aug 01, 2014 308   309 310 311 312 313 314 315  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 316   317 318 319  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 320   321 322  u_x_nodes = cache[strx] u_y_nodes = cache[stry]  Daniel Brown committed Aug 12, 2014 323   324 325 326 327 328 329 330  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 331 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 356 357 358 359 360  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: if cache == None: u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["x"].nodes) u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["y"].nodes) 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 401 402 403 404 405 406 407 408 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  return r def __S(m, _m, X, _X, F, _F, d=0): if m % 2 == 1: lim1 = (m-1)/2 else: lim1 = m/2 if _m % 2 == 1: lim2 = (_m-1)/2 else: lim2 = _m/2 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 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469  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)): if q1y == None: q1y = q1 if q2y == None: 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 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 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] hg1 = HG_beam(q, n=n, m=m) hg2 = HG_beam(q, n=_n, m=_m) kx = hg1.constant_x * hg2.constant_x.conjugate() ky = hg1.constant_y * hg2.constant_y.conjugate() print hg1.constant_x, hg2.constant_x, hg1.constant_y, hg2.constant_y 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 527 528 529  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 Jul 25, 2014 530 531 532 533 534 535 536  if q1y == None: q1y = q1 if q2y == None: q2y = q2 assert q1.wavelength == q2.wavelength and q1y.wavelength == q2y.wavelength and q1y.wavelength == q1.wavelength  Daniel Brown committed Jul 22, 2014 537 538 539 540 541 542 543 544  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 545 546 547 548 549 550 551 552 553 554 555 556  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)) if surface_map != None:  Daniel Brown committed Sep 04, 2014 557  Axy = surface_map.z_xy(wavelength=q1.wavelength)  Daniel Brown committed Aug 12, 2014 558 559 560  x = surface_map.x y = surface_map.y  Daniel Brown committed Jul 22, 2014 561 562 563 564 565 566  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 567 568 569  if profile: t0 = time.time()  Daniel Brown committed Aug 01, 2014 570   Daniel Brown committed Jul 25, 2014 571 572 573 574 575 576 577 578  if method == "romhom": if surface_map == None: raise BasePyKatException("Using 'romhom' method requires a surface map to be specified") weights = surface_map.ROMWeights if weights == None: raise BasePyKatException("The ROM weights need to be generated for this map before use.")  Daniel Brown committed Aug 01, 2014 579 580  cache = __gen_ROM_HG_knm_cache(weights, couplings, q1=q1, q2=q2, q1y=q1y, q2y=q2y)  Daniel Brown committed Jul 25, 2014 581   Daniel Brown committed Aug 01, 2014 582  elif method == "riemann":  Daniel Brown committed Aug 12, 2014 583 584 585 586  if surface_map == None: 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 587  else:  Daniel Brown committed Aug 01, 2014 588  cache = None  Daniel Brown committed Jul 25, 2014 589  weights = None  Daniel Brown committed Aug 01, 2014 590   Daniel Brown committed Aug 12, 2014 591 592 593 594  if profile: cache_time = time.time() - t0 Ktime = np.zeros((couplings.size/4,), dtype=np.float64)  Daniel Brown committed Aug 01, 2014 595 596 597  if verbose: p = ProgressBar(maxval=couplings.size, widgets=["Knm (%s): " % method, Percentage(), Bar(), ETA()])  Daniel Brown committed Jul 22, 2014 598 599  while not it.finished: try:  Daniel Brown committed Aug 12, 2014 600 601 602  if profile: t0 = time.time()  Daniel Brown committed Jul 22, 2014 603 604  mode_in = [int(it.next()), int(it.next())] mode_out = [int(it.next()), int(it.next())]  Daniel Brown committed Jul 25, 2014 605   Daniel Brown committed Sep 04, 2014 606   Daniel Brown committed Jul 22, 2014 607  if method == "riemann":  Daniel Brown committed Aug 12, 2014 608  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 609  elif method == "romhom":  Daniel Brown committed Aug 01, 2014 610  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 611 612  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 613 614  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 615  else:  Daniel Brown committed Jul 25, 2014 616  raise BasePyKatException("method value '%s' not accepted" % method)  Daniel Brown committed Aug 12, 2014 617 618 619 620  if profile: Ktime[i] = time.time() - t0  Daniel Brown committed Jul 22, 2014 621  i +=1  Daniel Brown committed Aug 01, 2014 622 623 624 625  if verbose: p.update(i*4)  Daniel Brown committed Jul 22, 2014 626 627 628 629  except StopIteration: break  Daniel Brown committed Aug 12, 2014 630 631 632 633  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 634 635   Daniel Brown committed Mar 24, 2015 636 637 638   Daniel Brown committed Aug 12, 2014 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 665 666 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