#!/usr/bin/env python # -*- coding: utf-8 -*- """ Animate the result of oblique incidence on slab Daniel Sjoberg, 2011-09-03 """ import matplotlib matplotlib.use('WxAgg') matplotlib.interactive(True) from pylab import * from numpy.lib.scimath import sqrt # Use complex sqrt by default from mayavi.mlab import * from tvtk.tools import visual na = 1.5; epsa = na**2 nb = 1.0; epsb = nb**2 d = 5 k0 = 2 kx = k0*nb*1.01 kza = -1j*sqrt(-k0**2*na**2 + kx**2) kzb = -1j*sqrt(-k0**2*nb**2 + kx**2) rhoaTE = (kza - kzb)/(kza + kzb) rhoaTM = (kzb*epsa - kza*epsb)/(kzb*epsa + kza*epsb) rhoa = rhoaTE P = exp(-1j*kzb*d) Gamma = rhoa*(1 - P**2)/(1 - rhoa**2*P**2) T = (1 - rhoa**2)*P/(1 - rhoa**2*P**2) x0 = 10 dx = 0.1 z0 = 10 - d/2 dz = 0.1 [x1, z1] = mgrid[-x0:x0:dx, -z0:0+dz:dz] [x2, z2] = mgrid[-x0:x0:dx, 0:d+dz:dz] [x3, z3] = mgrid[-x0:x0:dx, d:z0+d+dz:dz] wave1 = exp(-1j*(kza*z1 + kx*x1)) + Gamma*exp(-1j*(-kza*z1 + kx*x1)) wave2 = (1 + rhoa)/(1 - rhoa**2*P**2)*(exp(-1j*kzb*z2) - rhoa*P**2*exp(1j*kzb*z2))*exp(-1j*kx*x2) wave3 = T*exp(-1j*(kza*(z3-d) + kx*x3)) scale = amax(abs(wave2)) #fig = figure(size=(500,500)) fig = figure(size=(800,800)) visual.set_viewer(fig) s1 = surf(x1, z1, real(wave1)/scale, vmin=-1, vmax=1) s2 = surf(x2, z2, real(wave2)/scale, vmin=-1, vmax=1); outline() s3 = surf(x3, z3, real(wave3)/scale, vmin=-1, vmax=1) @show @animate(delay=20) def anim(): while True: for phi in linspace(0, 2*pi, 20): fig.scene.disable_render = True s1.mlab_source.scalars = real(wave1*exp(1j*phi))/scale s2.mlab_source.scalars = real(wave2*exp(1j*phi))/scale s3.mlab_source.scalars = real(wave3*exp(1j*phi))/scale fig.scene.disable_render = False yield anim()