(skip to the R code)
L Shells: the toric strata of the magnetosphere
In radiation belt physics, we often want to compute what's called the “L value,” or the “L shell.” The L value, or just L for short, refers to how far away a magnetic field line intersects Earth's equatorial plane. It's just a number, given in units of Earth radii, that gives magnetospheric scientists an idea of what doughnut-like shell in the near-Earth magnetosphere some phenomenon under investigation roughly takes place in. The L value is analagous to how an atmospheric scientist might refer to altitude, which is also a single number, to speak of an entire spherical shell where certain atmospheric activity takes place. Altitude corresponds to spherical shells because atmospheric behavior is largely organized in a radial, gravitational sense. L-shell corresponds to doughnut-like, tubular shells because magnetospheric behavior is largely organized by the geomagnetic field.
If I'm standing on a ledge and drop a rock, neglecting any winds strong enough to nudge the rock from its dominantly radial descent, it will rapidly fall earthward (“straight down”), predictably making contact with the ground according to some basic equations describing projectile motion. Mass behavior like this in the Earth's gravitational field is fairly intuitive and familiar. In a similar way, the motions and dynamics of the gaseous fluids that make up the lower strata of the atmosphere (the troposphere, stratosphere, and mesosphere) are largely organized by Earth's gravity (e.g., gravity waves). One can refer to vast regions of the atmosphere by referencing just the altitude. This single number encodes a host of characteristics (density, temperature, etc) that roughly describe an entire spherically-shaped shellular region of the atmosphere.
In the magnetosphere, instead of a gas, the region is dominated by a rarefied plasma—that is, the magnetosphere is mostly empty space, sparsely populated by charged particles (electrons, protons, oxygen ions, etc). Although these charged particles feel the radial pull of Earth's gravity due to their mass, they more strongly respond to the geomagnetic field due to their electrically charged nature. This ultimately amounts to mass organization that isn't approximated by spherical shells, but by torii (“doughnuts”), at least in the near-Earth magnetosphere (going further out, the shells are more-or-less toric, homeomorphically speaking). So, instead of using altitude to refer to spherical shells of gas organized by the 'geogravitic' field, magnetospheric scientists use the L value to refer to dougnut-like shells (L shells) of plasma organized by the geomagnetic field.
The L of Ground Locations
Given a ground location of a magnetometer, one can roughly assign an L shell to it that is useful for many practical purposes. In the real world, given the near-Earth space weather conditions, the L value associated with a ground location can vary, but not by so much that it isn't worth coming up with some 'nominal' L values for particular ground locations (that is, we can assign L values to ground locations that are true given average/typical conditions). Having nominal L values, we then have a good idea concerning which region of the magnetosphere corresponds to a particular magnetometer location (e.g., our magnetometers in Antarctica). This is helpful, for example, when trying to study magnetospheric phenomena using both ground-based instrumentation and instrumentation aboard spacecraft, like the Van Allen Probes (a.k.a. the “Radiation Belt Storm Probes”), which is exactly what the Terrestrial Lab at NJIT is doing.
It might be asked if a ground location's nominal L value changes over time as the geomagnetic field itself slowly shifts over the years (secular variation). Strictly, the answer to this question is “yes,” but for back-of-the-envelope purposes you wouldn't be remiss to assume the answer is “no.” For example, using VITMO (more on VITMO below), one can quickly compute how the nominal L-value for the ground location at (alt, glat, glon) = (0, 45.0, 270.0) changed between the years 1900 and 2000. The change in L was by about ~1.2%, slowly drifting from 4.19 to 4.24. So, approximating a ground location's nominal L value as a constant over 10 year intervals isn't too radical an assumption.
Computing L
Actually computiing the nominal L shell associated with a ground magnetometer is a fairly trivial computation (well, not exactly “trivial”; more on this below). Often, the challenge is obtaining data from a given magnetometer for a specific time of interest: the data might not exist or you simply might not have access to it. Thanks to efforts like Intermagnet (and, of course, having our own magnetometers in NJ and Antarctica), we are able to study and analyze geomagnetic data from many sources located all around the world in “near” real time (take note that “near” can sometimes mean a few weeks of lag!).
To compute the L value, one must first compute the radial distance of the ground location from the Earth's center (in units of Earth radii) and the location's magnetic field strength, B. For back-of-the-envelope estimates, you can safely assume any ground location is \( r=1[{R}_{E}] \). For B, you can do a number of things—for one, you can take lots of measurements over a couple of days and compute an average value to represents the ground location. Or, alternatively, if you know the ground location's geomagnetic latitude in ordinary geomagnetic dipole coordinates, then to find B, you can assume the Earth's field is basically dipolar and plug it into
If you don't know the site's geomagnetic latitude in ordinary dipole coordinates, then you can use a service like the one found at the World Data Center for Geomagnetism, which only demands you plug in the ground location's geographic coordinates. Note that the coordinates this calculator returns are different than simple dipole coordinates; they derive from a more advanced geomagnetic field model known as the International Geomagnetic Reference Field (IGRF). If you don't even know the geographic coordinates, fret not: you can get a good approximation by using Google. For example, type “geographic coordinates of jenny jump” into the Google search. With the geographic coordinates in hand, go ahead and plug them into the WDC coordinate calculator to get the IGRF latitude, which you can then plug into the above dipole equation (it should work ok!) to get an approximate, nominal value for B. For a better nominal value of B, you can employ more advanced geomagnetic field models to obtain its value, e.g., one of the Definitive/International Geomagnetic Reference Field (IGRF) models.
In relation to the L coordinate, there is also an equation for B (e.g., McIlwain, 1961; Walt, 1994,Ch. 4):
MatLab>> syms B B0 L r |
MatLab>> solve(B-(B0*L^3/r^3)*sqrt(4-3*r/L), L) |
ans =RootOf(X35^6 - (3*X35^5*r)/4 - (B^2*r^6)/(4*B0^2), X35) |
Don't ask me what X35 is. Just think of it as MatLab's way of saying, “I don't know!” Although algebraically hopeless, you can do it numerically and many computer programs already exist that will do it for you, free of charge. One such service is VITMO (Virtual Ionosphere-Thermosphere-Mesosphere Observatory). Technically, with VITMO you can do the entire process at once: plug in the geographic coordinates and output the L. (We could have skipped straight to this, but at least now you've seen some equations and understand a little about what VITMO is doing.)
Automating the Process with R
Using the programming language R (or the language of one's choice), the process can be automated if the geomagnetic coordinates and L shells of multiple ground locations must be computed. The list of magnetometers, along with their associated geographic coordinates, that we will use were obtained from INTERMAGNET. VITMO is accessed automatically via the curl routine available on UNIX (and related) systems.
Get the Ground Location Data
To make the process reproducible, I created a cleaned-up csv file (html view) of Intermagnet magnetometer location info, which can be downloaded and loaded into R fairly easily:
download.file("https://docs.google.com/spreadsheet/pub?key=0AgIi2HrVxV_ldFM2bU1TNzc3TGthS3o2UUZNS0h0Qmc&single=true&gid=0&output=csv",
destfile = "intermagnet_mags.csv", method = "curl")
list.of.mags = read.csv("intermagnet_mags.csv")
One has to do some backflips to make that file on Google Drive downloadable.
It's a little bit easier to access data sets stored in a Dropbox account. For example, I obtained a larger list of magnetometers from the UK Solar System Data Centre, cleaned it up for use in R, and stored it in my Dropbox. Using R, it can be accessed like so:
# csv file
download.file("https://dl.dropbox.com/s/ctj9cy5f2d4k3ee/ground_mags.csv", destfile = "list_of_mags.csv",
method = "curl")
# rda file
download.file("https://dl.dropbox.com/s/39vm23l8t9d6t05/ground_mags.rda", destfile = "list_of_mags.rda",
method = "curl")
In our conjunction studies between the VAP spacecraft and ground stations, we currently only focus on the 132 Intermagnet magnetometers and our own magnetometers at Jenny Jump State Forest, NJ, and the PENGUIn AGOs in Antarctica (however, most of the AGOs are much too high in geomagnetic latitude to be of service in these studies). The exclusion of magnetometers from the larger list, which has a total of 311 magnetometers (many overlapping with Intermagnet) is simply due to data availability and accessibility.
To see the Intermagnet magnetometers on a map, simply load the maps package into R and enter map('world') into the R console to bring up the world map. On this map, latitude is on its usual range, [-90,90]. Note, however, that instead of the usual azimuthal range of [0,360] in eastward degrees from the prime meridian, longitude takes on the range [-180,180], where a negative number indicates westward degrees from the prime meridian while a positive indicates eastward degrees.
library("maps")
# Plot the world map
map("world")
# Simplify var names
lat = list.of.mags$glat
lon = list.of.mags$glon
# Convert longitude range: [0,360) -> (-180,180]
long = lon
long[long > 180] = long[long > 180] - 360
# Plot magnetometer locations on world map
points(long, lat, pch = 19, col = "red")
VITMO, L, and R
The version of VITMO we will call upon here only computes outputs for ground locations with geographic latitudes whose absolute values are greater than \( \pm {20}^{\circ} \). Perhaps in a future post, we will find a way to include these dismissed ground locations using another tool, but for the current purposes of instruction, it's not a dire necessity. The excluded ground locations can be assumed to have an L of about 1 since most of them are on L shells way below VAP orbital altitudes. So the list of mags we want to compute geomagnetic paramaters for are those with latitudes whose magnitudes are greater than 20-deg.
One more thing: In the final output, besides L, other parameters associated with the ground location might be nice to include for later use, such as the geographic longitude and latitude of the ground location and it conjugate point, the corresponding corrected geomagnetic (CGM) longitudes and latitudes of both, and local magnetic midnight (given in Universal Time). More information on CGM coordinates and “magnetic midnight” can be found in my CGM post.
Without further ado… The R Code:
# Select year
year = "2013"
station.info = data.frame(matrix(nrow = 0, ncol = 10))
names(station.info) = c("glat", "glon", "cgmlat", "cgmlon", "magMid", "cp.glat",
"cp.glon", "cp.cgmlat", "cp.cgmglon", "L")
# Extract stations with abs(geo-lat) > 20
lat = lat[abs(lat) > 20]
lon = lon[abs(lat) > 20]
# Ask VITMO about L :-p
for (i in 1:length(lat)) {
dataSpecs = paste0("\"model=cgm&year=", year, "&height=0.&geo_flag=1&latitude=",
lat[i], "&longitude=", lon[i], "&profile=1&start=0.&stop=0.&step=50.&format=0&vars=02&vars=03&vars=04&vars=05&vars=11&vars=12&vars=13&vars=14&vars=15&vars=42&linestyle=solid&charsize=&symbol=0&symsize=&imagex=640&imagey=480\"")
query = paste0(dataSpecs, " http://omniweb.gsfc.nasa.gov/cgi/vitmo/vitmo_model.cgi > temp-4-L-shell-blog.txt")
system(paste0("curl -d ", query))
station.info[i, ] = as.numeric(strsplit(readLines("temp-4-L-shell-blog.txt")[25],
" +")[[1]][2:11])
}
write.table(station.info, file = "station-info.txt")
head(station.info)
## glat glon cgmlat cgmlon magMid cp.glat cp.glon cp.cgmlat cp.cgmglon
## 1 43.20 76.90 39.13 150.01 18.69 -27.80 82.42 -39.13 150.01
## 2 68.36 38.77 64.89 117.23 20.50 -56.69 71.53 -64.89 117.23
## 3 -65.25 72.87 -71.43 104.89 21.80 74.14 17.51 71.43 104.89
## 4 82.50 18.82 79.25 122.97 20.21 -68.38 95.83 -79.25 122.97
## 5 -37.80 295.75 -25.81 5.44 4.11 15.95 289.15 25.81 5.44
## 6 42.38 297.65 49.32 19.67 3.87 -61.82 314.20 -49.32 19.67
## L
## 1 1.66
## 2 5.55
## 3 9.86
## 4 999.99
## 5 1.23
## 6 2.35
Feel free to take a look at the data frame we just created. (Note that an L value of 999.99 basically means its not realistic to even assign a value because, e.g., the ground location is at such a high geomagnetic latitude that it often resides within the open-closed boundary.)
This document was created using RMarkdown in RStudio.
Last edited 2013-Mar-20 by Kevin Urban (Terrestrial Lab)
