C Appendix A Fortran routine called s_ab2h.f that applies the C bin-size, axis polynomial, orthogonality, and attitude coefficients C to basic TMGS analog and bin data. C C________________________________________________________________ C C SUBROUTINE S _ A B 2 H C________________________________________________________________ C C SUBROUTINE S_AB2H RETURNS THE MAGNETIC FIELD COMPONENTS AND A C GRADIENT TENSOR IN A FIXED REFERENCE FRAME USING DATA AND C COEFFICIENTS FROM THE USGS TMGS PROTOTYPE SYSTEM. C C THIS SUBROUTINE IS DESIGNED PRIMARILY TO SHOW HOW THE C COEFFICIENTS DERIVED FROM THE YUMA, 13MAR03 SPIN CALIBRATION C ARE TO BE APPLIED TO THE RAW TMGS DATA. THE PRELIMINARY C COEFFICIENTS MAY BE FOUND IN AN ASCII TEXT FILE CALLED: C C "tmgs_Coefs_n0430_yuma_prelim.txt". C C C INPUT ARGUMENTS (SIGNALS) C VANA - REAL*8. ARRAY OF DIMENSIONS (3,4) CONTAINING THE C ANALOG FIELD MEASUREMENTS FROM THE 3 AXES OF C THE 4 MAGNETOMETERS. ROWS 1,2,3 CORRESPOND TO THE C X,Y,Z AXES AND COLUMNS 1,2,3,4 CORRESPOND TO THE C MAGNETOMETER NUMBERS. IT IS EXPECTED THAT THESE C NUMBERS WILL HAVE BEEN CORRECTED FOR SPIKES AND C NOISE, AND THAT ANALOG FILTER RESPONSES WILL HAVE C BEEN DECONVOLVED (PARTICULARLY IMPORTANT IMMEDIATELY C AFTER A BIN STEP.) PRIOR TO CALLING THIS SUBROUTINE. C TYPICALLY, THESE SIGNALS WILL RANGE FROM +-5 BUT MAY C RAIL AT +-6.55 VOLTS. THESE SIGNALS TRANSLATE TO C MAGNETIC FIELD NOMINALLY AT ABOUT 100 nT/VOLT (SEE C ARGUMENT CANA BELOW). C VBIN - REAL*8. ARRAY OF DIMENSIONS (3,4) CONTAINING THE C BIN NUMBERS FROM THE 3 AXES OF THE 4 C MAGNETOMETERS. ROWS 1,2,3 CORRESPOND TO X,Y,Z AXES C AND COLUMNS 1,2,3,4 CORRESPOND TO THE MAGNETOMETER C NUMBERS. IT IS POSSIBLE TO HAVE CORRECTED THESE C NUMBERS FOR BIN-CURRENT VARIANCES AND THEREFORE THEY C ARE NOT NECESSARILY INTEGERS. TYPICALLY, BIN C NUMBERS WILL RANGE FROM +-150 BUT MAY RAIL AT +-255. C BIN NUMBERS INDI\CATE THE AMOUNT OF OFFSET IN THE C ANALOG FIELD MEASUREMENT VOLTAGE AND TRANSLATE TO C VOLTS NOMINALLY AT ABOUT 5.0 VOLTS/BIN# FOR MAGS 1 C AND 2 AND 3.3 VOLTS/BIN FOR MAGS 3 AND 4. (SEE C ARGUMENT CBIN BELOW). C C INPUT ARGUMENTS (COEFFICIENTS) C CBIN - REAL*8. ARRAY OF DIMENSIONS (3,4) CONTAINING THE C AVERAGE BIN SIZES (S) (VOLTS/BIN#). ROWS 1,2,3 C CORRESPOND TO X,Y,Z AXES AND COLUMNS 1,2,3,4 C CORRESPOND TO THE MAGNETOMETER NUMBERS. TYPICALLY, C FOR MAGS 1 AND 2, CBIN IS ABOUT 5.0 AND C FOR MAGS 3 AND 4 IT'S ABOUT 3.3 VOLTS/BIN#. C CANA - REAL*8. ARRAY OF DIMENSIONS (4,3,4) CONTAINING THE C POLYNOMIAL COEFFICIENTS (A,B,C,D) THAT TRANSLATE C ANALOG VOLTAGE TO MAGNETIC FIELD . ROWS 1,2,3,4 CORRESPOND TO C THE 3RD ORDER, 2ND ORDER, LINEAR, AND CONSTANT TERMS C RESPECTIVELY OF THE POLYNOMIAL. COLUMNS 1,2,3 C CORRESPOND TO X,Y,Z AXES AND LEVELS 1,2,3,4 C CORRESPOND TO THE MAGNETOMETER NUMBERS. TYPICALLY, C THE 3RD ORDER TERM IS EXTREMELY SMALL, THE 2ND C ORDER TERM IS VERY SMALL, THE LINEAR TERM IS NEAR C 100, AND THE CONSTANT TERM IS SMALL. THESE C COEFFICIENTS ARE AUGMENTED BY THE SCALE DIVISOR THAT C WAS APPLIED IN THEIR DERIVATION (SEE ARGUMENT C SCALE). C SCALE - REAL*8. SCALAR NUMBER (SCALE) THAT WAS USED DURING C DERIVATION OF THE 3RD-ORDER-POLYNOMIAL COEFFICIENTS C DESCRIBED IN ARGUMENT CANA. THIS SUBROUTINE WILL C CORRECTLY APPLY SCALE TO THE COEFFICIENTS. BUT IF C IT IS DESIRED TO APPLY IT BEFORE THE SUBROUTINE C CALL, SIMPLY DIVIDE THE 3RD-ORDER TERM BY SCALE^2, C DIVIDE THE SECOND ORDER TERM BY SCALE, AND MULTIPLY C THE CONSTANT TERM BY SCALE; THEN SET SCALE, ITSELF, C EQUAL TO 1. C CORT - REAL*8. ARRAY OF DIMENSIONS (3,4) CONTAINING THE C ORTHOGONALITY CORRECTION ANGLES (DEGREES). ROWS C 1,2,3 CORRESPOND TO ANGLES (ALP), (BET), (GAM), AND C COLUMNS 1,2,3,4 CORRESPOND TO THE MAGNETOMETER C NUMBERS. ALPHA IS THE ANGLE BETWEEN THE X AND Y C AXIS. BETA IS THE ANGLE BETWEEN THE Y AND Z AXIS. C AND, GAMMA IS THE ANGLE BETWEEN THE X AND Z AXIS. C TYPICALLY, EACH OF THESE ANGLES VARIES FROM 88 TO 92 C DEGREES. C CATT - REAL*8. ARRAY OF DIMENSIONS (3,4) CONTAINING THE C ATTITUDES (DEGREES) OF THE MAGS RELATIVE TO MAG1. C ROWS 1,2,3 CORRESPOND TO ANGLES (DEL), (EPS), (ZET), C AND COLUMNS 1,2,3,4 CORRESPOND TO THE MAGNETOMETER C NUMBERS. IN TRANSFORMING FROM EACH MAGNETOMETER'S C PHYSICAL ATTITUDE TO THE MAG1 ATTITUDE, THERE ARE 3 C ROTATIONS. THE POSITIVE SENSE IS ALWAYS CLOCKWISE C ABOUT THE ROTATION AXIS AS VIEWED FROM THE ORIGIN. C DELTA IS THE FIRST ROTATION ANGLE AND PROCEEDS ABOUT C THE Z''' AXIS. EPSILON IS THE SECOND ROTATION ANGLE C AND PROCEEDS ABOUT THE Y'' AXIS. ZETA IS THE THIRD C ROTATION ANGLE AND PROCEEDS ABOUT THE Z' AXIS. C THEREFORE, THE UNPRIMED SYSTEM IS THE MAG1 ATTITUDE. C TYPICALLY, EACH OF THESE ANGLES VARIES BY +-2 C DEGREES ABOUT THE IDEAL. THE IDEAL ANGLES ARE AS C FOLLOWS: C C ROTATION MAG1 MAG2 MAG3 MAG4 C DEL 0 0 +120 -120 C EPS 0 -109.471 -109.471 -109.471 C ZET 0 180 +60 -60 C C C OUTPUT ARGUMENTS: C HFLD - REAL*8. ARRAY OF DIMENSIONS (3,4) CONTAINING THE C MAGNETIC FIELD COMPONENTS (nT). ROWS 1,2,3 C CORRESPOND TO X,Y,Z AXES AND COLUMNS 1,2,3,4 C CORRESPOND TO THE MAGNETOMETER NUMBERS. THE FIELDS C ARE CALIBRATED, THE AXES ARE ORTHOGONAL, AND THE C MAGNETOMETERS ARE IN THE MAG1 COMMON REFERENCE C FRAME. THE ARRAY IS LAYED OUT PICTORIALLY AS C FOLLOWS: C X1 X2 X3 X4 C Y1 Y2 Y3 Y4 C Z1 Z2 Z3 Z4 C C GTSR - REAL*8. ARRAY OF DIMENSIONS (3,3) CONTAINING THE C MAGNETIC GRADIENT TENSOR (nT/m) DERIVED FROM HFLD, C ASSUMING A TETRAHEDRAL CONFIGURATION OF THE 4 C MAGNETOMETERS WITH 0.97 METERS DISTANCE BETWEEN C CENTERS. WITH 4 SENSORS, A DIRECT DERIVATION OF C EACH OF THE 9 COMPONENTS OF THE TENSOR IS AVAILABLE. C HOWEVER, IF A) THE FIELD IS STATIC, B) ALL OF THE C GRADIENTS ARE CONSTANT WITHIN THE VOLUME OF TESSA, C AND C) THERE IS NO MEASUREMENT NOISE, THEN THE C TENSOR WILL BE TRACELESS AND SYMMETRIC, HAVING ONLY C 5 INDEPENDENT COMPONENTS. THEREFORE, TO THE DEGREE C THAT IT IS NOT TRACELESS AND SYMMETRIC, ONE OR MORE C OF THESE CONDITIONS HAS NOT BEEN REACHED. THE ARRAY C IS LAYED OUT PICTORIALLY AS FOLLOWS: C C Gxx Gxy Gxz C Gyx Gyy Gyz C Gzx Gzy Gzz C C WHERE: C C Gij = dHi/dj = gtsr(i,j) C C C SUBROUTINE S_AB2H WRITTEN BY ROB BRACKEN, USGS. C FORTRAN 77, HP FORTRAN/9000, HP-UX RELEASE 11.0 C VERSION 1.0, 20030430 ( ORIGINAL CODE ). C C subroutine s_ab2h(vana,vbin, cbin,cana,scale,cort,catt & ,hfld,gtsr ) C C C DECLARATIONS C C INPUT ARGUMENTS (SIGNALS) real*8 vana(3,4),vbin(3,4) C C INPUT ARGUMENTS (COEFFICIENTS) real*8 cbin(3,4),cana(4,3,4),scale,cort(3,4),catt(3,4) C C OUTPUT ARGUMENTS (FIELDS AND TENSOR) real*8 hfld(3,4),gtsr(3,3) C C C INTERNAL VARIABLES C C AXIS VOLTAGE RECONSTRUCTION C real*8 va(3,4) C C FIELD CORRECTION real*8 hf(3,4) C C ORTHOGONALITY CORRECTION real*8 hfl(3,4) real*8 alp,bet,gam real*8 calp,salp, cbet, cgam real*8 coseta C C ATTITUDE ROTATION TO MAG1 REFERENCE FRAME real*8 del,eps,zet real*8 cdel,sdel, ceps,seps, czet,szet C C GRADIENT TENSOR CALCULATIONS real*8 edge, face, cent real*8 edge2,face3,cent4 C C C RECONSTRUCT AXIS VOLTAGE FROM ANALOG VOLTAGE AND BIN NUMBER C do k=1,4 do j=1,3 va(j,k)= & vana(j,k) & +vbin(j,k) * cbin(j,k) enddo enddo C C C CONVERT AXIS VOLTAGE TO CORRECTED MAGNETIC FIELD (nT) C do k=1,4 do j=1,3 hf(j,k)= & 1 * cana(4,j,k) * scale & +va(j,k) * cana(3,j,k) * 1 & +va(j,k)**2 * cana(2,j,k) / scale & +va(j,k)**3 * cana(1,j,k) / scale**2 enddo enddo C C C APPLY ORTHOGONALITY CORRECTION C do k=1,4 C C X TO Y ANGLE alp=cort(1,k) calp=dcosd(alp) salp=dsind(alp) C C Y TO Z ANGLE bet=cort(2,k) cbet=dcosd(bet) C C X TO Z ANGLE gam=cort(3,k) cgam=dcosd(gam) C C COSINE OF THE ANGLE BETWEEN Z AND Z-PRIME coseta=dsqrt(1.d0-cbet**2-cgam**2) C C X-AXIS CANNOT "MOVE" hfl(1,k)=hf(1,k) C C Y-AXIS CAN "MOVE" ONLY IN THE X-Y PLANE hfl(2,k)= & ( -hf(1,k)*( calp ) & +hf(2,k) ) / salp C C Z-AXIS CAN "MOVE" ANYWHERE hfl(3,k)= & ( hf(1,k)*(cbet*calp/salp-cgam) & -hf(2,k)*(cbet /salp ) & +hf(3,k) ) / coseta C enddo C C C APPLY ATTITUDE XFRM (ROTATION) TO COMMON (MAG1) REFERENCE FRAME C do k=1,4 C C ROTATION 1 IS ABOUT Z''' del=catt(1,k) cdel=dcosd(del) sdel=dsind(del) C C ROTATION 2 IS ABOUT Y'' eps=catt(2,k) ceps=dcosd(eps) seps=dsind(eps) C C ROTATION 3 IS ABOUT Z' zet=catt(3,k) czet=dcosd(zet) szet=dsind(zet) C C FIND THE ROTATED COMPONENTS hfld(1,k)= & hfl(1,k)*( czet*ceps*cdel-szet*sdel ) & +hfl(2,k)*( czet*ceps*sdel+szet*cdel ) & +hfl(3,k)*(-czet*seps ) hfld(2,k)= & hfl(1,k)*(-szet*ceps*cdel-czet*sdel ) & +hfl(2,k)*(-szet*ceps*sdel+czet*cdel ) & +hfl(3,k)*( szet*seps ) hfld(3,k)= & hfl(1,k)*( seps*cdel ) & +hfl(2,k)*( seps*sdel ) & +hfl(3,k)*( ceps ) enddo C C C FIND GRAD TENSOR WHERE: gtsr(i,j) = Gij = dHi/dj C C Note: With 4 sensors, a direct measurement of all 9 C components of the tensor is available. However, if C a) the field is static, b) all of the gradients are C constant within the volume of TESSA, and c) there is C no measurement noise, then the tensor will be C traceless and symmetric, having only 5 independent C components. Therefore, to the degree that it is not C traceless and symmetric, one or more of these C conditions has not been reached. C call s_h2g(hfld, gtsr) C C C EXIT PROCEDURE C 990 return end