SHARE
TWEET

Untitled

a guest Mar 8th, 2019 71 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.  
  33.     lat0, lat1 = rads*(90-50), rads*(90+50)
  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. OK, I Understand
 
Top