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

# MagBottleMat.py

"""  <p>From "COMPUTATIONHL PHYSICS" & "COMPUTER PROBLEMS in PHYSICS"
    by RH Landau, MJ Paez, and CC Bordeianu (deceased).
    Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, 
    C Bordeianu, Univ Bucharest, 2020. 
    Please respect copyright & acknowledge our work.
        <p>
In terms of elliptic integrals, the B field between two parallel coils with I in the same direction is:
     \begin{eqnarray}
B_x =B_0 \frac{1}{\pi \sqrt{Q}}\bigl[E(k) \frac{1-\alpha^2-\beta^2}{Q-4\alpha} +K(k)\bigr ]\\
B_r =B_0 \frac{\gamma}{\pi \sqrt{Q}}\bigl[E(k) \frac{1+\alpha^2+\beta^2}{Q-4\alpha} -K(k)\bigr ]\\
\alpha= \frac{r}{a}, \ \ \  \beta= \frac{x}{a}, \ \ \  \gamma= \frac{x}{r}, \ \ \  Q=(1+\alpha)^2 +\beta^2, \ \ \  k= \sqrt{\frac{4\alpha}{Q}}.\\
 K(k) = \int_0^{\pi/2} \, \frac{d\theta}{\sqrt{1-k^2\sin^2 \theta}}, \ \ 
 E(k)= \int_0^{\pi/2} \, \sqrt{1-k^2\sin^2 t}\,dt
 \end{eqnarray}
 The program computes the integrals numerically. """

# In[ ]:

get_ipython().run_line_magic('matplotlib', 'notebook')
from IPython.display import Image
import numpy as np
import matplotlib.pyplot as plt
ra=0.1   # coil radius
ra2=ra**2
nmax=100
dx=0.01                           #take  dx = dy
n=300
nm=int(nmax/2)
nplot=200
Bx=np.zeros((nmax,nmax),float)
By=np.zeros((nmax,nmax),float)
xx=np.zeros(nplot,float)
yy=np.zeros(nplot,float)
tt=np.zeros(nplot,float)
def funct(k,t):
    f=np.sqrt(1-k*np.sin(t)**2)
    return f
def funth(k,th):                       # integrand elliptic integral
    y=np.sqrt(1-k*np.sin(th)**2)
    return 1./y        
   
def trap(a,b,n,func,x):            # trapezoidal rule limits: a, b. n=points
    h=(b-a)/(n-1)                  # step, N-1 intervals
    sum=(func(x,a)+func(x,b))/2    # first, last mulp. by h/2
    for i in range(1, n -1):       #trap rule
        sum+=func(x,a+i*h)
    return h*sum
def Ellipints(k):      # returns 1st and 2nd kind elliptic int for given k
    a=0
    b=np.pi/2.
    Ek=trap(a,b,n,funct,k)
    Kk=trap(a,b,n,funth,k)
    return Kk,Ek                      # return integrals 1st and 2nd kinds
i=0
j=1
def Bfield(dx,x,y):         
    alpha = y/ra
    beta=x/ra
    gamma=x/y
    xx=1.-x
    betp=xx/ra
    gamp=xx/y
    term=(1+alpha)**2
    Q=term + beta**2
    Qp=term + betp**2
    Qsq=np.sqrt(Q)
    Qpsq=np.sqrt(Qp)
    k2=4*alpha/Q
    kp2=4*alpha/Qp
    Kk,Ek=Ellipints(k2)
    Kkp,Ekp=Ellipints(kp2)
    fac=Ek/(Q-4*alpha)
    albet=alpha*alpha+beta*beta
    facp=Ekp/(Qp-4*alpha)
    albetp=alpha*alpha+betp*betp
    Bxp=(facp*(1-albetp)+Kkp)/Qpsq
    Bx=(fac*(1-albet)+Kk)/Qsq+Bxp
    Byp=gamp*(facp*(1+albetp)-Kkp)/Qpsq
    By=gamma*(fac*(1+albet)-Kk)/Qsq-Byp 
    return 100*Bx,100*By
x=np.arange(0,1.0,0.01)
y=np.arange(-0.5,0.5,0.01) 
X,Y=np.meshgrid(x,y)   
 
fig = plt.figure()
ax = fig.add_subplot(111) 
Bx,By=Bfield(dx,X,Y)

ax.streamplot(x,y,Bx,By)
ax.set_aspect('equal')
ax.set_title('Magnetic field lines')
ax.set_xlabel('x')
ax.set_ylabel('y')

def Euler(dx):
    q=1
    m=10
    dt=0.1
    v=[0.1,0.1]
    r=[0.6,0.8]
    i=0
    for t in range(0,10):
        ii=int(r[0]/dx)    
        jj=int(r[1]/dx)
        bx=(Bx[ii,jj]+ Bx[ii+1,jj])/2
        by=(Bx[ii,jj]+ Bx[ii,jj+1])/2
        B=[bx,by]
        F=q*(np.cross(v,B))
        a=F/m
        v=v+a*dt
        r=r+v*dt
        print(r)
        xx[i]=r[0]
        yy[i]=r[1]
        print(xx[i],yy[i])
        i=i+1
    #F=
#l = plt.axvline(x=0, linewidth=2, color='g')# surface
Euler(dx)

#ax.plot(xx,yy,c='r')
print("finished")
plt.show()


# In[ ]:





# In[ ]: