C file CEA.FOR
C ----------------------------------------------------------------------!-------
      SUBROUTINE BETTERCEA(A,B,C,COSA,COSB,COSC,SINA,SINB,SINC,GAMMA)
C ----------------------------------------------------------------------!-------
C     ADDED 6 SEPT 2005 TO HANDLE ANGLES NEAR 0 AND PI
C
C     SPHERICAL TRIANGLE WITH SIDES A,B,C ARE OPPOSITE ANGLE ALPHA,
C     BETA,GAMMA.
C     ALPHA AND BETA ARE NOT USED
C
C     INPUT: A,B,COSA,COSB,SINA,SINB,GAMMA
C     OUTPUT: C,COSC,SINC
      implicit double precision (a-h,o-z)
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
      DATA K/0/
      SING=DSIN(GAMMA)
      COSG=DSIN(PID2-GAMMA)
      COSC=COSA*COSB+SINA*SINB*COSG
C
      IF (COSG.GT.0) THEN
         SINASB=DSIN(A-B)
         COSASB=DSIN(PID2-A+B)
         NTEST=2.5D0+DSIGN(1.D0,COSASB)+DSIGN(.5D0,COSC)
C         PRINT*,'NTEST = ',NTEST
         SELECT CASE (NTEST)
         CASE (1)
            SINC=DSQRT((1-COSC)*((SINASB**2/(1.D0-COSASB))
     1      -SINA*SINB*(SING**2/(1.D0+COSG))))
            PRINT*,'CASE #: ',NTEST
            PRINT*,'COSC = ',COSC
            PRINT*,'SINC = ',SINC
            ONE=SINC**2+COSC**2
            PRINT*,'ONE = ',ONE
            PRINT*,'SINA = ',SINA
            PRINT*,'SINB = ',SINB
            PRINT*,'COSA = ',COSA
            PRINT*,'COSB = ',COSB
            PRINT*,'COSASB = ',COSASB
            PRINT*,'A = ',A
            PRINT*,'B = ',B
         CASE (2)
            SINC=DSQRT((1+COSC)*(2.D0-((SINASB**2)/(1.D0-COSASB))
     1      +SINA*SINB*(SING**2/(1.D0+COSG))))
            X=1+COSC
            Y=(SINASB**2)/(1.D0-COSASB)
            Z=SINA*SINB*(SING**2/(1.D0+COSG))
            PRINT*,'CASE #: ',NTEST
            PRINT*,'COSASB = ',COSASB
            PRINT*,'COSC = ',COSC
            PRINT*,'SINC = ',SINC
            ONE=SINC**2+COSC**2
            PRINT*,'ONE = ',ONE
            PRINT*,'X = ',X
            PRINT*,'Y = ',Y
            PRINT*,'Z = ',Z
            PRINT*,'SINA = ',SINA
            PRINT*,'SINB = ',SINB
            PRINT*,'COSA = ',COSA
            PRINT*,'COSB = ',COSB
            PRINT*,'COSASB = ',COSASB
            PRINT*,'A = ',A
            PRINT*,'B = ',B
            ASB=A-B
            PRINT*,'ASB = ',ASB
         CASE (3)
            SINC=DSQRT((1-COSC)*(2.D0-((SINASB**2)/(1.D0+COSASB))
     1      -SINA*SINB*(SING**2/(1.D0+COSG))))
            PRINT*,'CASE #: ',NTEST
            PRINT*,'COSASB = ',COSASB
            PRINT*,'COSC = ',COSC
            PRINT*,'SINC = ',SINC
            ONE=SINC**2+COSC**2
            PRINT*,'ONE = ',ONE
         CASE (4)
            SINC=DSQRT((1+COSC)*((SINASB**2/(1.D0+COSASB))
     1      +SINA*SINB*(SING**2/(1.D0+COSG))))
C             K=K+1
C             IF (K.GT.1000) THEN
C             K=0
C             PRINT*,'CASE #: ',NTEST
C             PRINT 999,'COSC = ',COSC
C             PRINT 999,'SINC = ',SINC
C             ONE=SINC**2+COSC**2
C             PRINT*,'ONE = ',ONE
C             ENDIF
         ! THIS CASE IS FUNCTIONAL
999    FORMAT(A,E60.50)
         END SELECT
C
      ELSE
         SINAPB=DSIN(A+B)
         COSAPB=DSIN(PID2-A-B)
         NTEST=6.5D0+DSIGN(1.D0,COSAPB)+DSIGN(.5D0,COSC)
C         PRINT*,'NTEST = ',NTEST
         SELECT CASE (NTEST)
         CASE (5)
            SINC=DSQRT((1-COSC)*((SINAPB**2/(1.D0-COSAPB))
     1      +SINA*SINB*(SING**2/(1.D0-COSG))))
C            PRINT*,'CASE #: ',NTEST
C            PRINT*,'COSC = ',COSC
C            PRINT*,'SINC = ',SINC
C            ONE=SINC**2+COSC**2
C            PRINT*,'ONE = ',ONE
         ! THIS CASE IS FUNCTIONAL
         CASE (6)
            SINC=DSQRT((1+COSC)*(2.D0-((SINAPB**2)/(1.D0-COSAPB))
     1      -SINA*SINB*(SING**2/(1.D0-COSG))))
C            PRINT*,'CASE #: ',NTEST
C            PRINT*,'COSC = ',COSC
C            PRINT*,'SINC = ',SINC
C            ONE=SINC**2+COSC**2
C            PRINT*,'ONE = ',ONE
         ! THIS CASE IS FUNCTIONAL
         CASE (7)
            SINC=DSQRT((1-COSC)*(2.D0-((SINAPB**2)/(1.D0+COSAPB))
     1      +SINA*SINB*(SING**2/(1.D0-COSG))))
            PRINT*,'CASE #: ',NTEST
            PRINT*,'COSC = ',COSC
            PRINT*,'SINC = ',SINC
            ONE=SINC**2+COSC**2
            PRINT*,'ONE = ',ONE
         CASE (8)
            SINC=DSQRT((1+COSC)*((SINAPB**2/(1.D0+COSAPB))
     1      -SINA*SINB*(SING**2/(1.D0-COSG))))
C            PRINT*,'CASE #: ',NTEST
C            PRINT*,'COSC = ',COSC
C            PRINT*,'SINC = ',SINC
C            ONE=SINC**2+COSC**2
C            PRINT*,'ONE = ',ONE
         ! THIS CASE IS FUNCTIONAL
         END SELECT
      ENDIF
C
      C=DATAN2(SINC,COSC)
C
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE BETTERCEA
C ----------------------------------------------------------------------!-------
