Python Contour plot of Earth's magnetic field in Python

[jackthestudent]

The Earth's magnetic field is cylindrically symmetric and can be described using cylindrical coordinates, (ρ,z,φ). The components of the magnetic field are given by:

(B_φ) = 0

(B_ρ) = −∂A/∂z,

(B_z) =(1/ρ) * ∂(ρA)/∂ρ.

The vector potential for a dipole field is given by

(A_φ) = (B_e)(ρ) * (R_e)^3/r^3

where B_e and R_e are the equatorial surface magnetic field strength and the radius of the Earth, respectively, and r is the radial distance from the center of the planet. Magnetic field lines are given by contours of a flux function F given by

F = ρ(A_φ).

The Earth's surface magnetic field at the equator B_e is 3,120 nT and its radius R_e is 6,371 km.

Using this information, create a contour plot of the Earth's magnetic field along a bisecting plane covering ρ<20R_e (it's symmetric, so the horizontal coordinate goes from -20R_e to +20R_e and −10R_e < z< 10R_e. Use the following steps:

1) Start by creating the coordinate grids for ρ and z, and the radius r. A grid cell size of ~0.02 R_e should be fine.

2) Compute the vector potential A and the flux function F.

Show the magnetic field lines (which are contours of F, computed from A as shown above) in a plot using about 5 or so black lines roughly equally spaced at the equator.

-------------------------------------------------
I have created 2 arrays, one for rho and one for z and then I have turned them into a 2D array in the form of F, the flux function. This is what I have entered into the interpreter:

Re = 6371000
rho = np.linspace(-20*Re, 20*Re, 2000)
z1 = np.linspace(-10*Re, 10*Re, 2000)
z = z1[:, np.newaxis]
Be = 3120*10**(-9)
F = (rho**2)*Be*(Re/np.sqrt(rho**2 + z**2))**3
plt.contour(rho, z1, F)

It returns a contour map with one dot on zero and nothing else even though there should be several contours. I can't figure out what I've done wrong or what to try next? Any help would be greatly appreciated :)

 

[phyzguy]

I haven't gone through all of your code, but you haven't set up the plots properly. For one thing, the radius rho should never be negative. For another thing, what is the horizontal coordinate in your plots? You don't want it to be rho, you want to define an x-coordinate for plotting. A third thing I'm not clear on in your coordinate system is what is the difference between rho and r?

 

[jackthestudent]

Thank you for your help! It turns out the problem was that Matplotlib didn't automatically show any contours so I had to define the levels I wanted manually. This is the code I used (I have also changed rho and z to meshgrids):
Re = 6371000
rho1 = np.linspace(-20*Re, 20*Re, 2000)
z1 = np.linspace(-10*Re, 10*Re, 2000)
rho, z = np.meshgrid(rho1, z1)
Be = 3120*10**(-9)
r = np.sqrt(rho**2 + z**2)
F = rho**2*Be*(Re/r)**3
plt.contour(rho, z, F, [6E6, 7E6, 8E6, 9E6, 1E7], colors = 'k')

In answer to your questions:
rho is the horizontal coordinate, z is the vertical coordinate and r (the radius from the centre of the planet) is equal to sqrt(rho^2 + z^2). The cylindrical coordinates are basically Cartesian coordinates in this case because the B_φ = 0.

 

[phyzguy]

In answer to your questions:
rho is the horizontal coordinate, z is the vertical coordinate and r (the radius from the centre of the planet) is equal to sqrt(rho^2 + z^2). The cylindrical coordinates are basically Cartesian coordinates in this case because the B_φ = 0.

I think you're still confused on your coordinate system. In cylindrical coordinates, it is definitely not true that "The cylindrical coordinates are basically Cartesian coordinates in this case because the B_φ = 0" I suggest that you study these coordinate systems a little more.

 

[jackthestudent]

I think you're still confused on your coordinate system. In cylindrical coordinates, it is definitely not true that "The cylindrical coordinates are basically Cartesian coordinates in this case because the B_φ = 0" I suggest that you study these coordinate systems a little more.

Sorry I worded that badly- I meant that rho and z are perpendicular to each other and when φ = 0 they form the opposite and adjacent sides of a right angled triangle so I could have used a Cartesian coordinate system in the first place (not that the cylindrical and Cartesian coordinates are the same) if I'm making any sense?