# coding: utf-8

# # Kussmaul-Martensen quadrature (also often called "Kress quadrature")

# In[1]:

import numpy as np
import matplotlib.pyplot as pt


# In[2]:

t = np.linspace(0, 2*np.pi, 300,endpoint=False)


# ## Setup

# Let's make a curve and pick a target point:

# In[29]:

uncircleness = 1

path = np.array([
    np.cos(t) + uncircleness*0.2*np.sin(3*t),
    np.sin(t) + uncircleness*0.1*np.sin(3*t)
    ])

tgt_index = len(t)//2
tgt_t = t[tgt_index]
tgt = path[:, tgt_index]

pt.gca().set_aspect("equal")
pt.plot(path[0], path[1])
pt.plot(tgt[0], tgt[1], "o")


# Get some derivatives of the curve:

# In[28]:

import scipy.fftpack as fft

dpath_dt = np.array([
    fft.diff(path[0]),
    fft.diff(path[1]),
    ])

dpdt_squared = dpath_dt[0]**2 + dpath_dt[1]**2
pt.plot(dpdt_squared)


# Get normals to the curve:

# In[27]:

normals = np.array([
    dpath_dt[1],
    -dpath_dt[0]
    ]) / np.sqrt(dpdt_squared)

pt.plot(path[0], path[1])
pt.quiver(path[0], path[1], normals[0], normals[1])


# In[26]:

dist_vec = tgt[:, np.newaxis] - path

dist = np.sqrt(np.sum(dist_vec**2, axis=0))

pt.plot(dist)


### Single-layer potential

# Let's look at the integrand for the SLP:

# In[25]:

slp_integrand = np.log(1/dist)
pt.plot(slp_integrand)


# Even if this is integrable--Gaussian quadrature will do a terrible job. **Why?**

# In[24]:

near_sing_slice = slice(tgt_index-20, tgt_index+20)

log_sin_squared = -0.5*np.log(4*np.sin((tgt_t - t)/2)**2)
pt.plot(log_sin_squared[near_sing_slice])
pt.plot(slp_integrand[near_sing_slice])


# In[23]:

slp_subtracted = slp_integrand - log_sin_squared
pt.plot(slp_subtracted[near_sing_slice])


# In[17]:

pt.plot(slp_subtracted)


# **How does this help?**

#### Double-layer potential

# In[22]:

grad_slp = dist_vec/dist**2

dlp_integrand = np.sum(grad_slp * normals, axis=0)
pt.plot(dlp_integrand)


#### S'

# In[21]:

sp_integrand = np.sum(grad_slp * normals[:, tgt_index, np.newaxis], axis=0)
pt.plot(sp_integrand)


# ## Questions
# 
# * How would you apply this for **Helmholtz**?
# * Name aspects that make this rule slightly **impractical**
# * How would this apply to D'?