import plotly
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import plotly.graph_objects as go
from plotly.offline import plot

from IPython.core.display import display,HTML
import sympy as smp
from sympy.vector import cross

phi = np.linspace(0, 2*np.pi, 100)

def l(phi):
    return (1+3/4 * np.sin(3*phi)) * np.array([np.cos(phi), np.sin(phi), np.zeros(len(phi))])

lx, ly, lz = l(phi)

plt.figure(figsize=(7,7))
plt.plot(lx,ly)
plt.xlabel('$x/R$', fontsize=25)
plt.ylabel('$y/R$', fontsize=25)
plt.show()

solve integrand using sympy

t, x, y, z = smp.symbols('t, x, y, z') 

l = (1+(3/4)*smp.sin(3*t))*smp.Matrix([smp.cos(t), smp.sin(t), 0])
r = smp.Matrix([x,y,z])
sep=r-l

l
$\displaystyle \left[\begin{matrix}\left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)}\\\left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)}\\0\end{matrix}\right]$ 
r
$\displaystyle \left[\begin{matrix}x\\y\\z\end{matrix}\right]$ 
sep
$\displaystyle \left[\begin{matrix}x - \left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)}\\y - \left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)}\\z\end{matrix}\right]$

Define integrand

integrand = smp.diff(l,t).cross(sep) / sep.norm()**3

Get x, y z components of the integrand

integrand
$\displaystyle \left[\begin{matrix}\frac{z \left(\left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)} + 2.25 \sin{\left(t \right)} \cos{\left(3 t \right)}\right)}{\left(\left|{z}\right|^{2} + \left|{x - \left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)}}\right|^{2} + \left|{y - \left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)}}\right|^{2}\right)^{\frac{3}{2}}}\\- \frac{z \left(- \left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)} + 2.25 \cos{\left(t \right)} \cos{\left(3 t \right)}\right)}{\left(\left|{z}\right|^{2} + \left|{x - \left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)}}\right|^{2} + \left|{y - \left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)}}\right|^{2}\right)^{\frac{3}{2}}}\\\frac{- \left(x - \left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)}\right) \left(\left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)} + 2.25 \sin{\left(t \right)} \cos{\left(3 t \right)}\right) + \left(y - \left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)}\right) \left(- \left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)} + 2.25 \cos{\left(t \right)} \cos{\left(3 t \right)}\right)}{\left(\left|{z}\right|^{2} + \left|{x - \left(0.75 \sin{\left(3 t \right)} + 1\right) \cos{\left(t \right)}}\right|^{2} + \left|{y - \left(0.75 \sin{\left(3 t \right)} + 1\right) \sin{\left(t \right)}}\right|^{2}\right)^{\frac{3}{2}}}\end{matrix}\right]$ 
dBxdt = smp.lambdify([t,x,y,z], integrand[0])
dBydt = smp.lambdify([t,x,y,z], integrand[1])
dBzdt = smp.lambdify([t,x,y,z], integrand[2])

dBxdt(np.pi, 1, 1, 1)

-0.0680413817439772
quad(dBxdt, 0, 2*np.pi, args=(1,1,1)) # gives value and error

(0.367215052854198, 6.916483780662426e-09)

Get the magnetic field by performing the integral over each component

def B(x,y,z):
    return np.array([quad(dBxdt, 0, 2*np.pi, args=(x,y,z))[0],
                     quad(dBydt, 0, 2*np.pi, args=(x,y,z))[0],
                     quad(dBzdt, 0, 2*np.pi, args=(x,y,z))[0]])

B(0.5,0.5,0)

array([ 0.        ,  0.        , 10.87779227])
B(0.5,0.5,1)

array([0.19069963, 0.52786431, 1.52524645])

Set up a meshgrid to solve for the field in some 3D volume

x = np.linspace(-2,2,20)

xv,yv,zv = np.meshgrid(x,x,x)

B_field = np.vectorize(B, signature='(),(),()->(n)')(xv,yv,zv)
Bx = B_field[:,:,:,0]
By = B_field[:,:,:,1]
Bz = B_field[:,:,:,2]

use plotly for 3D intractive plot

xv.ravel()

array([-2., -2., -2., ...,  2.,  2.,  2.])
data = go.Cone(x=xv.ravel(), y=yv.ravel(), z=zv.ravel(),
               u=Bx.ravel(), v=By.ravel(), w=Bz.ravel(),
               colorscale='Inferno', colorbar=dict(title='$x^2$'),
               sizemode="absolute", sizeref=20)

layout = go.Layout(title=r'Biot Savard Law ',
                     scene=dict(xaxis_title=r'x',
                                yaxis_title=r'y',
                                zaxis_title=r'z',
                                aspectratio=dict(x=1, y=1, z=1),
                                camera_eye=dict(x=1.2, y=1.2, z=1.2)))

fig = go.Figure(data = data, layout=layout)
   
fig.add_scatter3d(x=lx, y=ly, z=lz, mode='lines',
                  line = dict(color='green', width=10))