MAGNETIC SOURCE DEPTH: CONSTANT PRISM DIMENSIONS
C***********************************************************************
C
C
C PROGRAM MAGSDQ (MAGNETIC SOURCE DEPTH: CONSTANT PRISM DIMENSIONS)
C (QUARTER DEGREE RESOLUTION)
C
C***********************************************************************
C
C
C PROGRAMMED BY: JOHN M. QUINN 10/1/88
C GEOMAGNETICS DIVISION CODE GMD
C U.S. NAVAL OCEANOGRAPHIC OFFICE
C STENNIS SPACE CENTER (SSC), MS 39522-5001
C
C
C***********************************************************************
C
C
C NOTE: MODIFIED BY JOHN M. QUINN 07/09/04
C
C THIS IS A MODIFICATION OF THE ORIGINAL MASDEP PROGRAM. IT
C ALLOWS THE MAGNETIC INPUT DATA TO BE GENERATED FROM A
C MAGNETIC FIELD MODEL RATHER THAN A SET OF LAT X LON GRID
C VALUES.
C
C
C***********************************************************************
C
C
C PURPOSE: THIS ROUTINE COMPUTES, USING THE MAGPOT SUBROUTINE,
C 3 VECTOR MAGNETIC FIELD COMPONENTS AND 6 GRADIENT
C MAGNETIC FIELD COMPONENTS GENERATED FROM A SPHERICAL
C HARMONIC POTENTIAL FIELD MODEL WHICH WAS ITSELF
C DERIVED FROM MAGNETIC SURVEY DATA. THESE DATA ARE
C INVERTED TO YIELD THE DEPTH TO THE MAGNETIC SOURCE AND
C THE SOURCE MAGNETIZATION VECTOR.
C
C
C***********************************************************************
C
C
C ASSUMPTION: THE OBSERVED MAGNETIC AND GRADIENT FIELD COMPONENTS
C AT A SPECIFIED (LAT, LON, ALT) OBSERVATION POINT, ARE
C ASSUMED TO BE DUE TO A MAGNETIZED RECTANGULAR PRISM
C LOCATED DIRECTLY BELOW THE OBSERVATION POINT, BELOW
C EARTH'S SURFACE. THE LATERAL DIMENSIONS OF THE PRISM
C ARE VARY WITH DEPTH (ROUGHLY BEING EQUAL TO THE
C DISTANCE BETWEEN THE SOURCE DEPTH AND THE OBSERVATION
C POINT ABOVE THE EARTH'S SURFACE (APPROX. 400 KM FOR
C MAGNETIC FIELD DATA DERIVED FROM THE CHAMP SATELLITE
C OBSERVATIONS).
C
C IT IS ASSUMED THAT THE PRISM'S CENTER IS DIRECTLY
C BELOW THE OBSERVATION POINT.
C
C THE POINT OF OBSERVATION IS ASSUMED TO BE ABOUVE THE
C EARTH'S SURFACE. DUE TO COMPUTER WORD SIZE
C LIMITATIONS IN DOUBLE PRECISION, IT IS NOT
C POSSIBLE TO PLACE THE OBSERVATION POINT TOO HIGH
C BECAUSE THE HIGHER HARMONICS OF THE MAGNETIC FIELD
C MODEL WILL NOT CONTRIBUTE TO THE OBSERVED MAGNETIC
C FIELD VALUES. THE MF3 MODEL OF STEFAN MAUS IS BASED
C
C ON CHAMP SATELLITE DATA COLLECTED AT 400 KM. BUT,
C DUE TO COMPUTER LIMITATIONS OF A TYPICAL PERSONAL
C COMPJTER (PC), ONE CAN NOT CHOOSE THE OBSERVATION
C POINT TO BE AT SATELLITE ALTITUDE AND EXPECT THAT
C ALL HARMONICS FROM 0 TO 90 WILL CONTRIBUTE TO THE
C OBSERVED MAGNETIC FIELD AND GRADIENT FIELD VALUES.
C GENERALLY, 50 KM IS HIGH ENOUGH TO MINIMZE NOISE
C ASSOCIATED WITH DOWNWARD CONTINUATION AND YET LOW
C ENOUGH TO HAVE ALL HARMONICS CONTRIBUTING TO THE
C FIELD VALUES. PC'S ARE 32 BIT MACHINES WITH 64 BIT
C DOUBLE PRECISION. MORE EXOTIC MACHINES LIKE THE
C SGI ORIGIN 2000 ARE 48 BIT MACHINES THAT HAVE 96 BITS
C DOUBLE PRECISION AND MAY THEREFORE ALLOW THE
C OBSERVATION POINT TO BE RAISED TO SATELLITE
C ALTITUDES.
C
C ONCE THE OBSERVATION ALTITUDE IS CHOSEN, THE LATERAL
C DIMENSIONS OF THE PRISM ARE DETERMINED BY NOTING
C THAT THE MAGNETIZATION SOURCE DEPTH MUST LIE
C BETWEEN THE EARTH'S SURFACE AND THE CURIE DEPTH,
C WHICH GLOBALLY AVERAGES ABOUT 25 KM. ACTUAL CURIE
C DEPTHS IN SOME LOCATIONS MAY BE ROUGHLY TWICE
C THE AVERAGE CURIE DEPTH. THUS THE SOURCE DEPTHS
C THAT THIS PROGRAM ATTEMPTS TO COMPUTE MUST LIE
C BETWEEN THE EARTH'S SURFACE AND TWICE THE AVERAGE
C CURIE DEPTH. ADJUSTING THE LATERAL DIMENSIONS
C OF THE PRISM WILL HAVE A DIRECT EFFECT ON THE
C RESULTING SOURCE DEPTH. IF THE LATERAL DIMENSIONS
C ARE TO LARGE, THE SOURCE DEPTH WILL BE FAR BELOW
C ANY REASONABLE CURIE DEPTH. IF THE LATERAL
C DIMENSIONS ARE TOO SMALL, THE SOURCE DEPTH WILL BE
C
C COMPUTED TO BE ABOVE THE EARTH'S SURFACE, WHICH IS
C PHYSICALLY IMPOSSIBLE. SO, FOR ANY GIVEN OBSERVATION
C ALTITUDE, THERE IS A NARROW RANGE OF LATERAL PRISM
C DIMENSIONS THAT WILL YIELD ACCEPTABLE SOURCE
C DEPTHS. FOR AN OBSERVATION ALTITUDE OF 50 KM, THE
C LATERAL DIMENSIONS SHOULD BE NEAR 450 KM. THIS IS
C DETERMINED BY RUNNING THE PROGRAM USING DIFFERENT
C LATERAL DIMENSIONS AND EXAMINING THE GLOBAL
C DISTRIBUTION OF MAGNETIC SOURCE DEPTHS. THE GLOBAL
C AVERAGE SOURCE DEPTH FOR A 50 KM OBSERVATION POINT,
C WITH A PRISM OF 450 KM LATERAL DIMENSIONS AND 2 KM
C THICKNESS IS ABOUT 9 KM.
C
C THE PRISM THICKNESS IS BASED ON OCEAN DRILLING
C WHICH HAVE A STRONGLY MAGNETIZED LAYER OF ABOUT
C 0.5 KM THICKNESS, AND A COUPLE OF ADJACENT LAYERS
C OF LESSER MAGNETIC STRENGTH TOTALLING ABOUT 1.5 KM
C THICK. THUS, 2 KM IS USED AS THE THICKNESS OF THE
C MAGNETIZED LAYER. TO THE EXTENT THAT THE THICKNESS
C IS NOT ENTIRELY CORRECT, THE RESULTING MAGNETIZATION
C CAN BE TREATED AS AN EQUIVALENT MAGNETIZATION.
C
C THE PRISM SOURCE DEPTH STARTS AT THE SURFACE AND IS
C INCREMENTALLY MOVED TOWARD THE EARTH'S CENTER. AT
C EACH INCREMENTAL STEP, THE GEOMETRICAL ARRAYS
C LAMDA AND DLAMDA ARE COMPUTED AND USED ALONG WITH
C KNOWN MAGNETIC FIELD AND GRADIENT MAGENTIC FIELD
C VALUES TO CONSTRUCT A CHI-SQUARE FUNCTION. THIS
C FUNCTION IS MINIMIZED AT THE SOURCE DEPTH. THE
C CHI-SQUARE FUNCTION IS COMPUTED FROM EQ. (5) BELOW.
C
C
C THIS PROCESS IS REPEATED AT EACH LAT X LON GRID
C POINT. THE LAT LON GRID CAN COVER EITHER ALL OR PART
C OF THE EARTH AT ANY DESIRED GRID SPACING. HOWEVER.
C FOR CONVENIENCE, THIS PROGRAM ASSUMES A 1 DEG. X
C 1 DEG. LAT X LON GRID SPACING.
C
C
C***********************************************************************
C
C
C METHOD: THE MAGNETIC FIELD IN A SURVEYED AREA IS DIVIDED
C INTO 3 DISTINCT PARTS, THE MAIN (CORE GENERATED)
C FIELD, THE REGIONAL (DEEP CRUST AND UPPER MANTLE
C GENERATED) FIELD AND THE LOCAL (MID TO UPPER
C CRUST GENERATED) FIELD.
C
C IN TERMS OF SPHERICAL HARMONICS, WE TAKE THE CORE FIELD
C TO BE ASSOCIATED WITH SPHERICAL HARMONIC DEGREES 1
C THROUGH 12, THE REGIONAL OR UPPER MANTLE/LOWER
C LITHOSPHERE IS ASSOCIATED WITH SPHERICAL HARMONIC DEGREES
C 16 THROUGH 30, AND THE LOCAL OR CURUSTAL/UPPERC LITHOSPHERIC FIELD WITH SPHERICAL HARMONIC DEGRES OF 31
C AND ABOVE. DEGREES 13, 14, AND 15 ARE CONSIDERED TO BE
C ASSOCIATED WITH THE CORE-MANTLE BOUNDARY TRANSITION.
C
C THE REGIONAL AND LOCAL (I.E., LOWER AND UPPER
C LITHOSPHERIC) SOURCE RGEIONS, BOTH OF WHICH ARE ASSUMED
C TO BE ABOVE THE CURRIE DEPTH, ARE INVERTED SEPARATELY,
C USING THE QUINN-SHIEL INVERSION TECHNIQUE, WHICH ASSUMES
C THAT THE OBSERVED FIELDS ARE DUE TO MAGNETIZED PRISMS.
C THE INVERSION ROUTINE SEEKS TO DETERMINE THE LOCATION AND
C
C DIMENSIONS OF EACH SOURCE PRISM, AS WELL AS THE MAGNITUDE
C OF THE X, Y AND Z COMPONENTS OF MAGNETIZATION WITHIN EACH
C PRISM.
C
C INVERSIONS MAY BE PERFORMED ON ANY GROUP OF USER DEFINED
C SPHERICAL HARMONICS (I USE 16-30, 31-60, AND 61-90).
C EACH GROUP IS COMPUTED INDEPENDENTLY OF THE OTHERS.
C ALTERNATIVELY, ONE CAN USE ALL OF THE SPHERICAL HARMONICS
C 16-90) AT ONE TIME.
C
C THE LOWER HARMONICS APPEAR TO BE ASSOCIATED WITH GREATER
C DEPTHS. EACH HARMONIC GROUP IS TREATED INDEPENDENTLY OF
C THE OTHERS. THAT IS, ONE PICKS ANY DESIRED HARMONIC GROUP
C AND DOES THE COMPUTATIONS ON THAT GROUP ONLY. WE ARE
C INTERESTED ONLY IN THE CENTER OF THE PRISM.
C
C THUS, DIRECTLY BELOW EACH OBSERVATION POINT THERE WILL BE
C THREE PRIMS. THESE PRISMS MAY OVERLAP IN THE VERTICAL
C DIRECTION, IN WHICH CASE, THAT PART WHICH OVERLAPS
C CONSISTS OF A MAGNETIZATION THAT IS THE VECTOR SUM OF THE
C TWO OR THREE OVERLAPPING PRISMS. THUS, ALONG VERTICAL
C CENTERLINE THROUGH EACH PRISM AND THE SURFACE GRID POINT,
C A MAGNETIZATION PROFILE IS CONSTRUCTED. IF THE ENTIRE
C IS PROCESSED IN THIS MANNER, THERE WILL BE 12 GRIDS, FOUR
C GRIDS FOR EACH OF THE 3 HARMONIC GROUPS ANALYZED. ONE
C GRID CORRESPONDS TO THE GLOBAL SOURCE DEPTH, WHILE THE
C THREE GRIDS CORRESPOND TO MAGNETIZATION (OR EQUIVALENT
C MAGNETIZATION IF THE SOURCE DEPTH IS BELOW THE CURIE DEPTH)
C OF EACH MAGNETIZATION VECTOR COMPONENT.
C
C IN THE EVENT THAT THAT THE SOURCE DEPTH IS DETERMINED TO
C
C BE BELOW THE CURIE DEPTH, THEN THE CONCLUSION MUST BE
C THAT AN ELECTRIC CURRENT EXISTS AT THAT DEPTH. THEN,
C GIVEN THE ELECTRICAL CONDUCTIVITY AT THAT DEPTH, ONE MAY
C ATTEMPT TO USE THE EQUIVALENT MAGNETIZATION VECTOR
C COMPONENTS TO DETERMINE THE ELECTRIC CURRENT VECTOR AT
C THAT DEPTH.
C
C BY CONTINUOUSLY SHIFTING FROM ONE GRID POINT TO THE NEXT
C AND REPEATING THE COMPUTATIONS, THE ENTIRE WORLD WILL
C BE ANALYZED.
C
C IT IS POSSIBLE TO SOLVE FOR THE LENGTH, WIDTH, AND
C THICKNESS OF EACH PRISM, AS WELL AS FOR THE DEPTH OF EACH
C PRISM. THIS GIVES 4 GEOMETRIC PARAMETERS. ADDITIONALLY,
C ONE CAN ALSO SOLVE FOR THE 3 MAGNETIZATION COMPONENTS.
C THUS, THERE ARE A TOTAL OF 7 POSSIBLE UNKNOWNS.
C
C HOWEVER, THE LATERAL DIMENSIONS OF THE PRISM ARE NOT
C CRITICAL, AS LONG AS THEY ARE SUFFICIENTLY LARGE (I.E.,
C ON THE ORDER OF THE DISTANCE FROM THE PRISM CENTER TO
C THE OBSERVATION POINT.
C
C AS MENTIONED ABOVE, THE PRISM THICKNESS IS ARBITRARY.
C HOWEVER, SEISMIC AND OTHER GEOPHYSICAL DATA MAY BE
C USEFUL TO CONSTRAIN THIS PARAMETER. FOR INSTANCE,
C IN OCEAN AREAS, THERE ARE KNOW TO BE THREE DISTINCT
C MAGNETIZED LAYERS TOTALING 2.5 KM OR SO IN THICKNESS.
C ONE LAYER IS STRONGLY MAGNETIZED AND THE OTHER TWO ARE
C OF LESSER STRENGTH. THIS ROUTNE ASSUMES A 2
C THISKNESS OF THE MAGNETIZATION LAYER, WORLDWIDE.
C
C
C THUS, IT IS POSSIBLE TO LIMIT THE NUMBER OF UNKNONS FOR
C EACH PRISM TO 4. BUT ONLY THE SOURCE DEPTH NEED BE
C COMPUTED VIA NONLINEAR INVERSION. ONCE THE
C GEOMETRIC PARAMETERS HAVE BEEN DETERMINED, THE
C MAGNEIZATION COMPONENTS MAY BE COMPUTED DIRECTLY, AND
C MORE OR LESS UNIQUELY BY SOLVING A SIMPLE LINEAR
C SYSTEM OF THREE EQUATIONS.
C
C THE FACT THAT THE MAGNETIC GRADIENTS ARE USED IN
C TO THE MAGNETIC FIELDS ALLOWS FOR THE SEPPARTION OF
C SOURCE DEPTH PROBLEM INTO A GEOMETRIC PART AND A
C GEOPHYSICAL PART. THS HELPS TO CIRCUMVENT THE
C NOTORIOUS NON-UNIQUENESS PROBLEM TO SOME EXTENT.
C
C THAT IS, IF ONE ONLY USED THE MAGNETIC FIELD DATA AND
C NOT THE GRADIENT FIELD DATA, THEN ONE COULD
C ARBITRARILY CHANGE THE SOURCE DEPTH AND ARRIVE AT AN]
C EQUALLY ACCEPTABLE MAGNETIZATION, AND VICE VERSA.
C
C USING BOTH THE MAGNETIC FIELD AND ITS GRADIENTS,
C CAUSES A SEPARATION OF THE TWO PARTS OF THE PROBLEM.
C THUS, ONCE THE GEOMETRIC PARAMETERS ARE DETERMINED,
C VIA THE NONLINEAR LEAST SQUARES PROCESS, THERE IS
C JUST ONE SET OF MAGNETIZATION PARAMETERS THAT WILL
C SATISFY THE GEOMETRIC PART OF THE PROBLEM.
C
C GENERALLY, ONCE A SET OF GEOMETRIC SET OF PARAMETERS
C ARE DETERMINED, USING THE GUIDELINES SET ABOVE, ANY
C CHANGE IN THE LATERAL DIMENSIONS OF THE PRRISM WILL
C NOT DRASTICALLY AFFECT EITHER THE SOURCE DEPTH OR THE
C MAGNETIZATION PARAMETERS. SO, THE SOLUTION TO THE
C
C SOURCE DEPTH PROBLEM GIVEN HERE
C 1) SOLVE THE NONLINEAR QUINN-SHIEL X**2 EQUATION FOR THE
C UNKNOWN GEOMETRIC PARAMETERS. SINCE THE PRISM IS
C ASSUMED TO BE DIRECTLY BELOW THE OBSERVATION POINT
C THE X AND Y POSITIONS ARE AT THE ORIGIN, SO THAT ONLY
C THE PRISM DEPTH IS UNKNOWN. THAT IS, 6 OF THE 7
C GEOMETRIC PARAMETERS ARE FIXED A PRIORI. THE
C 7'TH PARAMETER (PRISM DEPTH BELOW THE EARTH'S SURFACE)
C IS DETERMINED BY LOWERING THE PRISM IN INCREMENTAL
C STEPS (E.G., .1 KM) AND AT EACH STEP SOLVING A
C NONLINER LEAST-SQUARES PROBLEM. THE STEP HAVING
C THE LEAST RESIDUAL ERROR IS THE SOLUTION.
C
C 2) USING THE GEOMETRICAL PARAMETERS OBTAINED IN STEP (1)
C SOLVE A LINEAR INVERSION PROBLEM TO OBTAIN THE 3
C COMPONENTS OF THE MAGNTIZATION VECTOR (MX,MY,MZ) FOR
C EACH THE PRISM.
C
C 3) SHIFT THE OBSERVATION POINT AROUND THE WORLD IN
C INCREMENTAL LAT X LON STEPS UNTIL THE ENTIRE WORLD
C IS COVERED. A 1 DEGREE LAT X LON STEP SIZE MAY BE
C ADEQUATE. HOWEVER, ONE CAN CHOOSE SMALLER STEP
C LENGTHS IF DESIRED. THE SMALLER THE STEP SIZE, THE
C LONGER THE PROGRAM MUST RUN.
C
C
C***********************************************************************
C
C
C PROCEDURE: THE NONLINEAR INVERSION OF THE LOCAL OR REGIONAL
C MAGNETIC FIELD AND GRADIENT DATA DETERMINES THE
C GEOMETRIC PARAMETERS ASSOCIATED WITH THE PRISM
C
C AND IS CARRIED OUT VIA UNIFORM SAMPLING OF THE
C ENTIRE PARAMETER SPACE (A USEFUL METHOD FOR A SMALL
C NUMBER OF UNKNOWN PARAMETERS, WHICH IN THIS CASE
C IS JUST THE PRISM DEPTH BELOW THE EARTH'S SURFACE).
C
C THE SUBSEQUENT LINEAR INVERSION FOR THE PRISM
C MAGNETIZATIONS IS CARRIED OUT USING THE STANDARD
C LINEAR INVERSION TECHNIQUES.
C
C
C***********************************************************************
C
C
C NOTE: IT IS ASSUMED IN THIS PROGRAM THAT THE OBSERVED
C RESIDUAL MAGNETIC FIELD COMPONENTS AND GRADIENT
C FIELD COMPONENTS COLLECTED ABOVE THE EARTH'S SURFACE
C ARE DUE TO MAGNETIZED SOURCES LOCATED ABOVE THE
C CURIE DEPTH. THIS ASSUMPTION IS NOT ENTIRELY
C CORRECT, SINCE PART OF THOSE OBSERVED FIELD
C VALUES ARE ALSO LIKELY TO BE DUE TO CURRENTS
C ASSOCIATED WITH THE MOTION OF CONDUCTING FLUIDS
C (EG. VISCOUS FLUIDS LIKE MAGMA IN THE MANTLE AND
C LOWER CRUST OR SALT WATER FLOWING THROUGH CRACKS
C AND CREVASES OF THE UPPER CRUST) THROUGH THE EARTH'S
C MAIN MAGNETIC FIELD. THE COUPLING OF THE MOTION OF
C CONDUCTING FLUIDS TO THE EARTH'S MAIN MAGNETIC FIELD
C INDUCES CURRENTS, WHICH IN TURN PRODUCE SECONDARY
C MAGNETIC FIELDS. THIS COUPLING CAN OCCUR IN THE
C MANTLE OR THE CRUST, ABOVE OR BELOW THE CURIE DEPTH.
C ABOVE THE CURIE DEPTH IT IS NOT POSSIBLE TO SEPARATE
C MAGNETIZATION EFFECTS FROM INDUCTION EFFECTS. BELOW
C
C THE CURIE DEPTH, PRISMS IDENTIFIED AS SOURCES OF THE
C OBSERVED FIELDS CAN ONLY BE COMPOSED OF INDUCED
C CURRENTS. ELECTRIC CURRENTS CAN, VIA INDUCTION, ALSO BE
C GENERATED BY THE SECUALR VARIATION OF EARTH'S AMBIENT
C INTERNAL MAGNETIC FIELD AS WELL AS FROM FLUCTUATIONS OF
C EARTH'S EXTERNAL MAGNETIC FIELD, ESPECIALLY DURING
C MAGNETIC STORMS.
C
C CONSEQUENTLY, SOME CAUTION MUST BE EXCERCISED WHEN
C INTERPRETING THE INVERSION RESULTS OF THIS PROGRAM.
C
C
C
C***********************************************************************
C
C
C COORDINATE SYSTEMS:
C
C ALL INPUT DATA TO THIS ROUTINE HAS BEEN TRANSFORMED
C TO GEODETIC COORDINATES, DENOTED BY (ALT,LAT,LON), OR
C IN RECTANGULAR CORRDINATES, DENOTED BY (X,Y,Z).
C
C EITHER THE GEODETIC OR RECTANGULAR COORDIATE SYSTEM, THE
C X-COMPONENT IS POSITIVE IN THE NORTH DIRECTION, THE
C Y-COMPONENT IS POSITIVE IN THE EAST DIRECTION AND
C THE Z-COORDIATE IS POSITIVE IN THE DOWN (TOWARD THE
C CENTER OF THE EARTH) DIRECTION. THIS SYSTEM IS
C RIGHT HANDED. THE PARAMETER ZS IS THE SURVEY ALTITUDE
C IN THIS RECTANGULAR SYSTEM. CONSEQUENTLY, FOR OBSERVATIONS
C ABOVE THE EARTH'S SURFACE, ZS WILL BE NEGATIVE.
C
C
C EARTH CURVATUIRE EFFECTS ARE AUTOMATICALLY TAKEN INTO
C ACCOUNT, SINCE THE GEODESIC/RECTANGULAR SYSTEM BEING USED
C WILL ROTATE AUTOMATICALLY AS ONE TRAVERSES FROM ONE
C OBSERVATION POIT TO THE NEXT. THUS, THE SOURCE DEPTH AND
C THE COMPUTED MAGNETIZATIONS WILL BE AUTOMATICALLY
C REFERENCED TO THE GEODETIC COORDINATED SYSTEM. THAT IS,
C THERE IS NO NEED FOR EXPLICIT ROTATIONS.
C
C THE COORDINATE ORIGIN IS LOCATED AT THE EARTH'S SURFACE
C DIRECTLY BELOW THE OBSERVATION POINT, REGARDLESS OF WHICH
C OBSERVATION POINT IS BEING EXAMINED.
C
C AN OBSERVATION POINT CORRESPONDS TO THE LOCATION (ALT,LAT,LON)
C FOR WHICH A MAGNETIC FIELD VECTOR AND GRADIENT FIELD TENSOR
C ARE COMPUTED FROM THE SPHERICAL HARMONIC MODEL.
C
C
C**********************************************************************
C
C
C REFERENCES:
C
C QUINN, JOHN M. AND SHIEL, DONALD L.: FORWARD AND
C INVERSE GEOMAGNETIC AND GRAVITATIONAL POTENTIAL
C FIELD MODELING OF THE OCEANIC CRUST, U.S. NAVAL
C OCEANOGRAPHIC OFFICE, GEOMAGNETICS DIVISION
C TECNICAL REPORT (IN PRESS) 1987
C
C TARANTOLA, ALBERT: INVERSE PROBLEM THEORY - METHODS
C FOR DATA FITTING AND MODEL PARAMETER ESTIMATION:
C ELSEVIER (1987) PP 246, 285, 289 AND 371
C
C
C PRESS, WILLIAM H., FLANNERY, BRIAN P., TEUKOLSKY, SAUL A.
C AND VETTERLING, WILLIAM T.: NUMERICAL RECIPES,
C CAMBRIDGE UNIVERSITY PRESS (1986) PP 307-311
C
C POLLACK, HENRY N. AND CHAPMAN, DAVID S.: GLOBAL
C HEAT FLOW REVISITED, PRESENTED AT IUGG
C CONFERENCE, VANCOUVER, CANADA (1987)
C
C QUINN, JOHN M.: CURIE DEPTH DETERMINATION, U.S.
C NAVAL OCEANOGRAPHIC OFFICE, GEOMAGNETICS
C DIVISION MONOGRAPH (UNPUBLISHED) (1989)
C
C ADDITIONAL REFERENCES:
C
C JOSEPH C. CAIN ET AL.; A PROPOSED MODEL FOR THE INTERNATIONAL
C GEOMAGNETIC REFERENCE FIELD - 1965, JOURNAL OF
C GEOMAGNETISM AND GEOELECTRICITY, VOL. 19, NO. 4,
C PP. 335-355 (1967) (SEE APPENDIX)
C
C R. A. LANGEL, THE MAIN FIELD, CHAPTER 4, IN GEOMAGNETISM,
C Vol. 1, Ed. By J. A. JACOBS, ACADEMIC PRESS
C Inc., ORLANDO, FL, USA (1987)
C
C ALFRED J. ZMUDA, ED.; WORLD MAGNETIC SURVEY 1957-1969,
C INTERNATIONAL ASSOCIATION OF GEOMAGNETISM AND
C AERONOMY (IAGA) BULLETIN #28, PP. 186-188 (1971)
C
C DEFENSE MAPPING AGENCY (DMA); DEPARTMENT OF DEFENSE WORLD
C GEODETIC SYSTEM 1984, DMA TECHNICAL REPORT
C #8350.2, WASHINGTO, D.C. (1987)
C
C
C JOHN M. QUINN, RACHEL J. COLEMAN, DONAL L. SHIEL, AND JOHN
C M. NIGRO, THE JOINT US/UK 1995 EPOCH WORLD
C MAGNETIC MODEL, TECHNICAL REPORT #314, U. S.
C NAVAL OCEANOGRAPHIC OFFICE, STENNIS SPACE
C CENTER, MS (1995)
C
C
C***********************************************************************
C
C
C***********************************************************************
C
C
C PARAMETERS:
C
C LONMIN - MINIMUM GEODETIC LONGITUDE OF SURVEY AREA (DEG.)
C LONMAX - MAXIMUM GEODETIC LONGITUDE OF SURVEY AREA (DEG.)
C LATMIN - MINIMUM GEODETIC LATITUDE OF SURVEY AREA (DEG.)
C LATMAX - MAXIMUM GEODETIC LATITUDE OF SURVEY AREA (DEG.)
C LATGRD - GEODETIC SPACING BETWEEN COMPUTATIONAL POINTS (DEG.)
C LONGRD - GEODETIC SPACING BETWEEN COMPUTATIONAL POINTS (DEG.)
C HDR - THE HEADER RECORD OF MAGNETIZATION O/P FILE
C P - ARRAY OF GEOMETRIC PRAMETERS TO BE DETERMINED (KM)
C MAG - ARRAY OF MAGNETIZATIONS DETERMINED FOR EACH PRISM (NT)
C X0 - NORTH COORDINATE OF OBSERVATION POINT (KM)
C Y0 - EAST COORDINATE OF OBSERVATION POINT (KM)
C Z0 - VERT. COORD. (I.E., NEG. ALT) OF OBSERVA
C ALT - GEODETIC ALTITUDE OF OBSERVATION POINT WRT SURFACE (KM)
C EPOCH - DATE OF SURVEY (YRS)
C CDEP - CURRIE DEPTH (KM)
C Q - HEAT FLOW (MILIWATTS/M**2)
C Q0 - MAXIMUM HEATFLOW (MILIWATTS/M**2)
C Q1 - HEATFLOW BIAS (MILIWATTS/M**2)
C HETALT - ALTITUDE/DEPTH (+/-) OF HEATFLOW COMUTATION (KM)
C RIGDEP - AVERAGE DEPTH TO SPREADING ZONE RIDGE CREST (KM)
C HETCON - CURIE DEPTH/HEAT FLOW COUPLING CONSTANT (KM-MW/M**2)
C CNTAM - NANOTESLA TO AMPS/METER CONVERSION FACTOR (1/400*PI)
C LAMDA - 3X3 GEOMETRIC PRISM MATRIX (DIMENSIONLESS)
C LAMDAI - INVERSE OF LAMBDA MATRIX (DIMENSIONLESS)
C POSIN - POSITION VECTOR (LAT,LON,ALT)
C POSOUT - POSITION VECTOR (LAT,LON,ALT) - EQUALS POSIN
C LAMX - LENGTH OF PRISM IN NORTH DIRECTION (KM)
C LAMY - LENGTH OF PRISM IN EAST DIRECTION (KM)
C LAMZ - LENGTH OF PRISM IN VERTICAL DIRECTION (KM)
C MAG - INTERMEDIATE MAGNETIZATION VECTOR (nT)
C MAGX - FINAL GLOBAL ARRAY OF NORTH MAGNETIZATIONS (nT)
C MAGY - FINAL GLOBAL ARRAY OF EAST MAGNETIZATIONS (A/M)
C MAGZ - FINAL GLOBAL ARRAY OF VERT. DOEN MAGNETIZATIONS (nT)
C SDEP - FINAL GLOBAL ARRAY OF PRISM SOURCE DEPTHS (KM)
C TIM - DATE IN DECIMAL YEARS FOR MODEL COMPUTATION (YRS)
C VI - MAGNETIC POTENTIAL OF INTERNAL SH MODEL (nT*KM)
C BI - MAGNETIC FIELD VECTOR OF INTERNAL SH MODEL (nT)
C BFI - MAGNETIC TOTAL INTENSITY OF INTERNAL SH MODEL (nT)
C GI - MAG. FIELD GRADIENT TENSOR OF INTERNAL MODEL (nT)/KM)
C GFI - TOTAL MAG. GRADIENT INTENSITY OF INTERNAL MODEL (nT/KM)
C VE - MAGNETIC POTENTIAL OF EXTERNAL MODEL (nT*KM)
C BE - MAGNETIC FIELD VECTOR OF EXTERNAL SH MODEL (nT)
C BFE - MAGNETIC TOTAL INTENSITY OF EXTERNAL SH MODEL (nT)
168
C GE - MAGNETIC GRADIENT TENSOR OF EXTERNAL MODEL (nT/KM)
C GFE - TOTAL MAG. GRADIENT INTENSITY OF INTERNAL MODEL (nT/KM)
C GIP - MAGNETIC GRADIENT TENSOR DUE TO PRISM nT/KM)
C CNTAM - COMVERSION FACTOR NANOTESLAS (NT) TO AMP/METER (A/M)
C MAG - MAGNETIZATION VECTOR IN CHOSEN COORDINATE SYSTEM (A/M)
C MAGX - NORTH COMPONENT MAGNETIZATION GRID (A/M)
C MAGY - EAST COMPONENT MAGNETIZATION GRID (A/M)
C MAGZ - VERTICALLY DOWN COMPONENT MAGNETIZATION GRID (A/M)
C MX - NORTH COMPONENT MAGNETIZATION (A/M)
C MY - EAST COMPONENT MAGNETIZATION (A/M)
C MZ - VERTICALLY DOWN COMPONENT MAGNETIZATION (A/M)
C
C
C**********************************************************************
C
C
C BASIC MATHEMATICAL RELATIONS FOR RECTANGULAR PRISMS:
C
C B(m)=LAMDA(m,l)*M(l) (1)
C
C ===> M(l)=LAMDAI(l,k)*B(k) (2)
C
C
C G(m,n)=DLAMDA(m,l,n)*M(l) (3)
C
C ===> G(m,n)=DLAMDA(m,l,n)*LAMDAI(l.k)*B(k) (4)
C
C
C WHERE SUMMATION IS IMPLIED OVER REPEATED INDICES
C (i,j,k,l,m,n =1,2,3)
C
C
C
C THE CHI-SQUARE FUNTION IS FORMED FROM THE SQUARE OF THE
C DIFFERENCE BETWEEN EQ. (4), AND THE OBSERVED MAGNETIC FIELD
C GRADIENT TENSOR, SUMMING OVER ALL NINE GRADIENT MAGNETIC
C TENSOR COMPONENTS AS INDICATED IN EQ. (5) BELOW:
C
C
C X**2 = [Gp(m,n) - Go(m,n)]**2 (5)
C
C
C WHERE THE SUBSCRIPT P REFERS TO THE PRISM AND THE SUBSCRIPT o
C REFERS TO THE MAGNETIC FIELD OBSERVATIONS COMPUTED BY THE
C MAGPOT SUBROUTINE USING THE MF3 MODEL BASED ON CHAMP
C SATELLITE OBSERVATIONS.
C
C ACTUALLY, ONLY SEVEN COMPONENTS ARE USED. THIS IS BECAUSE
C G(1,2)=G(2,1)=0 FOR AN OBSERVATION POINT CENTERED DIRCTLY OVER
C A PRISM. SO THE GRADIENT DERIVED FROM THE MODEL IS SET TO
C ZERO SINCE IN SOME CASES, COMPUTER ROUNDOFF ERRORS WILL GIVE
C SPURIOUS RESULTS IF G(1,2),G(2,1) ARE EITHER TOO SMALL OR TOO
C NOISY.
C
C
C HERE. Bo AND Go ARE THE MAGNETIC FIELD VECTOR AND GRADIENT
C TENSOR COMPUTED FROM THE SPHERICAL HARMONIC MODEL, WHILE
C LAMDA, DLAMDA, AND LAMDAI ARE THE GEOMETRIC ARRAY: ITS
C GRADIENT ARRAY, AND ITS INVERSE ARRAY ASSOCIATED WITH THE
C PRISM, WHICH ARE USE TO COMPUTE THE PRISM GENERATED MAGNETIC
C FIELD, Bp, AND MAGNETIC GRADEIENT TENSOR FIELD, Gp.
C
C FROM THE ABOVE DISCUSSION, ONE SEES THAT THERE ARE FOUR
C
C BASIC STEPS NEEDED TO GENEARATE A SOLUTION FOR A GRIDDED
C SURVEY AREA IN GEODETIC OR SPHERICAL COORDINTAES:
C
C (1) DETERMINE THE PRISM DIMENSIONS USING THE SURFACE AND CURIE
C DEPTH AS LATERAL CONSTRAINTS, AND OCEAN DRILLING PROGRAM
C RESULTS AS A THICKNESS CONSTRAINT.
C
C (2) SOLVE EQ. (5), THE CHI-SQUARE EQUATION, FOR THE PRISM
C SOURC EDEPTH.
C
C (3) THEN SOLVE EQ.(2) FOR THE PRISM MAGNETIZATION VECTOR.
C
C (4) REPEAT THIS PROCESS FOR ALL GRID POINTS IN THE SELECTED AREA.
C
C
C***********************************************************************
C
C
C
C
C SUBROUTINES:
C
C MAGPOT - COMPUTES THE INTERNAL AND EXTERNAL PARTS OF THE
C MAGNETIC POTENTIAL, THE THREE MAGNETIC VECTOR
C COMPONENTS, AND THE NINE MAGNETIC GRADIENT
C COMPONENTS. ADDITIONALLY, THE PROGRAM COMPUTES
C THREE MAGNETIC SCALAR INVARIANTS, THAT ARE
C PARTICULARLY USEFUL IN THE ANALYSIS OF THE CRUSTAL
C MAGNETIC FIELD. THE USER MAY SPECIFIY THE LOWER AND
C UPPER SPHERICAL HARMONIC DEGREES OF THE MODEL TO BE
C INCLUDED IN THE CALCULATIONS FOR BOTH THE INTERNAL
C
C AND EXTERNAL PARTS OF THE FIELD, THEREBY PERMITTING
C ANALYSIS OF VRIOUS HARMONIC CONTRIBUTIONS TO THE
C MAGNETIC FIELD. THE INTERNAL AND EXTERNAL PARTS OF
C THE MAGNETIC FIELD ARE RETURNED SEPARATELY. THIS
C SUBROUTINE IS INITIALIZED ONCE THROUGH THE MAIN
C ENTRY POINT MAGPOT( ), WHILE ALL SUBSEQUENT
C CALLS ARE MADE THROUGH THE ENTRY POINT MAGPT1( )
C
C LAMDLA - COMPUTES THE GEOMETRIC ARRAYS: LAMDA AND ITS
C GRADIENT GENERATED BY A RECTANGULAR PRISM.
C SEE THE QUINN-SHIEL REFERENCE REPORT- A UNIFIED
C APPROACH .......AND ITS COMPANION REPORT APPLYING
C THE PRISM TECHNIQUE TO THE JUAN DE FUCA REGION
C NEAR VANCOUVER AND SEATTLE ...
C
C MATINV - COMPUTES THE INVERSE OF THE 3 X 3 LAMDA MATRIX
C
C
C***********************************************************************
C
C
C MAGPOT INPUT PARAMETERS:
C
C TIME - DECIMAL YEAR OF OUTPUT CHARTS (E.G., 2001.7321156)
C EPOCH - REFERENCE YEAR OF MAIN FIELD MODEL (E.G., 2000.000)
C ALT - ALTITUDE IN KILOMETERS (KM)
C NIMIN - MINUMUM INTERNAL SPHERICAL HARMONIC DEGREE OF MODEL
C NIMAX - MAXIMUM INTERNAL SPHERICAL HARMONIC DEGREE OF MODEL
C NEMIN - MINUMUM EXTERNAL SPHERICAL HARMONIC DEGREE OF MODEL
C NEMAX - MAXIMUM EXTERNAL SPHERICAL HARMONIC DEGREE OF
C ISLCT - INPUT/OUTPUT REFERENCE COORDINATES USED
C
C 1 = GEODETIC IN/GEODETIC OUT
C 2 = GEODETIC IN/SPHERICAL OUT
C 3 = SPHERICAL IN/SPHERICAL OUT
C JSLCT - DESIRED MAGNETIC FIELD TYPE
C 1 = INTERNAL
C 2 = EXTERNAL
C 3 = COMBINED (INTERNAL + EXTERNAL)
C
C
C OUTPUT PARAMETERS
C
C 1 = V MAGNETIC POTENTIAL
C 2 = BX NORTH MAGNETIC COMPONENT
C 3 = BY EAST MAGNETIC COMPONENT
C 4 = BF TOTAL MAGNETIC INTENSITY
C 5 = BZ VERTICALLY-DOWN MAGNETIC COMPONENT
C 6 = GXX NORTH_NORTH MAGNETIC GRADIENT COMPONENT
C 7 = GXY NORTH_EAST MAGNETIC GRADIENT COMPONENT
C 8 = GXZ NORTH_VERTICALLY-DOWN MAG. GRADIENT COMP.
C 9 = BYY EAST_EAST MAGNETIC GRADIENT COMPONENT
C 10 = BYZ EAST_VERTIALLY-DOWN MAG. GRADIENT COMP.
C 11 = BZZ VERT-DWN_VERT-DWN MAGNETIC GRADIENT COMP.
C 12 = GZ TOTAL MAGNETIC GRADIENT INTENSITY
C 13 = P MAGNETIC SCALAR INVARIANT (BF/V)
C 14 = Q MAGNETIC SCALAR INVARIANT (GF/BF)
C 15 = R MAGNETIC SCALAR INVARIANT (GF/V)
C
C
C***********************************************************************
C
C
C
C MAGPOT OUTPUT PARAMETERS:
C
C
C GEODETIC COORDINATES:
C INTERNAL
C
C VI - MAGNETIC POTENTIAL (NT-KM)
C BXI - X-NORTH MAGNETIC COMPONENT (NT)
C BYI - Y-EAST MAGETIC COMPONETN (NT)
C BZI - Z-VERTICALLY DOWN MAGNETIC COMPONENT (NT)
C GIXX - NORTH-NORTH MAGNETIC GRADIENT COMPONENT (NT/KM)
C GIXY - NORTH-EAST MAGNETIC GRADIENT COMPONENT
C GIXZ - NORTH-VERT. DOWN MAGNETIC GRADIENT COMPONENT (NT/KM)
C GIYY - EAST-EAST MAGNETIC GRADIENT COMPONENT (NT/KM)
C GIYZ - EAST-VERT. DOWN MAGNETIC GRADIENT COMPONENT (NT/KM)
C GIZZ - VERT. DOWN-VERT. DOWN MAG. GRADIENT COMP. (NT/KM)
C
C EXTERNAL
C
C VE - MAGNETIC POTENTIAL (NT-KM)
C BXE - X-NORTH MAGNETIC COMPONENT (NT)
C BYE - Y-EAST MAGETIC COMPONETN (NT)
C BZE - Z-VERTICALLY DOWN MAGNETIC COMPONENT (NT)
C GEXX - NORTH-NORTH MAGNETIC GRADIENT COMPONENT (NT/KM)
C GEXY - NORTH-EAST MAGNETIC GRADIENT COMPONENT
C GEXZ - NORTH-VERT. DOWN MAGNETIC GRADIENT COMPONENT (NT/KM)
C GEYY - EAST-EAST MAGNETIC GRADIENT COMPONENT (NT/KM)
C GEYZ - EAST-VERT. DOWN MAGNETIC GRADIENT COMPONENT (NT/KM)
C GEZZ - VERT. DOWN-VERT. DOWN MAG. GRADIENT COMP. (NT/KM)
C
C COMBINED (INTERNAL + EXTERNAL)
C
C
C V - MAGNETIC POTENTIAL (NT-KM)
C BX - X-NORTH MAGNETIC COMPONENT (NT)
C BY - Y-EAST MAGETIC COMPONETN (NT)
C BZ - Z-VERTICALLY DOWN MAGNETIC COMPONENT (NT)
C GXX - NORTH-NORTH MAGNETIC GRADIENT COMPONENT (NT/KM)
C GXY - NORTH-EAST MAGNETIC GRADIENT COMPONENT
C GXZ - NORTH-VERT. DOWN MAGNETIC GRADIENT COMPONENT (NT/KM)
C GYY - EAST-EAST MAGNETIC GRADIENT COMPONENT (NT/KM)
C GYZ - EAST-VERT. DOWN MAGNETIC GRADIENT COMPONENT (NT/KM)
C GZZ - VERT. DOWN-VERT. DOWN MAG. GRADIENT COMP. (NT/KM)
C
C
C SPHERICAL COORDINATES:
C INTERNAL
C
C VI - MAGNETIC POTENTIAL (NT-KM)
C BRI - RADIAL VECTOR MAG. COMP. (NT)
C BTI - THETA VECTOR MAG. COMP. (NT)
C BPI - PHI VECTOR MAG. COMP. (NT)
C GIRR - RADIAL-RADIAL GRADIENT MAG. COMP. (NT/KM)
C GIRT - RADIAL-THETA GRADIENT MAG. COMP. (NT/KM)
C GIRP - RADIAL-PHI GRADIENT MAG. COMP. (NT/KM)
C GITR - THETA-RADIAL GRADIENT MAG. COMP. (NT/KM)
C GITT - THETA-THETA GRADIENT MAG. COMP. (NT/KM)
C GITP - THETA-PHI GRADIENT MAG. COMP. (NT/KM)
C GIPR - PHI-RADIAL GRADIENT MAG. COMP. (NT/KM)
C GIPT - PHI-THETA GRADIENT MAG. COMP. (NT/KM)
C GIPP - PHI-PHI GRADIENT MAG. COMP. (NT/KM)
C
C EXTERNAL
C
C
C VE - MAGNETIC POTENTIAL (NT-KM)
C BRE - RADIAL VECTOR MAG. COMP. (NT)
C BTE - THETA VECTOR MAG. COMP. (NT)
C BPE - PHI VECTOR MAG. COMP. (NT)
C GERR - RADIAL-RADIAL GRADIENT MAG. COMP. (NT/KM)
C GERT - RADIAL-THETA GRADIENT MAG. COMP. (NT/KM)
C GERP - RADIAL-PHI GRADIENT MAG. COMP. (NT/KM)
C GETR - THETA-RADIAL GRADIENT MAG. COMP. (NT/KM)
C GETT - THETA-THETA GRADIENT MAG. COMP. (NT/KM)
C GETP - THETA-PHI GRADIENT MAG. COMP. (NT/KM)
C GEPR - PHI-RADIAL GRADIENT MAG. COMP. (NT/KM)
C GEPT - PHI-THETA GRADIENT MAG. COMP. (NT/KM)
C GEPP - PHI-PHI GRADIENT MAG. COMP. (NT/KM)
C
C COMBINED (INTERNAL + EXTERNAL)
C
C V - MAGNETIC POTENTIAL (NT-KM)
C BR - RADIAL VECTOR MAG. COMP. (NT)
C BT - THETA VECTOR MAG. COMP. (NT)
C BP - PHI VECTOR MAG. COMP. (NT)
C GRR - RADIAL-RADIAL GRAD. COMP. (NT/KM)
C GRT - RADIAL-THETA GRAD. COMP. (NT/KM)
C GRP - RADIAL-PHI GRAD. COMP. (NT/KM)
C GTR - THETA-RADIAL GRAD. COMP. (NT/KM)
C GTT - THETA-THETA GRAD. COMP. (NT/KM)
C GTP - THETA-PHI GRAD. COMP. (NT/KM)
C GPR - PHI-RADIAL GRAD. COMP. (NT/KM)
C GPT - PHI-THETA GRAD. COMP. (NT/KM)
C GPP - PHI-PHI GRAD. COMP. (NT/KM)
C
C
C
C SCALAR INVARIANT PARAMTERS
C
C INTERNAL
C
C VI - MAGNETIC POTENTIAL (NT-KM)
C BFI - MAGNETIC TOTAL INTENSITY (NT)
C GFI - TOTAL MAGNETIC GRADIENT INTENSITY (NT/KM)
C PI - SCALAR MAG. INVARIANT #1 - BFI/VI (1/KM)
C QI - SCALAR MAG. INVARIANT #2 - GFI/BFI (1/KM)
C RI - SCALAR MAG. INVARIANT #3 - GFI/VI (1/KM**2)
C
C EXTERNAL
C
C VE - MAGNETIC POTENTIAL (NT-KM)
C BFE - MAGNETIC TOTAL INTENSITY (NT)
C GFE - TOTAL MAGNETIC GRADIENT INTENSITY (NT/KM)
C PE - SCALAR MAG. INVARIANT #1 - BFE/VE (1/KM)
C QE - SCALAR MAG. INVARIANT #2 - GFE/BFE (1/KM)
C RE - SCALAR MAG. INVARIANT #3 - GFE/VE (1/KM**2)
C
C COMBINED INT./EXT.
C
C V - MAGNETIC POTENTIAL (NT-KM)
C BF - TOTAL MAGNETIC INTENSITY (NT)
C GF - TOTAL MAGNETIC GRADIENT (NT/KM)
C P - SCALAR MAG. INVARIANT #1 - BF/V (1/KM)
C Q - SCALAR MAG. INVARIANT #2 - GF/BF (1/KM)
C R - SCALAR MAG. INVARIANT #3 - GF/V (1/KM**2)
C
C
C
C***********************************************************************
C
C
C UNITS: GAUSSIAN UNTS ARE EMPLOYED IN THIS PROGRAM DUE TO THE
C SYSTEM OF UNITS USED TO DERIVE THE EQUATIONS. THIS
C SYSTEM IS BASICALLY A CGS SYSTEM OF UNITS EXCEPT THAT
C ELECTROMAGNETIC UNITS ARE EXPRESSED IN A MIXTURE OF
C ELECTROSTATIC UNITS (ESU) AND ELECTROMAGNETIC UNITS (EMU).
C IN PARTICULAR, ELECTRIC PARAMETERS USE ESU UNITS AND
C MAGNETIC PARAMETERS USE EMU UNITS. THUS, MAGNETIC FIELDS
C ARE EXPRESSEN IN UNITS OF GAUSS. HOWEVER, CONVENTIONAL
C MAGNETIC FIELD OBSERVATIONS ARE TYPICALLY GIVEN IN UNTIS
C OF nanoTeslas (nT). ON THE OTHER HAND, MAGNETIZATION
C UNITS IN THE GAUSSIAN SYSTEM ARE IN TERMS OF EMU/CM***3.
C
C ALSO, DUE TO THE SIZE OF THE SURVEY AREAS INVOLVED, THIS
C PROGRAM USES Kilometers (Km) RATHER THAN METERS (M)OR
C CENTIMETERS (CM) AS THE BASIC UNIT OF LENGTH.
C
C SO, ONE MUST BE CAREFUL AS ONE CONVERTS FROM OBSERVATIONAL
C UNITS (IE.: nT, Km) TO GAUSSIAN UNITS (I.E.: GAUSS, CM,
C EMU/CM**3) AND THEN BACK TO SI UNITS (I.E.: TESLAS,
C METERS, AMPS/METER)
C
C THE COMPUTER ALGORITHM USES THE FOLLOWING UNITS:
C
C MAGPOT COMPUTES THE MAGNETIC FIELD IN nanoTeslas (nT)
C MAGPOT COMPUTES THE MAGNETIC GRADIENTS IN (nT/Km)
C LAMDLA COMPUTES THE DIMENSIONLESS LAMDA MATRIX FOR A PRISM
C PRISM GENERATED MAGNETIC GRADIENTS HAVE UNITS OF (nT/Km)
C
C
C THUS EQUATIONS (1) AND (2) ABOVE INDICATE THAT THE
C MAGNETIZATION WILL BE COMPUTED IN TERMS OF nT. TO CONVERT
C THIS TO EMU/CM**3, NOTE THE FOLLOWING CONVERSIONS:
C
C OBSERVATIONAL UNITS TO GAUSSIAN UNITS
C
C 1 nT = 10**-5 GAUSS = 10**-5 EMU/CM**3
C
C WHERE:
C
C 1 GAUSS = 1 EMU/CM**3
C
C 1 EMU = 1 MAGNETIC MOMENT
C
C
C GAUSSIAN UNITS TO SI UNITS
C
C 1 EMU/CM**3 = 10**3 AMPS/M
C
C
C THUS:
C
C 1 nT = 10**-2 AMPS/M
C
C SO, 10**-2 IS THE CONVERSION FACTOR USED TO PRESENT
C THE MAGNETIZATION IN SI UNITS.
C
C
C TO VERIFY THESE RELATIONSHIPS, THE FOLLOWING
C REFERENCES ARE HANDY:
C
C
C
C RONALD T. MERRILL AND MICHAEL W. McELHINNY; THE EARTH'S
C MAGNETIC FIELD, TABLE A.1 PP. 354, ACADEMIC PRESS,
C LONDON & NEW YORK (1983)
C
C GEORGE BACKUS, ROBERT PARKER, AND CATHERINE CONSTABLE;
C FOUNDATIONS OF GEOMAGNETISM, TABLE ON PAGE 28,
C CAMBRIDGE UNIERSITY PRESS, CAMBRIDGE, ENGLAND
C (1996)
C
C J. D. JACKSON, CLASSICAL ELECTRODYNAMICS 2nd Ed., APPENDIX
C ON UNITS AND DIMENSIONS, TABLE 4, PP. 820, JOHN WILEY
C AND SONS, NEW YORK (1975)
C
C
C***********************************************************************
C
C
C NOTE: THE MAGNETIC GRADIENT IS SYMMETRIC WITH ZERO TRACE.
C
C
C***********************************************************************
C
C
C NOTE: MAGPOT DOES NOT ANALYTICALLY COMPUTE THE MAGNETIC FIELD
C PARMETERS FOR THE SPECIAL CASE OF THE NORTH AND SOUTH
C GEOGRAPHIC POLES. THIS IS BECAUSE YOU HAVE TO BE VERY
C CAREFUL ABOUT HOW THE LEGENDRE FUNCTIONS ARE COMPUTED
C SINCE THE SIN(THETA) IS ZERO AT THE POLES: THETA=0 AND
C THETA=180. THIS CAN BE CIRCUMVENTED BY SETTING THETA
C TO BE A VERY SMALL NUMBER CLOSE TO ZERO, IN LIEU OF
C
C PROGRAMMING THE SPECIAL CASES. SO, MAGPOT AS IT
C CURRENTLY STANDS IS VALID EVERY WHERE ON THE EARTH
C EXCEPT AT TWO POINTS: THE NORTH AND SOUTH POLES, WHERE
C IT IS ONLY APPROXIMATELY CORRECT.
C
C
C***********************************************************************
C
C
C NOTE: P(N,M), Q(N,M), AND R(N,M) ARE SCALAR INVARIANTS THAT HAVE
C INTERESTING PROPERTIES. IN PARTICULAR THEY DEFINE THE
C EXACT LOCATION AND LATERAL BOUNDARIES OF LITHOSPHERIC
C MAGNETIC ANOMALIES. THIS IS ESPECIALLY TRUE OF Q(N,M).
C
C
C***********************************************************************
C
C
C NOTE: THIS PROGRAM ASSUMES THAT THE MAIN FIELD (IE., INTERNAL)
C MODEL COEFFICIENTS ARE SEPARATED FROM THE EXTERNAL FIELD
C MODEL COEFFICIENTS BY A SET OF NINES (I.E., 999) IN THE
C DEGREE COLUMNS OF THE MODEL COEFFICIENT FILE.
C
C THE END OF THE EXTERNAL FIELD MODEL COEFFICIENTS IS DENOTED
C BY TWO CONSECUTIVE SETS OF NINES IN THE DEGREE COLUMNS.
C
C THE INPUT MODEL CONSISTS OF ONE RECORD CONTAINING THE MODEL
C REFERENCE DATE AND MODEL IDENTITY, FORMAT:
C
C 1X,F14.9,A60
C
C
C FOLLOWED BY MANY RECORDS CONTAINING THE MODEL COEFFICIENTS
C WHICH ARE EXPECTED TO BE IN THE FOLLOWING FIXED FORMAT:
C
C 2(1X,I3),4(2X,F14.7)
C
C CORRESPONDING TO THE INTERNAL PARAMETERS:
C
C N,M,G(N,M),H(N,M),GT(N,M),HT(N,M)
C
C AND TO THE EXTERNAL PARAMETERS:
C
C N,M,Q(N,M),S(N,M),QT(N,M),ST(N,M)
C
C WHERE GT, HT, QT, AND ST ARE SECULAR VARIATION PARAMETERS.
C
C GENERALLY, GT(N,M) AND HT(N,M) ARE ZERO ABOVE SPHERICAL
C HARMONIC DEGREE AND ORDER 12, WHILE QT(N,M), AND ST(N,M)
C ARE TYPICALLY SET TO ZERO FOR ALL DEGREES AND ORDERS. THAT
C IS, THEY ARE GENERALLY NOT COMPUTED AND ARE CONSIDERED
C SMALL.
C
C
C NOTE: THE SYMBOLS ABOVE ASSSOCIATED WITH THE MODEL ARE DIFFERENT
C FROM THE SCALAR INVARIANTS LISTED FURTHER ABOVE.
C
C
C***********************************************************************
C
C
C LAMDLA INPUT PARAMETERS:
C
C
C X0 - NORTH DISTANCE OF OBSERVATION POINT WRT ORIGIN (KM)
C Y0 - EAST DISTANCE OF OBSERVATION POINT WRT ORIGIN (KM)
C Z0 - VERTICALLY DOWN DISTANCE OF OBSERVATION POINT WRT
C ORIGIN ON ELLIPSOID SURFACE (I.E., NEG. ALT) (KM)
C XB - NORTH DISTANCE OF PRISM CENTER WRT ORIGIN (KM)
C YB - EAST DISTANCE OF PRISM CENTER WRT ORIGIN (KM)
C ZB - DEPTH (POSITIVE) OF PRISM CENTER WRT ORIGIN ON ELLIPSOID
C SURFACE (KM)
C LAMX - LENGTH OF PRISM IN NORTH-SOUTH DIRECTION (KM)
C LAMY - LENGTH OF PRISM IN EAST-WEST DIRECTION (KM)
C LAMZ - LENGTH OF PRISM IN VERTICAL DIRECTION (KM)
C
C
C
C NOTE: THE ORIGIN IN THIS PROGRAM IS LOCATED AT THE EARTH'S
C SURFACE DIRECTLY BELOW THE OBSERVATION POINT BEING
C CONSIDERED. THUS, XB=YB=X0=Y0=0, WHILE Z0 <= 0 AND
C ZB => 0, SINCE THE PRISM IS ALWAYS BELOW THE ORIGIN AND
C OBSERVATION POINT IS ALWAYS ABOVE THE ORIGIN, AND SINCE
C POSITIVE IS DOWN IN OUR COORDINATE SYSTEM.
C
C
C***********************************************************************
C
C
C NOTE: THIS PROGRAM WAS WRITTEN IN FORTRAN 77 USING MICROSOFT'S
C FORTRAN DEVELOPER STUDIO IN WINDOWS XP PROFESSIONAL.
C
C THIS CODE CAN BE TRANSFERRED DIRECTLY TO ANY UNIX MACHINE
C OPERATING FORTRAN 77. HOWEVER, YOU WILL NEED TO ADD
C THE FOLLOWING INCLUDE STATEMENTS, APPROPRIATLY REVISED
C
C TO INDICATE THE CORRECT PATH TO THE SUBROUTINES AND
C ALSO CHANGING THE .for TO .ftn OR SOMETHING LIKE THAT.
C FURTHERMORE, THE INCLUDE STATEMENTS ARE TO BE LOCATED
C WITHIN THIS PROGRAM JUST BELOW AS INDICATED (PATH CHANGES
C NECESSARY):
C
C
C***********************************************************************
C
C
c INCLUDE "C:\MSDEV\Projects\Source Depth\magpot.ftn"
c INCLUDE "C:\MSDEV\Projects\Source Depth\lamdla.ftn"
c INCLUDE "C:\MSDEV\Projects\Source Depth\matinv.ftn"
C
C
C***********************************************************************
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
CHARACTER*60 MODEL
C
DOUBLE PRECISION LAMDA(3,3),DLAMDA(3,3,3),LAMDAI(3,3)
DOUBLE PRECISION LAMX,LAMY,LAMZ
DOUBLE PRECISION LATGRD,LONGRD
DOUBLE PRECISION MAG(3)
DOUBLE PRECISION MAGX(-360:360,-720:720)
DOUBLE PRECISION MAGY(-360:360,-720:720)
DOUBLE PRECISION MAGZ(-360:360,-720:720)
DOUBLE PRECISION MAGT(-360:360,-720:720)
DOUBLE PRECISION MX,MY,MZ,MT
184
C
DIMENSION CI0(0:100,0:100),CIT0(0:100,0:100)
DIMENSION CE0(0:100,0:100),CET0(0:100,0:100)
DIMENSION POSIN(3),POSOUT(3),SDEP(-360:360,-720:720)
DIMENSION BI(3),GI(3,3),BE(3),GE(3,3),GIP(3,3)
C
COMMON /PRISM/ LAMX,LAMY,LAMZ,XB,YB,ZB,ALPHA,BETA,GAMA
COMMON /LIMITS/ NIMN,NIMX,NEMN,NEMX
COMMON /COEFFS/ CI0,CIT0,CE0,CET0
COMMON /MDLEPC/ MODEL,EPOC
C
C
C***********************************************************************
C
C
C INITIALIZE CONSTANTS
C
C
NORDER=3
XB=0.000D0
YB=0.000D0
X0=0.000D0
Y0=0.000D0
ALPHA=0.000D0
BETA=0.000D0
GAMA=0.000D0
NIMAX=0
NEMAX=0
EPS=1.0D9
R0=6371.2D0
CNTAM=1.0000D0/100.0000D0
C
C
C
PRINT *, ' '
PRINT *, ' ENTER NORTH, EAST, VERTICALLY DOWN PRISM DIMENSIONS'
PRINT *, ' IN DECIMAL KILOMETERS => 2*ALT'
PRINT *, ' LAMX,LAMY,LAMZ'
READ *, LAMX,LAMY,LAMZ
C
C
PRINT *, ' '
PRINT *, ' INTERMEDIATE RESULTS TO THE COMPUTER SCREEN?'
PRINT *, ' 1 = YES 0 = NO'
READ *, IPRINT
C
PRINT *, ' '
PRINT *, ' ENTER ALTITUDE OF OBSERVATION POINTS IN kILOMETERS'
PRINT *, ' RECOMMENDED VALUE 50 Km'
PRINT *, ' HIGHER ALTITUDES RUN INTO COMPUTER LIMITATIONS'
PRINT *, ' ESPECIALLY FOR THE HIGHER HARMONICS => 61'
PRINT *, ' THIS IS BECAUSE OF THE COMPUTER WORD SIZE'
PRINT *, ' HERE, ONLY A 32 BIT PC WAS USED'
PRINT *, ' '
PRINT *, ' '
READ *, ALT
Z0=-ALT
POSIN(1)=ALT
C
PRINT *, ' '
PRINT *, ' ENTER SHALLOWEST POSSIBLE SOURCE DEPTH: INTGER KM'
PRINT *, ' RANGE: => 0'
READ *, MINDEP
C
MINDEP=DFLOAT(MINDEP)+LAMZ/2.000D0
MINDEP=MINDEP*10
C
PRINT *, ' '
PRINT *, ' ENTER MAXIMUM POSSIBLE SOURCE DEPTH: INTEGER KM'
PRINT *, ' RANGE: <= 999'
READ *, MAXDEP
MAXDEP=MAXDEP*10
PRINT *, ' '
PRINT *, ' ENTER DEPTH STEP INTERVAL INTEGER: 1 - 10'
PRINT *, ' 1 = 0.1 KM'
PRINT *, ' 2 = 0.2 KM'
PRINT *, ' ETC.'
READ *, ISTEP
C
PRINT *, ' '
PRINT *, ' ENTER CHART LIMITS IN INTEGER DEGREES'
PRINT *, ' WITHIN THE FOLLOWING LAT X LON LIMITS:'
PRINT *, ' -90 <= LAT <= +90'
PRINT *, ' -180 <= LON <= +180'
PRINT *, ' LATMIN,LATMAX,LONMIN,LONMAX'
READ *, LATMIN,LATMAX,LONMIN,LONMAX
RLATMN=DFLOAT(LATMIN)
RLATMX=DFLOAT(LATMAX)
RLONMN=DFLOAT(LONMIN)
RLONMX=DFLOAT(LONMAX)
C
PRINT *, ' '
PRINT *, ' ENTER LAT/LON GRID INTERVALS IN INTEGER DEGREES'
PRINT *, ' MINIMUM GRID SPACING IS 1 DEGREE'
PRINT *, ' LATGRD,LONGRD'
187
READ *, LATGRD,LONGRD
RLTGRD=DFLOAT(LATGRD)
RLNGRD=DFLOAT(LONGRD)
C
C
C***********************************************************************
C
C
C READ SPHERICAL HARMONIC MAGNETIC FIELD MODEL COEFFICIENTS
C AND INITIALIZE THE MAGPOT SUBROUTINE
C
C
C***********************************************************************
C
C
C THIS SECTION READS A SET OF SPHERICAL HARMONIC
C MAGNETIC POTIENTIAL MODEL COEFFICIENTS IN SCHMIDT
C QUASI-NORMALIZED FORM AND PASSES THEM TO THE MAGPOT
C SUBROUTINE WHICH COMPUTES THE MAGNETIC POTENTIAL, THE
C MAGNETIC VECTOR COMPONENTS, AND THE MAGNETIC GRADIENT
C COMPONENTS AT ANY USER SPECIFIED TIME, ALTITUDE (DEPTH),
C AND LOCATION.
C
C
C***********************************************************************
C
C
OPEN(UNIT=7,FILE='C:\mf3\mf3_2001.cof',FORM='FORMATTED')
OPEN(UNIT=11,FILE='D:\src-depth.src',FORM='FORMATTED')
C
C
C
PRINT *, ' '
PRINT *, ' '
PRINT *, ' ENTER DESIRED INPUT/OUTPUT COORDINATE SYSTEMS'
PRINT *, ' 1 = GEODETIC IN/GEODETIC OUT'
PRINT *, ' 2 = GEODETIC IN/SPHERICAL OUT'
PRINT *, ' 3 = SPHERICAL IN/SPHERICAL OUT'
READ *, ISLCT
C
C
PRINT *, ' '
PRINT *, ' '
PRINT *, ' ENTER DESIRED PORTION OF THE FIELD'
PRINT *, ' 1 = INTERNAL'
PRINT *, ' 2 = EXTERNAL'
PRINT *, ' 3 = COMBINED INTERNAL/EXTERNAL'
READ *, JSLCT
C
C
PRINT *, ' '
PRINT *, ' '
PRINT *, ' MAGNETIC FIELD PARAMETER COMPUTED ARE:'
PRINT *, ' '
PRINT *, ' 1 = V POTENTIAL'
PRINT *, ' 2 = BX-NORTH OR BT-THETA (SOUTH)'
PRINT *, ' 3 = BY-EAST OR BP-PHI (EAST)'
PRINT *, ' 4 = BZ-VERT DOWN OR BR-RADIAL (OUTWARD)'
PRINT *, ' 5 = BF TOTAL VECTOR INTENSITY'
PRINT *, ' 6 = GXX NORTH_NORTH OR GTT-THETA_THETA'
PRINT *, ' 7 = GXY NORTH_EAST OR GTP-THETA_PHI'
PRINT *, ' 8 = GXZ NORTH_VERT. DOWN OR GTR-THETA_RADIAL'
PRINT *, ' 9 = GYY EAST_EAST OR GTT-THETA_THETA'
C
PRINT *, ' 10 = GYZ EAST_VERT. DOWN OR GPR-PHI_RADIAL'
PRINT *, ' 11 = GZZ VERT. DOWN_VERT. DOWN OR GRR-RADIAL_RADIAL'
PRINT *, ' 12 = GF TOTAL GRADIENT INTENSITY'
PRINT *, ' 13 = P SCALAR INVARIANT'
PRINT *, ' 14 = Q SCALAR INVARIANT'
PRINT *, ' 15 = R SCALAR INVARIANT'
C
C
PRINT *, ' '
PRINT *, ' '
PRINT *, ' ENTER DESIRED SPHER. HARMONIC LIMITS INT/EXT'
PRINT *, ' NIMN,NIMX,NEMN,NEMX'
READ *, NIMN,NIMX,NEMN,NEMX
C
C
PRINT *, ' '
PRINT *, ' '
PRINT *, ' ENTER DESIRED CHART DATE (E.G. 2001.72534122D0)'
READ *, TIME
TIM=TIME
C
C
C READ MODEL EPOCH AND COEFFICIENTS
C
C
READ(7,10) EPOCH,MODEL
10 FORMAT(1X,F14.9,A60)
PRINT *, ' '
PRINT *, ' '
PRINT 10, EPOCH,MODEL
PRINT *, ' '
C
PRINT *, ' '
EPOC=DBLE(EPOCH)
PRINT *, EPOC,MODEL
PAUSE
C
20 CONTINUE
READ(7,30) N,M,GNM,HNM,GTNM,HTNM
30 FORMAT(2(1X,I3),4(2X,F14.7))
C
IF (N .EQ. 999) GO TO 40
IF ( N .GT. 100 .OR. M .GT. 100) THEN
PRINT *, ' ERROR: N OR M > 100'
STOP 1
ENDIF
C
NIMAX=N
CI0(N,M)=DBLE(GNM)
CIT0(N,M)=DBLE(GTNM)
C
IF (M .NE. 0) THEN
CI0(M-1,N)=DBLE(HNM)
CIT0(M-1,N)=DBLE(HTNM)
ENDIF
C
GO TO 20
C
C
40 CONTINUE
C
READ(7,30) N,M,QNM,SNM,QTNM,STNM
C
191
IF (N .EQ. 999) GO TO 50
C
IF (N .GT. 10 .OR. M .GT. 10) THEN
PRINT *, ' ERROR: N OR M . 10'
STOP 2
ENDIF
C
NEMAX=N
CE0(N,M)=DBLE(QNM)
CET0(N,M)=DBLE(QTNM)
C
IF (M .NE. 0) THEN
CE0(M-1,N)=DBLE(SNM)
CET0(M-1,N)=DBLE(STNM)
ENDIF
C
GO TO 40
C
C
50 CONTINUE
REWIND 7
CLOSE(7)
PRINT *, ' '
PRINT *, ' '
PRINT *, ' ALL COEFICIENTS HAVE BEEN READ'
PRINT *, ' '
PAUSE
C
C
PRINT *, ' '
PRINT *, ' '
C
PRINT *, ' ISLCT = ',ISLCT
PRINT *, ' JSLCT = ',JSLCT
PRINT *, ' ALT = ',ALT,' KM'
PRINT *, ' NIMN = ',NIMN
PRINT *, ' NIMX = ',NIMX
PRINT *, ' NEMN = ',NEMN
PRINT *, ' NEMX = ',NEMX
PRINT *, ' NIMAX = ',NIMAX
PRINT *, ' NEMAX = ',NEMAX
PRINT *, ' '
PRINT *, ' '
PRINT *, ' RLATMN = ',RLATMN,' DEG.'
PRINT *, ' RLATMX = ',RLATMX,' DEG.'
PRINT *, ' RLONMN = ',RLONMN,' DEG.'
PRINT *, ' RLONMX = ',RLONMX,' DEG.'
PRINT *, ' '
PRINT *, ' '
PRINT *, ' RLTGRD = ',RLTGRD,' DEG.'
PRINT *, ' RLNGRD = ',RLNGRD,' DEG.'
PRINT *, ' '
PRINT *, ' '
C
C
PRINT *, ' '
PRINT *, ' '
PAUSE
C
C
C INITIALIZE MAGPOT SUBROUTINE
C
C
C
CALL MAGPOT(NIMAX,NEMAX,ISLCT)
C
C
C
C
PRINT *, ' '
PRINT *, ' '
PRINT *, ' MAGPOT INITIALIZED'
PRINT *, ' '
PRINT *, ' '
C
C
PAUSE
C
C
C***********************************************************************
C
C
C THE AREA OF THE WORLD SELECTED FOR GEOMAGNETIC INVERSION IS
C DETERMINED BY THE CHART LIMITS (LATMIN,LATMAX,LONMIN,LONMAX).
C
C THESE LIMITS COULD COVER THE ENTIRE EARTH OR ANY PORTION OF IT.
C
C
C THESE PARAMETERS AND THEIR CORRESPONDING GRID INTERVALS
C (LATGRD,LONGRD) ARE GIVEN IN INTEGER DEGREES: => 1
C
C
C***********************************************************************
C
C
C
DO 180 LAT=LATMIN,LATMAX,LATGRD
DO 179 JLAT=0,3
IF (LAT .EQ. 90 .AND. JLAT .GE. 1) GO TO 179
RLAT=DFLOAT(LAT)+0.25D0*DFLOAT(JLAT)
IF (IPRINT .NE. 1) PRINT *, ' LAT = ',SNGL(RLAT)
C
C
C AVOID DIVISION BY ZERO PROBLEMS AT THE POLES
C
C
IF (LAT .EQ. 90 .AND. JLAT .EQ. 0) RLAT=89.999998D0
IF (LAT .EQ. -90 .AND. JLAT .EQ. 0) RLAT=-89.999998D0
C
C
POSIN(2)=RLAT
DO 170 LON=LONMIN,LONMAX,LONGRD
DO 169 JLON=0,3
IF (LON .EQ. 180 .AND. JLON .GE. 1) GO TO 169
RLON=DFLOAT(LON)+0.25D0*DFLOAT(JLON)
POSIN(3)=RLON
IF (IPRINT .EQ. 1) THEN
PRINT *, ' LAT = ',SNGL(RLAT),' LON = ',SNGL(RLON)
ENDIF
C
C
C COMPUTE OBSERVED MAGNETIC FIELD AND MAGNETIC GRADIENT FIELD
C AT ONE ALT/LAT/LON OBSERVATION POINT.
C
C
CALL MAGPT1(POSIN,POSOUT,TIM,VI,BI,BFI,GI,GFI,VE,BE,BFE,GE,GFE)
C
C
C
IF (LAT .EQ. 90 .OR. LAT .EQ. -90) GI(1,2)=0.000D0
IF (LAT .EQ. 90 .OR. LAT .EQ. -90) GI(2,1)=0.000D0
C
C
C COMPUTE THE PRISM SOURCE DEPTH TO NEAREST 0.1 KILOMETERS
C BY MINIMIZING THE CHI-SQUARE FUNCTION. THE COORDINATE
C SYSTEM ORIGIN IS LOCATED AT THE EARTH'S SURFACE. THUS,
C ZB IS POSITIVE, WHILE Z0 IS NEGATIVE.
C
C
CHIMIN=1.000D10
DO 160 I=MINDEP,MAXDEP,ISTEP
ZB=DFLOAT(I)/10.000D0
C
C
C COMPUTE LAMDA MATRIX: ITS GRADIENT AND ITS INVERSE
C
C
CALL LAMDLA(X0,Y0,Z0,LAMDA,DLAMDA)
C
C
DO 70 J=1,3
DO 60 K=1,3
LAMDAI(J,K)=LAMDA(J,K)
60 CONTINUE
70 CONTINUE
C
C
CALL MATINV(LAMDAI,NORDER,DET)
C
C
C
DO 110 MU=1,3
DO 100 NU=1,3
GIP(MU,NU)=0.000D0
DO 90 KP=1,3
DO 80 LM=1,3
GIP(MU,NU)=GIP(MU,NU)+DLAMDA(MU,KP,NU)*LAMDAI(KP,LM)*BI(LM)
80 CONTINUE
90 CONTINUE
100 CONTINUE
110 CONTINUE
C
C
IF (IPRINT .EQ. 1) THEN
PRINT *, ' '
PRINT *, SNGL(ZB),SNGL(LAMX),SNGL(LAMY),SNGL(LAMZ)
PRINT *, ' '
PRINT *, SNGL(BI(1)),SNGL(BI(2)),SNGL(BI(3))
PRINT *, ' '
PRINT *, SNGL(GI(1,1)),SNGL(GI(1,2)),SNGL(GI(1,3))
PRINT *, SNGL(GI(2,1)),SNGL(GI(2,2)),SNGL(GI(2,3))
PRINT *, SNGL(GI(3,1)),SNGL(GI(3,2)),SNGL(GI(3,3))
PRINT *, ' '
PRINT *, SNGL(GIP(1,1)),SNGL(GIP(1,2)),SNGL(GIP(1,3))
PRINT *, SNGL(GIP(2,1)),SNGL(GIP(2,2)),SNGL(GIP(2,3))
PRINT *, SNGL(GIP(3,1)),SNGL(GIP(3,2)),SNGL(GIP(3,3))
PRINT *, ' '
ENDIF
C
CHISQ=0.000D0
197
DO 130 MU=1,3
DO 120 NU=1,3
CHISQ=CHISQ+(GI(MU,NU)-GIP(MU,NU))**2
120 CONTINUE
130 CONTINUE
C
IF (IPRINT .EQ. 1) THEN
PRINT *, CHISQ
PRINT *, ' '
PAUSE
ENDIF
C
IF (CHISQ .GT. CHIMIN) GO TO 165
C
IF (CHISQ .LT. CHIMIN) THEN
CHILST=CHIMIN
CHIMIN=CHISQ
SRCDEP=ZB
DELCHI=CHILST-CHIMIN
DO 150 K=1,3
MAG(K)=0.000D0
DO 140 L=1,3
MAG(K)=MAG(K)+LAMDAI(K,L)*BI(L)
140 CONTINUE
150 CONTINUE
MX=CNTAM*MAG(1)
MY=CNTAM*MAG(2)
MZ=CNTAM*MAG(3)
MT=DSQRT(MX**2+MY**2+MZ**2)
IF (IPRINT .EQ. 1) THEN
PRINT *, ' CHISQ = ',CHISQ
C
PRINT *, ' DEPTH = ',SRCDEP,' KM'
PRINT *, ' MAGX = ',MX,' A/M'
PRINT *, ' MAGY = ',MY,' A/M'
PRINT *, ' MAGZ = ',MZ,' A/M'
PRINT *, ' MAGT = ',MT,' A/M'
PRINT *, ' '
ENDIF
ENDIF
C
160 CONTINUE
165 CONTINUE
C
SDEP(4*LAT+JLAT,4*LON+JLON)=SRCDEP
MAGX(4*LAT+JLAT,4*LON+JLON)=MX
MAGY(4*LAT+JLAT,4*LON+JLON)=MY
MAGZ(4*LAT+JLAT,4*LON+JLON)=MZ
MAGT(4*LAT+JLAT,4*LON+JLON)=MT
C
WRITE(11,166) SNGL(RLAT),SNGL(RLON),SNGL(SRCDEP),
* SNGL(MX),SNGL(MY),SNGL(MZ),SNGL(MT)
166 FORMAT(3(2X,F12.4),4(2X,F16.11))
C
169 CONTINUE
170 CONTINUE
C
179 CONTINUE
180 CONTINUE
C
C
ENDFILE 11
REWIND 11
199
CLOSE(11)
C
C
PRINT *, ' END OF JOB'
PAUSE
STOP
END
C
C***********************************************************************
C
C
SUBROUTINE MAGPOT(NIMAX,NEMAX,ISLCT)
C
C
C***********************************************************************
C
C
C PROGRAMMED BY: JOHN M. QUINN 12/10/03
C SOLAR-TERRESTRIAL PHYSICS INSTITUTE
C 2732 SOUTH BRAUN WAY
C LAKEWOOD, CO 80228
C USA
C
C PHONE: (303) 989-6201
C EMAIL: p-j_quinn@comcast.net
C
C
C***********************************************************************
C
C
C PURPOSE:
C
C THIS ROUTINE COMPUTES THE INTERNAL AND EXTERNAL PARTS OF THE
C EARTH'S MAGNETIC POTENTIAL, MAGNETIC FIELD VECOTR, AND
C MAGNETIC GRADIENT TENSOR AT ANY ALTITUDE (DEPTH), AT ANY
C LOCATION (LATITUDE/LONGITUDE), AND AT ANY TIME USING A
C SPHERICAL HARMONIC MAGNETIC FIELD MODEL DETERMINED FROM
C LAND, SEA, AIR, AND SPACE-BORNE MAGNETIC FIELD SENSOR DATA.
C
C
C
C NOTE:
C
C THE USER SPECIFIES THE INPUT AND OUTPUT COORDINATE SYSTEMS
C ACCORDING TO THE VALUES OF THE PARAMETER ISLCT:
C
C IF ISLCT=1 THEN THE INPUT POSITIONS ARE REFERENCED
C TO GEODETIC COORDS. (ALT,GLAT,GLON), WHILE THE
C OUTPUT POSITIONS AND MAGNETIC FIELD PARAMETERS ARE
C REFERENCED TO GEODETIC COORDINATES (ALT,GLAT,GLON)
C
C IF ISLCT=2 THEN THE INPUT POSITIONS ARE REFERECED TO GEODETIC
C COORDS (GALT,GLAT,GLON) WHILE THE OUTPUT POSITIONS AND
C MAGNETIC PARMATERS ARE REFERENCED TO SPHERICAL COORDS.
C (R,THETA,PHI).
C
C IF ISLCT=3 THEN THE INPUT POSITIONS ARE REFERENCED TO SPHERICAL
C COORDS. (R,THETA,PHI), WHILE THE OUTPUT POSITIONS
C AND MAGNETIC FIELD PARAMETERS ARE REFERENCED TO
C SPHERICAL COORDS.
C
C
C***********************************************************************
C
C
C REFERENCES:
C
C LANGEL, R. A., THE MAIN FIELD, IN GEOMAGNETISM, CHAPTER 4,
C ED. BY: J. A. JACOBS, ACADEMIC PRESS, INC., ORLANDO,
C FLORIDA (1987)
C
C
C DEFENSE MAPPING AGENCY (DMA), DEPARTMENT OF DEFENSE WORLD GEODETIC
C SYSTEM 1984 - ITS DEFINITION AND RELATIONSHIPS WITH
C LOCAL COORDINATE SYSTEMS, DMA TECHNICAL REPORT #8350.2,
C WASHINGTON, D.C. (1987)
C
C NOTE: DMA CHANGED ITS NAME A NUMBER OF YEARS AGO TO NIMA
C (NATIONAL IMAGERY AND MAPPING AGENCY) AND MORE RECENTLY HAS
C CHANGED ITS NAME AGAIN TO GIA (GEOSPATIAL INFORMATION AGENCY)
C
C
C***********************************************************************
C
C
C INPUT PARAMETERS:
C
C ISLCT - COORDINATE SYSTEM INPUT/OUTPUT SELECTION MODE
C 1 = GEODETIC IN/GEODETIC OUT
C 2 = GEODETIC IN/SPHERICAL OUT
C 3 = SPHERICAL IN/SPHERICAL OUT
C NIMAX - INPUT MODEL MAXIMUM INTERNAL SPHERICAL HARMONIC DEG.
C NEMAX - INPUT MODEL MAXIMUM EXTERNAL SPHERICAL HARMONIC DEG.
C NIMN - MINIMUM DESIRED INTERNAL SPHERICAL HARMONIC DEG.
C NIMX - MAXIMUM DESIRED INTERNAL SPHERICAL HARMONIC DEG.
C NEMN - MINIMUM DESIRED EXTERNAL SPHERICAL HARMONIC DEG.
C NEMX - MAXIMUM DESIRED EXTERNAL SPHERICAL HARMONIC DEG.
C MODEL - NAME OF MAGNETIC FIELD MODEL (HOLLERETH)
C EPOC - REFERENCE YEAR OF MODEL
C POSIN - POSITION VECTOR (INPUT)
C IF ISLCT= 1 OR 2, THEN IF ISLCT= 3, THEN
C POS(1)=ALT (KM) POS(1)=R (KM)
C POS(2)=GLAT (DEG.) POS(2)=THETA (DEG.)
C
C POS(3)=GLON (DEG.) POS(3)=GLON (DEG.)
C TIME - TIME (YEARS)
C
C
C NOTE: GEODETIC OUTPUT IS REFERENCED TO THE WGS84 ELLIPSOID
C
C
C***********************************************************************
C
C
C NOTE:
C
C IF INPUT OR OUTPUT IS REFERENCED TO SPHERICAL COORDS.:
C
C R = RADIALLY OUTWARD
C THETA = SOUTH
C PHI = EAST
C
C
C IF INPUT OR OUTPUT IS REFERENCED TO GEODETIC COORDS.:
C
C Z = VERTICALLY DOWN
C X = NORTH
C Y = EAST
C
C
C***********************************************************************
C
C
C OUTPUT PARAMETERS:
C
C
C VI - MAGNETIC POTENTIAL INTERNAL PART (nT-Km)
C BI - MAGNETIC FIELD VECTOR INTERNAL PART (nT)
C BFI - TOTAL MAGNETIC INTENSITY INTERNAL PART (nT)
C GI - MAGNETIC GRADIENT TENSOR INTERNAL PART (nT/Km)
C GFI - MAGNETIC GRADIENT INTENSITY INTERNAL PART (nT/Km)
C VE - MAGNETIC POTENTIAL EXTERNAL PART (nT-Km)
C BE - MAGNETIC FIELD VECTOR EXTERNAL PART (nT)
C BFE - TOTAL MAGNETIC INTENSITY EXTERNAL PART (nT)
C GE - MAGNETIC GRADIENT TENSOR EXTERNAL PART (nT/Km)
C GFE - MAGNETIC GRADIENT INTENSITY EXTERNAL PART (nT/Km)
C POSOUT - POSITION VECTOR
C IF ISLCT = 1 OR 2 THEN IF ISLCT = 3 THEN
C POS(1)=ALT (KM) POS(1)=R (KM)
C POS(2)=GLAT (DEG.) POS(2)=THETA (DEG.)
C POS(3)=GLON (DEG.) POS(3)=PHI (DEG.)
C
C
C NOTE: GEODETIC PARAMETERS ARE REFERENCED TO THE WGS-84 ELLIPSOID
C
C
C***********************************************************************
C
C
C OTHER PARAMETERS:
C
C A - WGS84 ELLIPSOID SEMI-MAJOR AXIS LENGTH (KM)
C B - WGS84 ELLIPSOID SEMI-MINOR AXIS LENGTH (KM)
C F - WGS84 ELLIPSOID FLATTENING (A-B)/A
C FAC(N) - N FACTORIAL
C SNORM - LEGENDRE FUNCTION SCHMIDT QUASI-NORMALIZATION
C R - SPHERICAL COORD.GEOCENTRIC RADIAL RADIAL DISTANCE (KM)
C
C THETA - SPHERICAL COORD. CO-LATITUDE (RADIANS)
C PHI - SPHERICAL COORD. LONGITUDE (RADIANS)
C ST - SIN(THETA)
C CT - COSIN(THETA)
C CP(M) - COSIN(M*PHI)
C SP(M) - SIN(M*PHI)
C BETA - GEODETIC LATITUDE (RADIANS)
C SB - SIN(BETA)
C CB - COS(BETA)
C ALPHA - SPHERICAL TO GEODETIC COORD. ROTATION ANGLE (RADIANS)
C SA - SIN(ALPHA)
C CA - COS(ALPHA)
C P(N,M) - ASSOCIATED LEGENDRE POLYNOMIAL (DEG. N, ORD. M)
C DP(N,M)- 1'ST DERIVATIVE OF P(N,M) WRT. THETA
C D2P(N,M)- 2'ND DERIVATIVE OF P(N,M) WRT. THETA
C CI0(N,M)- INPUT MODEL INTERNAL MAIN FIELD GAUSS COEFFS. (nT)
C CE0(N,M)- INPUT MODEL EXTERNAL GAUSS COEFFS. (nT)
C CIT0(N,M)- INPUT MODEL INTERNAL SECULAR VARIATION COEFS. (nT/Yr)
C CET0(N,M)- INPUT MODEL EXTERNAL SECULAR VARIATION COEFS. (nT/Yr)
C CI(N,M)- TIME ADJUSTED INTERNAL GAUSS COEFS. UNNORMALIZED (nT)
C CE(N,M)- TIME ADJUSTED EXTERNAL GAUSS COEFS. UNNORMALIZED (nT)
C OTIM - TIME (YRS) ON PREVIOUS CALL TO SUBROUTINE
C OALT - ALTITUDE (RADIUS) (KM) ON PREVIOUS CALL TO SUBROUTINE
C OLAT - LATITIUDE (CO-LATITUDE) (DEG.) ON PREVIOUS CALL
C OLON - LONGITUDE (DEG.) ON PREVIOUS CALL TO SUBROUTINE
C
C
C***********************************************************************
C
C
C PROCEDURE:
C
C
C THIS ROUTINE ACCEPTS GEODETIC OR SPHERICAL POSITION AND TIME
C INFORMATION, CONVERTS THESE INFORMATION FROM GEODETIC TO
C SPHERICAL CORRDINATES IF NECESSARY, PERFORMS ALL MAGNETIC
C POTENTIAL, VECTOR FIELD, AND GRADIENT FIELD COMPUTATIONS IN
C SPHERICAL COORDS., AND FINALLY, IF NECESSARY, ROTATES THE
C COMPUTED VECTOR AND GRADIENT FIELD VALUES BACK INTO GEODETIC
C COORDINATES, SUCH THAT THE OUTPUT POSITIONS AND MAGNETIC FIELD
C PARAMETERS ARE REFERENCED TO EITHER SPHERICAL OR GEODETIC
C COORDINTES, ACCORDING TO THE USERS SELECTION.
C
C ALL OR PART OF THE MAGNETIC FIELD WILL BE COMPUTED DEPENDING
C ON THE USER INPUT VALUES FOR NIMN, NIMX, NEMN, AND NEMX. THIS
C PERMITS SEPARATE ANALYSES OF SHORT AND LONG WAVELENGTH
C CONTRIBUTIONS TO THE MAGNETIC FIELD.
C
C
C***********************************************************************
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
DOUBLE PRECISION K(0:100,0:100)
C
CHARACTER*60 MODEL
C
DIMENSION POSIN(3),POSOUT(3)
DIMENSION SP(0:100),CP(0:100)
DIMENSION BI(3),GI(3,3),BE(3),GE(3,3)
DIMENSION CI0(0:100,0:100),CIT0(0:100,0:100)
DIMENSION CE0(0:100,0:100),CET0(0:100,0:100)
207
DIMENSION CI(0:100,0:100),CE(0:100,0:100)
DIMENSION SNORM(0:100,0:100)
DIMENSION P(0:100,0:100),DP(0:100,0:100),D2P(0:100,0:100)
C
COMMON /LIMITS/ NIMN,NIMX,NEMN,NEMX
COMMON /COEFFS/ CI0,CIT0,CE0,CET0
COMMON /MDLEPC/ MODEL,EPOC
C
C
C***********************************************************************
C
C
C INITIALIZATION MODULE
C
C
C***********************************************************************
C
C
NMAX=NIMX
LSLCT=ISLCT
C
IF (NEMX .GT. NIMX) NMAX=NEMX
C
IF (NEMX .GT. NEMAX .OR. NIMX .GT. NIMAX) THEN
PRINT *, ' REQUESTED DEGREE EXCEEDS THAT OF INPUT MODEL'
STOP
ENDIF
C
NMAX2=2*NMAX
PI=3.141592654D0
DTR=PI/180.000000000D0
C
RE=6371.200000000D0
SP(0)=0.000000000D0
CP(0)=1.000000000D0
P(0,0)=1.000000000D0
DP(0,0)=0.000000000D0
D2P(0,0)=0.000000000D0
OTIM=-9999.0D0
OLAT=-9999.0D0
OLON=-9999.0D0
OALT=-9999.0D0
C
C INITIALIZE ELLIPSOID CONSTANTS
C
A=6378.13700D0
F=1.000000000D0/298.257223563D0
B=A*(1.0000000D0-F)
A2=A*A
B2=B*B
H2=A2-B2
A4=A2*A2
B4=B2*B2
H4=A4-B4
EPS=0.000000001D0
C
C INITIALIZE ASSOCIATED LEGENDRE FUNCTION RECURSION COEFS.
C
C FIRST COMPUTE: (N-M)!/(N+M)!
C
SNORM(0,0)=1.000000000D0
DO 20 N=1,NMAX
SNORM(N,0)=SNORM(N-1,0)*DFLOAT(2*N-1)/DFLOAT(N)
C
J=2
DO 10 M=0,N
K(N,M)=DFLOAT((N-1)**2-M**2)/DFLOAT((2*N-1)*(2*N-3))
IF (M .GT. 0) THEN
FLNMJ=DFLOAT((N-M+1)*J)/DFLOAT(N+M)
SNORM(N,M)=SNORM(N,M-1)*DSQRT(FLNMJ)
J=1
CI0(M-1,N)=SNORM(N,M)*CI0(M-1,N)
CIT0(M-1,N)=SNORM(N,M)*CIT0(M-1,N)
CE0(M-1,N)=SNORM(N,M)*CE0(M-1,N)
CET0(M-1,N)=SNORM(N,M)*CET0(M-1,N)
ENDIF
CI0(N,M)=SNORM(N,M)*CI0(N,M)
CIT0(N,M)=SNORM(N,M)*CIT0(N,M)
CE(N,M)=SNORM(N,M)*CE0(N,M)
CET0(N,M)=SNORM(N,M)*CET0(N,M)
C
C
10 CONTINUE
20 CONTINUE
K(1,1)=0.000000000D0
C
C
RETURN
C
C
C***********************************************************************
C
C
C COMPUTATION MODULE
C
C
C
C***********************************************************************
C
C
ENTRY MAGPT1(POSIN,POSOUT,TIM,VI,BI,BFI,GI,GFI,VE,BE,BFE,GE,GFE)
C
C
C***********************************************************************
C
C
DT=TIM-EPOC
C
C
C CONVERT FROM GEODETIC TO SPHERICAL COORDINATES
C
C
IF (LSLCT .EQ. 3) THEN
C
IF (POSIN(1) .NE. OALT) R=POSIN(1)
C
IF (POSIN(2) .NE. OLAT)THEN
THETA=POSIN(2)*DTR
CT=DCOS(THETA)
ST=DSIN(THETA)
ENDIF
C
IF (POSIN(3) .NE. OLON) THEN
GLON=POSIN(3)
PHI=GLON*DTR
SP(1)=DSIN(PHI)
CP(1)=DCOS(PHI)
211
ENDIF
C
ELSE
C
IF (POSIN(1) .NE. OALT) GALT=POSIN(1)
C
IF (POSIN(3) .NE. OLON) THEN
GLON=POSIN(3)
PHI=GLON*DTR
SP(1)=DSIN(PHI)
CP(1)=DCOS(PHI)
ENDIF
C
IF (POSIN(2) .NE. OLAT) THEN
GLAT=POSIN(2)
CB=DCOS(GLAT*DTR)
SB=DSIN(GLAT*DTR)
CB2=CB*CB
SB2=SB*SB
E1=A2-H2*SB2
E2=DSQRT(E1)
E3=A4-H4*SB2
E4=DSQRT(E3)
E5=(A2+GALT*E1)/(B2+GALT*E1)
E6=E5*E5
E7=DSQRT(A2*CB2+B2*SB2)
CT=SB/DSQRT(E6*CB2+SB2)
ST=DSQRT(1.000000000D0-CT**2)
R=DSQRT(GALT**2+2.00000000D0*GALT*E2+E3/E1)
THETA=ACOS(CT)/DTR
CA=(GALT+E7)/R
C
SA=(H2*SB*CB)/(R*E7)
ENDIF
ENDIF
C
IF (DABS(ST) .LT. 0.0000000001D0) ST=0.0000000001D0
C
C
C COMPUTE ASSOCIATED LEGENDRE FUNCTIONS & 1'ST/2'ND DERIVATIVES
C AND TIME ADJUST THE MODEL COEFFICIENTS
C
C
DO 40 N=1,NMAX
C
IF (POSIN(3) .NE. OLON) THEN
IF (N .GT. 1) SP(N)=SP(1)*CP(N-1)+SP(N-1)*CP(1)
IF (N .GT. 1) CP(N)=CP(1)*CP(N-1)-SP(1)*SP(N-1)
ENDIF
C
DO 30 M=0,N
C
IF (TIM .NE. OTIM) THEN
CI(N,M)=CI0(N,M)+DT*CIT0(N,M)
CE(N,M)=CE0(N,M)+DT*CET0(N,M)
IF (M .NE. 0) THEN
CI(M-1,N)=CI0(M-1,N)+DT*CIT0(M-1,N)
CE(M-1,N)=CE0(M-1,N)+DT*CET0(M-1,N)
ENDIF
ENDIF
C
C
IF (POSIN(2) .NE. OLAT) THEN
C
C
IF (N .EQ. M .AND. N .GT. 0) THEN
P(N,M)=ST*P(N-1,N-1)
DP(N,M)=CT*P(N-1,N-1)+ST*DP(N-1,N-1)
D2P(N,M)=-ST*P(N-1,N-1)+2.000000000D0*CT*DP(N-1,N-1)+
* ST*D2P(N-1,N-1)
ENDIF
C
IF (N .EQ. 1 .AND. M .EQ. 0) THEN
P(N,M)=CT*P(N-1,M)
DP(N,M)=-ST*P(N-1,M)+CT*DP(N-1,M)
D2P(N,M)=-CT*P(N-1,M)-2.000000000D0*ST*DP(N-1,M)+CT*D2P(N-1,M)
ENDIF
C
IF (N .GT. 1 .AND. N .NE. M) THEN
C
IF (M .GT. N-2) THEN
P(N-2,M)=0.000000000D0
DP(N-2,M)=0.000000000D0
D2P(N-2,M)=0.000000000D0
ENDIF
C
P(N,M)=CT*P(N-1,M)-K(N,M)*P(N-2,M)
DP(N,M)=-ST*P(N-1,M)+CT*DP(N-1,M)-K(N,M)*DP(N-2,M)
D2P(N,M)=-CT*P(N-1,M)-2.000000000D0*ST*DP(N-1,M)+
* CT*D2P(N-1,M)-K(N,M)*D2P(N-2,M)
ENDIF
ENDIF
C
C
30 CONTINUE
C
40 CONTINUE
C
C
C COMPUTE MAGNETIC POTENTIALS, VECTORS, AND GRADIENTS
C FROM SPHERICAL HARMONIC EXPANSIONS
C
C
VI=0.0D0
BRI=0.0D0
BTI=0.0D0
BPI=0.0D0
GIRR=0.0D0
GIRT=0.0D0
GIRP=0.0D0
GITR=0.0D0
GITT=0.0D0
GITP=0.0D0
GIPR=0.0D0
GIPT=0.0D0
GIPP=0.0D0
C
VE=0.0D0
BRE=0.0D0
BTE=0.0D0
BPE=0.0D0
GERR=0.0D0
GERT=0.0D0
GERP=0.0D0
GETR=0.0D0
GETT=0.0D0
GETP=0.0D0
C
GEPR=0.0D0
GEPT=0.0D0
GEPP=0.0D0
C
C
AR=RE/R
RA=R/RE
C
AR1=AR
AR2=AR*AR
AR3=AR2*AR
C
RA0=1.000000000D0
RA1=RA
RA2=RA*RA
C
C
C print *, r,re,ar,ra
C print *, ar1,ar2,ar3
C print *, ra0,ra1,ra2
C pause
C
C
DO 80 N=1,NMAX
FNM1=DFLOAT(N-1)
FN0=DFLOAT(N)
FN1=DFLOAT(N+1)
FN2=DFLOAT(N+2)
DO 70 M=0,N
FM0=DFLOAT(M)
FM2=FM0*FM0
C
AR1=AR1*AR
AR2=AR2*AR
AR3=AR3*AR
RA0=RA0*RA
RA1=RA1*RA
RA2=RA2*RA
IF (M .EQ. 0) THEN
UI1=CI(N,M)*CP(M)
UI2=CI(N,M)*SP(M)
UE1=CE(N,M)*CP(M)
UE2=CE(N,M)*SP(M)
ELSE
UI1=CI(N,M)*CP(M)+CI(M-1,N)*SP(M)
UI2=CI(N,M)*SP(M)-CI(M-1,N)*CP(M)
UE1=CE(N,M)*CP(M)+CE(M-1,N)*SP(M)
UE2=CE(N,M)*SP(M)-CE(M-1,N)*CP(M)
ENDIF
C
C
C print *, n,m
C print *, n,m,ui1,ui2
C print *, ci(n,m),ci(m-1,n)
C print *, sp(m),cp(m)
C pause
C
C
IF (N .LT. NIMN .OR. N .GT. NIMX) GO TO 50
VI=VI+AR1*UI1*P(N,M)
BRI=BRI+FN1*AR2*UI1*P(N,M)
BTI=BTI-AR2*UI1*DP(N,M)
BPI=BPI+FM0*AR2*UI2*P(N,M)/ST
C
GIRR=GIRR-FN1*FN2*AR3*UI1*P(N,M)
GIRT=GIRT+FN2*AR3*UI1*DP(N,M)
GIRP=GIRP-FN2*FM0*AR3*UI2*P(N,M)/ST
GITR=GITR+FN2*AR3*UI1*DP(N,M)
GITT=GITT+(FN1*P(N,M)-D2P(N,M))*AR3*UI1
GITP=GITP+FM0*(DP(N,M)-DCOTAN(THETA)*P(N,M))*AR3*UI2/ST
GIPR=GIPR-FN2*FM0*AR3*UI2*P(N,M)/ST
GIPT=GIPT+FM0*(DP(N,M)-DCOTAN(THETA)*P(N,M))*AR3*UI2/ST
GIPP=GIPP+(FM2*P(N,M)/ST+FN1*ST*P(N,M)-CT*DP(N,M))*AR3*UI1/ST
C
C
C print *, n,m
C print *, bri,bti,bpi
C print *, girr,girt,girp
C print *, gitr,gitt,gitp
C print *, gipr,gipt,gipp
C pause
C
C
50 CONTINUE
IF (N .LT. NEMN .OR. N .GT. NEMX) GO TO 60
VE=VE+RA0*UE1*P(N,M)
BRE=BRE-DFLOAT(N)*RA1*UE1*P(N,M)
BTE=BTE-RA1*UE1*DP(N,M)
BPE=BPE+DFLOAT(M)*RA1*UE2*P(N,M)/ST
GERR=GERR-FN0*FNM1*RA2*UE1*P(N,M)
GERT=GERT-FNM1*RA2*UE1*DP(N,M)
GERP=GERP+FNM1*FM0*RA2*UE2*P(N,M)/ST
GETR=GETR-FNM1*RA2*UE1*DP(N,M)
GETT=GETT-(D2P(N,M)+FN0*P(N,M))*RA2*UE1
GETP=GETP+FM0*(DP(N,M)-DCOTAN(THETA)*P(N,M))*RA2*UE2/ST
C
GEPR=GEPR+FM0*FNM1*RA2*UE2*P(N,M)/ST
GEPT=GEPT+FM0*(DP(N,M)-DCOTAN(THETA)*P(N,M))*RA2*UE2/ST
GEPP=GEPP+(FM2*P(N,M)/ST-FN0*ST*P(N,M)-CT*DP(N,M))*RA2*UE1/ST
C
C
C print *, n,m
C print *, bre,bte,bpe
C print *, gerr,gert,gerp
C print *, getr,gett,getp
C print *, gepr,gept,gepp
C pause
C
60 CONTINUE
C
C
70 CONTINUE
80 CONTINUE
C
C
VI=RE*VI
VE=RE*VE
C
C
C print *, vi,ve
C pause
C
C
GIRR=GIRR/RE
GIRT=GIRT/RE
GIRP=GIRP/RE
C
GITR=GITR/RE
GITT=GITT/RE
GITP=GITP/RE
GIPR=GIPR/RE
GIPT=GIPT/RE
GIPP=GIPP/RE
C
C
C print *, bri,bti,bpi
C print *, girr,girt,girp
C print *, gitr,gitt,gitp
C print *, gipr,gipt,gipp
C pause
C
C
GERR=GERR/RE
GERT=GERT/RE
GERP=GERP/RE
GETR=GETR/RE
GETT=GETT/RE
GETP=GETP/RE
GEPR=GEPR/RE
GEPT=GEPT/RE
GEPP=GEPP/RE
C
C
C print *, gerr,gert,gerp
C print *, getr,gett,getp
C print *, gepr,gept,gepp
C pause
C
C
C
BFI=DSQRT(BRI**2+BTI**2+BPI**2)
GFI=DSQRT(GIRR**2+GIRT**2+GIRP**2+GITR**2+GITT**2+GITP**2+
* GIPR**2+GIPT**2+GIPP**2)
C
BFE=DSQRT(BRE**2+BTE**2+BPE**2)
GFE=DSQRT(GERR**2+GERT**2+GERP**2+GETR**2+GETT**2+GETP**2+
* GEPR**2+GEPT**2+GEPP**2)
C
C
C
C
IF (LSLCT .EQ. 2 .OR. LSLCT .EQ. 3) THEN
C
C
C FILL OUTPUT ARRAYS
C
C
POSOUT(1)=R
POSOUT(2)=THETA/DTR
POSOUT(3)=PHI/DTR
C
BI(1)=BRI
BI(2)=BTI
BI(3)=BPI
C
GI(1,1)=GIRR
GI(1,2)=GIRT
GI(1,3)=GIRP
GI(2,1)=GITR
GI(2,2)=GITT
C
GI(2,3)=GITP
GI(3,1)=GIPR
GI(3,2)=GIPT
GI(3,3)=GIPP CC
BE(1)=BRE
BE(2)=BTE
BE(3)=BPE C
GE(1,1)=GERR
GE(1,2)=GERT
GE(1,3)=GERP
GE(2,1)=GETR
GE(2,2)=GETT
GE(2,3)=GETP
GE(3,1)=GEPR
GE(3,2)=GEPT
GE(3,3)=GEPP CC
OTIM=TIM
OALT=POSIN(1)
OLAT=POSIN(2)
OLON=POSIN(3) CC
RETURN
ENDIF CC
C
C***********************************************************************
C
C
C ROTATE MAGNETIC VECTORS AND GRADIENTS TO GEODETIC COORDINATES
C
C Bnew=(ROT)Bold
C
C Gnew=(ROT)Gold(ROTI)
C
C WHERE ROT IS THE ROTATION MATRIX AND ROTI IS ITS INVERSE
C AND WHERE B AND G ARE THE NEW AND OLD MAGNETIC FIELD VECTORS
C AND GRADIENTS
C
C
C FIRST ROTATE FROM (R,THETA,PHI) TO (-Z,-X,Y) COORDINATES,
C
C THEN INVERT AND REORDER AXES FROM (-Z,-X,Y) TO (X,Y,Z) COORDS.
C
C THUS:
C
C BX=-BTnew
C BY=+BPnew
C BZ=-BRnew
C
C GXX=+GTTnew
C GXY=-GTPnew
C GXZ=+GTRnew
C GYX=-GPTnew
C GYY=+GPPnew
C GYZ=-GPRnew
C GZX=+GRTnew
C
C GZY=-GRPnew
C GZZ=+GRRnew
C
C
C WHERE:
C
C
C | |
C | COS(ALPHA) -SIN(ALPHA) 0 |
C | |
C ROT = | SIN(ALPHA) COS(ALPHA) 0 |
C | |
C | 0 0 1 |
C | |
C
C
C | |
C | COS(ALPHA) SIN(ALPHA) 0 |
C | |
C ROTI = |-SIN(ALPHA) COS(ALPHA) 0 |
C | |
C | 0 0 1 |
C | |
C
C
C
C
C
C
C***********************************************************************
C
C
C
IF (LSLCT .EQ. 1) THEN
C
C
BXI=-(SA*BRI+CA*BTI)
BYI=BPI
BZI=-(CA*BRI-SA*BTI)
C
GIXX=SA*(SA*GIRR+CA*GIRT)+CA*(SA*GITR+CA*GITT)
GIXY=-(SA*GIRP+CA*GITP)
GIXZ=SA*(CA*GIRR-SA*GIRT)+CA*(CA*GITR-SA*GITT)
GIYX=-(SA*GIPR+CA*GIPT)
GIYY=GIPP
GIYZ=-(CA*GIPR-SA*GIPT)
GIZX=CA*(SA*GIRR+CA*GIRT)-SA*(SA*GITR+CA*GITT)
GIZY=-(CA*GIRP-SA*GITP)
GIZZ=CA*(CA*GIRR-SA*GIRT)-SA*(CA*GITR-SA*GITT)
C
C
BXE=-(SA*BRE+CA*BTE)
BYE=BPE
BZE=-(CA*BRE-SA*BTE)
C
GEXX=SA*(SA*GERR+CA*GERT)+CA*(SA*GETR+CA*GETT)
GEXY=-(SA*GERP+CA*GETP)
GEXZ=SA*(CA*GERR-SA*GERT)+CA*(CA*GETR-SA*GETT)
GEYX=-(SA*GEPR+CA*GEPT)
GEYY=GEPP
GEYZ=-(CA*GEPR-SA*GEPT)
GEZX=CA*(SA*GERR+CA*GERT)-SA*(SA*GETR+CA*GETT)
GEZY=-(CA*GERP-SA*GETP)
C
GEZZ=CA*(CA*GERR-SA*GERT)-SA*(CA*GETR-SA*GETT)
C
C
C FILL OUTPUT ARRAYS
C
C
POSOUT(1)=GALT
POSOUT(2)=GLAT
POSOUT(3)=GLON
C
C
BI(1)=BXI
BI(2)=BYI
BI(3)=BZI
C
GI(1,1)=GIXX
GI(1,2)=GIXY
GI(1,3)=GIXZ
GI(2,1)=GIYX
GI(2,2)=GIYY
GI(2,3)=GIYZ
GI(3,1)=GIZX
GI(3,2)=GIZY
GI(3,3)=GIZZ
C
C
BE(1)=BXE
BE(2)=BYE
BE(3)=BZE
C
GE(1,1)=GEXX
226
GE(1,2)=GEXY
GE(1,3)=GEXZ
GE(2,1)=GEYX
GE(2,2)=GEYY
GE(2,3)=GEYZ
GE(3,1)=GEZX
GE(3,2)=GEZY
GE(3,3)=GEZZ
C
C
OTIM=TIM
OALT=POSIN(1)
OLAT=POSIN(2)
OLON=POSIN(3)
C
C
RETURN
ENDIF
C
C
STOP
END
C
C***********************************************************************
C
C
C SUBROUTINE LAMDLA (LAMDA TRANSFORMATION MATRIX AND DERIVATES)
C
C
C***********************************************************************
C
C
C PROGRAMMED BY: JOHN M. QUINN AND DONALD L. SHIEL 8/8/86
C GEOMAGNETICS DIVISION CODE GM
C U.S. NAVAL OCEANOGRAPHIC OFFICE
C STENNIS SPACE CEENTER (SSC), MS 39522-5001
C
C
C***********************************************************************
C
C
C PURPOSE: THIS ROUTINE COMPUTES ELEMENTS OF THE MAGNETIC
C TRANSFORMATION MATRIX LAMDA AND ITS GRADIENT AT THE
C OBSERVATION POINT (X,Y,Z). THE MATRICES ARE ASSOCIATED
C WITH A RECTANGULAR PRISM LOCATED AT (XB,YB,ZB) RELATIVE
C TO SOME ORIGIN LOCATED AT THE OCEAN SURFACE AT THE
C LOWER LEFT-HAND CORNER OF THE SURVEY AREA. THE PRISM
C HAS DIMENSIONS (LAMX,LAMY,LAMZ) AND ITS ORIENTATION IS
C DESCRIBED BY EULER ANGLES (ALPHA,BETA,GAMA)
C CORRESPONDING TO YAW, PITCH AND ROLL ACCORDING TO THE
C 3-2-1 CONVENTION, RELATIVE TO THE USUAL GEODETIC
C COORDINATES FOR WHICH X=NORTH, Y=EAST AND Z=DOWN.
C
C
C
C***********************************************************************
C
C
C REFERENCE: FORWARD AND INVERSE GEOMAGNETIC AND GRAVITATIONAL
C FIELD MODELING OF THE OCEANIC CRUST, BY: JOHN M.
C QUINN AND DONALD L. SHIEL (IN PRESS)
C
C
C***********************************************************************
C
C
C PARAMETER DESCRIPTIONS:
C
C X - INERTIAL X (NORTH) COORDINATE OF OBSERVATION POINT (KM)
C Y - INERTIAL Y (EAST) COORDINATE OF OBSERVATION POINT (KM)
C Z - INERTIAL Z (DOWN) COORDINATE OF OBSERVATION POINT (KM)
C XB - INERTIAL X (NORTH) COORDINATE OF CENTER OF PRISM (KM)
C YB - INERTIAL Y (EAST) COORDINATE OF CENTER OF PRISM (KM)
C ZB - INERTIAL Z (DOWN) COORDINATE OF CENTER OF PRISM (KM)
C XP - PRISM FIXED X-AXIS COORD. OF OBSERVATION POINT (KM)
C YP - PRISM FIXED Y-AXIS COORD. OF OBSERVATION POINT (KM)
C ZP - PRISM FIXED Z-AXIS COORD. OF OBSERVATION POINT (KM)
C LAMX - DIMENSION OF PRISM ALONG X-AXIS (KM)
C LAMY - DIMENSION OF PRISM ALONG Y-AXIS (KM)
C LAMZ - DIMENSION OF PRISM ALONG Z-AXIS (KM)
C
C
C************************************************************************
C
C
C NOTE: THE ORIGINAL INTENTION WAS TO ALLOW THE PRISM TO BE ROTATED.
C
C HOWEVER, THE ROTATION OPTION WAS UNNECESSARY AND NEVER USED
C (I.E., ALL OF ROTATION ANGLES (ALPHA BETA, AND GAMMA) WERE SET
C TO ZERO. LATER IT WAS DETERMINED THAT THE ROTATION
C ALGORITHM WAS INCORRECT. IT HAS NEVER BEEN FIXED, SINCE IT
C WAS NEVER NEEDED. SO THE ROTATION ANGLES ARE ALL JUST SET TO
C ZERO.
C ALPHA - YAW ROTATION ANGLE ABOUT Z-XIS OF PRISM (DEG.)
C BETA - PITCH ROTATION ANGLE ABOUT Y-AXIS OF PRISM (DEG.)
C GAMA - ROLL RATATION ANGLE ABOUT X-AXIS OF PRISM (DEG.)
C ROT - ROTATION MATRIX FROM INERTIAL TO PRISM FIXED COORDS.
C ROTI - INVERSE ROT. MATRIX: PRISM-FIXED TO INERTIAL COORDS.
C LAMDA - MAGNETIC FIELD TRANSFORMATION MATRIX (DIMENSIONLESS)
C DLAMDA - GRADIENT OF LAMDA MATRIX [KM**(-1)]
C
C
C NOTE:
C
C * IS A MULTIPLICATION OPERATOR
C ** IS A RAISE TO THE POWER OPERATOR
C
C
C************************************************************************
C
C
C NOTE: THE PRISM IS ROTATED THROUGH EULER ANGLES ALPHA,
C BETA AND GAMA. THESE ANGLES DEFINE A NET
C ROTATION R(GAMA,BETA,ALPHA). THESE ROTATION ANGLES
C ARE DEFINED IN ACCORDANCE WITH THE 3-2-1 CONVENTION
C THAT IS IN GENERAL USE BY BRITISH AND AMERICAN
C AERODYNAMICISTS. IN THIS CONVENTION THE ANGLE ALPHA
C CORRESPONDS TO A COUNTERCLOCKWISE ROTATION ABOUT THE
C
C POSITIVE Z-AXIS, THE ANGLE BETA CORRESPONDS TO A
C CLOCKWISE ROTATION ABOUT THE NEW Y-AXIS AND THE ANGLE
C GAMMA CORRESPONDS TO A COUNTERCLOCKWISE ROTATION ABOUT THE
C FINAL X AXIS. THE CONSECUTIVE ROTATIONS MUST BE PERFORMED
C IN THE ABOVE ORDER. THE INVERSE ROTATION MUST BE
C PERFORMED IN THE REVERSE ORDER.
C
C THE INERTIAL COORDINATE SYSTEM IS REFERENCED TO AN
C ORIGIN AT THE LOWER LEFT-HAND CORNER OF THE SURVEY AREA
C
C THE PRISM FIXED COORDINATE SYSTEM IS REFERENCED TO AN
C ORIGIN THAT IS AT THE CENTER OF THE PRISM. THE PRISM
C FIXED COORDINATES ROTATE WITH THE PRISM RELATIVE TO
C THE INERTIAL (I.E. GEODETIC) COORDINATES.
C
C
C NOTE: The prism rotation part of tis routine is incorrect,
C so just set the rotation angles, ALPHA, BETA,
C and GAMA to Zero. Contact Donald Shiel at the
C U. S. Naval Oceanogaphic Office for a corrected
C version. I have a corrected copy somewhere, but it is
C burried somewhere in one of my archive files (JMQ).
C
C*********************************************************************
C
C
SUBROUTINE LAMDLA(X,Y,Z,LAMDA,DLAMDA)
C
C
C*********************************************************************
C
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION LAMX,LAMY,LAMZ,LAMDA(3,3),DLAMDA(3,3,3)
DOUBLE PRECISION LAMDAP(3,3),DLAMDP(3,3,3)
DIMENSION ROT(3,3),ROTI(3,3)
COMMON /PRISM/ LAMX,LAMY,LAMZ,XB,YB,ZB,ALPHA,BETA,GAMA
C
C
C*********************************************************************
C
C
PI=3.141592654D0
ALPH=ALPHA*PI/180.D0
BET=BETA*PI/180.D0
GAM=GAMA*PI/180.D0
C
C SET THE ROTATION ANGLES TO ZERO SINCE THE ROTATION ALGORITHM IS
C INCORRECT (I.E., DO NOT ROTATE THE PRISM). IT CAN BE EASILY FIXED
C IF NECESSARY.
C
ALPH=0.0D0
BET=0.0D0
GAMA=0.0D0
C
C COMPUTE FORWARD AND INVERSE ROTATION MATRICES
C
ROT(1,1)=+DCOS(ALPH)*DCOS(BET)
ROT(1,2)=+DSIN(ALPH)*DCOS(BET)
ROT(1,3)=-DSIN(BET)
ROT(2,1)=-DSIN(ALPH)*DCOS(GAM)+DCOS(ALPH)*DSIN(BET)*DSIN(GAM)
ROT(2,2)=+DCOS(ALPH)*DCOS(GAM)+DSIN(ALPH)*DSIN(BET)*DSIN(GAM)
C
ROT(2,3)=+DCOS(BET)*DSIN(GAM)
ROT(3,1)=+DSIN(ALPH)*DSIN(GAM)+DCOS(ALPH)*DSIN(BET)*DCOS(GAM)
ROT(3,2)=-DCOS(ALPH)*DSIN(GAM)+DSIN(ALPH)*DSIN(BET)*DCOS(GAM)
ROT(3,3)=+DCOS(BET)*DCOS(GAM)
C
ROTI(1,1)=+DCOS(ALPH)*DCOS(BET)
ROTI(1,2)=+DCOS(ALPH)*DSIN(BET)*DSIN(GAM)-DSIN(ALPH)*DCOS(GAM)
ROTI(1,3)=+DCOS(ALPH)*DSIN(BET)*DCOS(GAM)+DSIN(ALPH)*DSIN(GAM)
ROTI(2,1)=+DSIN(ALPH)*DCOS(BET)
ROTI(2,2)=+DSIN(ALPH)*DSIN(BET)*DSIN(GAM)+DCOS(ALPH)*DCOS(GAM)
ROTI(2,3)=+DSIN(ALPH)*DSIN(BET)*DCOS(GAM)-DCOS(ALPH)*DSIN(GAM)
ROTI(3,1)=-DSIN(BET)
ROTI(3,2)=+DCOS(BET)*DSIN(GAM)
ROTI(3,3)=+DCOS(BET)*DCOS(GAM)
C
C ROTATE FROM INERTIAL TO PRISM-FIXED COORDINATES
C
XP=ROT(1,1)*(X-XB)+ROT(1,2)*(Y-YB)+ROT(1,3)*(Z-ZB)
YP=ROT(2,1)*(X-XB)+ROT(2,2)*(Y-YB)+ROT(2,3)*(Z-ZB)
ZP=ROT(3,1)*(X-XB)+ROT(3,2)*(Y-YB)+ROT(3,3)*(Z-ZB)
C
C COMPUTE THE LAMDA TRANSFORMATION MATRIX AND ITS
C GRADIENT IN PRISM FIXED COORDINATES
C
C
C CLEAR ARRAYS
C
DO 30 I=1,3
DO 20 J=1,3
LAMDAP(I,J)=0.0D0
DO 10 K=1,3
233
DLAMDP(I,J,K)=0.0D0
10 CONTINUE
20 CONTINUE
30 CONTINUE
DO 60 JB=1,2
IF (JB .EQ. 1) XL=+LAMX/2.0D0
IF (JB .EQ. 2) XL=-LAMX/2.0D0
IF (JB .EQ. 1) SX=+1.0D0
IF (JB .EQ. 2) SX=-1.0D0
XD=XP-XL
DO 50 KB=1,2
IF (KB .EQ. 1) YL=+LAMY/2.0D0
IF (KB .EQ. 2) YL=-LAMY/2.0D0
IF (KB .EQ. 1) SY=+1.0D0
IF (KB .EQ. 2) SY=-1.0D0
YD=YP-YL
DO 40 LB=1,2
IF (LB .EQ. 1) ZL=+LAMZ/2.0D0
IF (LB .EQ. 2) ZL=-LAMZ/2.0D0
IF (LB .EQ. 1) SZ=+1.0D0
IF (LB .EQ. 2) SZ=-1.0D0
ZD=ZP-ZL
S=SX*SY*SZ
R2=XD**2+YD**2+ZD**2
R=DSQRT(R2)
IF (R+ZD .EQ. 0.0D0) PRINT *, ' A'
IF ((R-YD)/(R+YD) .EQ. 0.0D0) PRINT *, ' B'
IF ((R-XD)/(R+XD) .EQ. 0.0D0) PRINT *, ' C'
IF (R+XD .EQ. 0.0D0) PRINT *, ' D'
LAMDAP(1,1)=LAMDAP(1,1)+S*DATAN(YD*ZD/(XD*R))
LAMDAP(1,2)=LAMDAP(1,2)-S*DLOG(ABS(R+ZD))
C
LAMDAP(1,3)=LAMDAP(1,3)-S*DLOG(ABS(R+YD))
LAMDAP(2,2)=LAMDAP(2,2)+S*DATAN(XD*ZD/(YD*R))
LAMDAP(2,3)=LAMDAP(2,3)-S*DLOG(ABS(R+XD))
LAMDAP(3,3)=LAMDAP(3,3)+S*DATAN(XD*YD/(ZD*R))
C
C
DLAMDP(1,1,2)=DLAMDP(1,1,2)-(S/R)*(XD/(R+ZD))
DLAMDP(1,1,3)=DLAMDP(1,1,3)-(S/R)*(XD/(R+YD))
DLAMDP(1,2,2)=DLAMDP(1,2,2)-(S/R)*(YD/(R+ZD))
DLAMDP(1,2,3)=DLAMDP(1,2,3)-(S/R)
DLAMDP(1,3,3)=DLAMDP(1,3,3)-(S/R)*(ZD/(R+YD))
DLAMDP(2,2,3)=DLAMDP(2,2,3)-(S/R)*(YD/(R+XD))
DLAMDP(2,3,3)=DLAMDP(2,3,3)-(S/R)*(ZD/(R+XD))
40 CONTINUE
50 CONTINUE
60 CONTINUE
LAMDAP(2,1)=LAMDAP(1,2)
LAMDAP(3,1)=LAMDAP(1,3)
LAMDAP(3,2)=LAMDAP(2,3)
C
C
DLAMDP(1,2,1)=DLAMDP(1,1,2)
DLAMDP(1,3,1)=DLAMDP(1,1,3)
DLAMDP(1,3,2)=DLAMDP(1,2,3)
DLAMDP(2,1,1)=DLAMDP(1,1,2)
DLAMDP(2,1,2)=DLAMDP(1,2,2)
DLAMDP(2,1,3)=DLAMDP(1,2,3)
DLAMDP(2,2,1)=DLAMDP(1,2,2)
DLAMDP(2,3,1)=DLAMDP(1,2,3)
DLAMDP(2,3,2)=DLAMDP(2,2,3)
DLAMDP(3,1,1)=DLAMDP(1,1,3)
C
DLAMDP(3,1,2)=DLAMDP(1,2,3)
DLAMDP(3,1,3)=DLAMDP(1,3,3)
DLAMDP(3,2,1)=DLAMDP(1,2,3)
DLAMDP(3,2,2)=DLAMDP(2,2,3)
DLAMDP(3,2,3)=DLAMDP(2,3,3)
DLAMDP(3,3,1)=DLAMDP(1,3,3)
DLAMDP(3,3,2)=DLAMDP(2,3,3)
DLAMDP(1,1,1)=-DLAMDP(1,3,3)-DLAMDP(1,2,2)
DLAMDP(2,2,2)=-DLAMDP(2,3,3)-DLAMDP(1,1,2)
DLAMDP(3,3,3)=-DLAMDP(2,2,3)-DLAMDP(1,1,3)
C
C
C ROTATE MAGNETIC FIELD COMPONENTS FROM PRISM FIXED TO
C INERTIAL COORDINATES
C
C
C
C
DO 110 I=1,3
DO 100 J=1,3
LAMDA(I,J)=0.0D0
DO 70 K=1,3
LAMDA(I,J)=LAMDA(I,J)+ROTI(I,K)*LAMDAP(K,J)
70 CONTINUE
DO 90 L=1,3
DLAMDA(I,J,L)=0.0D0
DO 80 K=1,3
DLAMDA(I,J,L)=DLAMDA(I,J,L)+ROTI(I,K)*DLAMDP(K,J,L)
80 CONTINUE
90 CONTINUE
100 CONTINUE
C
110 CONTINUE
RETURN
END
237
C***********************************************************************
C
C
C SUBROUTINE MATINV (MATRIX INVERSION)
C
C
C***********************************************************************
C
C
C REFERENCE:
C DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES
C BY PHILLIP R. BEVINGTON, MC GRAW-HILL 1969
C APPENDIX B
C
C
C PURPOSE:
C INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT
C
C
C USAGE:
C CALL MATINV(ARRAY,NORDER,DET)
C
C
C DESCRIPTION OF PARAMETERS:
C ARRAY - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE
C NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT)
C DET - DETERMINANT OF INPUT MATRIX
C
C
C COMMENTS:
C DIMENSION STATEMENT VALID FOR NORDER UP TO 500
C
C
C
C***********************************************************************
C
C
SUBROUTINE MATINV(ARRAY,NORDER,DET)
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION ARRAY(NORDER,NORDER),IK(500),JK(500)
10 DET=1.0D0
11 DO 100 K=1,NORDER
C
C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
AMAX=0.0D0
21 DO 30 I=K,NORDER
DO 30 J=K,NORDER
23 IF (DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30
24 AMAX=ARRAY(I,J)
IK(K)=I
JK(K)=J
30 CONTINUE
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
31 IF (AMAX) 41,32,41
32 DET=0.0D0
GO TO 140
41 I=IK(K)
IF (I-K) 21,51,43
C
43 DO 50 J=1,NORDER
SAVE=ARRAY(K,J)
ARRAY(K,J)=ARRAY(I,J)
50 ARRAY(I,J)=-SAVE
51 J=JK(K)
IF (J-K) 21,61,53
53 DO 60 I=1,NORDER
SAVE=ARRAY(I,K)
ARRAY(I,K)=ARRAY(I,J)
60 ARRAY(I,J)=-SAVE
C
C ACCUMULATE ELEMENTS OF INVERSE MATRIX
C
61 DO 70 I=1,NORDER
IF (I-K) 63,70,63
63 ARRAY(I,K)=-ARRAY(I,K)/AMAX
70 CONTINUE
71 DO 80 I=1,NORDER
DO 80 J=1,NORDER
IF (I-K) 74,80,74
74 IF (J-K) 75,80,75
75 ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
80 CONTINUE
81 DO 90 J=1,NORDER
IF (J-K) 83,90,83
83 ARRAY(K,J)=ARRAY(K,J)/AMAX
90 CONTINUE
ARRAY(K,K)=1.0D0/AMAX
100 DET=DET*AMAX
C
C RESTORE ORDERING OF MATRIX
240
C
101 DO 130 L=1,NORDER
K=NORDER-L+1
J=IK(K)
IF (J-K) 111,111,105
105 DO 110 I=1,NORDER
SAVE=ARRAY(I,K)
ARRAY(I,K)=-ARRAY(I,J)
110 ARRAY(I,J)=SAVE
111 I=JK(K)
IF (I-K) 130,130,113
113 DO 120 J=1,NORDER
SAVE=ARRAY(K,J)
ARRAY(K,J)=-ARRAY(I,J)
120 ARRAY(I,J)=SAVE
130 CONTINUE
140 RETURN
END