• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Mar 8th, 2019 79 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1.     def dotem(a, b):
2.         return (a*b).sum(axis=0)
3.
4.     def get_B2D(x):
5.
6.         r    = (x - x0)/R0
7.         rabs = np.sqrt((r**2).sum(axis=0)) # [..., None]
8.         rhat = r / rabs
9.
10.         B = B0 * (3*dotem(p, rhat) * rhat - p) / rabs**3
11.
12.         return B
13.
14.     import numpy as np
15.     import matplotlib.pyplot as plt
16.
17.     halfpi, pi, twopi, fourpi = [f*np.pi for f in (0.5, 1, 2, 4)]
18.     degs, rads                = 180/pi, pi/180
19.
20.     mu0  = fourpi * 1E-07
21.     R0   = 6371.     # km
22.     B0   = 3.12E-5   # T
23.     inc  = rads * 11.0
24.
25.     p    = np.array([f(inc) for f in (np.sin, np.zeros_like, np.cos)], dtype=float)[:, None, None]
26.     x0   = np.array([-500, 0, 0], dtype=float)[:, None, None]  # offset (manifests as South Atlantic Anomaly)
27.
28.     pp   = p[:, 0, 0]
29.     xx0  = x0[:, 0, 0]
30.
31.     R    = 6778.
32.
34.     lat      = np.linspace( lat0,  lat1, 301) # 0 to 180 degrees (spherical coordinates)
35.     lon      = np.linspace(-pi, pi, 601)
36.     LAT, LON = np.meshgrid(lat, lon, indexing='ij')
37.     (clat, clon), (slat, slon) = [[f(x) for x in (LAT, LON)] for f in (np.cos, np.sin)]
38.
39.     XYZ = np.stack([R * x for x in (slat*clon, slat*slon, clat)])
40.
41.     B = get_B2D(XYZ)
42.     Bmag = np.sqrt((B**2).sum(axis=0))
43.
44.     if True:
45.         plt.figure()
46.         plt.imshow(Bmag)
47.         plt.colorbar()
48.         plt.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.
Not a member of Pastebin yet?