PROGRAM CONFND
C PROGRAM TO FIND CONSTELLATION NAME  FROM A POSITION
C ---> Modified FO, CDS Strasbourg, DARSIN is NOT a standard F77 function...
C	... but DASIN !!!
C
C THE FIRST RECORD IN THE DATA SET MUST BE THE EQUINOX FOR THE POSITIONS
C  IN THE FORMAT XXXX.X
C THE REMAINING RECORDS MUST BE THE POSITION IN HOURS AND DECIMALS OF
C  AN HOUR FOLLOWED BY THE DECLINATION IN DEGREES AND FRACTIONS OF A
C  DEGREE IN THE FORMAT: F7.4,F8.4
C THE OUTPUT WILL REPEAT THE POSITION ENTERED AND GIVE THE THREE LETTER
C  ABBREVIATION OF THE CONSTELLATION IN WHICH IT IS LOCATED
C
      CHARACTER * 3 CON
      DOUBLEPRECISION ARAD,DRAD,A,D,PI4,E75,E,CONVH,CONVD
      DATA CONVH/0.2617993878D0/,CONVD/0.1745329251994D-01/,
     1PI4/6.28318530717948D0/,E75/1875.D0/
      OPEN(2,FILE='data.dat',READONLY,STATUS='OLD')
      WRITE(6,1000)
 1000 FORMAT('   CONSTALLATION derived from POSITION    ')
      WRITE(6,1001)
 1001 FORMAT('Type the Equinox, format' /'YYYY.y')
      READ(5,FMT='(D6.1)') E
    5 WRITE (6,1002)
 1002 FORMAT('Enter the position '/ 'HH.hhhh+DD.dddd')
      READ(5,FMT='(F7.4,F8.4)',END=200,ERR=210) RAH,DECD
C PRECESS POSITION TO 1875.0 EQUINOX
      ARAD = CONVH * RAH
      DRAD = CONVD * DECD
      CALL HGTPRC(ARAD,DRAD,E,E75,A,D)
      IF(A.LT.0.D0)A=A+PI4
      IF(A.GE.PI4)A=A-PI4
      RA= A/CONVH
      DEC=D/CONVD
C FIND CONSTELLATION SUCH THAT THE DECLINATION ENTERED IS HIGHER THAN
C  THE LOWER BOUNDARY OF THE CONSTELLATION WHEN THE UPPER AND LOWER
C  RIGHT ASCENSIONS FOR THE CONSTELLATION BOUND THE ENTERED RIGHT
C  ASCENSION
   30 READ (2,120,END=220,ERR=230) RAL,RAU,DECL,CON
  120 FORMAT(2F8.4,F9.4,1X,A3)
      IF(DECL.GT.DEC) GOTO 30
   10 IF (RAU.LE.RA) THEN
          READ (2,120,END=220,ERR=230) RAL,RAU,DECL,CON
          GO TO 10
      ELSE
      END IF
   20 IF(RAL.GT.RA) THEN
          READ (2,120,END=220,ERR=230) RAL,RAU,DECL,CON
          GO TO 20
      ELSE
      END IF
C IF CONSTELLATION HAS BEEN FOUND, WRITE RESULT AND RERUN PROGRAM FOR
C  NEXT ENTRY.  OTHERWISE, CONTINUE THE SEARCH BY RETURNING TO RAU
      IF (RA.GE.RAL.AND.RA.LT.RAU.AND.DECL.LE.DEC) THEN
          WRITE(6,100) RAH, DECD,CON
  100 FORMAT(' RA =',F8.4,' DEC =',F9.4,'  is in Constellation ',A3)
      ELSE IF(RAU.LT.RA) THEN
          GOTO 10
      ELSE
          WRITE(6,110) RAH,DECD
  110 FORMAT(' Constellation NOT FOUND for: RA =',F8.4,'  DEC =',F9.4)
      END IF
      REWIND 2
      GOTO 5
  200 WRITE(6,300) RAH,DECD,E
  300 FORMAT(' End of input positions after: RA = ',F7.4,3X,'DEC = ',
     1 F8.4 / ' The Equinox for these positions is ',F6.1)
      GOTO 400
  210 WRITE(6,310) RA,DEC
  310 FORMAT(' INPUT ERROR IN POSITION LIST AFTER: RA = ',F7.4,'  DEC =
     1 ',F7.4)
      GOTO 400
  220 WRITE(6,320) RAL,RAU,DECL,CON
  320 FORMAT(' END OF CONSTELLATION BOUNDARIES AFTER ',2F8.4,F9.4,1X,A3)
      GOTO 400
  230 WRITE(6,330) RAL,RAU,DECL,CON
  330 FORMAT(' INPUT ERROR AFTER RECORD: ',2F8.4,F9.4,1X,A3)
  400 STOP
      END
      SUBROUTINE HGTPRC (RA1,DEC1,EPOCH1,EPOCH2,RA2,DEC2)
C      HERGET PRECESSION, SEE P. 9 OF PUBL. CINCINNATI OBS. NO. 24
C INPUT=  RA1 AND DEC1 MEAN PLACE, IN RADIANS, FOR EPOCH1, IN YEARS A.D.
C OUTPUT= RA2 AND DEC2 MEAN PLACE, IN RADIANS, FOR EPOCH2, IN YEARS A.D.
      DOUBLE PRECISION CDR,RA1,DEC1,EPOCH1,EPOCH2,RA2,DEC2,X1(3),X2(3),
     DR(3,3),T,ST,A,B,C,EP1,EP2,CSR
      DATA CDR/0.17453292519943D-01/
      DATA EP1,EP2/0.D0,0.D0/
C      COMPUTE INPUT DIRECTION COSINES
      A=DCOS(DEC1)
      X1(1)=A*DCOS(RA1)
      X1(2)=A*DSIN(RA1)
      X1(3)=DSIN(DEC1)
C      SET UP ROTATION MATRIX (R)
      IF(EP1.EQ.EPOCH1.AND.EP2.EQ.EPOCH2) GO TO 10
      CSR=CDR/3600.D0
      T=0.001D0*(EPOCH2-EPOCH1)
      ST=0.001D0*(EPOCH1-1900.D0)
      A=CSR*T*(23042.53D0+ST*(139.75D0+0.06D0*ST)+T*(30.23D0-0.27D0*ST+
     $18.0D0*T))
      B=CSR*T*T*(79.27D0+0.66D0*ST+0.32D0*T)+A
      C=CSR*T*(20046.85D0-ST*(85.33D0+0.37D0*ST)+T*(-42.67D0-0.37D0*ST
     $-41.8D0*T))
      SINA=DSIN(A)
      SINB=DSIN(B)
      SINC=DSIN(C)
      COSA=DCOS(A)
      COSB=DCOS(B)
      COSC=DCOS(C)
      R(1,1)=COSA*COSB*COSC-SINA*SINB
      R(1,2)=-COSA*SINB-SINA*COSB*COSC
      R(1,3)=-COSB*SINC
      R(2,1)=SINA*COSB+COSA*SINB*COSC
      R(2,2)=COSA*COSB-SINA*SINB*COSC
      R(2,3)=-SINB*SINC
      R(3,1)=COSA*SINC
      R(3,2)=-SINA*SINC
      R(3,3)=COSC
C      PERFORM THE ROTATION TO GET THE DIRECTION COSINES AT EPOCH2
   10 DO 20 I=1,3
      X2(I)=0.D0
      DO 20 J=1,3
   20 X2(I)=X2(I)+R(I,J)*X1(J)
      RA2=DATAN2(X2(2),X2(1))
      IF(RA2.LT.0) RA2 = 6.28318530717948D0 + RA2
      DEC2=DASIN(X2(3))
      RETURN
      END