A.3 Useful tricks

In this section, we have assembled assorted mathematical and plotting techniques that come in handy through out this book.

A.3.1 Spherical trigonometry

Spherical trigonometry has widespread applications throughout the book. It is used in the transformations of observed directions to virtual poles (Chapter  2) and transformation of coordinate sytems, to name a few. Here we summarize the two most useful relationships: the Law of Sines and the Law of Cosines.


EPSfiles/strig.eps

Figure A.1: Rules of spherical trigonometry. a,b,c are all great circle tracks on a sphere which form a triangle with apices A,B,C. The lengths of a,b,c on a unit sphere are equal to the angles subtended by radii that intersect the globe at the apices, as shown in the inset. α,β,γ are the angles between the great circles.


In Figure A.1, α,β and γ are the angles between the great circles labelled a, b, and c. On a unit sphere, a,b and c are also the angles subtended by radii that intersect the globe at the apices A, B, and C (see inset on Figure A.1). Two formulae from spherical trigonometry come in handy in paleomagnetism, the Law of Sines:

sin α    sin β    sinγ
-----= ----- = -----,
sin a    sin b    sin c
(A.10)

and the Law of Cosines:

cosa =  cosbcosc+  sin bsin ccosα.
(A.11)

A.3.2 Vector addition


EPSfiles/vectors.eps

Figure A.2: Vectors A and B, their components Ax,y, Bx,y and the angles between them and the X axis, α and β. The angle between the two vectors is α -β = Δ. Unit vectors in the directions of the axes are ˆx and ŷ respectively.


To add the two vectors (see Figure A.2) A and B, we break each vector into components Ax,y and Bx,y. For example, Ax = |A|cosα,Ay = |A|sinα where |A| is the length of the vector A. The components of the resultant vector C are: Cx = Ax + Bx,Cy = Ay + By. These can be converted back to polar coordinates of magnitude and angles if desired, whereby:

      ∘ --------             C
|C| =   C2x + C2yand γ = cos-1 -x-.
                             |C |

A.3.3 Vector subtraction

To subtract two vectors, compute the components as in addition, but the components of the vector difference C are: Cx = Ax - Bx,Cy = Ay - By.

A.3.4 Vector multiplication

There are two ways to multiply vectors. The first is the dot product whereby A B = AxBx + AyBy. This is a scalar and is actually the cosine of the angle between the two vectors if the A and B are taken as unit vectors (assume a magnitude of unity in the component calculation.


EPSfiles/cross.eps

Figure A.3: Illustration of cross product of vectors A and B separated by angle θ to get the orthogonal vector C.


The other way to perform vector multiplication is the cross product (see Figure A.3), which produces a vector orthogonal to both A and B and whose components are given by:

       |             |
       || xˆ   ˆy    ˆz ||
C = det| Ax  Ay   Az |.
       || B   B    B  ||
          x    y    z

To calculate the determinant, we follow these rules:

Cx = AyBz  - AzBy, Cy = AzBx  - AxBz, Cz = AxBy - AyBx.

or

Ci = AjBk -  AkBj    i ⁄= j ⁄= k.

A.3.5 Tricks with tensors

Vectors belong to a more general concept called tensors. While a vector describes a magnitude of something in a given direction, tensors allow calculation of magnitudes as a function of orientation. Velocity is a vector relating speed to direction, but speed may change depending on direction, so we might need a tensor to calculate speed as a function of direction. Many properties in Earth science require tensors, like the indicatrix in mineralogy which relates the speed of light to crystallographic direction, or the relationship between stress and strain. Tensors in paleomagnetism are used, for example, to transform coordinate systems and to characterize the anisotropy of magnetic properties such as susceptibility. We will cover transformation of coordinate systems in the following.

Direction cosines

We use direction cosines in paleomagnetism in a variety of applications, from mineralogy to transformation from specimen to geographic or stratigraphic coordinate systems. Direction cosines are the cosines of the angles between different axes in given coordinate systems, here X and Xrespectively (see, e.g., Figure A.4a.) The direction cosine a12 is the cosine of the angle between the X1 and the X2, α12 axes. We can define four of these direction cosines to fully describe the relationship between the two coordinate systems:

a11 = cosα11,a21 = cosα21,
a12 = cosα12,a22 = cosα22.

The first subscript always refers to the X system and the second refers to the X.


EPSfiles/dircosines.eps

Figure A.4: Definition of direction cosines in two dimensions. a) Definition of vector in one set of coordinates, x1,x2. b) Definition of angles relating X axes to X.


Changing coordinate systems

One application of using direction cosines is the transformation of coordinates systems from one set (X) to a new set X. To find new coordinates x1,x2,.. from the old (x1,x2,...), we have:

x′ = a  x + a  x ,
 1′    11 1    12 2
x2 = a21x2 + a22x2.

In three dimensions we have:

x ′= a  x  + a  x + a  x ,
  1′   11 1    12 2   13 3
x 2′ = a21x2 + a22x2 + a23x3,
x 3 = a31x3 + a32x2 + a33x3,

which can also be written as:

(    )    (              )(    )
  x ′1       a11  a12  a13     x1
( x ′) =  ( a21  a22  a23 )(  x2),
  x 2′       a    a   a       x
    3        31   32    33      3
(A.12)

with a short cut notation as: xi = aijxj. However we write this, it means that for each axis i, just sum through the j’s for all the dimensions. The matrix aij is an example of a 3 x 3 tensor and equations of the form Ai = BijCj relating two vectors with a tensor will be used throughout the book. A more common notation is with bold-faced variables which indicate vectors or tensors, e.g., A = B C.


EPSfiles/transform.eps

Figure A.5: a) Sample coordinate system. b) Trigonometric relations between two cartesian coordinate systems, Xi and Xi. λ,ϕ,ψ are all known and the angles between the various axes can be calculated using spherical trigonometry. For example, the angle α between X1 and X1forms one side of the triangle shown by dash-dot lines. Thus, cosα = cosλcosϕ + sinλsinϕcosψ. [Figure from Tauxe, 1998.]


Now we would like to apply this to changing coordinate systems for a paleomagnetic specimen in the most general case. The specimen coordinate system is defined by a right-hand rule where the thumb (X1) is directed parallel to an arrow marked on the sample, the index finger (X2) is in the same plane but at right angles and clockwise to X1 and the middle finger (X3) is perpendicular to the other two (Figure A.5a). The transformation of coordinates (xi) from the Xi axes to the coordinates in the desired Xcoordinate system requires the determination of the direction cosines as described in Appendix A.3.5. The various aij can be calculated using spherical trigonometry as in Appendix A.3.1. For example, a11 for the general case depicted in Figure A.5 is cosα, which is given by the Law of Cosines (see Appendix A.3.1) by using appropriate values, or:

cos α = cosλ cosϕ + sin λsinϕ cosψ.

The other aij can be calculated in a similar manner. In the case of most coordinate system rotations used in paleomagnetism, X2 is in the same plane as X1 and X2 (and is horizontal) so ψ = 90. This problem is much simpler. The directions cosines for the case where ψ = 90 are:

    (                                )
       cosλcosϕ   - sin ϕ  - sin λcos ϕ
a = (  cosλsinϕ    cosϕ   - sin λ sinϕ ) .
         sin λ       0        cos λ
(A.13)

The new coordinates can be obtained from Equation A.12, as follows:

 x′=  a x  + a  x +  a x
  1′    11  1   12 2    13  3
 x2′ = a21x1 + a22x2 + a23x3
x 3 = a31x1 + a32x2 + a33x3.
(A.14)

The declination and inclination can be calculated by inserting these values in the equations in Chapter  2.

Method for rotating points on a globe using finite rotation poles

Given the coordinates of the point on the globe Pp with latitude λp, longitude ϕp the finite rotation pole Pf with latitude λf, longitude ϕf, the way to transform coordinates is as follows (you should also review Appendix A.3.5).

  1. Convert the latitudes and longitudes to cartesian coordinates by:
    P1 =  cosϕcosλ, P2 = sin ϕ cos λ,P3 = sinλ
    where P is the point of interest.
  2. Set up the rotation matrix R as:

      R   = P  P  (1 - cosΩ )+ cos Ω
   11    f1 f1
R12 = Pf1Pf2 (1-  cosΩ) - Pf3sinΩ
R13 = Pf1Pf3 (1-  cosΩ) + Pf2sinΩ
R21 = Pf2Pf1 (1-  cosΩ) + Pf3sinΩ
  R22 = Pf2Pf2(1 - cosΩ )+ cos Ω  .
R23 = Pf2Pf3 (1-  cosΩ) - Pf1sinΩ
R31 = Pf3Pf1 (1-  cosΩ) - Pf2sinΩ
R32 = Pf3Pf2 (1-  cosΩ) + Pf1sinΩ
  R   = P  P  (1 - cosΩ )+ cos Ω
   33    f3 f3
  3. The coordinates of the transformed pole (Pt) are:
    Pt1 = R11Pp1 + R12Pp2 + R13Pp3
Pt2 = R21Pp1 + R22Pp2 + R23Pp3,
Pt3 = R31Pp1 + R32Pp2 + R33Pp3
    which can be converted back into latitude and longitude in the usual way (see Chapter  2).

Table A.4: Finite rotations for selected Gondwana continents.
Rotates continent to South African fixed coordinates.
AUS: Australia, ANT: East Antarctica; IND: India; SAM:
South American Craton. [Rotations of Torsvik et al. (2008)
- see for additional data.]

Table A.5: Finite rotations for selected Laurentian continents. Rotates continent to South African fixed coordinates. EUR: Europe; NAM: North America; GRN: Greenland. [Rotations of Torsvik et al. (2008) - see for additional data.]










Age EUR NAM GRN










Ma λϕ Ωλϕ Ωλϕ Ω










EUR NAM GRN
5 17.9-27.10.6 80.922.8 1.3 80.922.8 1.3
1018.4-26.31.2 80.922.9 2.6 80.922.9 2.6
1518.9-24.61.8 80.923.2 4.1 80.923.2 4.1
2017.2-22.72.4 80.624.4 5.5 80.624.4 5.5
2520.7-19 3 79.528.1 6.8 79.528.1 6.8
3024.9-19.54.3 77.312.5 8.6 77.312.5 8.6
3527.2-19.35.8 75.43.5 10.574.87.2 10.2
4028.7-18.57.5 74.5-1.1 12.672.69.5 11.5
4530.3-18.29 74.3-4.3 14.671.411.4 12.7
5030.8-16.710 75.9-3.5 16.271 20.7 14.2
5532.7-15.411.379.84.1 17.671.829.6 16.8
6034.8-15.712.681.65.1 19.171.930.5 17.5
6536 -15.813.682.63.2 20.771.332.9 17.6
7035.4-16.114.981.6-6.5 22.469.829 17.9
7535.5-15.715.580.4-13.1 24.669 26.6 18.5
8036.1-15.216.978.2-18.8 27.567.621 19.8
8537 -14.218.876.2-21.3 30.566.316.4 21.5
9039.6-13.721.974.6-23 33.865.911.5 24.2
9539.8-13.725.272 -24.7 36.964.25.5 26.9
10040.2-12.528.570 -24 40.262.72.830.1
10541.6-11.231.769.1-23.343.662.41.633.3
11042.6-9.8 34.568.3-22.647 62.10.936.5
11543.4-8.5 37.367.6-21.850.461.80.539.7
12044.5-6.9 40.367.1-20.453.961.80.843.1
12545.3-6.3 42 67 -19.755.661.91 44.9
13045.9-5.7 43 67 -19.156.762.21.346
13546.6-5.3 44 67.1-18.757.962.41.647.1
14047.3-4.9 45.267.2-18.459.262.71.648.4
14547.8-4.8 46.467.1-18.360.562.91.349.7
15048.6-4 47.967.3-17.662.263.21.851.4
15549.8-2.2 50 67.6-15.564.663.73.653.8
16050.6-1.2 52.167.6-14.566.864.14.256
16551.4-0.3 54.267.7-13.669.164.44.858.3
17052.10.6 56.367.8-12.871.464.75.360.6
17552.91.9 59.667.7-11.574.864.86 64.1
18053 2 60 67.7-11.575.364.96 64.5
18553 2 60.467.7-11.575.764.95.964.9
19053.12.1 60.867.7-11.576.165 5.965.4
19553.22.2 61.167.7-11.576.665 5.865.8
20053.32.261.567.7-11.577 65.15.866.2
20553.22.659.767.7-11.577.465.15.766.7
21053.12.957.867.7-11.577.965.25.767.1
21553.13.355.967.7-11.578.365.25.667.5
22052.93.653.667.7-11.578.365.25.667.5
22552.74 51.467.7-11.578.365.25.667.5
23052.44.449.167.7-11.578.365.25.667.5
23552.24.846.867.7-11.578.365.25.667.5
24051.95.344.567.7-11.578.365.25.667.5
24551.95.344.567.7-11.578.365.25.667.5
25051.95.344.567.7-11.578.365.25.667.5
25551.95.344.567.7-11.578.365.25.667.5
26051.95.344.567.7-11.578.365.25.667.5
26551.95.344.567.7-11.578.365.25.667.5
27051.95.344.567.7-11.578.365.25.667.5
27551.95.344.567.7-11.578.365.25.667.5
28051.95.344.567.7-11.578.365.25.667.5
28551.95.344.567.7-11.578.365.25.667.5
29051.95.344.567.7-11.578.365.25.667.5
29551.95.344.567.7-11.578.365.25.667.5
30051.95.344.567.7-11.578.365.25.667.5
30551.95.344.567.7-11.578.365.25.667.5
31051.95.344.567.7-11.578.365.25.667.5
31551.95.344.567.7-11.578.365.25.667.5
32051.95.344.567.7-11.578.365.25.667.5






















Table A.6: Finite rotations for South Africa to the paleomagnetic reference frame. All rotation pole latitudes are assumed to be zero. Ages are in Ma and rotation angles are in degree.[Rotations of Torsvik et al. (2008) - see for additional options.]












Ageϕ Ω Ageϕ Ω Ageϕ Ω Ageϕ Ω












5 56 2.2 85 138.119.3165 157 30.7245 137.436.5
10 57.6 2.5 90 142.919.6170 159.532.5250 143.139.6
15 53.9 2.5 95 144.720.5175 167.628.8255 145.440.4
20 66.5 3 100 144.320.8180 167.827.7260 145.641.8
25 75.5 4.7 105 150.822.3185 167.425.9265 144.841.9
30 84.1 6.8 110 160.226.9190 168.421.6270 141.647.1
35 95.8 7.9 115 169.232.1195 158.818.2275 140.346.8
40 98.8 8.7 120 170.335.6200 147.917.8280 138.251.1
45 107.59.2 125 171.336.2205 144.419.2285 138.651.6
50 110.910.3130 172.137.5210 137.420.7290 136.551.8
55 111.613.2135 170 39.4215 133.623.1295 135.852.8
60 115.713.9140 172.642.1220 129.926.4300 136.853.5
65 123.515.7145 163.140.8225 127.227.2305 136.955.4
70 127.817.5150 155.238.1230 128 29.4310 138.956.3
75 137.217.5155 155 34.8235 130 31.4315 139.959.5
80 140.319.2160 155 33.2240 133.635.3320 138.960.8













The orientation tensor and eigenvectors

The orientation tensor T (Scheidegger, 1965) (also known as the matrix of sums of squares and products), is extremely useful in paleomagnetism. This is found as follows:

  1. Convert the D, I, and M for a set of data points (e.g., a sequence of demagnetization data, or a set of geomagnetic vectors or unit vectors where M = 1) to corresponding xi values (see Chapter  2).
  2. Calculate the coordinates of the “center of mass” (x) of the data points:
          1  N∑              1 ∑N              1 ∑N
x1 =  N-(   x1i);  x2 = N- (   x2i);  x3 = N-(    x3i),
          1                 1                1
    (A.15)

    where N is the number of data points involved. Note that for unit vectors, the center of mass is the same as the Fisher mean (Chapter  11).

  3. Transform the origin of the data cluster to the center of mass:
    x′ = x  - x ;  x ′ = x  - x ;  x′ =  x  - x ,
 1i   1i    1    2i    2i   2    3i    3i   3
    (A.16)

    where xi are the transformed coordinates.

  4. The orientation matrix is defined as:
        ( ∑         ∑         ∑       )
      ∑  x′1ix ′1i  ∑  x′1ix2i  ∑  x′1ix′3i
T = ( ∑  x′1ix ′2i  ∑  x′2ix ′2i  ∑  x′2ix′3i) .
         x′1ix ′3i     x′2ix ′3i     x′3ix′3i
    (A.17)

T is a 3 x 3 matrix, where only six of the nine elements are independent. It is constructed in some coordinate system, such as the geographic or sample coordinate system. Usually, none of the six independent elements are zero. There exists, however, a coordinate system along which the “off-axis” terms are zero and the axes of this coordinate system are called the eigenvectors of the matrix. The three elements of T in the eigenvector coordinate system are called eigenvalues. In terms of linear algebra, this idea can be expressed as:

TV  = τV,
(A.18)

where V is the matrix containing three eigenvectors and τ is the diagonal matrix containing three eigenvalues. Equation A.18 is only true if:

det|T - τ| = 0.
(A.19)

If we expand equation A.19, we have a third degree polynomial whose roots (τ) are the eigenvalues:

(T11 - τ)[(T22 - τ)(T33 - τ)- T 223]-
T  [T  (T   - τ)- T  T  ]+ T  [T  T  - T  (T  - τ)] = 0.
 12  12  33        13 23    13  13 23   13  22

The three possible values of τ (τ123) can be found with iteration and determination. In practice, there are many programs for calculating τ. My personal favorite is the Numpy Module for Python (see many free websites, especially Scientific Python (SciPy) for hints). Please note that the conventions adopted here are to scale the τ’s such that they sum to one; the largest eigenvalue is termed τ1 and corresponds to the eigenvector V1.

Inserting the values for the transformed components calculated in equation  A.16 into T gives the covariance matrix for the demagnetization data. The direction of the axis associated with the greatest scatter in the data (the principal eigenvector V1) corresponds to a best-fit line through the data. This is usually taken to be the direction of the component in question. This direction also corresponds to the axis around which the “moment of inertia” is least. The eigenvalues of T are the variances associated with each eigenvector. Thus the standard deviations are σi = √ τi.

A.3.6 Upside down triangles,

Gradient

We often wish to differentiate a function along three orthogonal axes. For example, imagine we know the topography of a ski area (see Figure A.6). For every location (in say, X and Y coordinates), we know the height above sea level. This is a scalar function. Now imagine we want to build a ski resort, so we need to know the direction of steepest descent and the slope (red arrows in Figure A.6).


EPSfiles/ski.eps

Figure A.6: Illustration of the relationship between a vector field (direction and magnitude of steepest slope at every point, e.g., red arrows) and the scalar field (height) of a ski slope.


To convert the scalar field (height versus position) to a vector field (direction and magnitude of greatest slope) mathematically, we would simply differentiate the topography function. Let’s say we had a very weird two dimensional, sinusoidal topography such that z = f(x) = sinx with z the height and x is the distance from some marker. The slope in the x direction (ˆx), then would be ˆx d _ dxf(x). If f(x,y,z) were a three dimentional topography then the gradient of the topography function would be:

  ∂--     ∂--    -∂-
(xˆ∂x f + ˆy ∂yf + ˆz∂z f).

For short hand, we define a “vector differential operator” to be a vector whose components are

∇ = (ˆx ∂-,ˆy ∂-,ˆz-∂-).
       ∂x   ∂y  ∂z

This can also be written in polar coordinates:

     -∂- -∂------∂---
∇  = ∂r, r∂θ,r sin θ∂ϕ.


EPSfiles/div.eps

Figure A.7: Example of a vector field with a non-zero divergence.


Divergence

The divergence of a vector function (e.g. H) is written as:

∇ ⋅H.

The trick here is to treat as a vector and use the rules for dot products described in Appendix A.3.2. In cartesian coordinates, this is:

         ∂Hx-    ∂Hy-    ∂Hz-
∇ ⋅H  = ˆx ∂x  + ˆy ∂y  + ˆz ∂z  .

Like all dot products, the divergence of a vector function is a scalar.


EPSfiles/divzero.eps

Figure A.8: Example of a vector field with zero divergence.



EPSfiles/curl.eps

Figure A.9: Example of a vector field with non-zero curl.


The name divergence is well chosen because ∇⋅ H is a measure of how much the vector field “spreads out” (diverges) from the point in question. In fact, what divergence quantifies is the balance between vectors coming in to a particular region versus those that go out. The example in Figure A.7 depicts a vector function whereby the magnitude of the vector increases linearly with distance away from the central point. An example of such a function would be v(r) = r. The divergence of this function is:

        ∂
∇  ⋅v = --r = 1.
        ∂r
(a scalar). There are no arrows returning in to the dashed box, only vectors going out and the non-zero divergence quantifies this net flux out of the box.

Now consider Figure A.8, which depicts a vector function that is constant over space, i.e. v(r) = k. The divergence of this function is:

∇  ⋅v = ∂-k = 0.
        ∂r

The zero divergence means that for every vector leaving the box, there is an equal and opposite vector coming in. Put another way, no net flux results in a zero divergence. The fact that the divergence of the magnetic field is zero means that there are no point sources (monopoles), as opposed to electrical fields that have divergence related to the presence of electrons or protons.

Curl

The curl of the vector function B is defined as ∇× B. In cartesian coordinates we have

           ∂--     ∂--       -∂-     -∂-        ∂--     ∂--
∇ × B =  ˆx(∂yBz -  ∂zBy )+ ˆy(∂z Bx - ∂x Bz) + ˆz(∂xBy -  ∂yBx ).

Curl is a measure of how much the vector function “curls” around a given point. The function describing the velocity of water in a whirlpool has a significant curl, while that of a smoothly flowing stream does not.

Consider Figure A.9 which depicts a vector function v = -yˆx + xŷ. The curl of this function is:

            |            |
            || ˆx    ˆy   ˆz ||
∇ ×  v = det| ∂-  -∂   ∂-|,
            ||-∂xy  ∂xy   ∂z0 ||

or

ˆx ( ∂-0 - ∂-x )+ ˆy(-∂-0- ∂-(- y)) + ˆz(-∂-x-  ∂-(- y)).
   ∂y    ∂z       ∂x     ∂z          ∂x     ∂y
= 0xˆ+ 0ˆy + 2ˆz
So there is a positive curl in this function and the curl is a vector in the direction.

The magnetic field has a non-zero curl in the presence of currents or changing electric fields. In free space, away from currents (lightning!!), the magnetic field has zero curl.

A.3.7 The statistical bootstrap

Sometimes things just are not normal. Statistically that is. When you can not assume that your data follow some known distribution, like the normal distribution, or the Fisher distribution, what do you do? In this section, we outline a technique called the bootstrap, which allows us to make statistical inferences when parametric assumptions fail. The reader should also refer to Efron and Tibshirani (1993) for a more complete discussion.

In Figure A.10, we illustrate the essentials of the statistical bootstrap. We will develop the technique using data drawn from a normal distribution. First, we generate a synthetic data set by drawing 500 data points from a normal distribution with a mean x of 10 and a standard deviation σ of 2. The synthetic data are plotted as a histogram in Figure A.10a. In Figure A.10b we plot the data as a Q-Q plot (see Appendix B.1.5) against the zi expected for a normal distribution.


EPSfiles/bootstrap.eps

Figure A.10: Bootstrapping applied to a normal distribution. a) 500 data points are drawn from a Gaussian distribution with mean of 10 and a standard deviation of 2. b) Q-Q plot of data in a). The 95% confidence interval for the mean is given by Gauss statistics as 0.17. 10,000 new (para) data sets are generated by randomly drawing N data points from the original data set shown in a). c) A histogram of the means from all the para-data sets. 95% of the means fall within the interval 10.06-0.16+0.16, hence the bootstrap confidence interval is similar to that calculated with Gaussian statistics. [Figure from Tauxe, 1998.]


The data in Figure A.10a plot in a line on the Q-Q plot (Figure A.10b). The value for D is 0.0306. Because N = 500, the critical value of D, Dc at the 95% confidence level is 0.0396. Happily, our normal distribution simulation program has produced a set of 500 numbers for which the null hypothesis of a normal distribution has not been rejected. The mean of the synthetic dataset is about 10 and the standard deviation is 1.9. The usual Gaussian statistics allow us to estimate a 95% confidence interval for the mean as 1.96σ∕√ --
  N or 0.17.


EPSfiles/sundefs.eps

Figure A.11: Calculation of the azimuth of the shadow direction (β) relative to true North, using a sun compass. L is the site location (at λLL), S is the position on the Earth where the sun is directly overhead (λSS). [Figure from Tauxe, 1998.]


In order to estimate a confidence interval for the mean using the bootstrap, we first randomly draw a list of N data by selecting data points from the original data set. This list is called a pseudo-sample of the data. Some data points will be used more than once and others will not be used at all. We then calculate the mean of the pseudo-sample. We repeat the procedure of drawing pseudo-samples and calculating the mean many times (say 10,000 times). A histogram of the “bootstrapped” means is plotted in Figure A.10c. If these are sorted such that the first mean is the lowest and the last mean is the highest, the 95% of the means are between the 250th and the 9,750th mean. These therefore are the 95% confidence bounds because we are approximately 95% confident that the true mean lies between these limits. The 95% confidence interval calculated for the data in Figure A.10 by bootstrap is about 0.16 which is nearly the same as that calculated the Gaussian way. However, the bootstrap required orders of magnitude more calculations than the Gaussian method, hence it is ill-advised to perform a bootstrap calculation when a parametric one will do. Nonetheless, if the data are not Gaussian, the bootstrap provides a means of calculating confidence intervals when there is no quick and easy way. Furthermore, with a modern computer, the time required to calculate the bootstrap illustrated in Figure A.10 was virtually imperceptible.

A.3.8 Directions using a sun compass

In a sun compass problem, we have the direction of the sun’s shadow and an angle between that and the desired direction (α). The declination of the shadow itself is 180 from the direction toward the sun. In Figure A.11, the problem of calculating declination from sun compass information is set up as a spherical trigonometry problem, similar to those introduced in Chapter  2 and Appendix A.3.1. The declination of the shadow direction β, is given by 180 - β. We also know the latitude of the sampling location L (λL). We need to calculate the latitude of S (the point on the Earth’s surface where the sun is directly overhead), and the local hour angle H.

Knowing the time of observation (in Universal Time), the position of S (λs = δ,ϕs in Figure A.11) can be calculated with reasonable precision (to within 0.01) for the period of time between 1950 and 2050 using the procedure recommended in the 1996 Astronomical Almanac:

  1. First, calculate the Julian Day J. Then, calculate the fraction of the day in Universal Time U. Finally, calculate the parameter d which is the number of days from J2000 by:

    d = J - 2451545 + U.
  2. The mean longitude of the sun (ϕs), corrected for aberration, can be estimated in degrees by:

    ϕ  = 280.461 + 0.9856474d.
 s
  3. The mean anomaly g = 357.528 + 0.9856003d (in degrees).
  4. Put ϕs and g in the range 0 360.
  5. The longitude of the ecliptic is given by ϕE = ϕs + 1.915sing + 0.020sin2g (in degrees).
  6. The obliquity of the ecliptic is given by ϵ = 23.439 - 0.0000004d.
  7. Calculate the right ascension (A) by:

                              2
A = ϕE - ft sin 2ϕE + (f∕2)t sin 4ϕE,
    where f = 180∕π and t =tan2ϵ∕2.
  8. The so-called “declination” of the sun (δ in Figure A.11 which should not be confused with the magnetic declination D), which we will use as the latitude λs, is given by:

    δ = sin -1(sinϵsinϕe).
  9. Finally, the equation of time in degrees is given by E = 4(ϕs - A).

We can now calculate the Greenwich Hour Angle GHA from the Universal Time U (in minutes) by GHA = (U + E)4 + 180. The local hour angle (H in Figure A.11) is GHA + ϕL. We calculate β using the laws of spherical trigonometry (see Appendix A.3.1). First we calculate θ by the Law of Cosines (remembering that the cosine of the colatitude equals the sine of the latitude):

cosθ = sinλ  sin λ  + cosλ  cos,λ cosH
            L    s       L      s

and finally using the Law of Sines:

sin β = (cosλssinH )∕sinθ.

If λs < λL, then the required angle is the shadow direction β, given by: β= 180 - β. The azimuth of the desired direction is βplus the measured shadow angle α.