## Problems in Chapter 2

### Problem 1a: write a script to convert declination, inclination, intensity data to North, East and Down. First we need to import numpy, the module with lots of math functions and pandas with nice data manipulation functions



In [2]:

import numpy as np
import pandas as pd



Let's write a little function to do the conversion.



In [3]:

def dir2cart(data):
X=ints*np.cos(decs)*np.cos(incs)
Y=ints*np.sin(decs)*np.cos(incs)
Z=ints*np.sin(incs)
cart=np.array([X,Y,Z]).transpose()
return cart



Now let's read in a data file with some geomagnetic field vectors in it.



In [4]:

# read in the data and transpose it to rows of dec, inc, int
print dir2cart(data)




[[ 21352.55524831   2093.63634727  23332.08409238]
[  4504.44337072   -259.7245706   -1225.86288284]
[ 23546.1300489    3141.72451736  33426.255268  ]
[ 14629.0911691    1022.96570709  21021.51776849]
[ 23150.99484809   2965.71083857  30861.24994328]
[ 14767.09147922   2127.97038951    651.40495181]
[ 18929.94924879   -231.28446662  25961.37752135]
[  8342.98700429    759.27129675   8495.26107758]
[ 10858.76521357   2606.95887762  19818.79867013]
[ 30243.76260383   1532.08292009  41375.84902637]]



### Problem 1b: Read in locations from 10 random spots on Earth and calculate the IGRF vectors at each place.

First we have to understand how the function pmag.get_unf() works. To do this, we need to tell the notebook where the pmag module lives, import it and print out the doc string for get_unf():



In [5]:

import sys
sys.path.insert(0,'/Users/ltauxe/PmagPy')
import pmag
print pmag.get_unf.__doc__




Called with get_unf(N).
subroutine to retrieve N uniformly distributed directions
using the way described in Fisher et al. (1987).



Now we can use that function to generate a list of random points on the Earth's surface.



In [6]:

places=pmag.get_unf(10)
print places




[[ 105.31380297  -42.28827184]
[ 104.43893471  -16.26947532]
[ 225.43412517  -23.91520594]
[ 319.15116615   -2.68760951]
[ 354.27977495  -47.15291859]
[  86.28188506  -75.6932067 ]
[ 114.57001234  -38.57553035]
[ 312.18516165  -33.87893672]
[ 336.91620331  -27.52984633]
[ 313.74866005   66.61986125]]



Now let's find out about ipmag.igrf



In [7]:

import ipmag
print ipmag.igrf.__doc__




prints out Declination, Inclination, Intensity data from an input with format: Date, Altitude, Latitude, Longitude



And now we can ship the data in places to doigrf.



In [8]:

for place in places:
print ipmag.igrf([2006,0,place[1],place[0]])




[   338.62643377    -75.45243189  61051.29186633]
[   338.62643377    -75.45243189  61051.29186633]
[  3.59061529e+02  -4.93586702e+01   5.01630409e+04]
[  3.59061529e+02  -4.93586702e+01   5.01630409e+04]
[  1.47198793e+01  -3.73405456e+01   3.63248144e+04]
[  1.47198793e+01  -3.73405456e+01   3.63248144e+04]
[  3.38838448e+02  -9.29221666e+00   2.63942741e+04]
[  3.38838448e+02  -9.29221666e+00   2.63942741e+04]
[   339.69790061    -62.08266839  25520.99603321]
[   339.69790061    -62.08266839  25520.99603321]
[   265.30243821    -73.81812537  56762.24184327]
[   265.30243821    -73.81812537  56762.24184327]
[   353.14746029    -72.66799533  61283.79766406]
[   353.14746029    -72.66799533  61283.79766406]
[   344.31586161    -42.73959706  23232.50846823]
[   344.31586161    -42.73959706  23232.50846823]
[   334.5593929     -54.80640697  24901.30712016]
[   334.5593929     -54.80640697  24901.30712016]
[   328.06679306     79.2464479   54944.13402417]
[   328.06679306     79.2464479   54944.13402417]



## Problem 1c: Take the output from 1b and call dir2cart''.



In [9]:

data=[] # make a blank list
for place in places:
Dir=ipmag.igrf([2006,0,place[1],place[0]])
data.append(Dir) # append to the data list
data=np.array(data).transpose() # dir2cart takes arrays of data
print dir2cart(data)




[   338.62643377    -75.45243189  61051.29186633]
[  3.59061529e+02  -4.93586702e+01   5.01630409e+04]
[  1.47198793e+01  -3.73405456e+01   3.63248144e+04]
[  3.38838448e+02  -9.29221666e+00   2.63942741e+04]
[   339.69790061    -62.08266839  25520.99603321]
[   265.30243821    -73.81812537  56762.24184327]
[   353.14746029    -72.66799533  61283.79766406]
[   344.31586161    -42.73959706  23232.50846823]
[   334.5593929     -54.80640697  24901.30712016]
[   328.06679306     79.2464479   54944.13402417]
[[ 14280.40373426  -5588.83019226 -59093.95303089]
[ 32667.89651857   -535.12897139 -38063.79966558]
[ 27931.99561184   7338.18109798 -22032.85862437]
[ 24291.41065134  -9403.26937449  -4261.87710302]
[ 11206.54886205  -4145.89511602 -22550.96571805]
[ -1295.5074958  -15765.7784904  -54513.42942793]
[ 18126.52425584  -2178.31720519 -58501.17869484]
[ 16427.70320872  -4612.71237644 -15767.14627143]
[ 12959.97921918  -6165.11041094 -20349.58104118]
[  8700.30268401  -5422.46006446  53979.25086292]]



## Problem 2: Take the output from Problem 1c and plot as an equal area projection (first by hand and then with ipmag functions). The ipmag functions call pmagplotlib and use matplotlib, so these will have to be imported as well.



In [10]:

import pmagplotlib
import matplotlib.pyplot as plt
%matplotlib inline




In [11]:

ipmag.plot_net(1) # make an equal angle net
ipmag.plot_di(data[0],data[1]) # put on the dots






### Problem 3: Use the dipole formula ($\tan (I) = 2 \tan (\lambda)$ where $I$ is inclination and $\lambda$ is latitude and calculate the GAD field at 36 $^{\circ}$N. Note that declination is always zero for a GAD field.



In [12]:

inc = np.arctan(2.*np.tan(lat)) /deg2rad # and back!
print '%7.1f'%(inc) # and print it out




55.5



Let's use the pmag function dia_vgp. First let's figure out what it does:



In [13]:

print pmag.dia_vgp.__doc__




converts declination, inclination, alpha95 to VGP, dp, dm
takes input as (Decs, Incs, a95, Site latitudes, Site Longitudes).
These can be lists or individual values.



Now we can use it to convert our directions to VGPs. Note that alpha95 is require but is not given here so we supply a zero in its place. Note also that westward longitudes are indicated by minus signs



In [14]:

vgp_lat,vgp_lon,dp,dp= pmag.dia_vgp(345,47,0.,36,-112)
print '%7.1f %7.1f'%(vgp_lat,vgp_lon)




130.6    75.1