#!/usr/bin/env python
# coding: utf-8

# ## Vortex and antivortex in the Kosterlitz–Thouless phase transition<br> Nobel Prize in Physics, 2016 ##

# The basic structures in the Kosterlitz-Thouless phase transition are the vortex 
# and antivortex spin configurations.
# 
# These configurations are  represented by  unit vector fields exhibiting an isolated equilibrium point.

# In the following we define and plot such vector fields, and moreover animate the rotation of spins.

# A vortex is represented as the unit vector field at the points of a lattice, associated to the Hamiltonian vector field, V, of Hamilton function 
# $H(x,y)=\pm(x^2/2+y^2/2)$, i.e.  $V=(-\partial H/\partial y, \partial H/\partial x)$.
# 
# An antivortex is defined similarly, by any of the two  Hamilton functions $H(x,y)=\pm(x^2/2-y^2/2)$.

# In[1]:


import numpy as np

from plotly.offline import download_plotlyjs, init_notebook_mode,  iplot
init_notebook_mode(connected=True)


# Define a  function that computes the position and direction of arrows  representing spins located at the vertices of a  lattice. The vectors are rotated  about their starting points:

# In[2]:


def rot_matrix(theta):
    return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])

def get_arrows(start, end, arrow_angle, headsize=0.8):#function for drawing quivers
    #start is the numpy array of x,y-coordinates of the starting points of arrows; shape (2, m, n)
    #start[0,:,:] contains x-coordinates
    #end -numpy array with the same shape as start; contains the x and y-coords of ending points of the arrow  
    #the arrowhead is an isosceles triangle with the equal sides forming an angle of 2*arrow_angle radians 
    #headsize of 0.8, 0.75 is the most appropriate choice; 
    
    start=np.asarray(start)
    end=np.asarray(end)
    arr_dir=headsize*(start-end)/2.5 #
    v=np.einsum('ji, imn -> jmn', rot_matrix(arrow_angle), arr_dir)#Einstein summation to apply the rotation to 
                                                                   #all vectors at a time
    w=np.einsum('ji, imn -> jmn', rot_matrix(-arrow_angle), arr_dir)
    p=end+v
    q=end+w
    m=start[0,:,:].shape[0]
    n=start[0,:,:].shape[1]
    suppx=np.stack((start[0,:,:], end[0,:,:], np.nan*np.ones((m,n)) ))#supp is the line segment as support
                                                                      #for arrow; suppx, suppy are x,y coords
    suppy=np.stack((start[1,:,:], end[1,:,:], np.nan*np.ones((m,n)) ))

    x=suppx.flatten('F')#Fortran type flattening
    y=suppy.flatten('F')
    
    #headx, heady the x, y coordinates of triangle vertices
    headx=np.stack((end[0,:,:], p[0,:,:], q[0,:,:], end[0,:,:], np.nan*np.ones((m,n))))
    heady=np.stack((end[1,:,:], p[1,:,:], q[1,:,:], end[1,:,:], np.nan*np.ones((m,n))))
    xx=headx.flatten('F')
    yy=heady.flatten('F')
    x=map(lambda u: None if np.isnan(u) else u, x)
    y=map(lambda u: None if np.isnan(u) else u, y)
    xx=map(lambda u: None if np.isnan(u) else u, xx)
    yy=map(lambda u: None if np.isnan(u) else u, yy)
    return x,y, xx, yy
    


# Define data and frames for a vortex/antivortex animation:

# In[3]:


def anim_setup(start, end,  theta, traces):


    xs, ys, xx, yy = get_arrows(start, end, theta)
    
    data=[dict(x=xs, 
               y=ys,                        
               mode='lines', 
               line=dict(color='rgb(10,10,10)', width=1),
               name='',
               type='scatter',
              ),
          dict(x=xx, 
               y=yy,
               mode='lines', 
               line=dict(color='red', width=1), 
               name='', 
               type='scatter',
              )      
          ]

    frames=[]
    fin=2*np.pi + np.pi/18
    arrow_dir=end-start
    for phi in np.arange(0, fin, np.pi/18):
        v=np.einsum('ji, imn -> jmn', rot_matrix(phi), arrow_dir)
        end=start+v
    
        X, Y, XX, YY=get_arrows(start, end, theta)
        frames.append(dict(data=[dict(x=X, y=Y),
                                 dict(x=XX, y=YY) 
                                ],
                           traces=traces
                         )
                  )
    return data, frames


# In[4]:


def set_layout(title='',  width=900, height=500, xpos=0.5, ypos=0):
    return dict(title=title,
                font=dict(family='Balto'),
                autosize=False,
                width=width, 
                height=height,   
                showlegend=False, 
                hovermode='closest',
                
                updatemenus=[dict(type='buttons',
                                  showactive=False,
                                  y=ypos,
                                  x=xpos,
                                  xanchor='left',
                                  yanchor='top',
                                  pad=dict(t=1),
                                  buttons=[dict(label='Play',
                                                method='animate',
                                                args=[None, dict(frame=dict(duration=1, redraw=False), 
                                                                 transition=dict(duration=0),
                                                                 #easing=['bounce-in-out'],
                                                                 fromcurrent=True,
                                                                 mode='immediate'
                                                                 ),
                                                            
                                                     ]
                                                )
                                          ]
                                  )
                         ]
                )

#Easing functions specify the rate of change of a parameter over time.


# Define a lattice,  the numpy array of the spin starting points and directions:

# In[5]:


n=8
x=np.arange(-n, n)
y=np.arange(-n, n)
x,y=np.meshgrid(x,y)
start=np.stack((x,y))#(2,n,n) start[0,:,:] is x
norm=np.sqrt(x**2+y**2)
with np.errstate(over='ignore', divide='ignore', invalid='ignore'):# avoid warning messages pointing invalid 
                                                                   #values at the equilibria
    X=x/norm
    Y=y/norm
    
scale=0.6
theta=np.pi/12


# The unit vector fields of coordinates  $(-Y,X),  (Y, -X)$, represent two vortices, whereas 
# $(-Y,-X), (Y,X)$, are antivortices. We plot $(-Y,X)$ and $(-Y,-X)$.

# In[6]:


end1=start+scale*np.stack((-Y,X))
data1, framesV=anim_setup(start, end1,  theta, traces=[0,1])
for k in [0,1]:
    data1[k].update(xaxis='x1', yaxis='y1')
end2=start+scale*np.stack((-Y,-X))
data2, framesAV=anim_setup(start, end2,  theta, traces=[2,3])
for k in [0,1]:
    data2[k].update(xaxis='x2', yaxis='y1')


# In[7]:


data=data1+data2
frames=framesV+framesAV


# In[8]:


layout=set_layout(title='Gauge transformation of a vortex and antivortex')
axis=dict(showline=False, zeroline=False, mirror=False, showgrid=False, ticks='',showticklabels=False)
layout.update(xaxis1=dict(axis, **{'domain':[0,0.485]}),
              yaxis1=dict(axis),
              xaxis2=dict(axis, **{'domain':[0.515,1]})
               )
fig=dict(data=data, layout=layout, frames=frames)


# In[9]:


iplot(fig)


# In[10]:


from IPython.core.display import HTML
def  css_styling():
    styles = open("./custom.css", "r").read()
    return HTML(styles)
css_styling()