"""
(C) Steven Byrnes, 2013. This code is released under the MIT license
http://opensource.org/licenses/MIT

This code runs in Python 2.7 or 3.3. It requires imagemagick to be installed;
that's how it assembles images into animated GIFs.

Creates an animation of the electric field vectors of a
surface-plasmon-polariton wave. The metal-dielectric interface is at z=0,
with metal at z<0 and dielectric (such as air) at z>0.
"""
 
# Use Python 3 style division and print
# a/b is real division, a//b is integer division
from __future__ import division, print_function

import numpy as np
from cmath import exp, pi

import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches

import subprocess, os
directory_now = os.path.dirname(os.path.realpath(__file__))

# --- PART 1 OF 2: PHYSICS --- #

def Efield(x_times_kvac, z_times_kvac, t_times_omega, eps_m, eps_d, print_wave_properties=False):
    """"
    Inputs:
      * x_times_kvac and z_times_kvac are the coordinates (x,z) multiplied
          by the vacuum angular wavenumber kvac = omega / c.
      * t_times_nu is the time t multiplied by frequency nu.
      * eps_m and eps_d are the complex permittivities of the metal and dielectric.
    
    Output: The electric field vector (E_x,E_z). It is scaled so that
      E_z(0, 0, tw, e1, e2) = cos(tw) on the dielectric side of the interface.
    """
    # Calculate components of the angular wavevector in the dielectric and
    # metal, as a multiple of kvac = omega / c:
    eps_m = complex(eps_m) #cast to complex, so square roots won't raise errors.
    kx_over_kvac = (eps_m * eps_d / (eps_m + eps_d))**(1/2)
    kzd_over_kvac = (eps_d - kx_over_kvac**2)**(1/2)
    kzm_over_kvac = (eps_m - kx_over_kvac**2)**(1/2)
    
    # Pick the correct square-roots, so that e^(i*k*z) decays away from interface
    if kzd_over_kvac.imag < 0:
        kzd_over_kvac *= -1
    if kzm_over_kvac.imag > 0:
        kzm_over_kvac *= -1
    
    #double-check the boundary condition
    almost_equal = lambda a,b,tolerance : (abs(a-b) / (abs(a) + abs(b))) < tolerance
    if not almost_equal(kzd_over_kvac * eps_m, kzm_over_kvac * eps_d, 1e-10):
        raise ValueError('Something is wrong! Boundary condition fails!')
    
    if print_wave_properties:
        print('kx / kvac = ', kx_over_kvac)
        print('kzd / kvac = ', kzd_over_kvac)
        print('kzm / kvac = ', kzm_over_kvac)
        print('Wavelength / (Vacuum wavelength) = ', 1/(kx_over_kvac.real))
        if kx_over_kvac.imag != 0:
            print('(Decay length) / (Vacuum wavelength) = ', 1/(kx_over_kvac.imag))
        else:
            print('Wave does not decay, it propagates forever.')
        print('(Decay length into dielectric) / (Vacuum wavelength) = ', 1/(kzd_over_kvac.imag))
        print('(Decay length into metal) / (Vacuum wavelength) = ', -1/(kzm_over_kvac.imag))
   
    if z_times_kvac > 0:
        # dielectric
        Ez = exp(1j * kx_over_kvac * x_times_kvac + 1j * kzd_over_kvac * z_times_kvac - 1j * t_times_omega)
        Ex = -Ez * kzd_over_kvac / kx_over_kvac
    else:
        # metal
        Ez = (kzd_over_kvac / kzm_over_kvac) * exp(
                 1j * kx_over_kvac * x_times_kvac + 1j * kzm_over_kvac * z_times_kvac - 1j * t_times_omega)
        Ex = -Ez * kzm_over_kvac / kx_over_kvac
    return (Ex.real, Ez.real)

# --- PART 2 OF 2: DRAWING --- #

def draw_box(x1, x2, y1, y2, color):
    """
    Code to draw a rectangular box in matplotlib ... used to color the metal area.
    See http://matplotlib.org/users/path_tutorial.html
    """
    vertices = [ (x1, y1), (x1, y2), (x2, y2), (x2, y1), (x1, y1)]
    codes = [Path.MOVETO,Path.LINETO,Path.LINETO,Path.LINETO,Path.CLOSEPOLY]
    path = Path(vertices, codes)
    #lw=0 means no outline; zorder=-1 means it shouldn't block the arrows.
    return patches.PathPatch(path, facecolor=color, lw=0, zorder=-1)

def draw_frame(t_times_omega, eps_m, eps_d, aspect_ratio=1.5,
               frac_metal=0.5, fig_width_px=480, x_range_times_kvac=2,
               img_filename=None):
    """
    Draw one frame of the animation.
    
    Inputs:
      * t_times_omega is time multiplied by angular frequency,
          it goes 0 --> 2pi each cycle.
      * eps_m and eps_d are the dielectric constants of the metal and dielectric
      * aspect_ratio is width over height
      * frac_metal is how much of the image is taken up by the metal.
      * fig_width_px is figure width in pixels
      * x_range_times_kvac is the width of the image as a multiple of kvac.
      * "img_filename" is what to save the frame as (or None to not save it).
    """
    # Figure geometry...
    fig_height_px = fig_width_px // aspect_ratio
    dpi = 80 #This number doesn't affect the final animation...
    fig_width_inches = fig_width_px / dpi
    fig_height_inches = fig_height_px / dpi
    
    # Coordinate limits in figure. All are implicitly multiplied by kvac.
    z_range_times_kvac = x_range_times_kvac / aspect_ratio
    xmin = 0
    xmax = x_range_times_kvac
    zmin = -z_range_times_kvac * frac_metal
    zmax = z_range_times_kvac * (1-frac_metal)
    
    # How many arrows to draw?
    num_arrows_x = 15
    num_arrows_z = num_arrows_x // aspect_ratio
    
    # Pick arrow coordinates...
    arrow_x_list, spacing = np.linspace(xmin, xmax, num=num_arrows_x,
                                     endpoint=False, retstep=True)
    arrow_x_list += spacing / 2
    
    arrow_z_list, spacing = np.linspace(zmin, zmax, num=num_arrows_z,
                                     endpoint=False, retstep=True)
    arrow_z_list += spacing / 2

    X,Z = np.meshgrid(arrow_x_list, arrow_z_list)
    
    # Arrow length scale: Larger number = smaller arrows
    arrow_len_scale = 15
    
    # Calculate the length of each arrow
    Ex_func = np.vectorize(lambda x,z : Efield(x,z,t_times_omega,eps_m,eps_d)[0])
    Ex_array = Ex_func(X, Z)
    Ez_func = np.vectorize(lambda x,z : Efield(x,z,t_times_omega,eps_m,eps_d)[1])
    Ez_array = Ez_func(X, Z)
    
    # Open a new figure with correct aspect ratio and pixel count and white background
    fig = plt.figure(figsize = (fig_width_inches,fig_height_inches), dpi=dpi, facecolor='w')
    
    # Draw a new set of axes that fill the entire figure area
    ax=fig.add_axes((0,0,1,1),axisbg='w')
    ax.set_axis_off()
    
    # Color the metal part
    metal_color = '#dddddd' #light gray
    ax.add_patch(draw_box(xmin,xmax,zmin,0,metal_color))
    
    # Draw the arrows
    ax.quiver(X, Z, Ex_array , Ez_array , scale=arrow_len_scale, scale_units='width', pivot='middle')
    
    ax.set_xlim(xmin, xmax)
    ax.set_ylim(zmin, zmax)
    if img_filename is not None:
        fig.savefig(os.path.join(directory_now, img_filename), dpi=dpi)
 
def draw_anim(eps_m, eps_d, anim_filename='anim.gif', frames_in_anim=30,
              total_anim_time_in_sec=2, keep_frame_images=False, **kwargs):
    """
    Create an animated GIF. **kwargs are all the keyword arguments to
    draw_frame()
    
    keep_frame_images=True to save the individual frame image files that make
    up the animation; otherwise they are created and immediately deleted.
    """
    filename_list = ['temp' + str(n) + '.png' for n in range(frames_in_anim)]
    
    for n in range(frames_in_anim):
        draw_frame(2*pi*n/frames_in_anim, eps_m, eps_d,
                   img_filename=filename_list[n], **kwargs)
    
    seconds_per_frame = total_anim_time_in_sec / frames_in_anim
    frame_delay = str(seconds_per_frame * 100)
    command_list = ['convert', '-delay', frame_delay, '-loop', '0'] + filename_list + [anim_filename]
    # Use the "convert" command (part of ImageMagick) to build the animation
    subprocess.call(command_list, cwd=directory_now)
    if keep_frame_images is False:
        for filename in filename_list:
            os.remove(os.path.join(directory_now, filename))

###################################################################

if True:
    # Silver-air interface at 9.919um = 0.125eV, using data from Palik (p357)
    eps_m = (13.11+53.7j)**2
    eps_d = 1
    # Print diagnostics
    Efield(0, 0, 0, eps_m, eps_d, print_wave_properties=True)
    # Create animation
    draw_anim(eps_m, eps_d, anim_filename='plas_Silver_10um_Palik.gif',
              frac_metal=.2, x_range_times_kvac=4)