In [1]:
############################################
# This code adds davitpy to your python path
# Eventually, this won't be necessary
import sys
sys.path.append('/davitpy')
############################################
In [2]:
from datetime import datetime as dt
from models import *
import utils
In [3]:
# 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)
In [11]:
# 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')
In [12]:
# 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
In [28]:
# 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)
In [29]:
# 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})
In [30]:
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()
In [31]:
# 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)
In [32]:
# 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])
Out[32]:
In [34]:
# 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.
In [35]:
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]
In [36]:
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()
Input arguments:
Output argument:
In [39]:
w = hwm.hwm07(11001, 0., 200., 40., -80., 0, 0, 0, [0, 0])
In [40]:
print w
models.aacgm.aacgmConv(lat,lon,alt,flg)
convert between geographic coords and aacgm
Input arguments:
Outputs:
In [41]:
#geo to aacgm
glat,glon,r = aacgm.aacgmConv(42.0,-71.4,300.,2000,0)
print glat, glon, r
In [42]:
#aacgm to geo
glat,glon,r = aacgm.aacgmConv(52.7,6.6,300.,2000,1)
print glat, glon, r
models.aacgm.aacgmConvArr(lat,lon,alt,flg)
convert between geographic coords and aacgm (array form)
Input arguments:
Outputs:
In [43]:
#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
models.aacgm.mltFromEpoch(epoch,mlon)
calculate magnetic local time from epoch time and mag lon
Input arguments:
Outputs:
In [44]:
import datetime as dt
myDate = dt.datetime(2012,7,10)
epoch = utils.timeUtils.datetimeToEpoch(myDate)
mlt = aacgm.mltFromEpoch(epoch,52.7)
print mlt
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:
In [45]:
mlt = aacgm.mltFromYmdhms(2012,7,10,0,0,0,52.7)
print mlt
models.aacgm.mltFromYrsec(yr,yrsec,mlon)
calculate magnetic local time from seconds elapsed in the year and mag lon
Input arguments:
Outputs:
In [46]:
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
In [ ]: