text:       SUBROUTINE AXOA
        IMPLICIT None
  C
  C 1.    AXOA
  C
  C 1.1   AXOA PROGRAM SPECIFICATION
  C
  C 1.1.1 AXOA adds entries to the table of contents for the Axis Offset 
  C       text message, the axis offset flow control message, the axis
  C       offset partial derivatives and contributions, the feedbox 
  C       rotation routine test message, the feedbox rotation angles, and
  C       the group and phase delay rate corrections. 
  C
  C 1.2   AXOA PROGRAM INTERFACE
  C
  C 1.2.4 DATA BASE ACCESS -
  C           ACCESS CODES ADDED:
  C              1.  'AXO MESS'  -  THE DATA BASE ACCESS CODE FOR THE 
  C                                 AXIS OFFSET MODULE TEXT MESSAGE.
  C              2.  'AXO CFLG'  -  THE DATA BASE ACCESS CODE FOR THE AXIS
  C                                 OFFSET CONTROL FLAG MESSAGE.
  C              3.  'AXO PART'  -  THE DATA BASE ACCESS CODE FOR THE AXIS 
  C                                 OFFSET MODULE PARTIAL DERIVATIVES ARRAY.
  C                                 (Using new, simple axis offset model)
  C              4.  'AXO CONT'  -  The database access code for the new axis
  C                                 offset module contributions array.
  C                                 (Using the new, simple axis offset model.)
  C              5.  'AXO2CONT'  -  The database access code for the correction
  C                                 to convert to the more complicated new axis
  C                                 offset model.
  C              6.  'AXIS OLD'  -  The database access code for the correction
  C                                 to convert to the old (Calc 7.6 and earlier)
  C                                 axis offset model.
  C              7.  'PAN MESS'  -  The database access code for the
  C                                 feedhorn rotation routine text message.
  C              8.  'PARANGLE'  -  The database access code for the
  C                                 feedhorn rotation angles.
  C              9.  'FEED.COR'  -  The database access code for the feedbox 
  C                                 rotation corrections for group delay and
  C                                 phase delay rate.
  C 
  C 1.2.6 SUBROUTINE INTERFACE -
  C             CALLER SUBROUTINES: TOCUP
  C             CALLED SUBROUTINES: ADDA, ADDR
  C
  C 1.2.9 PROGRAMMER - DALE MARKHAM  01/13/77
  C                    PETER DENATALE 07/11/77
  C                    BRUCE SCHUPLER 09/16/77
  C                    SAVITA GOEL    06/03/87 (CDS FOR A900)
  C                    Jim Ryan 86.29 Character stings used.
  C                    Jim Ryan 89.12.12 UNIX-like database interface
  C                    implemented.
  C                    David Gordon 94.03.25 AXO2CONT and AXIS_OLD Lcodes added.
  C                    David Gordon 94.07.22 Changed 'AXIS_OLD' to 'AXIS OLD' to 
  C                                 keep Chopo happy.
  C                    David Gordon 98.11.12 Merged in code from subroutine
  C                                 PANA for feedbox rotation stuff. 'AXIS OLD'
  C                                 contribution removed. 
  C                    David Gordon 98.12.30 FEED.COR changed from (2,2) to
  C                                 (2,1) dimensions.
  C
  C 1.3   AXOA PROGRAM STRUCTURE
  C
  C   ADD for axis offset module text message.
        CALL ADDA (1,'AXO MESS','Axis Offset Message Definition  ',
       *     40, 1, 1 )
  C
  C   ADD for axis offset module control flag message.
        CALL ADDA (1,'AXO CFLG','Axis Offset Control flag mes def',
       *     40,1,1)
  C
  C   ADD for axis offset partial derivatives.
        CALL ADDR (2,'AXO PART','Axis Offset partial deriv. def. ',
       *     2, 2, 1 )
  C
  C   ADD for new axis offset contributions (one-term formula).
        CALL ADDR (2,'AXO CONT','New Axis Offset Contributions   ',
       *     2, 2, 1 )
  C
  C   ADD for correction to convert to the alternate new axis offset
  C       (three-term formula).
        CALL ADDR (2,'AXO2CONT','Other new axis offset correction',
       *     2, 1, 1 )
  C
  C   Delete for old axis offset correction.
        CALL DELR (2,'AXIS OLD')
  C   ADD for correction to convert back to the old (pre Calc 8.0) axis offset.
  c     CALL ADDR (2,'AXIS OLD','Old axis offset correction      ',
  c    *     2, 1, 1 )
  C
  C     ADDA for feedbox rotation (paralactic angle) text message.
        CALL ADDA (1,'PAN MESS','Feedhorn rot. angle mod. ident. ',
       *     40, 1, 1 )
  C
  C     ADD for feedbox rotation angle (paralactic angle).
        CALL ADDR (2,'PARANGLE','Feedhorn rot. angle, ref and rem',
       *     2, 1, 1 )
  C
  C     ADD for feedbox corrections for standard observables.
        CALL ADDR (2,'FEED.COR','Feedhorn corr. in CORFIL scheme ',
       *     2, 1, 1 )
  C
        RETURN
        END
  C
  C******************************************************************************
        SUBROUTINE AXOI
        IMPLICIT None
  C
  C 2.    AXOI
  C
  C 2.1   AXOI PROGRAM SPECIFICATION
  C
  C 2.1.1 AXOI is the Axis Offset and Feedbox Rotation Module input and 
  C       initialization section.
  C
  C 2.2   AXOI PROGRAM INTERFACE
  C
  C 2.2.2 COMMON BLOCKS USED -
  C
         INCLUDE 'ccon.i'
  C      VARIABLES 'FROM'
  C        1. KAXOC - THE AXIS OFFSET MODULE FLOW CONTROL FLAG
  C                   = 0 ==> Axis offset corrections applied to theoreticals.
  C                   = 1 ==> Axis offset corrections NOT applied to theoreticals.
  C        1.  KPANC  -  The Feedbox Rotation flow control flag.
  C
  C 2.2.3 PROGRAM SPECIFICATIONS -
        INTEGER*2       LAXOM(40),      LON(40),    LOFF(40)
        CHARACTER*40  C_LAXOM(2),     C_LON(2),   C_LOFF(2)
        EQUIVALENCE(  C_LAXOM,LAXOM),(C_LON,LON),(C_LOFF,LOFF)
        INTEGER*2      LPANG_ON(40),  LPANG_OF(40)
        CHARACTER*40 C_LPANG_ON(2), C_LPANG_OF(2)
        EQUIVALENCE (  LPANG_ON, C_LPANG_ON)
        EQUIVALENCE (  LPANG_OF, C_LPANG_OF)
  C
        DATA C_LAXOM /
       .'Axis Offset Module - Last modified 98DEC',
       .'30, D. Gordon/GSFC.                    '/
        DATA C_LON  /
       .'Axis Offset Module is turned ON in CALC.',
       .'                                        '/
        DATA C_LOFF /
       .'Axis Offset Module is turned OFF in CALC',
       .'.                                       '/
  C
        DATA C_LPANG_ON  /
       .'Feedbox Rotation Angle Module, Last Modi',
       .'fied 98NOV12, D. Gordon/GSFC.           '/
        DATA C_LPANG_OF  /
       .'Feedbox Rotation angle zeroed out.      ',
       .'                                        '/
  C
  C 2.2.4 DATA BASE ACCESS -
  C           'PUT' VARIABLES:
  C              1. LAXOM(40)    - THE AXIS OFFSET MODULE TEXT MESSAGE.
  C              2. LON(40)      - THE AXIS OFFSET MODULE TURNED ON MESSAGE
  C              3. LOFF(40)     - THE AXIS OFFSET MODULE TURNED OFF MESSAGE.
  C              4. LPANG_ON(40) - The feedbox rotation (parallactic) angle
  C                                text message.
  C              5. LPANG_OF(40) - The feedbox rotation (parallactic) angle
  C                                turned OFF message.
  C           ACCESS CODES PUT:
  C              1. 'AXO MESS'   -  THE DATA BASE ACCESS CODE FOR THE
  C                                 AXIS OFFSET MODULE TEXT MESSAGE.
  C              2. 'AXO CFLG'   -  THE DATA BASE ACCESS CODE FOR THE AXIS
  C                                 OFFSET MODULE FLOW CONTROL MESSAGE.
  C              3. 'PAN MESS'   -  The data base access code for the feedbox 
  C                                 rotation text message.
  C 
  C 2.2.6 SUBROUTINE INTERFACE -
  C             CALLER SUBROUTINES: INITL
  C             CALLED SUBROUTINES: PUTA
  C
  C 2.2.9 PROGRAMMER - DALE MARKHAM  01/13/77
  C                    PETER DENATALE 07/11/77
  C                    BRUCE SCHUPLER 03/07/78
  C                    Jim Ryan 89.06.29 Character stings used.
  C                    Jim Ryan 89.12.12 UNIX-like database interface
  C                             implemented.
  C                    D. Gordon 98.11.12 Merged in code from subroutine PANI,
  C                             the feedbox rotation (parallactic) angle
  C                             initialization routine.
  C
  C 2.3   AXOI PROGRAM STRUCTURE
  C
  C     Put Axis Offset Module text message in database
        CALL PUTA ('AXO MESS      ', LAXOM, 40, 1, 1 )
  C
  C     Put Axis Offset Module flow control message into database header.
        IF (KAXOC .NE. 1) CALL PUTA('AXO CFLG      ',LON,40,1,1)
        IF (KAXOC .EQ. 1) CALL PUTA('AXO CFLG      ',LOFF,40,1,1)
  C
  C     PUT for the feedbox rotation (paralactic) angle routine text message.
        IF (KPANC .NE. 1) CALL PUTA ('PAN MESS      ', LPANG_ON,40,1,1)
        IF (KPANC .EQ. 1) CALL PUTA ('PAN MESS      ', LPANG_OF,40,1,1)
  C
  C     Normal conclusion.
        RETURN
        END
  C
  C******************************************************************************
        SUBROUTINE AXOG (KAXIS, R2000, SITLAT, STAR, TCTOCF, SITEV, AXOFF,
       *           EARTH, AZ, ELEV, STAR_ABERRATED,
       *           SITHEIGHT, axis2000, daxis2000 )
        IMPLICIT None
  C
  C 3.    AXOG
  C
  C 3.1   AXOG PROGRAM SPECIFICATION
  C
  C 3.1.1 AXOG is the geometry section of the Axis Offset Module. AXOG computes
  C       the 3 dimensional axis offset vector and its time derivative corrected
  C       for aberration and atmospheric refraction at each site.
  C       Beginning with Calc 9.0, AXOG also computes the feedbox rotation
  C       angle and the group delay and phase delay corrections, using a new
  C       algorithm that replaces subroutine PANG.
  C
  C 3.1.2 REFERENCES - NINER, EDGAR P., 'MEMORANDUM FOR RECORD', 25 OCTOBER
  C                    1967 (REFRACTION CORRECTION).
  C                    SMART, W.M., "TEXTBOOK ON SPHERICAL ASTRONOMY", 1965,
  C                    P. 68.
  C                    Gordon, David, "Axis offset VLBI delay correction", 
  C                    1994 March 7 memo, updated 1994 May 5.
  C
  C 3.2   AXOG PROGRAM INTERFACE
  C 
  C 3.2.1 CALLING SEQUENCE -
  C          INPUT VARIABLES:
  C             1. KAXIS(2)       -  THE ANTENNA AXIS TYPES AT EACH OBSERVATION
  C                                  SITE. (UNITLESS)
  C             2. R2000(3,3,3)   -  THE ROTATION MATRIX WHICH ROTATES THE CRUST
  C                                  FIXED REFERENCE FRAME TO THE J2000.0 FRAME 
  C                                  AND ITS FIRST TWO CT TIME DERIVATIVES. 
  C                                  (UNITLESS, 1/SEC, 1/SEC**2)
  C             3. SITLAT(2)      -  THE GEODETIC LATITUDES AT EACH OBSERVATION
  C                                  SITE. (RAD) 
  C             4. STAR(3)        -  THE J2000.0 SOURCE UNIT VECTOR. (UNITLESS) 
  C             5. TCTOCF(3,3,2)  -  THE ROTATION MATRIX WHICH ROTATES THE
  C                                  TOPOCENTRIC REFERENCE SYSTEM TO THE CRUST
  C                                  FIXED GEOCENTRIC REFERENCE SYSTEM AT EACH 
  C                                  OBSERVATION SITE. (UNITLESS)
  C             6. SITEV(3,2)     -  THE J2000.0 GEOCENTRIC VELOCITY VECTORS OF
  C                                  EACH OBSERVATION SITE. (M/SEC)
  C             7. AXOFF(2)       -  THE ANTENNA AXIS OFFSETS AT EACH OBSERVATION
  C                                  SITE.  (M)
  C             8. EARTH(3,3)     -  THE SSBC EARTH POSITION, VELOCITY, AND 
  C                                  ACCELERATION VECTORS. (M, M/SEC, M/SEC**2)
  C             9. ELEV(2,2)      -  The elevation angle of the source corrrected
  C                                  for aberration and its CT time derivative at
  C                                  each site (rad,rad/sec)
  C            10. AZ(2,2)        -  The azimuth angle of the source corrrected
  C                                  for aberration and its CT time derivative
  C                                  at each site (rad,rad/sec)
  C            11. STAR_ABERRATED(3,2) - The J2000 source unit vector corrected
  C                                  for annual and diurnal aberration at each 
  C                                  site. (unitless)
  C            12. SITHEIGHT(2)   -  Station heights above the geoid. (meters)
  C          OUTPUT VARIABLES:
  C             1. axis2000(3,2) -  Vector axis offset of antenna in the J2000.0
  C                                 frame (effect on baseline). First index is
  C                                 X,Y,Z (meters), second runs over sites.
  C             2. daxis2000(3,2) - Time derivative of axis2000, rate of change
  C                                 of vector axis offset of antenna in the
  C                                 J2000.0 frame (effect on baseline). First
  C                                 index is velocity, second runs over sites.
  C
  C 3.2.2 COMMON BLOCKS USED -
  C 
        Real*8  DAXOP(2,2), DCOMP(2,2), SURREF, RICHM(2), DCOMP_new(2,2),
       .        udel_h(2,2),DAXOP_new(2,2), Daxop_newer(2,2)
        COMMON / AXOCM / DAXOP, DCOMP, SURREF, RICHM, DCOMP_new, 
       .        udel_h, DAXOP_new, Daxop_newer
  C         VARIABLES 'TO':
  C            1. DCOMP(2,2) - THE PARTIAL DERIVATIVES OF THE COMPONENT OF THE
  C                            ANTENNA AXIS OFFSET IN THE APPARENT DIRECTION OF
  C                            THE SOURCE AND THE CT TIME DERIVATIVE OF THAT
  C                            COMPONENT WITH RESPECT TO THE ANTENNA AXIS
  C                            OFFSETS AT EACH OBSERVATION SITE. THE FIRST
  C                            INDEX RUNS OVER THE OBSERVATION SITES AND THE
  C                            SECOND INDEX RUNS OVER THE DELAY AND THE DELAY
  C                            RATE.  (UNITLESS, 1/SEC) 
  C            2. SURREF     - THE SURFACE REFRACTIVITY AT RADIO WAVELENGTHS
  C                            IN N UNITS (N = (1 - REFRACTIVITY) * 10**6)
  C                            (SEE REFERENCES) 
  C            3. RICHM(2)   - THE AXES ORIENTATION FOR THE RICHMOND ANTENNA
  C                            SYSTEM. IT IS ORIENTED AS FOR AN EQUATORIAL
  C                            MOUNT AT LATITUDE+39.06 DEGREES AND ROTATED 0.12
  C                            DEGREES WEST OF NORTH.
  C            4. udel_h(2,2)- Up component of the topocentric unit vector axis
  C                            offset and its derivative. First index runs
  C                            over the sites, second runs over height and 
  C                            its time derivative. (meters, meters/sec)
  C            5. DCOMP_new(2,2)-The first two parts of the partial derivative of
  C                            the new axis offset delay and rate with respect
  C                            to AXOFF. First index runs over the sites, second
  C                            runs over delay and rate. (sec/m, sec/sec/m)
  C
        INCLUDE 'cphys.i'               
  C           VARIABLES 'FROM':
  C              1. VLIGHT  -  The velocity of light in a vacuum (m/sec).
  C 
        INCLUDE 'ccon.i'               
  C       VARIABLES 'FROM':
  C         KAXOC - THE AXIS OFFSET MODULE FLOW CONTROL FLAG.
  C                  = 0 ==> Axis offset corrections applied to theoreticals.
  C                  = 1 ==> Axis offset corrections NOT applied to theoreticals.
  C         KAXOD - THE AXIS OFFSET MODULE DEBUG OUTPUT FLAG.
  C         KPANC - The feedhorn rotation module flow control flag.
  C         KPAND - The feedhorn rotation module debug flag (not used).
  C
        Real*8         PI,TWOPI,HALFPI,CONVD,CONVDS,CONVHS,SECDAY 
        COMMON /CMATH/ PI,TWOPI,HALFPI,CONVD,CONVDS,CONVHS,SECDAY 
        EXTERNAL CMATHB
  C           VARIABLES 'FROM':
  C              1. HALFPI - THE VALUE OF PI/2
  C              2. CONVD  - THE CONVERSION FACTOR FROM DEGREES TO RADIANS
  C                          (RAD/DEG)
  C
        INCLUDE 'cobsn.i'
  C          Variables from:
  C            1. Nzero  -  Set to 1 or 2 if station 1 or 2 is at the geocenter,
  C                         otherwise equals zero. For correlator usage.
  C       DATA BASE ACCESS -
  C            'PUT' VARIABLES:
  C               1. PANGL(2) - Feedbox rotation angles at site 1 and 2 in degrees
  C                             of phase.
  C             ACCESS CODES:
  C               1. 'PARANGLE' - The database access code for the feedbox 
  C                               rotation angles.
  C               2. 'FEED.COR' - The database access code for the feedbox
  C                               rotation correction to the group delay and 
  C                               the phase delay rate.
  C 
  C 3.2.3 PROGRAM SPECIFICATIONS -
        DATA  SURREF / 313.0D0/ 
        DATA RICHM /39.06D0, 0.12D0/  
        Integer*2 KAXIS(2), NDO(3), KERR
        Integer*4 L, LL, kx  
        Real*8 CFSTAR(3,2), CFTOTC(3,3), CTCSTR(3), SITEV(3,2), 
       1       R2000(3,3,3), SITLAT(2), STAR(3), TCAXIS(3), 
       2       TCSTAR(3,2), TCTOCF(3,3,2), TR2000(3,3,2),
       3       AZ(2,2), ELEV(2,2), STAR_ABERRATED(3,2), SITHEIGHT(2),
       4       ZC, ZEN, TANEL, AZMUTH, ZENCOR, REF_FREQ
        Real*8 axistc(3,2), daxistc(3,2), axiscf(3,2), daxiscf(3,2), 
       *       AXOFF(2), xnormal, ccfstr(3), dccfstr(3), vectr(3),
       *       daxocnew(2,2), EARTH(3,3), c2000str(3), axis2000(3,2),
       *       daxis2000(3,2), dctcstr(3), vec1(3), dvec1(3), absvec1,
       *       dabsvec1, vec2(3), vec3(3), uaxis2000(3,2),duaxis2000(3,2)
        Real*8 VECMG, DOTP, AMAG
        Real*8 zen_K, tanel_K, azmuth_K, zencor_K
        Real*8 term1, dterm1, sinEZ, dsinEZ, DotDot1, dDot1Dot1,
       *       term2, dterm2
        Real*8 El_Rad, Sithit, Temp_K, X, Press_Hg, Humid_F, sbend
        Real*8 TNCP(3), ZA, DZA, DAZM, RANGL, DRANGL, PANGL(2), 
       *       FEED_COR(2,2), FCONT(2)
  C 
  C 3.2.6 SUBROUTINE INTERFACE -
  C          CALLER SUBROUTINES: DRIVG 
  C          CALLED SUBROUTINES: DCOS, DSIN, DTAN, MTRAN, VECRT, VECMG, VECAD,
  C                              DOTP, CROSP, VECAD, FBOX
  C 
  C 3.2.8 PROGRAM VARIABLES - 
  C           1. TR2000(3,3,2)  -  THE COMPLETE J2000.0 TO CRUST FIXED
  C                                ROTATION MATRIX AND ITS FIRST CT
  C                                TIME DERIVATIVE. (UNITLESS, 1/SEC)
  C           2. CFTOTC(3,3)    -  THE 3x3 ROTATION MATRIX WHICH ROTATES 
  C                                THE GEOCENTRIC CRUST FIXED REFERENCE
  C                                SYSTEM TO THE TOPOCENTRIC REFERENCE 
  C                                SYSTEM. (UNITLESS) 
  C           3. CFSTAR(3,2)    -  THE SOURCE UNIT VECTOR IN THE CRUST FIXED
  C                                GEOCENTRIC REFERENCE SYSTEM CORRECTED FOR
  C                                ABERRATION AND ITS CT TIME DERIVATIVE.
  C                                (UNITLESS, 1/SEC)
  C           4. TCSTAR(3,2)    -  THE SOURCE UNIT VECTOR IN THE TOPOCENTRIC
  C                                SYSTEM CORRECTED FOR ABERRATION AND ITS CT 
  C                                TIME DERIVATIVE. (UNITLESS, 1/SEC)
  C           5. CTCSTR(3)      -  THE SOURCE UNIT VECTOR IN THE TOPOCENTRIC
  C                                REFERENCE SYSTEM CORRECTED FOR ABERRATION
  C                                AND ATMOSPHERIC REFRACTION. (UNITLESS)
  C           6. dctcstr(3)     -  Derivative of CTCSTR (above).
  C           7. ccfstr(3)      -  Aberrated/refracted source unit vector in the
  C                                crust fixed frame.
  C           8. dccfstr(3)     -  Time derivative of ccfstr (above).
  C           9. c2000str(3)    -  Aberrated/refracted source unit vector in the
  C                                J2000.0 frame.
  C          10. AZMUTH         -  THE AZIMUTH OF THE TOPOCENTRIC SOURCE UNIT
  C                                VECTOR CORRECTED FOR ABERRATION. (RAD) 
  C          11. TANEL          -  THE TANGENT OF THE ELEVATION ANGLE,
  C                                CORRECTED FOR ABERRATION. (UNITLESS)
  C          12. ZEN            -  THE ZENITH ANGLE OF THE SOURCE UNIT VECTOR.
  C                                (RADIANS)
  C          13. ZENCOR         -  THE CORRECTION ANGLE TO THE SOURCE UNIT VECTOR
  C                                DUE TO ATMOSPHERIC REFRACTION. (RADIANS)
  C          14. TCAXIS(3)      -  THE UNIT VECTOR REPRESENTING THE ANTENNA FIXED
  C                                AXIS IN THE TOPOCENTRIC SYSTEM. (UNITLESS)
  C          15. XI             -  THE ANGLE BETWEEN THE FIXED ANTENNA AXIS AND
  C                                THE SOURCE UNIT VECTOR CORRECTED FOR 
  C                                ABERRATION AND REFRACTION. (RADIANS)
  C          16. axistc(3,2)    -  Vector axis offset of the antenna in the
  C                                topocentric frame (effect on baseline). First
  C                                index is X,Y,Z position (meters); second runs
  C                                over the sites.
  C          17. daxistc(3,2)   -  Time derivative of axistc, rate of change of
  C                                vector axis offset of antenna in the
  C                                topocentric frame (effect on baseline). First
  C                                index is velocity, second runs over sites.
  C          18. axiscf(3,2)    -  Vector axis offset of the antenna in the
  C                                geocentric frame (effect on baseline). First
  C                                index is X,Y,Z position (meters); second runs
  C                                over the sites.
  C          19. daxiscf(3,2)   -  Time derivative of axiscf, rate of change of
  C                                vector axis offset of the antenna in the
  C                                geocentric frame (effect on baseline). First
  C                                index is velocity, second runs over sites.
  C          20. uaxis2000(3,2) -  Unit vector axis offset of antenna in the 
  C                                J2000.0 frame. First index is X,Y,Z, second
  C                                runs over sites. (unitless) [See axis2000]
  C          21. duaxis2000(3,2) - Time derivative of uaxis2000, rate of change
  C                                of unit vector axis offset in the J2000.0
  C                                frame. First index is velocity, second runs
  C                                over sites. (unitless) [See daxis2000]
  C          22. El_Rad          - Source elevation for use by function Sbend to
  C                                compute the atmospheric bending. (radians
  C          23. Sithit          - Station height above the geoid, for use in 
  C                                computing atmospheric bending. (meters) 
  C          24. Temp_K          - A priori value for absolute temperature, for
  C                                use in computing atmospheric bending. (Kelvins)
  C          25. Press_Hg        - A priori value for atmospheric pressure, for 
  C                                use in computing atmospheric bending. (mm of 
  C                                mercury)
  C          26. Humid_F         - A priori value for station humidity (.5).
  C          27. X               - Variable related to station height, for use in
  C                                computing atmospheric bending.
  C          28. PANGL(2)        - The clock-wise feed box rotation angles for
  C                                stations 1 and 2. (Degrees) 
  C          29. FEED_COR(2,2)   - The corrections to group delay and phase
  C                                delay rate due to feedbox rotation. First
  C                                index runs over the sites, second index runs
  C                                over delay and rate. Sec, Sec/sec)
  C
  C 3.2.9 PROGRAMMER - DALE MARKHAM  01/13/77
  C                    PETER DENATALE 07/11/77
  C                    DOUG ROBERTSON ??/??/84 (RICHMOND AXIS CODING)
  C                    DAVID GORDON   08/14/84 (RICHMOND)
  C                    Jim Ryan 89.06.29 Character stings used.
  C                    David Gordon  June - Aug, 1993 Modified to compute 3-D
  C                                  axis offset vector; Aberrated star vector,
  C                                  and altitude and azimuth from subroutine 
  C                                  ATMG added; new computation of DCOMP.
  C                    David Gordon  93.12.23 Corrected division by zero problem
  C                                  in DCOMP.
  C                    David Gordon  94.03.08 New axis offset delay computation.
  C                    David Gordon  94.03.25 Second (simpler) new axis offset 
  C                                  delay computation.
  C                    David Gordon  94.09.21 Added function Sbend to compute
  C                                  atmospheric bending. Sbend comes from the 
  C                                  Mark III Field System. Changed computation 
  C                                  of second new axis offset contribution to
  C                                  combine aberration and bending. SITHEIGHT(2)
  C                                  added to subroutine input arguments.
  C                    David Gordon  94.10.20 Changed computation of second axis
  C                                  offset correction to apply aberration to
  C                                  vacuum part, but not to aberrated/refracted 
  C                                  part.
  C                    David Gordon  94.10.24 Code cleanup.
  C                    David Gordon  98.08.05 Mods for geocenter station.
  C                    David Gordon  98.11.12 Code added to call the new feedbox
  C                                  rotation subroutine, FBOX, and to do PUT's 
  C                                  of the rotation angles and the group delay
  C                                  and phase delay rate corrections. Code to
  C                                  GET Reference frequency added. 
  C                    David Gordon  98.12.30 FEED.COR access code changed from
  C                                  a (2,2) array to a (2,1) array - sites 
  C                                  1 and 2 combined.
  C
  C 3.3   AXOG PROGRAM STRUCTURE
  C
  C   Get the reference frequency for use in the phase delay rate corrections. 
  C    Convert from MHz to Hz.
        CALL GET4('REF FREQ      ',REF_FREQ,1,1,1,NDO,KERR)
        IF(KERR.NE.0) then
          write(6,'("AXOG: Failure to obtain ref frequency.")')
          CALL CKILL(6HAXOG  ,1,KERR)
        Endif
        REF_FREQ = REF_FREQ*1.D6
  C
  C
  C    Loop through twice for the calculation of the axis offset vectors and
  C    their time derivatives and the partial derivatives of the component of
  C    the antenna axis offsets in the apparent direction of the source and the
  C    time derivative of that component with respect to the antenna axis
  C    offsets at each site.
  C
        DO 800  L = 1,2
             LL = L
  C
  C  Check for geocenter station
        IF (L. eq. Nzero) Go to 750
  C
  C    Identify the unit vector representing the antenna fixed axis in a
  C    topocentric reference system. The topocentric reference system sits on the
  C    Earth's surface at the observation site with axes pointing 1) radially Up
  C    (X-AXIS), 2) East (Y-axis), and 3) North (Z-axis).
  C
  C  Identify the antenna axis type in order to obtain the correct algorithm:
  C
  C    Equatorial mount - fixed axis points at the North Celestial Pole.
           If (KAXIS(L) .EQ. 1) Then
             TCAXIS(1) = DSIN ( SITLAT(L) )
             TCAXIS(2) = 0.D0
             TCAXIS(3) = DCOS ( SITLAT(L) )
             GO TO 310
           Endif
  C
  C    X-Y mount - antenna fixed axis points North-South. (i.e. the antenna fixed
  C      axis is in the plane of the horizon pointed North-South.)
           If (KAXIS(L) .EQ. 2) Then
             TCAXIS(1) = 0.D0
             TCAXIS(2) = 0.D0
             TCAXIS(3) = 1.D0
             GO TO 310
           Endif
  C
  C    Alt-Az mount - the antenna fixed axis points at the zenith.
           If (KAXIS(L) .EQ. 3) Then
             TCAXIS(1) = 1.D0
             TCAXIS(2) = 0.D0
             TCAXIS(3) = 0.D0
             GO TO 310
           Endif
  C
  C    X-Y mount -  antenna fixed axis pointing East-West. (i.e. The antenna
  C      fixed axis is in the plane of the horizon pointed East-West.)
           If (KAXIS(L) .EQ. 4) Then
             TCAXIS(1) = 0.D0
             TCAXIS(2) = 1.D0
             TCAXIS(3) = 0.D0
             GO TO 310
           Endif
  C
  C    Old Richmond antenna system - axes oriented as for an equatorial mount at
  C      latitude +39.06 degrees and rotated 0.12 degrees West of North.
           If (KAXIS(L) .EQ. 5) Then
             ZC        =     DCOS(RICHM(1)*CONVD)
             TCAXIS(1) =     DSIN(RICHM(1)*CONVD)
             TCAXIS(2) = -ZC*DSIN(RICHM(2)*CONVD)
             TCAXIS(3) =  ZC*DCOS(RICHM(2)*CONVD)
             GO TO 310
           Endif
  C
  C    Antenna axis type invalid, a message is written and the program stops.
           IF ((KAXIS(L) .le. 0 ) .OR. (KAXIS(L) .gt. 5)) Then
             WRITE ( 6, 9300 )  LL, KAXIS(LL)
             CALL CKILL (6HAXOG  , 0, 0)
           Endif
   9300 FORMAT (1X, "THE PROGRAM HAS TERMINATED IN SUBROUTINE AXOG.  ",
       1            'KAXIS (', I2, ') = ', I2, '.' )
  C
  C******************************************************************************
  C    Rotate the J2000.0 source unit vector to the topocentric reference system.
  C
  C        Compute the rotation matrix which rotates from the geocentric crust
  C        fixed reference system to the topocentric reference system.
    310      CALL MTRAN ( TCTOCF(1,1,L), CFTOTC )
  C
  C        Compute the rotation matrix which rotates from the J2000.0
  C        reference system to the geocentric crust fixed system.
             CALL MTRAN ( R2000(1,1,1), TR2000(1,1,1) )
  C
  C        Rotate the aberrated J2000.0 source unit vector to the crust fixed
  C        system and then to the topocentric system.
             CALL VECRT ( TR2000(1,1,1),STAR_ABERRATED(1,L),CFSTAR(1,1))
             CALL VECRT ( CFTOTC, CFSTAR(1,1), TCSTAR(1,1) )
  C
  C*****************************************************************************
  C    Correct the aberrated topocentric source unit vector for atmospheric
  C    refraction. 
  C
  C   Compute the zenith angle of the aberrated source.
          ZEN =  HALFPI - ELEV(L,1) 
  C   Compute the azimuth of the aberrated source.
          AZMUTH = AZ(L,1) 
  C   Compute values needed for function sbend.
          El_Rad = elev(l,1)
          Sithit = Sitheight(l)
          Temp_K = 293.15 - (6.5D-3)*sithit
          X = 1.D0 - (6.5D-3)*sithit/293.15D0
          Press_Hg = 760.D0 * (X**5.26D0)
          Humid_F = .5D0
  C   Compute atmospheric bending.
          zencor = sbend(El_rad,Temp_K,Humid_F,Press_Hg)
  C
  C *************************************************************************
  C    Compute the corrected (aberrated/refracted) topocentric star unit
  C    vector. (NOTE: This corrects the vector so that it is pointing nearer 
  C    to the zenith).
            CTCSTR(1) = DCOS ( ZEN - ZENCOR )
            CTCSTR(2) = DSIN ( ZEN - ZENCOR ) * DSIN ( AZMUTH )
            CTCSTR(3) = DSIN ( ZEN - ZENCOR ) * DCOS ( AZMUTH )
  C       Make sure we still have a unit vector
            amag = vecmg(CTCSTR)
            do kx = 1,3
             CTCSTR(kx) = CTCSTR(kx)/amag
            enddo
  C
  C     Compute the CT time derivative of the aberrated source unit vector
  C      in the topocentric reference system.
              CALL MTRAN ( R2000(1,1,2), TR2000(1,1,2) )
              CALL VECRT ( TR2000(1,1,2),STAR_ABERRATED(1,L),CFSTAR(1,2))
              CALL VECRT ( CFTOTC, CFSTAR(1,2), TCSTAR(1,2) )
  C
  C .............................................................................
  C   New code for computing the three dimensional components of the axis offset
  C    vector.  D. Gordon 93MAY13 thru 93JUL27
  C
  C  Compute topocentric axis offset vector (change in effective site position
  C   due to axis offset orientation) and its time derivative:
  C
  C   DELAY COMPUTATION:
  C   ------------------
  C    Take triple vector product to get direction of axis offset vector: 
          call crosp(TCAXIS,CTCSTR,vectr)
          call crosp(vectr,TCAXIS,vec1)
  C    Compute magnitude of vector so we can normalize it
         absvec1 = vecmg(vec1)
  C    Normalize to get topocentric unit axis offset vector
         do kx=1,3
          axistc(kx,L) = vec1(kx) / absvec1
         enddo
  C
  C    Rotate axis offset vector to the crust fixed frame
          call vecrt(TCTOCF(1,1,L),axistc(1,L),axiscf(1,L))
  C    Then rotate axis offset vector to the J2000 frame
          call vecrt(R2000(1,1,1),axiscf(1,L),uaxis2000(1,L))
  C
  C   Derivative of CTCSTR: (aberrated/refracted topocentric source unit vector)
  C   --------------------
  C     Rotate atmosphere corrected aberrated source unit vector to the
  C     crust fixed frame:
          call vecrt(TCTOCF(1,1,L),CTCSTR,ccfstr)
  C    Then rotate to the J2000 frame:
          call vecrt(R2000(1,1,1),ccfstr,c2000str)
  C
  C    Take time derivative and rotate back to crust fixed frame
          call vecrt(TR2000(1,1,2),c2000str,dccfstr)
  C
  C    Rotate source vector time derivative back to topocentric frame
          call vecrt(CFTOTC,dccfstr,dctcstr)
  C
  C   RATE COMPUTATION:
  C   -----------------
  C    Take derivative of triple vector product (TCAXIS is a constant)
          call crosp(TCAXIS,dctcstr,vectr)
          call crosp(vectr,TCAXIS,dvec1)
  C    Derivative of magnitude of vec1
          dabsvec1 = Dotp(vec1,dvec1) / absvec1
  C    Compute derivative of unit axis offset vector in topocentric frame
         do kx=1,3
          daxistc(kx,L) = dvec1(kx)/absvec1 -
       .   vec1(kx) * dabsvec1 / absvec1**2 
         enddo
  C
  C    Rotate axis offset time derivative to crust fixed frame
          call vecrt(TCTOCF(1,1,L),daxistc(1,L),daxiscf(1,L))
  C
  C    Then rotate to J2000 frame:
          call vecrt(R2000(1,1,1),daxiscf(1,L),vec2)
          call vecrt(R2000(1,1,2), axiscf(1,L),vec3)
          call vecad(vec2,vec3,duaxis2000(1,L))
  C
  C    Compute partial derivative of the components of the antenna axis offset
  C     in the apparent direction of the source with respect to the antenna axis
  C     offset.
          DCOMP(L,1) = Dotp(c2000str,uaxis2000(1,L))
          DCOMP(L,2) = Dotp(c2000str,duaxis2000(1,L))
  C
  C   Convert unit axis offset vector and its derivative to actual axis offset
  C     vector and its derivative.
         do kx=1,3
           axis2000(kx,L) =  uaxis2000(kx,L) * AXOFF(L)
          daxis2000(kx,L) = duaxis2000(kx,L) * AXOFF(L)
         enddo
  C
  C   For use later in AXOP, save the topocentric UP (change in height) component
  C    of the unit axis offset vector and its time derivative
         udel_h(L,1) =  axistc(1,L)    ! unit UP component
         udel_h(L,2) = daxistc(1,L)    ! unit UP time derivative
  C
  C .............................................................................
  C  Now compute new axis offset partial with respect to AXOFF (3-part formula,
  C    see D. Gordon reference)
  C
          sinEZ = DSIN(ELEV(L,1)+zencor)
         dsinEZ = DCOS(ELEV(L,1)+zencor) * ELEV(L,2)
  C
  C   Vacuum part 1:
  C!      term1 = -Dotp( uaxis2000(1,L),Star)
  C!     dterm1 = -Dotp(duaxis2000(1,L),Star)
  C   Vacuum part 1 with aberration correction:
          term1 = -Dotp( uaxis2000(1,L),Star_aberrated(1,L))
         dterm1 = -Dotp(duaxis2000(1,L),Star_aberrated(1,L))
  C
  C   Vacuum part 2 with aberration correction and refraction/geometric part:
  c       DotDot1 = Dotp(Star,c2000str) - 1.D0
  C      dDotDot1 = 0.D0  ! approximation
          DotDot1 = Dotp(Star_aberrated(1,L),c2000str) - 1.D0
          term2 =  axistc(L,1) * DotDot1 / sinEZ
         dterm2 = daxistc(L,1) * DotDot1 / sinEZ - 
       *           axistc(L,1) * DotDot1 * dsinEZ / sinEZ**2
  C
  C   Now combine everything to get the first two parts of the new axis 
  C    offset delay and rate partials with respect to AXOFF. (The last part, the
  C    additional atmosphere correction to the reference height is done in AXOP.)
         DCOMP_new(L,1) = (term1 + term2) / VLIGHT 
         DCOMP_new(L,2) = (dterm1 + dterm2) / VLIGHT
  C******************************************************************************
  C
  C          Now do feedbox rotation computations
  C
          IF(KPANC.EQ.1) THEN
            PANGL(L)      = 0.0D0      
            FEED_COR(L,1) = 0.0D0
            FEED_COR(L,2) = 0.0D0 
            GO TO 156
           ENDIF  
  C
  C  Define vector TNCP(3), a unit vector in the direction of the North 
  C   Celestial pole.
             TNCP(1) = DSIN(SITLAT(L))
             TNCP(2) = 0.D0
             TNCP(3) = DCOS(SITLAT(L))
  C
          ZA = ZEN - ZENCOR
          DZA = -ELEV(L,2)
          DAZM = AZ(l,2)
         If (KAXIS(L) .NE. 1) Then
          CALL FBOX (TNCP, TCAXIS, CTCSTR, AZMUTH, DAZM, ZA, DZA, 
       *            Rangl, dRangl)
         Else
  C           Equatorial Mount, no feedbox rotation
          Rangl = 0.0D0
          dRangl = 0.0D0
         Endif
  C
  C  Convert the feedbox rotation angle to degrees for this site. 
         PANGL(L) = RANGL/CONVD
  C  Compute the feedbox rotation corrections for group delay and phase delay
  C   rate for this site. The units are seconds and sec/sec.
         FEED_COR(L,1) = 0.D0
         FEED_COR(L,2) = DRANGL/TWOPI/REF_FREQ 
  C
  C    To correct the phases, either ADD the quantity [PANGL(2) - PANGL(1)]
  C     to the total or residual observed phases, or SUBTRACT it from the 
  C     theoretical (model) phases.
  C    Note that the group delay corrections are always zero.
  C    To correct the phase delay rates, either ADD the quantity 
  C     [FEED_COR(2,2) - FEEDCOR(1,2)] to the observed rates, or SUBTRACT
  C     it from the theoretical (model) rates.
  C
   156   CONTINUE
  C******************************************************************************
  C     See if debug is needed.
        IF ( KAXOD .ne. 0 )  Then
  C
        If (L.eq.1) then 
          WRITE ( 6, 9100 )
        endif
  C
        If (L.eq.2) then 
          WRITE ( 6, 9101 )
          WRITE ( 6, 9200 )  KAXIS, R2000, SITLAT, STAR, TCTOCF
          WRITE(6,8)' TR2000  ',TR2000
          WRITE(6,8)' STAR_ABERRATED ',STAR_ABERRATED
          WRITE(6,8)' AXOFF   ', AXOFF
          WRITE(6,8)' PANGL ', PANGL
          WRITE(6,8)' FEED_COR  ', FEED_COR
        endif
  C
   9100 FORMAT (20X, "Debug output for subroutine AXOG." )
   9200   FORMAT (1X, "KAXIS  = ", 2 ( I2, 10X ), /, 1X,
       1   "R2000  = ", 9 ( 3 ( D30.16, 10X ), /, 1X ), /, 1X,
       2   "SITLAT = ", 2 ( D30.16, 10X ), /, 1X,
       3   "STAR   = ", 3 ( D30.16, 4X ), /, 1X,
       4   "TCTOCF = ", 6 ( 3 ( D30.16, 10X ), /, 1X ) )
   9101 FORMAT (1X, "Dump for both stations:          " )
      8 FORMAT(A,3D25.16,/,5(7X,3D25.16))
  C
        write(6,'(/,10x," ******** DUMP FOR SITE",i2," ******** ")') L
        WRITE(6,8)' TCAXIS  ',TCAXIS
        WRITE(6,8)' CFTOTC  ',CFTOTC
        WRITE(6,8)' CFSTAR  ',CFSTAR
        WRITE(6,8)' TCSTAR  ',TCSTAR
        WRITE(6,8)' SURREF  ',SURREF
        WRITE(6,8)' ZEN     ',ZEN
        WRITE(6,8)' TANEL   ',TANEL
        WRITE(6,8)' AZMUTH  ',AZMUTH
        WRITE(6,8)' ZENCOR  ',ZENCOR
        write(6,'("vec1:  ",3D23.14)') vec1
        write(6,'("dvec1: ",3D23.14)') dvec1
        write(6,'("absvec1: ",3D23.14)') absvec1
        write(6,'("dabsvec1: ",3D23.14)') dabsvec1
        write(6,'("axistc:  ",3D23.14)') axistc(1,L),axistc(2,L)
       *                                ,axistc(3,L)
        write(6,'("daxistc: ",3D23.14)') daxistc(1,L), daxistc(2,L), 
       *                                 daxistc(3,L) 
        write(6,'("axiscf:  ",3D23.14)') axiscf(1,L),axiscf(2,L)
       *                                ,axiscf(3,L)
        write(6,'("daxiscf: ",3D23.14)') daxiscf(1,L), daxiscf(2,L)
       *                                ,daxiscf(3,L)
        write(6,'("c2000str:  ",3D23.14)') c2000str
        write(6,'("ccfstr:  ",3D23.14)') ccfstr
        write(6,'("dccfstr:  ",3D23.14)') dccfstr
        write(6,'("CTCSTR:  ",3D23.14)') CTCSTR
        write(6,'("dctcstr:  ",3D23.14)') dctcstr
        write(6,'("uaxis2000:",3D23.14)') uaxis2000(1,L),uaxis2000(2,L)
       *                                ,uaxis2000(3,L) 
        write(6,'("duaxis2000:",3D23.14)') duaxis2000(1,L),duaxis2000(2,L)
       *                                 ,duaxis2000(3,L)
        write(6,'("axis2000:",3D23.14)') axis2000(1,L),axis2000(2,L)
       *                                ,axis2000(3,L) 
        write(6,'("daxis2000:",3D23.14)') daxis2000(1,L),daxis2000(2,L)
       *                                 ,daxis2000(3,L)
        write(6,'("DCOMP:   ",2D23.14)') DCOMP(L,1),DCOMP(L,2)
        write(6,'("DCOMP_new",2D23.14)') DCOMP_new(L,1),DCOMP_new(L,2)
        WRITE(6,8)' SinEZ   ', SinEZ
        WRITE(6,8)' term1, dterm1 ', term1, dterm1 
        WRITE(6,8)' DotDot1 ', DotDot1
        WRITE(6,8)' term2, dterm2 ', term2, dterm2 
  C
        ENDIF
         Go to 800
  C
    750 CONTINUE
  C   Handling of geocenter station, set axis offset quantities to zero.
         do kx=1,3
           axis2000(kx,L) = 0.0D0
          daxis2000(kx,L) = 0.0D0
         enddo
          DCOMP(L,1) = 0.0D0 
          DCOMP(L,2) = 0.0D0 
          DCOMP_new(L,1) = 0.0D0
          DCOMP_new(L,2) = 0.0D0 
          udel_h(L,1) = 0.0D0 
          udel_h(L,2) = 0.0D0 
  C    Ditto for the feedbox rotation quantities.
         PANGL(L)      = 0.0D0      
         FEED_COR(L,1) = 0.0D0
         FEED_COR(L,2) = 0.0D0 
  C
  C     Close loop running over the observing sites.
    800 CONTINUE
  C
  C  PUT the feedbox rotation angles and delay/rate corrections in the database.
  C       FCONT(1) = FEED_COR(1,1) - FEEDCOR(2,1)
          FCONT(1) = 0.D0 
          FCONT(2) = FEED_COR(1,2) - FEED_COR(2,2)
        CALL PUT4 ('PARANGLE      ', PANGL, 2, 1, 1 )
  C     CALL PUT4 ('FEED.COR      ', FEED_COR,2,2,1)
        CALL PUT4 ('FEED.COR      ', FCONT,2,1,1)
  C
        IF ( KAXOD .ne. 0 )  Then
         WRITE(6,'(" KAXIS, PANGL ",2I3,2F10.2)')  KAXIS, PANGL
         WRITE(6,'(" FEED_COR ",4D20.10)')  FEED_COR
         WRITE(6,'(" FCONT    ",2D20.10)')  FEED_COR
        ENDIF
  C
  C     Normal conclusion.
        RETURN
        END
  C
  C******************************************************************************
        SUBROUTINE AXOP ( AXOFF, dATMCdh )
        IMPLICIT None 
  C
  C 5.    AXOP
  C
  C 5.1   AXOP PROGRAM SPECIFICATION
  C
  C 5.1.1 AXOP is the partial derivatives section of the axis offset
  C       module. AXOP computes the partial derivatives of the delay
  C       and rate with respect to the antenna axis offsets at the sites.
  C
  C 5.2   AXOP PROGRAM INTERFACE
  C
  C 5.2.1 CALLING SEQUENCE -
  C       INPUT VARIABLES:
  C          1. AXOFF(2)     - THE ANTENNA AXIS OFFSETS AT EACH SITE. (M)
  C          2. dATMCdh(2,2) - Partial derivatives of the Niell dry 
  C                            atmosphere contribution with respect to
  C                            station height. First index runs over the
  C                            sites and the second over delay and rate.
  C                            (sec/meter, sec/sec/meter)
  C
        Real*8  AXOFF(2), dATMCdh(2,2)
  C
  C 5.2.2 COMMON BLOCKS USED -
  C
        INCLUDE 'cphys.i'
  C         VARIABLES 'FROM':
  C            1.  VLIGHT - THE VELOCITY OF LIGHT IN VACUUM. (M/SEC)
  C
        Real*8  DAXOP(2,2), DCOMP(2,2), SURREF, RICHM(2), DCOMP_new(2,2),
       .        udel_h(2,2),DAXOP_new(2,2), Daxop_newer(2,2)
        COMMON / AXOCM / DAXOP, DCOMP, SURREF, RICHM, DCOMP_new, 
       .        udel_h, DAXOP_new, Daxop_newer
  C         VARIABLES 'FROM': 
  C            1. DCOMP(2,2)    -  THE PARTIAL DERIVATIVES OF THE COMPONENT OF
  C                                THE ANTENNA AXIS OFFSET IN THE APPARENT
  C                                DIRECTION OF THE SOURCE AND THE CT TIME
  C                                DERIVATIVE OF THAT COMPONENT WITH RESPECT TO
  C                                THE ANTENNA AXIS OFFSETS AT EACH SITE. THE
  C                                FIRST INDEX RUNS OVER THE OBSERVATION SITES
  C                                AND THE SECOND INDEX RUNS OVER THE DELAY AND
  C                                THE DELAY RATE. (UNITLESS, 1/SEC) 
  C            2. DCOMP_new(2,2) - The first two parts of the partial derivative 
  C                                of the new axis offset delay and rate with 
  C                                respect to AXOFF. First index runs over the
  C                                sites, second runs over delay and rate.
  C                                (sec/m, sec/sec/m)
  C            3. udel_h(2,2)    - Up component of the topocentric unit vector 
  C                                axis offset and its derivative. First index 
  C                                runs over the sites, second runs over height 
  C                                and its time derivative. (meters, meters/sec)
  C         VARIABLES 'TO': 
  C            1. DAXOP(2,2)     - THE PARTIAL DERIVATIVES OF THE DELAY AND THE
  C                                DELAY RATE WITH RESPECT TO THE ANTENNA AXIS 
  C                                OFFSETS AT EACH SITE (old model). FIRST INDEX
  C                                RUNS OVER THE OBSERVATION SITES AND THE
  C                                SECOND RUNS OVER THE DELAY AND DELAY RATE.
  C                                (SEC/M, SEC/SEC-M) 
  C            2. DAXOP_new(2,2) - The partial derivatives of the delay and delay
  C                                rate with respect to the antenna axis offsets
  C                                at each site using the new 3-term model. First
  C                                index runs over sites, second over delay and 
  C                                rate. (sec/m, sec/sec/m)
  C            3. DAXOP_newer(2,2)-The partial derivatives of the delay and delay
  C                                rate with respect to the antenna axis offsets
  C                                at each site using the new 1-term model. First
  C                                index runs over sites, second over delay and 
  C                                rate. (sec/m, sec/sec/m)
  C
        INCLUDE 'ccon.i'             
  C         VARIABLES 'FROM': 
  C            1. KAXOC  -  THE AXIS OFFSET MODULE FLOW CONTROL FLAG. 
  C            2. KAXOD  -  THE AXIS OFFSET MODULE DEBUG OUTPUT FLAG. 
  C
        Common /OBS/N_air, DAXOC_new, DAXOC_newer
        Real*8 N_air(2), DAXOC_new(2,2), DAXOC_newer(2,2)
  C 
  C 5.2.3 PROGRAM SPECIFICATIONS -
  C 
  C 5.2.4 DATA BASE ACCESS -
  C       'PUT' VARIABLES:
  C          1. Daxop_newer(2,2)-THE PARTIAL DERIVATIVES OF THE DELAY AND THE
  C                              DELAY RATE WITH RESPECT TO THE ANTENNA AXIS
  C                              OFFSETS AT EACH SITE (modified for new model).
  C                              THE FIRST INDEX RUNS OVER THE OBSERVATION SITES
  C                              AND THE SECOND RUNS OVER THE DELAY AND DELAY
  C                              RATE. (SEC/M, SEC/SEC-M) 
  C       ACCESS CODES: 
  C          1.  'AXO PART'  -  THE DATA BASE ACCESS CODE FOR THE AXIS OFFSET
  C                             MODULE PARTIAL DERIVATIVES ARRAY. 
  C 
  C 5.2.6 SUBROUTINE INTERFACE -
  C          CALLER SUBROUTINES: DRIVP
  C          CALLED SUBROUTINES: PUT4
  C
  C 5.2.7 CONSTANTS USED - VLIGHT
  C
  C 5.2.8 PROGRAM VARIABLES 
        Real*8 uatm1, duatm1, uatm2, duatm2 
        Integer*4 K
  C
  C 5.2.9 PROGRAMMER - DALE MARKHAM  01/13/77
  C                    PETER DENATALE 07/11/77
  C                    BRUCE SCHUPLER 03/07/78
  C                    Jim Ryan 89.06.29 Character stings used.
  C                    Jim Ryan 89.10.05 CPHYS common now an include file.
  C                    Jim Ryan 89.12.12 UNIX-like database interface
  C                    implimented.
  C                    D. Gordon 94.03.25 Old partials multiplied by index of
  C                              refraction and PUT into database (new, simple 
  C                              axis offset model). Added complete partials for
  C                              more complex (3 terms) axis offset model. 
  C                    D. Gordon 94.08.03 Corrected old axis multiplication by
  C                              N_air(site #). Was incorrrectly multiplying both
  C                              delays by site 1 and both rates by site 2 
  C                              indices of refraction.
  C
  C 5.3   AXOP PROGRAM STRUCTURE
  C
  C     Compute the partial derivatives of the delay and rate with
  C     respect to the antenna axis offsets at each site. Old compuation now
  C     modified by index of refraction in air to correct for actual 
  C     propagation speed. 94MAR25 -DG- 
  C     (NOTE: The minus sign appears due to the definition of the
  C     baseline vector being the site#1 vector minus the site#2 vector.)
  C
  C   Old axis offset partial compution
        DO   K = 1,2                       ! looping over delay and rate
             DAXOP(1,K) = + DCOMP(1,K) / VLIGHT      ! site 1
             DAXOP(2,K) = - DCOMP(2,K) / VLIGHT      ! site 2
        ENDDO
  C
  C   New axis offset partial (modified to use velocity of light in air)
        DO   K = 1,2                       ! looping over delay and rate
             Daxop_newer(1,K) = DAXOP(1,K) * N_air(1)      ! site 1
             Daxop_newer(2,K) = DAXOP(2,K) * N_air(2)      ! site 2
        ENDDO
  C
  C     PUT the new (simple formula) axis offset partial derivatives array.
        CALL PUT4 ('AXO PART      ', Daxop_newer, 2, 2, 1 )
  C
  C Compute axis offset delay and rate new partials using the alternate new
  C  model (D. Gordon 3-term formula)
  C   Station 1:
         uatm1  =  dATMCdh(1,1)*udel_h(1,1)
        duatm1  =  dATMCdh(1,1)*udel_h(1,2) + dATMCdh(1,2)*udel_h(1,1)
        DAXOP_new(1,1) = - DCOMP_new(1,1) + uatm1
        DAXOP_new(1,2) = - DCOMP_new(1,2) + duatm1
  C
  C   Station 2:
         uatm2 =  dATMCdh(2,1)*udel_h(2,1) 
        duatm2 =  dATMCdh(2,1)*udel_h(2,2) + dATMCdh(2,2)*udel_h(2,1)
        DAXOP_new(2,1) = + DCOMP_new(2,1) + uatm2
        DAXOP_new(2,2) = + DCOMP_new(2,2) + duatm2
  C  The partials are so similar we only use the first one (Daxop_newer)
  C
  C   See if debug is neccessary.
        IF (KAXOD .ne. 0) Then
         WRITE ( 6, 9100 )
   9100  FORMAT (1X, "Debug output for subroutine AXOP." )
     8   FORMAT(A,5D25.16/(9X,5D25.16))
         WRITE(6,8)' AXOFF  ',AXOFF 
         WRITE(6,8)' DCOMP      ',DCOMP
         WRITE(6,8)' DCOMP_new  ',DCOMP_new
         WRITE(6,8)' DAXOP      ',DAXOP
         WRITE(6,8)' Daxop_newer',Daxop_newer
         WRITE(6,8)' Daxop_new  ',Daxop_new
         WRITE(6,8)' VLIGHT ',VLIGHT
         WRITE(6,'(" N_air: ",2f12.8)')  N_air 
         write(6,8)' dATMCdh ', dATMCdh 
         write(6,8)' udel_h  ', udel_h
         write(6,8)' uatm1, duatm1 ', uatm1,duatm1
         write(6,8)' uatm2, duatm2 ', uatm2,duatm2
        Endif
  C
        RETURN
        END
  C
  C******************************************************************************
        SUBROUTINE AXOC ( AXOFF, DAXOC )
        IMPLICIT None
  C
  C 6.    AXOC
  C
  C 6.1   AXOC PROGRAM SPECIFICATION
  C
  C 6.1.1 AXOC is the contributions section of the axis offset module. It
  C       computes the contributions to the delay and rate due to the antenna
  C       axis offsets at each site. We now compute three sets of contributions:
  C       1) The old model contributions, 2) Contributions from a new 3-term 
  C       model, and 3) Contributions from a new but simpler model which
  C       multiplies the old contributions by the indexes of refraction at each 
  C       site. Models 2 and 3 ahould give the same answer but they don't quite
  C       because of different assumptions in estimating a priori the indexes of 
  C       refraction versus estimating the a priori Niell mapping function. The
  C       contribution that will be added to the theoretical delays and rates
  C       is from the simple new formula. Corrections are computed to convert
  C       to the more complicated new formula or to the old formula in program
  C       SOLVE 'on the fly'.
  C
  C 6.2   AXOC PROGRAM INTERFACE
  C
  C 6.2.1 CALLING SEQUENCE -
  C          INPUT VARIABLES:
  C             1. AXOFF(2)  -  THE ANTENNA AXIS OFFSETS AT EACH SITE. (M) 
  C          OUTPUT VARIABLES: 
  C             1. DAXOC(2,2) - THE CONTRIBUTIONS TO THE DELAY AND TO THE DELAY
  C                             RATE DUE TO ANTENNA AXIS OFFSETS AT EACH SITE.
  C                             -Using the new 1-part axis offset delay formula.
  C                             THE FIRST INDEX RUNS OVER THE OBSERVATION SITES
  C                             AND THE SECOND RUNS OVER THE DELAY AND THE DELAY
  C                             RATE. (SEC, SEC/SEC) 
  C 
  C 6.2.2 COMMON BLOCKS USED -
  C 
        Real*8  DAXOP(2,2), DCOMP(2,2), SURREF, RICHM(2), DCOMP_new(2,2),
       .        udel_h(2,2),DAXOP_new(2,2), Daxop_newer(2,2)
        COMMON / AXOCM / DAXOP, DCOMP, SURREF, RICHM, DCOMP_new, 
       .        udel_h, DAXOP_new, Daxop_newer
  C          VARIABLES 'FROM':
  C            1. DAXOP(2,2)    -  THE PARTIAL DERIVATIVES OF THE DELAY AND THE
  C                                DELAY RATE WITH RESPECT TO THE ANTENNA AXIS
  C                                OFFSETS AT EACH SITE. THE FIRST INDEX RUNS
  C                                OVER THE OBSERVATION SITES AND THE SECOND RUNS
  C                                OVER THE DELAY AND THE DELAY RATE. (Old model)
  C                                (SEC/M, SEC/SEC-M) 
  C            2. DAXOP_new(2,2) - The partial derivatives of the delay and delay
  C                                rate with respect to the antenna axis offsets
  C                                at each site using the new 3-term model. First
  C                                index runs over sites, second over delay and 
  C                                rate. (sec/m, sec/sec/m)
  C            3. DAXOP_newer(2,2)-The partial derivatives of the delay and delay
  C                                rate with respect to the antenna axis offsets
  C                                at each site using the new 1-term model. First
  C                                index runs over sites, second over delay and 
  C                                rate. (sec/m, sec/sec/m)
  C 
         INCLUDE 'ccon.i'
  C         VARIABLES 'FROM':
  C            1. KAXOC - THE AXIS OFFSET MODULE FLOW CONTROL FLAG.
  C            2. KAXOD - THE AXIS OFFSET MODULE DEBUG OUTPUT FLAG.
  C 
  C     Common /OBS/N_air, DAXOC_new, DAXOC_newer
  C     Real*8 N_air(2), DAXOC_new(2,2), DAXOC_newer(2,2)
  C
  C 6.2.3 PROGRAM SPECIFICATIONS -
        Real*8  AXOFF(2), DAXOC(2,2)
        Real*8  del_axis_new(2), del_axis_old(2)
        Real*8  DAXOC_new(2,2), DAXOC_newer(2,2), DAXOC_old(2,2)
        Integer*4 K, L
  C 
  C 6.2.4 DATA BASE ACCESS -
  C         'PUT' VARIABLES: 
  C            1. DAXOC_newer(2,2)-THE CONTRIBUTIONS TO THE DELAY AND THE DELAY
  C                              RATE DUE TO ANTENNA AXIS OFFSETS AT EACH SITE,
  C                              using the new, simple axis offset model.
  C                              THE FIRST INDEX RUNS OVER THE OBSERVATION
  C                              SITES AND THE SECOND RUNS OVER THE DELAY AND
  C                              THE DELAY RATE. (SEC, SEC/SEC)
  C         ACCESS CODES:
  C            1. 'AXO CONT' -  The database access code for the new, simple
  C                             version, axis offset contributions array.
  C            2. 'AXO2CONT' -  The database access code for the correction to
  C                             use 'on-the-fly' in SOLVE to convert to the new,
  C                             more complex version of the axis offset model.
  C            3. 'AXIS OLD' -  The database access code for the correction to
  C                             use 'on-the-fly' in SOLVE to convert to the old 
  C                             (pre-Calc8.0) axis offset model.
  C 
  C 6.2.6 SUBROUTINE INTERFACE -
  C          CALLER SUBROUTINES: DRIVC
  C          CALLED SUBROUTINES: PUT4
  C
  C 6.2.7 CONSTANTS USED - NONE
  C
  C 6.2.8 PROGRAM VARIABLES -
  C            1. DAXOC_new(2,2)   - The contributions to the delay and rate due
  C                                  to the antenna axis offset using the new,
  C                                  complex model. First index runs over sites,
  C                                  second runs over delay and rate.
  C            2. DAXOC_newer(2,2) - The contributions to the delay and rate due
  C                                  to the antenna axis offset using the new,
  C                                  simple model. First index runs over sites,
  C                                  second runs over delay and rate.
  C            3. DAXOC_old(2,2)   - The contributions to the delay and rate due
  C                                  to the antenna axis offset using the old axis
  C                                  offset model. First index runs over sites,
  C                                  second runs over delay and rate.
  C            4. del_axis_new(2)  - Delay and rate differences between the new
  C                                  simple and the new complex axis offset
  C                                  models. Add to theoreticals in SOLVE. 
  C                                  (sec, sec/sec)
  C            5. del_axis_old(2)  - Delay and rate differences between the new
  C                                  simple and the old axis offset models. Add to
  C                                  theoreticals in SOLVE. (sec, sec/sec)
  C
  C 6.2.9 PROGRAMMER - DALE MARKHAM  01/13/77
  C                    PETER DENATALE 07/11/77
  C                    BRUCE SCHUPLER 03/07/78
  C                    Jim Ryan 89.06.29 Character stings used.
  C                    Jim Ryan 89.12.12 UNIX-like database interface
  C                    implimented.
  C                    D. Gordon 94.03.25 Added computations for both new versions
  C                              of the axis offset contributions. Made the
  C                              simpler version the default (will be added to 
  C                              theoreticals in THERY). Computing corrections
  C                              to convert to alternate new model and to old
  C                              model ('on-the-fly' in SOLVE), and PUTting them
  C                              into the database. 
  C                    D. Gordon 94.05.23 Changed $Include to Include.
  C                    D. Gordon 94.08.04 Error corrected. Old axis offset was
  C                              still the default. Corrected to make new/simple
  C                              offset the default. Added variable DAXOC_old.
  C                    D. Gordon 98.11.12 Removed PUT of old axis offset 
  C                              differential contribution.
  C
  C 6.3   AXOC PROGRAM STRUCTURE
  C
  C     Compute the three (old, new complex, and new simple) contributions.
        DO 120  K = 1,2
             DO 110  L = 1,2
  C            DAXOC_old(L,K)   = DAXOP(L,K)       * AXOFF(L)
               DAXOC_new(L,K)   = DAXOP_new(L,K)   * AXOFF(L)
               DAXOC_newer(L,K) = DAXOP_newer(L,K) * AXOFF(L)
    110      CONTINUE
    120 CONTINUE
  C
  C   Set default axis offset to the new simple formula
         do  k = 1,2
             do  l = 1,2
                  DAXOC(l,k) = DAXOC_newer(l,k)
             enddo   
         enddo   
  C
  C     Check to determine if the axis offset module is to be turned off.
        IF (KAXOC .eq. 1) Then
         DO 320  K = 1,2
             DO 310  L = 1,2
                  DAXOC(L,K) = 0.D0
    310      CONTINUE
    320  CONTINUE
        Endif
  C
  C     PUT  the new/simple axis offset contributions.
    400 CALL PUT4 ('AXO CONT      ', DAXOC_newer, 2, 2, 1 )
  C
  C  Compute difference corrections to get the equivalent of the other new
  C   formula and the old formula.
         do 340 k=1,2
              del_axis_new(k) = (DAXOC_new(1,k) + DAXOC_new(2,k)) -
       .                  (DAXOC_newer(1,k) + DAXOC_newer(2,k))
  C           del_axis_old(k) = (DAXOC_old(1,k) + DAXOC_old(2,k)) -
  C    .                  (DAXOC_newer(1,k) + DAXOC_newer(2,k))
    340  continue
  C
        CALL PUT4 ('AXO2CONT      ', del_axis_new, 2, 1, 1 )
  C     CALL PUT4 ('AXIS OLD      ', del_axis_old, 2, 1, 1 )
  C
  C     Check for debug output.
        IF (KAXOD .ne. 0)  Then
         WRITE ( 6, 9100 )
   9100  FORMAT (1X, "Debug output for subroutine AXOC." )
      8  FORMAT(A,5D25.16/(9X,5D25.16))
         WRITE(6,8)' AXOFF ',AXOFF
         WRITE(6,8)' DAXOP       ',DAXOP
         WRITE(6,8)' DAXOP_new   ',DAXOP_new
         WRITE(6,8)' DAXOP_newer ',DAXOP_newer
         WRITE(6,8)' DAXOC       ',DAXOC
  C      WRITE(6,8)' DAXOC_old   ',DAXOC_old
         WRITE(6,8)' DAXOC_new   ',DAXOC_new
         WRITE(6,8)' DAXOC_newer ',DAXOC_newer
         WRITE(6,8)'del_axis_new  ',del_axis_new
  C      WRITE(6,8)'del_axis_old',del_axis_old
        Endif
  C
    600 RETURN
        END
  C
  C******************************************************************************
        double precision function sbend(El_rad,Temp_K,Humid_F,Press_Hg)
        IMPLICIT None
  C
  C input:
  C   El_rad   -- elevation angle in radians 
  C   Press_Hg -- Pressure in mm of Mercury (Hg)
  C   Temp_K   -- Temperature in Kelvins
  C   Humid_F  -- relative humidity (percent)
  C
  C output   --
  C   Sbend  -- bending angle in radians.
  C
        Real*8 El_rad, Temp_K, Humid_F, Press_Hg  
        Real*8 e(12),wp1(4),d3
        Real*8 fp,ft,fw,u,x,ad1,ad2,bd1,bd2,zd2,r,delta
        Real*8 a1,a2,b1,b2,c1,c2,e1,e2,e3,e4,e5,e6,e7,e8,e9
        Real*8 e10,e11,e12,p1,p2,t1,t2,z1,z2,w0,w1,w2,w3
        Real*8 conv
        Integer*4 I
  C
        equivalence (e( 1), e1),(e( 2), e2),(e( 3), e3),(e( 4), e4),
       &            (e( 5), e5),(e( 6), e6),(e( 7), e7),(e( 8), e8),
       &            (e( 9), e9),(e(10),e10),(e(11),e11),(e(12),e12)
        equivalence (wp1(1),w0),(wp1(2),w1),(wp1(3),w2),(wp1(4),w3)
  C 
        data a1, a2 /     0.40816d0, 112.30d0  / 
        data b1, b2 /     0.12820d0, 142.88d0  / 
        data c1, c2 /     0.80000d0,  99.344d0 / 
        data e   /    46.625d0  ,  45.375d0 ,     4.1572d0,  1.4468d0  ,
       &               0.25391d0,   2.2716d0,    -1.3465d0, -4.3877d0  ,
       &               3.1484d0 ,   4.5201d0,    -1.8982d0,  0.89000d0 /
        data p1 /   760.0d0 / 
        data t1 /   273.0d0 / 
        data wp1 / 22000.0d0    ,  17.149d0 ,  4684.1d0,    38.450d0   /
        data z1 /  91.870d0 / 
        data conv/57.295779512d0/ 
  C 
        Real*8         PI,TWOPI,HALFPI,CONVD,CONVDS,CONVHS,SECDAY 
        COMMON /CMATH/ PI,TWOPI,HALFPI,CONVD,CONVDS,CONVHS,SECDAY 
        EXTERNAL CMATHB
  C           VARIABLES 'FROM':
  C              1. HALFPI - THE VALUE OF PI/2
  C              2. CONVD  - THE CONVERSION FACTOR FROM DEGREES TO RADIANS
  C                          (RAD/DEG)
  C 
  C STATEMENT FUCNTION
        delta(ad1,ad2,bd1,bd2,zd2)=(ad2-ad1)*dexp(bd1*(zd2-bd2))
  C 
  C CONVERT UNITS 
  C  Zenith angle in degrees
        z2 = 90.0d0 - El_Rad/CONVD
  C  Temperature in Kelvins
        t2 = Temp_K
  C  Fractional humidity (0.0 -> 1.0)
        r = Humid_F 
  C  Pressure in mm of Hg
        p2 = Press_Hg
  C
  C      WRITE(6,9956) R,T2,P2,Z2                                          DEBUG
  C9956  FORMAT(" R,T2,P2,Z2",4F10.4)                                      DEBUG
  C 
  C CALCULATE CORRECTIONS FOR PRES, TEMP, AND WETNESS 
  C      WRITE(6,9980) Z1,Z2,C1,C2                                         DEBUG
  C9980  FORMAT(" Z1,Z2,C1,C2",4D15.6)                                     DEBUG
        d3=1.0d0+delta(z1,z2,c1,c2,z2)
        fp=(p2/p1)*(1.0d0-delta(p1,p2,a1,a2,z2)/d3) 
        ft=(t1/t2)*(1.0d0-delta(t1,t2,b1,b2,z2)/d3) 
        fw=1.0d0+(w0*r*dexp((w1*t2-w2)/(t2-w3))/(t2*p2))
  C      WRITE(6,9957) D3,FP,FT,FW                                         DEBUG
  C9957  FORMAT(" D3,FP,FT,FW ",4D12.6)                                    DEBUG
  C 
  C  CALCULATE OPTICAL REFRACTION 
        u=(z2-e1)/e2
        x=e11 
        do 10 i=1,8 
          x=e(11-i)+u*x 
  10    continue
  C 
  C  COMBINE FACTORS AND FINISH OPTICAL FACTOR
        sbend=ft*fp*fw*(dexp(x/d3)-e12) 
  C      WRITE(6,9958) SBEND,X,U                                           DEBUG
  C9958  FORMAT(" SBEND,X,U ",3D20.10)                                     DEBUG
  C 
  C BACK TO RADIANS FROM ARC SECONDS
        sbend=(sbend/3600.0d0)*CONVD
  C      WRITE(6,9959) SBEND                                               DEBUG
   9959  FORMAT(" SBEND (RADIANS) ",D20.10)                                DEBUG
        return
        end 
  C
  C-----------------------------------------------------------------------------
        SUBROUTINE FBOX ( TNCP, TCAXIS, CTCSTR, AZ, DAZ, ZA, DZA, 
       *            RANGL, DRANGL )
  C
  C     FBOX computes the feedbox rotation angle, relative to the source. 
  C     It computes the angle, looking up at the sky and measured at the 
  C     source, between the North Celestial Pole and the antenna fixed axis. 
  C     An increase in the feedbox rotation angle means a clockwise 
  C     rotation of the feedbox, for an observer looking up at the sky. 
  C     The angle itself is arbitrary. The important quantity is the 
  C     difference in the feedbox rotations at the two antennas, and how it
  C     changes from one observation to the next.
  C     
  C     It's not clear what to do with these angle in all case. For RCP data
  C     in geodesy, we correct the phase by adding the remote station 
  C     correction and subtracting the reference station correction. 
  C
  C     The algorithm makes two coordinate rotations to put the star vector
  C     along the X-axis. The azimuths of the NCP and the fixed axis in the
  C     transformed Y-Z plane are then easily found and differenced.  
  C
  C     This subroutine replaces subroutine PANG in Calc 9.0, and will handle
  C     all azis types. 
  C
  C     Input Variables:
  C         1) CTCSTR(3)  - Topocentric unit vector in the direction of the 
  C                         source. [X = Up = COS(Zenith Angle). Y = East =
  C                         SIN(Zenith Angle) * SIN(Azimuth). Z = North =
  C                         SIN(Zenith Angle) * COS(Azimuth). For greatest
  C                         accuracy, the effects of atmospheric refraction 
  C                         and aberration should be included.]
  C         2) TCAXIS(3)  - Topocentric unit vector in the direction of the
  C                         antenna's fixed axis. [Altazimuth: X = 1.0,
  C                         Y = 0.0, Z = 0.0.  X/Y North (Gilcreek): 
  C                         X = 0., Y = 0., Z = 1.0.  X/Y East (HARTRAO):
  C                         X = 0., Y = 1.0., Z = 0.] 
  C         3) TNCP(3)    - Topocentric unit vector in the direction of the 
  C                         North Celestial Pole. [X = SIN(Geodetic Latitude),
  C                         Y = 0.0, Z = COS(Geodetic Latitude)] 
  C         4) AZ         - Azimuth angle of source, in radians. 
  C         5) DAZ        - Time derivative of Azimuth angle. (radians/sec)
  C                         Zeros can be used here if the user does not care
  C                         about the rate of change of the rotation angle, or
  C                         about the correction to the phase delay rate.
  C         6) ZA         - Zenith angle of source, in radians. Include
  C                         refraction and aberration for greatest accuracy.
  C         7) DZA        - Time derivative of zenith angle. (radians/sec)
  C                         Zeros can be used here if the user does not care
  C                         about the rate of change of the rotation angle, or
  C                         about the correction to the phase delay rate.
  C     Output Variables:
  C         1) RANGL      - The clockwise rotation angle of the feedbox. (rad)
  C         2) DRANGL     - The time derivative of RANGL. (rad/sec)
  C
  C 
        Real*8         PI,TWOPI,HALFPI,CONVD,CONVDS,CONVHS,SECDAY 
        COMMON /CMATH/ PI,TWOPI,HALFPI,CONVD,CONVDS,CONVHS,SECDAY 
        EXTERNAL CMATHB
  C           VARIABLES 'FROM':
  C              1. HALFPI - THE VALUE OF PI/2
  C              2. CONVD  - THE CONVERSION FACTOR FROM DEGREES TO RADIANS
  C                          (RAD/DEG)
  C
  C    Program History:
  C       D. Gordon 98.11.10 Subroutine written at GSFC. Replaces subroutine 
  C                 PANG. 
  C
  C
        Real*8 TNCP(3), TCAXIS(3), CTCSTR(3), AZ, DAZ, 
       *       ZA, DZA, Rangl, dRangl
        REAL*8 P1AZ(3,3), P2ZA(3,3), P12(3,3), DP1AZ(3,3), DP2ZA(3,3), 
       *       DP12(3,3), TT1(3,3), TT2(3,3)
        REAL*8 NCP1(3), AXIS1(3), STAR1(3), DNCP1(3), DAXIS1(3)
        REAL*8 AZNCP, AZAXIS, DAZNCP, DAZAXIS
  C
  C  Construct rotation matrix to rotate about the X-axis by the angle -AZ
  C    and its derivative.
        CALL ROTAT (-AZ, 1, P1AZ)
         CALL DROTT (-AZ, -DAZ, 1, DP1AZ)
  C  Construct rotation matrix to rotate about the Y-axis by the angle -ZA
  C    and its derivative.
        CALL ROTAT (-ZA, 2, P2ZA)
         CALL DROTT (-ZA, -DZA, 2, DP2ZA)
  C  Combine the above rotation matrices and the derivative.
        CALL MMUL2 (P2ZA, P1AZ, P12)
         CALL MMUL2 (DP2ZA, P1AZ,  TT1)
         CALL MMUL2 (P2ZA, DP1AZ,  TT2)
         CALL MADD2 (TT1,    TT2, DP12)
  C  Rotate the three topocentric vectors and take derivatives of first two. 
        CALL VECRT (P12,   TNCP,  NCP1)
        CALL VECRT (P12, TCAXIS, AXIS1)
        CALL VECRT (P12, CTCSTR, STAR1)
         CALL VECRT (DP12,   TNCP,  DNCP1)
         CALL VECRT (DP12, TCAXIS, DAXIS1)
  C  The STAR1 vector should now be along the X-axis.
  C       print *,' STAR1  ', STAR1
  C
  C  Find the azimuth angle of the NCP and its time derivative in the new frame
         AZNCP = -ATAN2(NCP1(3),NCP1(2)) + HALFPI
         IF (AZNCP .LT. 0.D0) AZNCP = AZNCP + TWOPI
          IF(DABS(NCP1(2)) .GT. 1.D-16) THEN
            dAZNCP = -1.D0/(1.D0+((NCP1(3)/NCP1(2))**2)) * 
       *      ( 1.D0/NCP1(2) * DNCP1(3) - NCP1(3)/NCP1(2)**2 * DNCP1(2) ) 
          ELSE
            dAZNCP = 0.D0 
          ENDIF 
  C
  C  Find the azimuth angle of the fixed axis and its time derivative 
         AZAXIS = -ATAN2(AXIS1(3),AXIS1(2)) + HALFPI
         IF (AZAXIS .LT. 0.D0) AZAXIS = AZAXIS + TWOPI
          IF(DABS(AXIS1(2)) .GT. 1.D-16) THEN
            dAZAXIS = -1.D0/(1.D0+((AXIS1(3)/AXIS1(2))**2)) *
       *      ( 1.D0/AXIS1(2) * DAXIS1(3) - AXIS1(3)/AXIS1(2)**2 *
       *                                           DAXIS1(2) ) 
          ELSE
            dAZAXIS = 0.D0 
          ENDIF 
  C  Take difference (NCP to AXIS) and express as a clock-wise rotation, as
  C   seen from the ground looking up
         RANGL = AZNCP - AZAXIS 
          DRANGL = DAZNCP - DAZAXIS 
          IF (RANGL .LT. -PI) RANGL = RANGL + PI
          IF (RANGL .GT.  PI) RANGL = RANGL - PI
  C
  C  That's all!
        RETURN
        END
  