#!/usr/bin/python
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

def mmag(y,z):
    x=0.
    r=np.sqrt(x**2+y**2+z**2)
    if r >= 1.:
     bx=3*x*z/r**5
     by=3*y*z/r**5
     bz=(2.*z**2-x**2-y**2)/r**5
    else:
     bx=by=bz=0.
    return by,bz
    
mmag=np.vectorize(mmag)
y=np.linspace(-6,6,num=70)
z=np.linspace(-6,6,num=70)
Y,Z=np.meshgrid(y,z)
By,Bz=mmag(Y,Z)
#*** visualization ******
fig=plt.figure(figsize=(10,10))
p1=plt.subplot()
p1.set_aspect("equal")
maa=plt.Circle((0,0),1.,color="green")
air=plt.Circle((0,0),1.016,fill=False,color="blue")
p1.add_artist(maa)
p1.add_artist(air)
p1.streamplot(Y,Z,By,Bz,color="red",arrowsize=1.5,density=1)

#p2=plt.subplot(122)
#p2.set_aspect("equal")
#maa=plt.Circle((0,0),1.,color="green")
#air=plt.Circle((0,0),1.016,fill=False,color="blue")
#cs=p2.contour(Z,Y,bmaa,10)
#p2.clabel(cs,fontsize=10)
#p2.add_artist(maa)
#p2.add_artist(air)

plt.show()