#!/usr/bin/env python
# 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[3]:


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[4]:


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[5]:


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[6]:


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[7]:


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


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

# In[8]:


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[9]:


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


# In[10]:


pt.plot(slp_subtracted)


# **How does this help?**

# ### Double-layer potential

# In[11]:


grad_slp = dist_vec/dist**2

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


# ### S'

# In[12]:


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'?