subroutine PRECESSION_APPROX(JEpoch1,JEpoch2,LAMBDA,BETA)
            double precision    JEpoch1
            double precision    JEpoch2
            real                LAMBDA
            real                BETA

        double precision        T
        double precision        A
        double precision        B
        double precision        C

        ! Needed for g77 when using goniod.f
        ! What happens if these are available as intrinsic functions??

        double precision        dsind
        double precision        dcosd

        if (abs(BETA) .gt. 90.0) stop 'PRECESSION_APPROX: invalid BETA'

        T = 0.01*(JEpoch2-2000.)                                !Time interval in Julian centuries
        A = (1.396971 + 0.000309*T)*T                           !Form. (9.38), p. 224
        B = (0.013056 - 0.000010*T)*T
        C = 5.1236 + 0.2416*T
        if (abs(BETA) .eq. 90.0) then
            LAMBDA = 174.8764+0.9137*T-sign(1.,BETA)*90.0       !Form. (9.11a), p. 214
            DUM    = (0.47-0.0007*T)/36.0*T                     !Form. (9.11b), p. 214
            BETA   = BETA+DUM*sign(1.0,BETA)
        else
            DUM    = LAMBDA + A-B*tand(BETA)*dcosd(LAMBDA+C)
            BETA   = BETA   + B*dsind(LAMBDA+C)                 !J2000.0 ---> JEpoch2
            LAMBDA = DUM
        end if
        LAMBDA = amod(amod(LAMBDA,360.0)+360.0,360.0)

        return
        end