-
Notifications
You must be signed in to change notification settings - Fork 5
/
plot_J_surface.py
102 lines (80 loc) · 2.91 KB
/
plot_J_surface.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
"""
Current surfaces
----------------
This example demonstrates the resulting current surfaces
from the curl-free and divergence-free SEC poles.
"""
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
from pysecs import SECS
# Radius of Earth
R_E = 6378e3
# Pole of the current system at the North Pole
sec_loc = np.array([90., 0., R_E + 100e3])
# Set up the system with a single SEC
system_df = SECS(sec_df_loc=sec_loc)
system_cf = SECS(sec_cf_loc=sec_loc)
# Fit unit currents since we aren't fitting to any data
system_df.fit_unit_currents()
system_cf.fit_unit_currents()
# Scale the system corresponding to Figure 2 of
# Amm and Viljannen (1999) doi:10.1186/BF03352247
# 10 kA
total_current = 10000
system_df.sec_amps *= total_current
system_cf.sec_amps *= total_current
# Create prediction surface
delta = 1
lat_min, lat_max = -90, 90
lon_min, lon_max = -180, 180
lats = np.arange(lat_min, lat_max, delta)
lons = np.arange(lon_min, lon_max, delta)
lats2 = np.arange(lat_min-delta/2, lat_max+delta/2, delta)
lons2 = np.arange(lon_min-delta/2, lon_max+delta/2, delta)
nlat = len(lats)
nlon = len(lons)
xx, yy = np.meshgrid(lons, lats)
xx2, yy2 = np.meshgrid(lons2, lats2)
points = np.zeros((nlat*nlon, 3))
points[:, 0] = yy.ravel()
points[:, 1] = xx.ravel()
points[:, 2] = R_E
shell_points = points.copy()
shell_points[:, 2] = sec_loc[2]
# Predict the currents
J_df = system_df.predict(shell_points, J=True)
J_cf = system_cf.predict(shell_points, J=True)
J_df_x = J_df[:, 0].reshape(nlat, nlon)
J_df_y = J_df[:, 1].reshape(nlat, nlon)
J_df_phi = np.linalg.norm(J_df[:, :2], axis=1).reshape(nlat, nlon)
J_df_r = J_df[:, 2].reshape(nlat, nlon)
J_cf_x = J_cf[:, 0].reshape(nlat, nlon)
J_cf_y = J_cf[:, 1].reshape(nlat, nlon)
J_cf_theta = np.linalg.norm(J_cf[:, :2], axis=1).reshape(nlat, nlon)
J_cf_r = J_cf[:, 2].reshape(nlat, nlon)
fig, (ax_cf, ax_df) = plt.subplots(figsize=(8, 3), ncols=2,
constrained_layout=True)
fig.suptitle('Amm and Viljanen (1999) Figure 1')
# Limit quiver to only a few points
nstep = 30
norm = LogNorm(vmin=1e-5, vmax=1e-2)
cax = ax_cf.pcolormesh(xx2, yy2, J_cf_theta, norm=norm)
ax_cf.quiver(xx[nstep//2::nstep, nstep//2::nstep*2],
yy[nstep//2::nstep, nstep//2::nstep*2],
J_cf_y[nstep//2::nstep, nstep//2::nstep*2],
J_cf_x[nstep//2::nstep, nstep//2::nstep*2],
color='w', zorder=9)
ax_cf.set_ylabel('Latitude (deg)')
ax_cf.set_xlabel('Longitude (deg)')
cax = ax_df.pcolormesh(xx2, yy2, J_df_phi, norm=norm)
ax_df.quiver(xx[nstep//2::nstep, nstep//2::nstep*2],
yy[nstep//2::nstep, nstep//2::nstep*2],
J_df_y[nstep//2::nstep, nstep//2::nstep*2],
J_df_x[nstep//2::nstep, nstep//2::nstep*2],
color='w', zorder=9)
ax_df.set_ylabel('Latitude (deg)')
ax_df.set_xlabel('Longitude (deg)')
ax_cf.set_title('Curl Free')
ax_df.set_title('Divergence Free')
plt.show()