H:\EXCERC\GEOD98e.EX4.wpd

Department of Geophysics, Juliane Maries Vej 30, 2100 København Ø.

Geodesy Course, Exercise 4.

The purpose of th exercise is to demonstrate how geoidheights and gravity anomalies can be calculated using a spherical harmonic expansion. The nexessary data and programs are found on the Unix-net, /disk1/cct/geod/*4*, see also exercise 2. Please use your personal directory, so that the output does not owerwrite the output of other students.

The coefficients for the spherical harmonic series are found in the file /disk1/cct/cctf/osu91a1f . It consist of fully normalized coefficients up to and inclusive degree 360. After log-in go to directory geod and write head osu91a1f in order to view the first lines of the file. Each line (record) consist of:

degree (n), order (m), C(n,m), S(n,m),

and has the format (2I3,2D19.12) ( which is FORTRAN and means thet the record holds 2 integers with 3 digits and 2 double precision numbers with 12 digits).

Which coefficients are zero, and why are they zero ?

For the calculation on SGI computers (geb, seth, mani) use the program /disk1/cct/dgravsoft/geocol12. We want the geoid computed for Denmark in a grid of points with a distance of 0.1 degrees in latitude and 0.2 degrees in longitude. However, initially execute the calculation using the coefficients up to degree and order 10 with a grid spacing of 1 degree and 2 degrees. (This will be executed very fast). Then up to degree 360 with the smaller grid distance.

The selection of maximal order and of gridspacing is done by changing input-values to the program.

The program is run interactively or

by using the a prepared job-file: Type

ex4.inp

and the job will be done. You can change this file in order to vary the input parameters. Below is shown an interactive run using maximal degree 10 and the coarse grid-spacing. That geoid-heights are to be calculated is defined using the data-type input parameter equal to 11:

/disk1/cct/dgravsoft/geocol12

GEODETIC COLLOCATION, VERSION 1 FEB 1994, RELEASE 6 (UNIX)

NOTE THAT THE FUNCTIONALS ARE IN SPHERICAL APPROXIMATION

MEAN RADIUS = RE = 6371 KM AND MEAN GRAVITY 981 KGAL USED.

MAX NUMBER OF OBS= 3200, MAX NUMBER OF PARAMETERS=239

MAX NUMBER OF OBS IN GIVEN REF. FRAME = 3200

SIZE OF NORMAL EQ. BLOCKS=19800, SIZE OF POT.COFF. BLOCK= 130322

INTERACTIVE INPUT (T/F) t

INPUT: LTRAN, TRUE IF NON-STANDARD REF. SYSTEM IS USED

LPOT, TRUE IF SPHERICAL HARMONIC EXPANSION IS USED

LTEST, TRUE IF TEST-OUTPUT IS NEEDED

LLEG, TRUE IF LEGEND IS TO BE OUTPUT

LPARAM,TRUE IF PARAMETERS ARE TO BE DETERMINED

LNCOL, TRUE IF COLLOCATION IS NOT USED

LIOSOL,TRUE IF SOLUTION IS STORED OR RECOVERED

f t f f f t f

ARE ALL PARAMETERS OK ? t

INPUT CODE FOR BASIC REFERENCE SYSTEM:

0: USER DEFINED, 1: ED50 NORTH SEA, 2: ED50/EDOC,

3: NAD1927 /NEW MEXICO, 4: GRS67, 5: GRS80, 6: NWL9D,

7: BEST CURRENT, 8: BEST CUR. FAROE ISL, 9: ED50 FOR SF,

10: IAG-75, 11: KRASSOWSKY, DDR, 12: GERMAN DHDN, BESS. 5

REFERENCE SYSTEM: GRS1980.

A = 6378137.00 M

1/F = 298.2572221

GM= 0.3986005000E+15

REF.GRAVITY AT EQUATOR = 978032.6772 MGAL

POTENTIAL AT REF.ELL. = 62636860.8500 M**2/SEC**2

INPUT NAME OF POT.COEFF. SET OSU91A TO DEGREE 10.

SOURCE OF THE POTENTIAL COEFFICIENTS USED:

OSU91A TO DEGREE 10

INPUT: GM, SEMI-MAJOR AXIS (M), C(2,0), MAX. DEGREE

LFM, TRUE IF COEFF. IN INPUT STREEM AND *1.0D6

LBIN, TRUE IF ON BINARY FORM

LFORM, TRUE IF FORMAT IS INPUT

LINT, TRUE IF STORED AS INTEGERS

3.986005D14 6378136.2 0.0 10 f f t f

GM A COFF(5) MAX.DEGREE

0.39860050E+15 6378136.2 0.0000 10

INPUT FORMAT (2I4,2D18.0) F.EX. (2I3,2D19.2)

INPUT NAME OF FILE HOLDING COEFF. osu91a1f

NAME OF FILE HOLDING COEFFICIENTS: osu91a1f

COEFFICIENTS UP TO N=5

2 0 -0.484165533E-03 0.000000000E+00

2 1 0.857179552E-12 0.289607376E-11

2 2 0.243815798E-05 -0.139990175E-05

3 0 0.957139401E-06 0.000000000E+00

3 1 0.202968777E-05 0.249431310E-06

3 2 0.904648671E-06 -0.620437817E-06

3 3 0.720295507E-06 0.141470959E-05

4 0 0.540441630E-06 0.000000000E+00

4 1 -0.535373285E-06 -0.474065010E-06

4 2 0.350729847E-06 0.663967363E-06

4 3 0.991080200E-06 -0.202148896E-06

4 4 -0.190576532E-06 0.309704029E-06

INPUT: LGRID - TRUE IF COMPUTATIONS IN A GRID

LERR - TRUE IF ERROR ESTIMATES ARE TO BE COMPUTED

OR REPRODUCED IN OUTPUT

LCOMP- TRUE IF COMPUTED VALUES ARE SUBTRACTED FROM OBSERVED

t f f

INPUT GRID SPECIFICATION

MIN, MAX LATITUDE, MIN, MAX LONGITUDE, STEP IN LAT AND LONG

FUNCTIONAL TYPE (CODE), NEG. VALUE THEN SPH.EXP. SUBTRACTED

COORD.SYSTEM CODE (-1 THEN GLOBAL SYSTEM)

HEIGHT OF GRID POINTS (M)

LMAP - PRIMITIVE MAP OUTPUT

LPUNCH- OUTPUT TO UNIT 17

LMEAN - MEAN VALUES OUTPUT

LMAP&LPUNCH SIMULTANEOUS TRUE, ONLY MAP OUTPUT TO UNIT 17

GRID CONSIST OF 4 * 4 POINTS

54.5 57.5 7.0 13.0 1.0 2.0 11 -1 0.0 f t f (Tallet 11 ændres til 13 for Tyngder !!!)

INPUT NAME OF FILE TO HOLD RESULT grid11.ex4

SIMULTANEOUS OUTPUT TO FILE: grid11.ex4

ALL SPECIFICATIONS OK ? t

SELECTED GEOCENTRIC SYSTEM USED.

NO LATITUDE LONGITUDE H ZETA (M)

DEGREES DEGREES M

ERR POT

1 57.500000 7.000000 0.0 0.00 41.47

2 57.500000 9.000000 0.0 0.00 39.45

3 57.500000 11.000000 0.0 0.00 37.50

4 57.500000 13.000000 0.0 0.00 35.63

5 56.500000 7.000000 0.0 0.00 41.80

6 56.500000 9.000000 0.0 0.00 39.82

7 56.500000 11.000000 0.0 0.00 37.92

8 56.500000 13.000000 0.0 0.00 36.10

9 55.500000 7.000000 0.0 0.00 42.16

10 55.500000 9.000000 0.0 0.00 40.24

11 55.500000 11.000000 0.0 0.00 38.40

12 55.500000 13.000000 0.0 0.00 36.63

13 54.500000 7.000000 0.0 0.00 42.55

14 54.500000 9.000000 0.0 0.00 40.69

15 54.500000 11.000000 0.0 0.00 38.91

16 54.500000 13.000000 0.0 0.00 37.20

STOP ? t

After this contour the data using /disk1/cct/dgeoplot/geoplot21, see exercise 2. As data-file use the file grid11.ex4, if you use the names given in the example above. Remember first to type uniras .

/disk1/cct/dgeoplot/geoplot21

GEOPLOT-UNIRAS, VERSION 1.4, 2 NOV. 1995.

INPUT FIGURE TEXT (MAX. 120 CHAR).

Geoid heights from OSU91A to grad 10, c.i. 1 m. <Your name >.

SELECT PROJECTION:

0: LAMBERT CONFORMAL CONIC (2 STD. PARAL.) AUTO.

1: U T M, 2: TRANSVERSE MERCATOR, 3: MERCATOR

4: LAMBERT CONFORMAL CONIC (2 STD. PARAL.)

5: POLAR STEREOGRAPHIC (AZIMUTHAL-), 0

INPUT AREA BOUNDARIES (LATMIN,LATMAX,LONMIN,LONMAX)

54.5 57.5 7.0 13.0

INPUT NUMBER OF DIGITS IN GRID ANNOTATION 0

INPUT NAME OF FILE WITH COASTLINE COORDINATES ciakyst

DO YOU WANT CONTOUR MAP TO BE DRAWN ? (T/F) t

INPUT NAME OF FILE WITH DATA TO BE CONTOURED grid11.ex4

NUMBER OF DATA ELEMENTS, AND ELEMENT NUMBER USED

(IF DATA ON GRID FORM INPUT 0 1) 3 2

ARE DATA GRIDDED ? (T/F) t

INPUT NUMBER OF POINTS IN E AND S DIR. 4 4

INPUT CONTOUR INTERVAL 1.0

ANNOTATED CONTOURS ? (T/F) f

CONTOURING WITHIN A REGION ? (T/F) f

2 OR 3-D PLOT ? 2

OUTPUT TO UNIPICT FILE ? (T/F) f

COLOR CONTOURING ? (T/F) t

INPUT COLOR SCALE INDEX 3

IS COLOR SCALE DEFINED BY MIN, MAX OF DATA (T/F) t

INPUT NUMBER OF DEC. 0

DO YOU WANT POINTS TO BE PLOTTED ? (T/F) f

GROUTE> select MPOST; EXIT Eller: select lx11; EXIT

MIN,MAX LON,LAT 6.75 13.24 54.46 57.53

NUMBER OF POINTS INPUT 16, MIN CONT.= 35.000, MAX= 42.550

number of contours 8

PLOT SIZE 258.1 MM 202.2 MM

POINTS ON COASTLINE PLOTTED 335

After this the plot can be output to a post-script printer,

lpr -Pp243 POST, where the output is directed to the printer in room 243. Screen-plot is made, if the driver lx11 is used instead of MPOST.

After this run the jobs once more, now using data-code 13 in order to obtain gravity anomalies. Obviously the legend on the plot must be changed and it may be necessary to select another contouring interval. Then repeat the calculations by executing the calculations at 10000 m altitude (otherwise it is done at height 0). Then calculate gravity south of Ceylon in altitude 0 and 10000 m.The gravity anomalies are very large here.

It is possible to make post-script color plots see

the file /disk1/cct/instructions. You type:

more /disk1/cct/instructions .

It is maybe difficult to execute the runs fast, if many students are using the same computer. Type top to see who are running.