# kahansum.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

def kahansum(xs):
    s = 0.; e = 0.
    for x in xs:
        temp = s
        y = x + e
        s = temp + y
        e = (temp - s) + y
    return s

if __name__ == '__main__':
    xs = [0.7, 0.1, 0.3]
    print(sum(xs), kahansum(xs))

# naiveval.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

from math import sqrt

def naiveval(x):
    return 1/(sqrt(x**2 + 1) - x)

xs = [10**i for i in range(4,8)]
ys = [naiveval(x) for x in xs]
for x, y in zip(xs, ys):
    print(x, y)

# compexp.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

from math import exp

def compexp(x):
    n = 0
    oldsum, newsum, term = 0., 1., 1.
    while newsum != oldsum:
        oldsum = newsum
        n += 1
        term *= x/n
        newsum += term
        print(n, newsum, term)
    return newsum

for x in (0.1, 20., -20.):
    print("x, library exp(x):", x, exp(x))
    val = compexp(x)

# recforw.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

from math import exp

def forward(nmax=22):
    oldint = 1 - exp(-1)
    for n in range(1,nmax):
        print(n-1, oldint)
        newint = n*oldint - exp(-1)
        oldint = newint

print("n = 20 answer is 0.0183504676972562")
print("n, f[n]")
forward()

# recback.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

from math import exp

def backward(nmax=31):
    oldint = 0.01
    for n in reversed(range(20,nmax)):
        print(n, oldint)
        newint = (oldint + exp(-1))/n
        oldint = newint

print("n = 20 answer is 0.0183504676972562")
print("n, f[n]")
backward()

# cancel.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

from math import exp, log

def f(x):
    return (exp(x) - 1)/x

def g(x):
    w = exp(x)
    if w==0.:
        val = -1/x
    elif w==1.:
        val = 1.
    else:
        val = (w-1)/log(w)
    return val

xs = [10**(-i) for i in (14, 15, 16)]
xs += [-10**(-i) for i in (15, 16, 17)]

fvals = [f(x) for x in xs]
gvals = [g(x) for x in xs]

for x, fval, gval in zip(xs, fvals, gvals):
    print(x, fval, gval)

# chargearray.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

from kahansum import kahansum
from math import sqrt

def chargearray(nvals):
    vals = [-0.5 + i/(nvals-1) for i in range(nvals)]
    qtopos = {}
    for i,posx in enumerate(vals):
        for j,posy in enumerate(vals):
            count = j + nvals*i + 1
            key = 1.02*count if (i+j)%2==0 else -count
            qtopos[key] = posx, posy
    return qtopos

def vecmag(rs):
    sq = [r**2 for r in rs]
    return sqrt(kahansum(sq))

def fullpot(qtopos,rs):
    potvals = []
    for q,pos in qtopos.items():
        diffs = [r - po for r,po in zip(rs,pos)]
        R = vecmag(diffs)
        potvals.append(q/R)
    return kahansum(potvals)

if __name__ == '__main__':
    qtopos = chargearray(6)
    for y in 1,-1:
        rs = [0.,y]
        potval = fullpot(qtopos,rs)
        print(rs, potval)

# legendre.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

import matplotlib.pyplot as plt

def legendre(n,x):
    if n==0:
        val2 = 1.
        dval2 = 0.
    elif n==1:
        val2 = x
        dval2 = 1.
    else:
        val0 = 1.; val1 = x
        for j in range(1,n):
            val2 = ((2*j+1)*x*val1 - j*val0)/(j+1)
            val0, val1 = val1, val2
        dval2 = n*(val0-x*val1)/(1.-x**2)
    return val2, dval2

def plotlegendre(der,nsteps):
    plt.xlabel('$x$', fontsize=20)

    dertostr = {0: "$P_n(x)$", 1: "$P_n'(x)$"}
    plt.ylabel(dertostr[der], fontsize=20)

    ntomarker = {1: 'k-', 2: 'r--', 3: 'b-.', 4: 'g:', 5: 'c^'}
    xs = [i/nsteps for i in range (-nsteps+1,nsteps)]
    for n,marker in ntomarker.items():
        ys = [legendre(n,x)[der] for x in xs]
        labstr = 'n={0}'.format(n)
        plt.plot(xs, ys, marker, label=labstr, linewidth=3)

    plt.ylim(-3*der-1, 3*der+1)
    plt.legend(loc=4)
    plt.show()

if __name__ == '__main__':
    nsteps = 200
    plotlegendre(0,nsteps)
    plotlegendre(1,nsteps)

# multipole.py
# Author: Alex Gezerlis
# Numerical Methods in Physics with Python (CUP, 2020)

from kahansum import kahansum
from chargearray import chargearray, vecmag
from legendre import legendre

def decomp(rs,ris):
    rmag = vecmag(rs); rimag = vecmag(ris)
    prs = [r*ri for r,ri in zip(rs,ris)]
    vecdot = kahansum(prs)
    costheta = vecdot/(rmag*rimag)
    return rmag, rimag, costheta

def multicoes(rs,qtopos,nmax=60):
    coes = [0. for n in range(nmax+1)]
    for n in range(nmax+1):
        for q,pos in qtopos.items():
            rmag, rimag, costheta = decomp(rs,pos)
            val = q*(rimag**n)*legendre(n,costheta)[0]
            coes[n] += val
    return coes

def multifullpot(rs,qtopos):
    coes = multicoes(rs,qtopos)
    rmag = vecmag(rs)
    contribs = [coe/rmag**(n+1) for n,coe in enumerate(coes)]
    return kahansum(contribs)

if __name__ == '__main__':
    qtopos = chargearray(6)
    for y in 1,-1:
        rs = [0.,y]
        potval = multifullpot(rs,qtopos)
        print(rs, potval)