This notebook introduces useful space science models included in davitpy.
Currently we have ported/wrapped the following models to python:
%pylab inline
from datetime import datetime as dt
from davitpy.models import *
from davitpy import utils
Populating the interactive namespace from numpy and matplotlib
# INPUTS
itype = 1 # Geodetic coordinates
pyDate = dt(2006,2,23)
date = utils.dateToDecYear(pyDate) # decimal year
alt = 300. # altitude
stp = 5.
xlti, xltf, xltd = -90.,90.,stp # latitude start, stop, step
xlni, xlnf, xlnd = -180.,180.,stp # longitude start, stop, step
ifl = 0 # Main field
# Call fortran subroutine
lat,lon,d,s,h,x,y,z,f = igrf.igrf11(itype,date,alt,ifl,xlti,xltf,xltd,xlni,xlnf,xlnd)
# Check that it worked by plotting magnetic dip angle contours on a map
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import meshgrid
# Set figure
fig = figure(figsize=(10,5))
ax = fig.add_subplot(111)
rcParams.update({'font.size': 14})
# Set-up the map background
map = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90,\
llcrnrlon=-180,urcrnrlon=180,resolution='c')
map.drawmapboundary()
map.drawcoastlines(color='0.5')
# draw parallels and meridians.
map.drawparallels(np.arange(-80.,81.,20.))
map.drawmeridians(np.arange(-180.,181.,20.))
# The igrf output needs to be reshaped to be plotted
dip = s.reshape((180./stp+1,360./stp+1))
dec = d.reshape((180./stp+1,360./stp+1))
lo = lon[0:(360./stp+1)]
la = lat[0::(360./stp+1)]
x,y = meshgrid(lo,la)
v = arange(0,90,20)
# Plot dip angle contours and labels
cs = map.contour(x, y, abs(dip), v, latlon=True, linewidths=1.5, colors='k')
labs = plt.clabel(cs, inline=1, fontsize=10)
# Plot declination and colorbar
im = map.pcolormesh(x, y, dec, vmin=-40, vmax=40, cmap='coolwarm')
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", "3%", pad="3%")
colorbar(im, cax=cax)
cax.set_ylabel('Magnetic field declination')
cticks = cax.get_yticklabels()
cticks = [t.__dict__['_text'] for t in cticks]
cticks[0], cticks[-1] = 'W', 'E'
_ = cax.set_yticklabels(cticks)
savefig('dipdec.png')
/usr/local/lib/python2.7/dist-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal if rotation in ('horizontal', None): /usr/local/lib/python2.7/dist-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal elif rotation == 'vertical':
# Check that it worked by plotting magnetic dip angle contours on a map
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import meshgrid
# Set figure
fig = figure(figsize=(10,5))
ax = fig.add_subplot(111)
rcParams.update({'font.size': 14})
# Set-up the map background
map = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90,\
llcrnrlon=-180,urcrnrlon=180,resolution='c')
map.drawmapboundary()
map.drawcoastlines(color='0.5')
# draw parallels and meridians.
map.drawparallels(np.arange(-80.,81.,20.))
map.drawmeridians(np.arange(-180.,181.,20.))
# The igrf output needs to be reshaped to be plotted
babs = f.reshape((180./stp+1,360./stp+1))
lo = lon[0:(360./stp+1)]
la = lat[0::(360./stp+1)]
x,y = meshgrid(lo,la)
v = arange(0,90,20)
# Plot declination and colorbar
im = map.pcolormesh(x, y, babs, cmap='jet')
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", "3%", pad="3%")
colorbar(im, cax=cax)
cax.set_ylabel('Magnetic field intensity [nT]')
savefig('babs.png')
JF switches to turn off/on (True/False) several options
[0] : True
[1] : True
[2] : True
[3] : False
[4] : False
[5] : False
[6] : True
[7] : True
[8] : True
[9] : True
[10] : True
[11] : True
[12] : True
[13] : True
[14] : True
[15] : True
[16] : True
[17] : True
[18] : True
[19] : True
[20] : False
[21] : True
[22] : False
[23] : True
[24] : True
[25] : True
[26] : True
[27] : False
[28] : False
[29] : False
[28,29]:
[30] : True
[31] : True
[32] : False
[33] : True
[34] : False
[..] : ....
[50] : ....
# Inputs
jf = [True]*50
jf[2:6] = [False]*4
jf[20] = False
jf[22] = False
jf[27:30] = [False]*3
jf[32] = False
jf[34] = False
jmag = 0.
alati = 40.
along = -80.
iyyyy = 2012
mmdd = 806
dhour = 12.
heibeg, heiend, heistp = 80., 500., 10.
oarr = np.zeros(100)
# Call fortran subroutine
outf,oarr = iri.iri_sub(jf,jmag,alati,along,iyyyy,mmdd,dhour,heibeg,heiend,heistp,oarr)
# Check that it worked by plotting vertical electron density profile
figure(figsize=(5,8))
alt = np.arange(heibeg,heiend,heistp)
ax = plot(outf[0,0:len(alt)],alt)
xlabel(r'Electron density [m$^{-3}$]')
ylabel('Altitude [km]')
grid(True)
rcParams.update({'font.size': 12})
lats = range(10, 90, 10)
lons = zeros(len(lats))
rhos = 6372.*ones(len(lats))
trace = tsyganenko.tsygTrace(lats, lons, rhos)
print trace
ax = trace.plot()
/usr/local/lib/python2.7/dist-packages/davitpy-0.2-py2.7-linux-x86_64.egg/davitpy/models/tsyganenko/__init__.py:80: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future. assert (None not in [lat, lon, rho]) or filename, 'You must provide either (lat, lon, rho) or a filename to read from' /usr/local/lib/python2.7/dist-packages/davitpy-0.2-py2.7-linux-x86_64.egg/davitpy/models/tsyganenko/__init__.py:82: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future. if None not in [lat, lon, rho]: /usr/local/lib/python2.7/dist-packages/davitpy-0.2-py2.7-linux-x86_64.egg/davitpy/models/tsyganenko/__init__.py:431: RuntimeWarning: invalid value encountered in less indMask = sign*self.yTrace[ip,0:self.l[ip]] < 0
vswgse=[ -400, 0, 0] [m/s] pdyn= 2 [nPa] dst= -5 [nT] byimf= 0 [nT] bzimf= -5 [nT] Coords: geo (latitude [degrees], longitude [degrees], distance from center of the Earth [km]) (10.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(13.613, 359.855, 6371.131) --> SH( 9.924, 0.004, 6371.164) (20.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(20.019, 360.000, 6371.191) --> SH( 3.264, 0.744, 6371.132) (30.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(30.008, 360.000, 6371.198) --> SH(-8.075, 2.517, 6371.173) (40.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(40.005, 360.000, 6371.199) --> SH(-21.727, 6.465, 6371.188) (50.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(50.003, 360.000, 6371.200) --> SH(-39.383, 15.079, 6371.187) (60.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(60.002, 360.000, 6371.200) --> SH(-56.706, 31.776, 6371.182) (70.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(70.001, 360.000, 6371.200) --> SH(-61.419, 86.864, 6371.196) (80.000, 0.000, 6372.000) @ 05:47 UT (12-Mar-15) --> NH(80.001, 360.000, 6371.200) --> SH( nan, nan, nan)
# Inputs
# Date and time
year = 2000
doy = 1
hr = 1
mn = 0
sc = 0
# Solar wind speed
vxgse = -400.
vygse = 0.
vzgse = 0.
# Execution parameters
lmax = 5000
rlim = 60.
r0 = 1.
dsmax = .01
err = .000001
# Direction of the tracing
mapto = 1
# Magnetic activity [SW pressure (nPa), Dst, ByIMF, BzIMF]
parmod = np.zeros(10)
parmod[0:4] = [2, -8, -2, -5]
# Start point (rh in Re)
lat = 50.
lon = 0.
rh = 0.
# This has to be called first
tsyganenko.tsygFort.recalc_08(year,doy,hr,mn,sc,vxgse,vygse,vzgse)
# Convert lat,lon to geographic cartesian and then gsw
r,theta,phi, xgeo, ygeo, zgeo = tsyganenko.tsygFort.sphcar_08(1., np.radians(90.-lat), np.radians(lon), 0., 0., 0., 1)
xgeo,ygeo,zgeo,xgsw,ygsw,zgsw = tsyganenko.tsygFort.geogsw_08(xgeo, ygeo, zgeo,0,0,0,1)
# Trace field line
xfgsw,yfgsw,zfgsw,xarr,yarr,zarr,l = tsyganenko.tsygFort.trace_08(xgsw,ygsw,zgsw,mapto,dsmax,err,
rlim,r0,0,parmod,'T96_01','IGRF_GSW_08',lmax)
# Convert back to spherical geographic coords
xfgeo,yfgeo,zfgeo,xfgsw,yfgsw,zfgsw = tsyganenko.tsygFort.geogsw_08(0,0,0,xfgsw,yfgsw,zfgsw,-1)
gcR, gdcolat, gdlon, xgeo, ygeo, zgeo = tsyganenko.tsygFort.sphcar_08(0., 0., 0., xfgeo, yfgeo, zfgeo, -1)
print '** START: {:6.3f}, {:6.3f}, {:6.3f}'.format(lat, lon, 1.)
print '** STOP: {:6.3f}, {:6.3f}, {:6.3f}'.format(90.-np.degrees(gdcolat), np.degrees(gdlon), gcR)
** START: 50.000, 0.000, 1.000 ** STOP: -40.341, 17.675, 1.000
# A quick checking plot
from mpl_toolkits.mplot3d import proj3d
import numpy as np
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
# Plot coordinate system
ax.plot3D([0,1],[0,0],[0,0],'b')
ax.plot3D([0,0],[0,1],[0,0],'g')
ax.plot3D([0,0],[0,0],[0,1],'r')
# First plot a nice sphere for the Earth
u = np.linspace(0, 2 * np.pi, 179)
v = np.linspace(0, np.pi, 179)
tx = np.outer(np.cos(u), np.sin(v))
ty = np.outer(np.sin(u), np.sin(v))
tz = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(tx,ty,tz,rstride=10, cstride=10, color='grey', alpha=.5, zorder=2, linewidth=0.5)
# Then plot the traced field line
latarr = [10.,20.,30.,40.,50.,60.,70.,80.]
lonarr = [0., 180.]
rh = 0.
for lon in lonarr:
for lat in latarr:
r,theta,phi, xgeo, ygeo, zgeo = tsyganenko.tsygFort.sphcar_08(1., np.radians(90.-lat), np.radians(lon), 0., 0., 0., 1)
xgeo,ygeo,zgeo,xgsw,ygsw,zgsw = tsyganenko.tsygFort.geogsw_08(xgeo, ygeo, zgeo,0,0,0,1)
xfgsw,yfgsw,zfgsw,xarr,yarr,zarr,l = tsyganenko.tsygFort.trace_08(xgsw,ygsw,zgsw,mapto,dsmax,err,
rlim,r0,0,parmod,'T96_01','IGRF_GSW_08',lmax)
for i in xrange(l):
xgeo,ygeo,zgeo,dum,dum,dum = tsyganenko.tsygFort.geogsw_08(0,0,0,xarr[i],yarr[i],zarr[i],-1)
xarr[i],yarr[i],zarr[i] = xgeo,ygeo,zgeo
ax.plot3D(xarr[0:l],yarr[0:l],zarr[0:l], zorder=3, linewidth=2, color='y')
# Set plot limits
lim=4
ax.set_xlim3d([-lim,lim])
ax.set_ylim3d([-lim,lim])
ax.set_zlim3d([-lim,lim])
(-4, 4)
# Inputs
import datetime as dt
myDate = dt.datetime(2012, 7, 5, 12, 35)
glat = 40.
glon = -80.
mass = 48
# First, MSIS needs a bunch of input which can be obtained from tabulated values
# This function was written to access these values (not provided with MSIS by default)
solInput = msis.getF107Ap(myDate)
# Also, to switch to SI units:
msis.meters(True)
# Other input conversion
iyd = (myDate.year - myDate.year/100*100)*100 + myDate.timetuple().tm_yday
sec = myDate.hour*24. + myDate.minute*60.
stl = sec/3600. + glon/15.
altitude = linspace(0., 500., 100)
temp = zeros(shape(altitude))
dens = zeros(shape(altitude))
N2dens = zeros(shape(altitude))
O2dens = zeros(shape(altitude))
Odens = zeros(shape(altitude))
Ndens = zeros(shape(altitude))
Ardens = zeros(shape(altitude))
Hdens = zeros(shape(altitude))
Hedens = zeros(shape(altitude))
for ia,alt in enumerate(altitude):
d,t = msis.gtd7(iyd, sec, alt, glat, glon, stl, solInput['f107a'], solInput['f107'], solInput['ap'], mass)
temp[ia] = t[1]
dens[ia] = d[5]
N2dens[ia] = d[2]
O2dens[ia] = d[3]
Ndens[ia] = d[7]
Odens[ia] = d[1]
Hdens[ia] = d[6]
Hedens[ia] = d[0]
Ardens[ia] = d[4]
figure(figsize=(16,8))
#rcParams.update({'font.size': 12})
subplot(131)
plot(temp, altitude)
gca().set_xscale('log')
xlabel('Temp. [K]')
ylabel('Altitude [km]')
subplot(132)
plot(dens, altitude)
gca().set_xscale('log')
gca().set_yticklabels([])
xlabel(r'Mass dens. [kg/m$^3$]')
subplot(133)
plot(Odens, altitude, 'r-',
O2dens, altitude, 'r--',
Ndens, altitude, 'g-',
N2dens, altitude, 'g--',
Hdens, altitude, 'b-',
Hedens, altitude, 'y-',
Ardens, altitude, 'm-')
gca().set_xscale('log')
gca().set_yticklabels([])
xlabel(r'Density [m$^3$]')
leg = legend( (r'O',
r'O$_2$',
r'N',
r'N$_2$',
r'H',
r'He',
r'Ar',),
'upper right')
tight_layout()
/usr/local/lib/python2.7/dist-packages/matplotlib/cbook.py:137: MatplotlibDeprecationWarning: The "loc" positional argument to legend is deprecated. Please use the "loc" keyword instead. warnings.warn(message, mplDeprecation, stacklevel=1)
Input arguments:
Output argument:
w = hwm.hwm07(11001, 0., 200., 40., -80., 0, 0, 0, [0, 0])
print w
[ -4.81379271 67.79786682]
models.aacgm.aacgmConv(lat,lon,alt,flg)
convert between geographic coords and aacgm
Input arguments:
Outputs:
#geo to aacgm
glat,glon,r = aacgm.aacgmConv(42.0,-71.4,300.,2000,0)
print glat, glon, r
53.0477701663 6.31063730422 1.0
#aacgm to geo
glat,glon,r = aacgm.aacgmConv(52.7,6.6,300.,2000,1)
print glat, glon, r
41.4692793038 -71.090203565 1.0
models.aacgm.aacgmConvArr(lat,lon,alt,flg)
convert between geographic coords and aacgm (array form)
Input arguments:
Outputs:
#geo to aacgm
olat,olon,r = aacgm.aacgmConvArr([10.,20.,30.,40.],[80.,90.,100.,110.],[100.,150.,200.,250.],2000,0)
print olat
print olon
print r
[7.4318388600758, 15.446694598871858, 25.709559426821173, 35.92397286391025] [151.9994260976915, 162.04520124338288, 172.05073106755245, -177.5742593251219] [1.0, 1.0, 1.0, 1.0]
models.aacgm.mltFromEpoch(epoch,mlon)
calculate magnetic local time from epoch time and mag lon
Input arguments:
Outputs:
import datetime as dt
myDate = dt.datetime(2012,7,10)
epoch = utils.timeUtils.datetimeToEpoch(myDate)
mlt = aacgm.mltFromEpoch(epoch,52.7)
print mlt
22.864638597
models.aacgm.mltFromYmdhms(yr,mo,dy,hr,mt,sc,mlon)
calculate magnetic local time from year, month, day, hour, minute, second and mag lon
Input arguments:
Outputs:
mlt = aacgm.mltFromYmdhms(2012,7,10,0,0,0,52.7)
print mlt
22.864638597
models.aacgm.mltFromYrsec(yr,yrsec,mlon)
calculate magnetic local time from seconds elapsed in the year and mag lon
Input arguments:
Outputs:
yrsec = int(utils.timeUtils.datetimeToEpoch(dt.datetime(2012,7,10)) - utils.timeUtils.datetimeToEpoch(dt.datetime(2012,1,1)))
print yrsec
mlt = aacgm.mltFromYrsec(2013,yrsec,52.7)
print mlt
16502400 22.8620094835