# Presentation: 371575

## Geodesy for Neutrino Physicistsby Wes Smart, Fermilab 

- Based on: “GPS Satellite Surveying”                By Alfred Leick, Wiley (1990)
- Geodesy : a branch of applied mathematics that determines the exact positions of points and figures and areas of large portions of the earth’s surface, the shape and size of the earth, and the variations of terrestrial gravity and magnetism.
- Or: what’s needed beyond the Flat Earth Society

## Outline 

-  Ellipsoid model of the earth
-  Three geodetic coordinate systems and the         .      .   transformations between them
-  Method of calculation
-  Excel spreadsheet to do these transformations    _[http://](file:///calcs.xls)__[home.fnal.gov/~smart/geodesy/calcs.xls](file:///calcs.xls)_
-  Examples (in excel): Chicago – Barcelona, NuMI
-  Height above sea level, geoid, geoid height
-  Summary

## Earth Modeled by Reference EllipseSpin Causes Larger Diameter at Equator than at Poles

|  |  |  |  |
| --- | --- | --- | --- |

- a=semi-major axis=6378137 m
- b=semi-minor axis=6356752.3141
- f=flattening= 1/298.25722210
- e=eccentricity=(0.00669438)**0.5**
- f=(a-b)/a
- e**2**=2f-f**2**=1-(b/a)**2**
- a-b= 21385 m

- b

- a

- GRS 80 (Geodetic Reference System) = Ellipse parameters in
- NAD 83 (North American Datum)

## The Geodetic and Geocentric Cartesian Coordinate Systems

- 

- N

- h

- P

- xy

- z

- 

- x

- y

**Surface Normal**

**Meridian**

- Looking from above Equator

- Looking from above North Pole

- z is the spin axis
-  is latitude
-  is longitude

- x=(N+h)coscos
- y=(N+h)cossin
- z=[(N(1-e**2**)+h]sin

- P

- N=a/(1-e**2**sin**2**)**0.5**
- e**2**=1-(b/a)**2**

**Greenwich**

- +

- -

**North   pole**

- +

- -

**North**

**South**

**East**

**West**

- +

**East**

**West**

- -

**(Not Origin)**

## Local Geodetic Coordinates

**up**

- P**1**(,h)

- xy

- z

**Normal to**
**Ellipsoid**

- Looking from above Equator

- z is the spin axis
-  is latitude
-  is longitude

-         A second point P**2** relative to P**1** is given by:
**n**=-(x**2**-x**1**)sincos-(y**2**-y**1**)sinsin+(z**2**-z**1**)cos
**e**=-(x**2**-x**1**)sin+(y**2**-y**1**)cos
**u**=(x**2**-x**1**)coscos+(y**2**-y**1**)cossin(z**2**-z**1**)sin

**North   pole**

**north**

**east **

- Specified for a point P**1**, Cartesian
**up** is along the normal to Ellipsoid
**north** is the intersection of the plane perpendicular to the normal containing P**1** and the plane containing the z (spin) axis and P**1**
**east** = the cross product: **north x up**

- Into screen

- h

## Compare Coordinate Systems

- System		Coordinates	Range		Cartesian/	Familiarity
- _ __	__	__	__	__	__	__Easy Calcs ?                          ._
- Geodetic	Latitude		global		no		medium
- 		Longitude
- 		Ellipsoidal ht.
- Geocentric	x, y, z		global		yes		low
- Cartesian
- Local		north,		local		yes		high
- Geodetic	east, up

## Calculation Method

- Get Geodetic coordinates of points: may need to find ellipsoidal heights from elevations
- Use Spreadsheet to find Geocentric Cartesian coordinates
- Do desired calculations in the Geocentric Cartesian coordinate system (which you already know how to do)
- If needed, use the inverse transformation to calculate Geodetic coordinates of results

## Azimuth Example Chicago to Barcelona

**up**

- xy

- z

**Normal to**
**Ellipsoid**

- Looking from above Equator

**North   pole**

**north**

**east **

- Into screen

- Looking from above North Pole
- Dashed lines are not in the plane

- y

- x

- Chicago

- n**c**

- e**c**

- n**b**

- e**b**

- Barcelona

- Plane of right plot

- These 2 cities are both at 42**o** N Latitude and 90**o **apart in Longitude.  Beam must leave Chicago north of east and would arrive in Barcelona from north of west.  These directions are not 180**o** apart because east is a different direction in each city.  (This is also true for north and up.)  This applies as well for an airplane on the great circle route between the two cities.

## Spreadsheet Results; Chicago to Barcelona

| Part A |  | Geodetic | Coordinates |  | Geocentric | Cartesian | Coordinates |
| --- | --- | --- | --- | --- | --- | --- | --- |
| A1 |  | Ellipsoid height | Latitude   0 | Longitude  0 | x | y | z |
| "Chicago" | NAD 83 | 0.0000 | 42.0000000 | -88.0000000 | 165667.8748 | -4744107.2248 | 4245603.8360 |
| A2 |  |  |  |  |  |  |  |
| "Barcelona" | NAD83 | 0.0000 | 42.0000000 | 2.0000000 | 4744107.2248 | 165667.8748 | 4245603.8360 |
|  |             Angles A1 to A2 |  |  | d | dx | dy | dz |
|  | azimuth | vertical |  | 6713270.3216 | 4578439.3500 | 4909775.0997 | 0.0000 |
|  | 0.98108986 | -0.55328241 | rad |  | Local Geodetic Coordinates of A2; ref A1 |  |  |
|  | 56.21230819 | -31.7007471 | deg |  | dn | de | du |
|  |  |  |  |  | 3176362.2981 | 4746998.9683 | -3527707.7199 |
| Part B |  | Geodetic | Coordinates |  | Geocentric | Cartesian | Coordinates |
| B1 |  | Ellipsoid height | Latitude   0 | Longitude  0 | x | y | z |
| "Barcelona" | NAD83 | 0.0000 | 42.0000000 | 2.0000000 | 4744107.2248 | 165667.8748 | 4245603.8360 |
| B2 |  |  |  |  |  |  |  |
| "Chicago" | NAD 83 | 0.0000 | 42.0000000 | -88.0000000 | 165667.8748 | -4744107.2248 | 4245603.8360 |
|  |             Angles B1 to B2 |  |  | d | dx | dy | dz |
|  | azimuth | vertical |  | 6713270.3216 | -4578439.3500 | -4909775.0997 | 0.0000 |
|  | -0.98108986 | -0.55328241 | rad |  | Local Geodetic Coordinates of B2; ref B1 |  |  |
|  | 303.7876918 | -31.7007471 | deg |  | dn | de | du |
| All lengths in meters |  |  |  |  | 3176362.2981 | -4746998.9683 | -3527707.7199 |

## Spreadsheet Results; NuMI Target to Far

| Part A |  | Geodetic | Coordinates |  | Geocentric | Cartesian | Coordinates |
| --- | --- | --- | --- | --- | --- | --- | --- |
| A1 |  | Ellipsoid height | Latitude   0 | Longitude  0 | x | y | z |
| ACTRN1 | NAD 83 | 153.9683 | 41.8320235 | -88.2655587 | 144058.5523 | -4757396.9876 | 4231823.0745 |
| A2 |  |  |  |  |  |  |  |
| FARctr 2001 | (NAD83) | -248.3992 | 47.8202665 | -92.2414120 | -167796.9924 | -4287098.7216 | 4703296.8722 |
|  |             Angles A1 to A2 |  |  | d | dx | dy | dz |
|  | azimuth | vertical |  | 735337.9371 | -311855.5447 | 470298.2661 | 471473.7977 |
|  | -0.41723537 | -0.05829776 | rad |  | Local Geodetic Coordinates of A2; ref A1 |  |  |
|  | 336.0941743 | -3.3402156 | deg |  | dn | de | du |
|  |  |  |  |  | 671113.2727 | -297478.1082 | -42844.2768 |
| Part B |  | Geodetic | Coordinates |  | Geocentric | Cartesian | Coordinates |
| B1 |  | Ellipsoid height | Latitude   0 | Longitude  0 | x | y | z |
| FARctr 2001 | (NAD83) | -248.3992 | 47.8202665 | -92.2414120 | -167796.9924 | -4287098.7216 | 4703296.8722 |
| B2 |  |  |  |  |  |  |  |
| ACTRN1 | NAD 83 | 153.9683 | 41.8320235 | -88.2655587 | 144058.5523 | -4757396.9876 | 4231823.0745 |
|  |             Angles B1 to B2 |  |  | d | dx | dy | dz |
|  | azimuth | vertical |  | 735337.9371 | 311855.5447 | -470298.2661 | -471473.7977 |
|  | 2.67536239 | -0.05718496 | rad |  | Local Geodetic Coordinates of B2; ref B1 |  |  |
|  | 153.2869734 | -3.2764567 | deg |  | dn | de | du |
|  |  |  |  |  | -655781.0335 | 330010.3321 | -42027.3543 |

## Spreadsheet Results; MINOS Near to Far

| Part A |  | Geodetic | Coordinates |  | Geocentric | Cartesian | Coordinates |
| --- | --- | --- | --- | --- | --- | --- | --- |
| A1 |  | Ellipsoid height | Latitude   0 | Longitude  0 | x | y | z |
| Near | NAD 83 | 93.4971 | 41.8405633 | -88.2706209 | 143617.7785 | -4756732.2723 | 4232489.4513 |
| A2 |  |  |  |  |  |  |  |
| FARctr 2001 | (NAD83) | -248.3992 | 47.8202665 | -92.2414120 | -167796.9924 | -4287098.7216 | 4703296.8722 |
|  |             Angles A1 to A2 |  |  | d | dx | dy | dz |
|  | azimuth | vertical |  | 734298.6171 | -311414.7709 | 469633.5508 | 470807.4209 |
|  | -0.41729429 | -0.05813482 | rad |  | Local Geodetic Coordinates of A2; ref A1 |  |  |
|  | 336.0907986 | -3.3308801 | deg |  | dn | de | du |
|  |  |  |  |  | 670153.5882 | -297099.9605 | -42664.2796 |
| Part B |  | Geodetic | Coordinates |  | Geocentric | Cartesian | Coordinates |
| B1 |  | Ellipsoid height | Latitude   0 | Longitude  0 | x | y | z |
| FARctr 2001 | (NAD83) | -248.3992 | 47.8202665 | -92.2414120 | -167796.9924 | -4287098.7216 | 4703296.8722 |
| B2 |  |  |  |  |  |  |  |
| Near | NAD 83 | 93.4971 | 41.8405633 | -88.2706209 | 143617.7785 | -4756732.2723 | 4232489.4513 |
|  |             Angles B1 to B2 |  |  | d | dx | dy | dz |
|  | azimuth | vertical |  | 734298.6171 | 311414.7709 | -469633.5508 | -470807.4209 |
|  | 2.67536239 | -0.05718496 | rad |  | Local Geodetic Coordinates of B2; ref B1 |  |  |
|  | 153.2869734 | -3.2764567 | deg |  | dn | de | du |
| All lengths in meters |  |  |  |  | -654854.1585 | 329543.8985 | -41967.9533 |

## Spreadsheet “Subroutines”

| Convert degrees, minutes, seconds to degrees |  |  |  | Convert degrees to degrees, minutes, seconds |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 0 | ' | " | 0 | 0 | 0 | ' | " |
| 41 | 49 | 55.2846 | 41.83202351 | 41.83202351 | 41 | 49 | 55.284625 |
| Enter absolute value input data only in cells with black borders using paste special (value) |  |  |  |  |  |  |   |
| For negative input data, enter absolute value and append minus sign to result. |  |  |  |  |  |   |   |
|  |  |  |  |  |  |  |  |
| Find ellipsoid height, h, in meters |  |  |  |   | Linear Interpolation |  |  |
| from elevation, H, in feet. h=H+N |  |  |  |   | Parameter | Result |  |
| H (feet) | N (m) | h (m) |  | 1 | 0 | 0 |  |
| 1210 | -32.0 | 336.8080 |  | 2 | -10 | 100 |  |
|  |  |  |  |   | 7 | -70 |  |
|  |  |  |  |  |  |  |  |
|  | Find h, Latitude, Longitude from x,y,z |  |  |   | Only change cells with black borders. |  |  |
|  |   |  |  Lat k-1    0 |  |  |  |   |
|  |   | h       ( m) | 41.83202351 | Longitude    0 | x input | y input | z input |
|  |   | 153.9683 | 41.83202351 | -88.26555873 | 144058.5523 | -4757396.9876 | 4231823.0745 |
|  |   | Lat k  | 41.83202351 | 0 |  |  |   |
|  |   | Lat 1 | 41.83202813 | 0   3. Iterate Latitude equation by copying red border cell to  |  |  |   |
|  | 1. Enter x,y,z using paste special (value) |  |  |        black border cell, using paste special (value), |  |  |  |
|  | 2. Copy "Lat 1" or red border cell to black |  |  |        until the red and black bordered cells agree. |  |  |  |
|  |      border Lat k-1cell using paste special |  |  |    4. Copy h, Lat, Long to desired locations, using |  |  |  |
|  |      (value) |   |   |        paste special (value), |  |   |   |

- Linear Interpolation

- Use to find the speadsheet input parameter which gives the desired result for an output value.

- All data input should be by typing or paste special value.
- Input only into cells marked for input.
- Select the input parameter and output result you wish to use, put desired value of result into the answer line of the “subroutine”
- Guess a value for the parameter, put in spreadsheet, copy parameter and result into line 1 of the “subroutine”
- Repeat for line 2
- Put answer parameter value in spreadsheet, copy it and result into line 1 or 2 (pick the line which has its result further from the desired value).
- Repeat last step until the speadsheet result has the desired value.

|   | Linear Interpolation |  |
| --- | --- | --- |
|   | Parameter | Result |
| 1 | 0 | 0 |
| 2 | -10 | 100 |
| answer | 7 | -70 |

## Spreadsheet Results; Offaxis Detector

|  |  | Geodetic | Coordinates |  | Geocentric | Cartesian | Coordinates |
| --- | --- | --- | --- | --- | --- | --- | --- |
| A1 |  | Ellipsoid height | Latitude   0 | Longitude  0 | x | y | z |
| ACTRN1 | NAD 83 | 153.9683 | 41.8320235 | -88.2655587 | 144058.5523 | -4757396.9876 | 4231823.0745 |
| A2 |  |  |  |  |  |  |  |
| 10 mrad | NAD83 | 305.0000 | 47.1500000 | -91.8152727 | -137659.0444 | -4343503.4468 | 4653345.3173 |
|  |             Angles A1 to A2 |  |  | d | dx | dy | dz |
|  | azimuth | vertical |  | 654487.3327 | -281717.5967 | 413893.5409 | 421522.2428 |
|  | -0.42425381 | -0.05116410 | rad |  |  |  |  |
|  | 335.6920471 | -2.9314871 | deg |  | Local Geodetic Coordinates of A2; ref A1 |  |  |
|  |  |  |  |  | dn | de | du |
|  |  |  |  |  | 595683.9779 | -269061.1748 | -33471.6487 |
|  | along y' (m) | total transverse | angle (rad) |  |  |  | angle (deg) |
| A2 Relative to NuMI Beam |  |  |  |  |  |  |  |
| 10 mrad | 654454.6086 | 6544.7643 | 0.010000 | 500635.6119 | -6054.0536 | 2486.4382 | 0.5730 |
|  | Longitude  0 | Angle to B (rad) |  |  |  |  |  |
|  | -91.8152727 | 0.010000 |  |  |  |  |  |
|   | Parameter | Result |  |  |  |  |  |
| 1 | -91.81527 | 0.010000 |  |  |  |  |  |
| 2 | -91.81500 | 0.009981 |  |  |  |  |  |
| answer | -91.81527266 | 0.010000 |  |  | All lengths in meters |  |  |

- Find latitude, longitude, and ellipsoidal height
- from geocentric Cartesian coordinates x,y,z
- First approximate solution for    tan**1**=z/[(1-e**2**)(x**2**+y**2**)**0.5**]
- Then find  by iteration                 tan=[z+ae**2**sin/(1-e**2**sin**2**)**0.5**]/(x**2**+y**2**)**0.5**
- Finally                  tan=y/x   and   h=[(x**2**+y**2**)**0.5**)/cos]-N

- Inverse Transformation

## Heights

- Geoid

- Ellipsoid

- H

**N**

- h

- P

**H**, Orthometric height, is above “sea level”, ie elevation
**h** is the ellipsoidal height, GPS measures in h directly
**N****, **the geoid height, is about -32 m at Soudan and Fermilab
- To calculate **N:**  http://www.ngs.noaa.gov/GEOID/GEOID03/download.html

_**Geoid**_* *is the equipotential surface with gravity potential chosen such that on average it coincides with the global ocean surface.
_**N**_ accounts for the difference between the real earth and the ideal reference ellipsoid used for calculation. **N** varies with latitude and longitude.

- h=H+**N**

## Geoid Heights for North America

## Slide 19

## Summary

- Earth is modeled well by ellipsoid
- 3 geodetic coordinate systems
- Geodetic: Latitude, Longitude, Ellipsoidal height
- Geocentric Cartesian: x, y, z
- Local Geodetic: north, east, up
- Transformations between them with Excel
- Transform points to Geocentric Cartesian where calculations are easy and familiar
- If desired, transform answers back to Geodetic Coordinates