* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* radec.f
*
*
      SUBROUTINE RADEC(ALONG,ELAT,ALST,TELLAT,RADATE,DEDATE,DRA,DDEC)
C******************************************************************
C   THIS ROUTINE TAKES AS INPUT THE AZ AND EL OFFSET RELATIVE     *
C   TO THE FIELD CENTRE, THE L.S.T., THE OBSERVATORY LATITUDE     *
C   AND THE CURRENT EPOCH EQUATORIAL COORDINATES OF THE FIELD     *
C   CENTRE AND RETURNS THE RA AND DEC OFFSETS RELATIVE TO THE     *
C   FIELD CENTRE.                                                 *
C   ALONG , ELAT = OFFSETS RELATIVE TO FIELD CENTRE IN AZ-EL      *
C   ALST = L.S.T. OF DATA POINT                                   *
C   TELLAT = OBSERVATORY LATITUDE                                 *
C   RADATE , DEDATE = CURRENT EPOCH RA AND DEC OF FIELD CENTRE    *
C   DRA , DDEC = RA AND DEC OFFSETS RELATIVE TO FIELD CENTRE      *
C
C Input:
C   ALONG,ALAT : Offset to the field centre (arcsec)
C   ALST       : Local sideral time (sec)
C   TELLAT     : Telescope latitude (degrees)
C   RADATE,DEDATE : RA and DEC of the day (degrees)
C Output:
C   DRA,DDEC   : RA and DEC offsets relative to the field centre (degrees)
C******************************************************************
*
      INCLUDE 'radec.inc'
      INCLUDE 'const.inc'
*
      REAL*8   ALONG
      REAL*8   ALST,ALST1
      REAL*8   ELAT,TELLAT
      REAL*8   RADATE,DEDATE,RACUR
      REAL*8   DRA,DDEC
      REAL*8   DLO,DLA,HA
      REAL*8   SINR,COSR,SINE,COSE,SINPA,COSPA,SINLO,COSLO,SINLA,COSLA
      REAL*8   SINDD,COSDD,COSDRA,SINDRA
      REAL*8   PA
*
*
      ALST1 = ALST / 240. * DEG_TO_RAD
      RACUR = RADATE * DEG_TO_RAD
*
      DLO = ALONG / 3600. * DEG_TO_RAD
      DLA = ELAT /3600. * DEG_TO_RAD
      HA = RACUR - ALST1                 ! hour angle
*
C	SINDEC = DSIN(DECCUR)                ! computed in RDCON
C	COSDEC = DSQRT(1.D0 - SINDEC * SINDEC)
C	SINOL = DSIN(OBSLAT)
C	COSOL = DSQRT(1.D0 - SINOL * SINOL)
*
      SINR = DSIN(HA)
      COSR = DCOS(HA)
      SINE = SINDEC * SINOL + COSDEC * COSOL * COSR
C	E = DASIN(SINE)
      COSE = DSQRT(1.D0 - SINE * SINE)
      SINPA = COSOL * SINR / COSE
      COSPA = (SINOL - SINDEC *SINE) / (COSDEC * COSE)
      PA = DATAN2(SINPA,COSPA)
      SINLO = DSIN(DLO)
      COSLO = DCOS(DLO)
      SINLA = DSIN(DLA)
      COSLA = DCOS(DLA)
      SINDD = COSPA * SINLA - SINPA * COSLA *SINLO
      DDEC = DASIN(SINDD)
      COSDD = DCOS(DDEC)
      COSDRA = COSLA * COSLO / COSDD
      IF(ABS(SINPA).LT.1.D-6) THEN
         DRA = -DLO
      ELSE
         SINDRA = -(SINLA - COSPA * SINDD) / (SINPA * COSDD)
         DRA = DATAN2(SINDRA,COSDRA)
      ENDIF
      DRA = DRA / DEG_TO_RAD             ! DRA  in degrees
      DDEC = DDEC / DEG_TO_RAD           ! DDEC in degrees
      RETURN
      END
*
*
      SUBROUTINE RADEC_TO_AZEL(DRA,DDEC,LST,SITE_LATITUDE,RA_DAY,
     $DEC_DAY,DAZ,DEL)
C******************************************************************
C   THIS ROUTINE TAKES AS INPUT THE RA AND DEC OFFSET RELATIVE     *
C   TO THE FIELD CENTRE, THE L.S.T., THE OBSERVATORY LATITUDE     *
C   AND THE CURRENT EPOCH EQUATORIAL COORDINATES OF THE FIELD     *
C   CENTRE AND RETURNS THE AZ AND EL OFFSETS RELATIVE TO THE     *
C   FIELD CENTRE.                                                 *
C   DRA,DDEC  = OFFSETS RELATIVE TO FIELD CENTRE IN RA-DEC
C   LST = L.S.T. OF DATA POINT                                   *
C   SITE_LATITUDE = OBSERVATORY LATITUDE                                 *
C   RA_DAY , DEC_DAYDATE = CURRENT EPOCH RA AND DEC OF FIELD CENTRE    *
C   DAZ , DEL = AZ AND EL OFFSETS RELATIVE TO FIELD CENTRE      *
C
C Input:
C   DRA,DDEC   : Offset to the field centre (arcsec)
C   LST        : Local sideral time (sec)
C   SITE_LATITUDE     : Telescope latitude (degrees)
C   RA_DAY,DEC_DAY : RA and DEC of the day (degrees)
C Output:
C   DAZ,DEL   : AZ and EL offsets relative to the field centre (ARCSEC)
C******************************************************************
*
      INCLUDE 'radec.inc'
      INCLUDE 'const.inc'
*
      REAL*8   DRA,DDEC
      REAL*8   LST
      REAL*8   SITE_LATITUDE
      REAL*8   RA_DAY,DEC_DAY
      REAL*8   DAZ,DEL
      REAL*8   HA
      REAL*8   SIN_DEC, COS_DEC
      REAL*8   SIN_DRA,COS_DRA,SIN_DDEC,COS_DDEC
      REAL*8   SIN_EL, COS_EL, SIN_PA,  COS_PA
      REAL*8   SIN_DAZ,COS_DAZ, SIN_DEL,  COS_DEL
      REAL*8   SIN_LAT,COS_LAT
      REAL*8   SIN_HA ,COS_HA
      REAL*8   PA,EL
*
      HA = -(LST/240.D0 - RA_DAY)*DEG_TO_RAD   ! hour angle
*
      SIN_LAT = DSIN(SITE_LATITUDE*DEG_TO_RAD)
      COS_LAT = DCOS(SITE_LATITUDE*DEG_TO_RAD)
      SIN_DEC = DSIN(DEC_DAY*DEG_TO_RAD)
      COS_DEC = DCOS(DEC_DAY*DEG_TO_RAD)
      SIN_HA  = DSIN(HA)
      COS_HA  = DCOS(HA)
* Compute Az, El corresponding to RA, DEC, and LST
      SIN_EL  = SIN_DEC*SIN_LAT + COS_DEC*COS_HA*COS_LAT
      EL      = DASIN(SIN_EL)
      COS_EL  = SQRT(1.D0-SIN_EL**2)
      SIN_PA  = COS_LAT*SIN_HA/COS_EL
      COS_PA  = (SIN_LAT - SIN_DEC*SIN_EL)/(COS_EL*COS_DEC)
      PA      = ATAN2(SIN_PA,COS_PA)
* Compute Az and El offsets (Daz and Del)
      SIN_DRA  = DSIN(DRA*SEC_TO_RAD)
      COS_DRA  = DCOS(DRA*SEC_TO_RAD)
      SIN_DDEC = DSIN(DDEC*SEC_TO_RAD)
      COS_DDEC = DCOS(DDEC*SEC_TO_RAD)
      SIN_DEL = -SIN_DRA*SIN_PA*COS_DDEC + COS_PA*SIN_DDEC
      DEL = DASIN(SIN_DEL)               ! del in ]-pi/2,pi/2]  (del~0)
      COS_DEL = DCOS(DEL)
      COS_DAZ = COS_DRA*COS_DDEC/COS_DEL
      IF (ABS(PA).GT.1.D-6) THEN
         SIN_DAZ = (-SIN_DDEC + COS_PA*SIN_DEL)/(SIN_PA*COS_DEL)
         DAZ = DATAN2(SIN_DAZ,COS_DAZ)
      ELSE
         DAZ = -DRA*SEC_TO_RAD
      ENDIF
* Conversion in arcsec
      DAZ = DAZ/SEC_TO_RAD
      DEL = DEL/SEC_TO_RAD
      RETURN
      END
*
      SUBROUTINE PREC(ISTYLE,EEPOCH,RACEN,DECCEN,DRA,DDEC,RA50,DEC50)
C------------------------------------------------------------------------
C    THIS SUBROUTINE PERFORM PRECESSION TO EPOCH 1950.
C    ISTYLE = 1 MEANS THE OUTPUT WILL BE RACEN, DECCEN TRANSFORMED  TO
C               EPOCH 1950. AND PLACED IN RA50, DEC50 IN DECIMAL DEGREES.
C                * input  : RACEN,DECCEN (DRA=DDEC=0)
C                * output : RA50,DEC50
C
C    ISTYLE = 2 MEANS THAT THE OFFSETS FROM FIELD CENTRE DRA,DDEC WILL BE
C               TRANSFORMED TO OFFSETS IN EPOCH 1950. USING THE COORDS
C               OF THE FIELD CENTRE IN CURRENT EPOCH IN RACEN, DECCEN
C               AND IN EPOCH 1950. IN RA50,DEC50.
C                * input  : DRA,DDEC,RACEN,DECCEN
C                * output : RA50,DEC50
C
C    ISTYLE = 3  * input  : RA50,DEC50 (DRA=DDEC=0)
C                * output : RACEN,DECCEN
C
C    EPOCH = CURRENT EPOCH IN DECIMAL YEARS A.D.
C------------------------------------------------------------------------
      INCLUDE 'radec.inc'
      INCLUDE 'const.inc'
*
      INTEGER  ISTYLE
      REAL*8   EEPOCH,RACEN,DECCEN,DRA,DDEC,RA50,DEC50
*
      REAL*8   DELR,DELD,DC
      REAL*8   RACUR,DECCUR,RAH,DD,RAP50,DECP50
*
C	PI=DACOS(-1.0D0)          ! initialized in routine SETUP
C	CORR1 = PI / 180.D0
C	CORR2 = 180.D0 / PI
C       NYRI = 1950
C       NYRF = INT(EPOCH)
C       MO = 1
C       NDA = INT((EPOCH-DFLOAT(NYRF))*365.2422)
*
*
      IF (ISTYLE.EQ.3) THEN
         RACUR  = RA50 * DEG_TO_RAD
         DECCUR = DEC50 * DEG_TO_RAD
      ELSE
         RACUR  = RACEN * DEG_TO_RAD
         DECCUR = DECCEN * DEG_TO_RAD
      ENDIF
      DRA = DRA * DEG_TO_RAD
      DDEC = DDEC * DEG_TO_RAD
      DECCUR = DECCUR + DDEC
      RACUR = RACUR + DRA / DCOS(DECCUR)
      RAH = RACUR
      DD = DECCUR
*
      CALL MOVE(NYRI,NYRF,MO,NDA,RAH,DD,DELR,DELD,DC)
*
      IF (ISTYLE.EQ.3) THEN
         RAH = RAH + 0.5D0 * DELR
         DD = DD + 0.5D0 * DELD
      ELSE
         RAH = RAH - DELR * 0.5D0
         DD = DD - DELD * 0.5D0
      ENDIF
*
* Precession routine MOVE is called a second time with the mean position
      CALL MOVE(NYRI,NYRF,MO,NDA,RAH,DD,DELR,DELD,DC)
*
      IF (ISTYLE .EQ. 1) THEN
         RA50 = RACUR - DELR
         DEC50 = DECCUR - DELD
         RA50 = RA50 / DEG_TO_RAD
         DEC50 = DEC50 / DEG_TO_RAD
      ELSEIF (ISTYLE .EQ. 2) THEN
         RAP50 = RACUR - DELR
         DECP50 = DECCUR - DELD
         DRA = (RAP50 - RA50 * DEG_TO_RAD) * DCOS(DECP50)
         DDEC = DECP50 - DEC50 * DEG_TO_RAD
         DRA = DRA / DEG_TO_RAD
         DDEC = DDEC / DEG_TO_RAD
      ELSE
         RACEN = RACUR + DELR
         DECCEN = DECCUR + DELD
         RACEN = RACEN / DEG_TO_RAD
         DECCEN = DECCEN / DEG_TO_RAD
      ENDIF
      RETURN
      END
*
      SUBROUTINE MOVE(NNYRI,NNYRF,MMO,NNDA,RA,D,DELR,DELD,DC)
C----------------------------------------------------------------------
C THIS IS A SIMPLE PRECESSION ROUTINE.
C    NYRI = INITIAL EPOCH IN YRS A.D. (I.E. 1950).
C    NYRF = THE YEAR TO WHICH WE ARE PRECESSING.
C    MO = MONTH I!N YEAR NYRF (1 TO 12)
C    NDA = DAY IN MONTH MO.
C         IF DAY NO. IS KNOWN PUT IN NDA WITH MO SET TO 1.
C    RA = RIGHT ASCENSION OF SOURCE IN RADIANS.
C    D = DECLINATION OF SOURCE IN RADIANS.
C         BEST RESULTS IF RA & D ARE FOR THE MEAN EPOCH.
C    DELR = CORRN TO BE ADDED TO MEAN R.A. FOR YEAR NYRI.
C    DELD = CORRN TO BE ADDED TO MEAN DEC FOR YEAR NYRI.
C----------------------------------------------------------------------
      INCLUDE 'radec.inc'
      INCLUDE 'const.inc'
*
      INTEGER  NNYRI,NNYRF,MMO,NNDA
      REAL*8   RA,D,DELR,DELD,DC
      REAL*8   SND,CSD,TND,CSR,SNR
*
      SND = DSIN(D)
      CSD = DCOS(D)
      TND = SND/CSD
      CSR = DCOS(RA)
      SNR = DSIN(RA)
C	AL = 30 * (MO - 1) + NDA
C	TO = DFLOAT(NYRI - 1900)/100.0D0
C	T = (DFLOAT(NYRF - NYRI) + AL / 365.2421988D0) / 100.0D0
C	ZETAO=(2304.25D0+1.396D0*TO)*T+0.302D0*T**2+0.018D0*T**3
C	Z=ZETAO+0.791D0*T**2
C	THETA=(2004.682D0-0.853D0*TO)*T-0.426D0*T**2-0.042D0*T**3
C	AM=(ZETAO+Z)*4.848136811D-6
C	AN=THETA*4.848136811D-6
C	ALAM=(0.985647D0*AL+278.5D0)*0.0174532925D0
C	SNL=DSIN(ALAM)
C	CSL=DCOS(ALAM)
      DELR=-9.92413605D-5*(SNL*SNR+F3*CSR)/CSD
     $+AM+AN*SNR*TND
      DELD=-9.92413605D-5*(SNL*CSR*SND-F3*SNR*SND
     $+F4*CSD)+AN*CSR
C	OMEGA=259.183275D0-1934.142D0*(TO+T)
C	ARG=OMEGA*0.0174532925D0
C	DLONG=-8.3597D-5*DSIN(ARG)
C	DOBLQ=4.4678D-5*DCOS(ARG)
      DELR=DELR+DLONG*(0.91745051D0+0.39784993D0*SNR*TND)-CSR*TND*DOBLQ
      DELD=DELD+F5*CSR+SNR*DOBLQ
      RETURN
      END
*
      SUBROUTINE DELCEN(RA,DEC,RA2,DEC2,ARAARA,DELDEL)
C*****************************************************************
C THIS SUBROUTINE COMPUTES THE R.A. AND DEC DIFFERENCE BETWEEN   *
C MAP CENTRES FOR THE MAP IN AMAP AND THE MAP IN CMAP.           *
C THE DIFFERENCES IN DECIMAL DEGREES ARE OUTPUT IN ARAARA,DELDEL.*
c
c ra ,dec  : amap (input map)
c ra2,dec2 : cmap (resulting map)
C*****************************************************************
*
      INCLUDE 'const.inc'
      REAL*8   RA,DEC
      REAL*8   DEC1,RA1,DEC2,RA2,ARAARA,DELDEL
      REAL*8   DELAPH,SINDEL,COSDEL,SINARA,COSARA
*
      RA1  = RA   * DEG_TO_RAD
      DEC1 = DEC  * DEG_TO_RAD
      RA2  = RA2  * DEG_TO_RAD
      DEC2 = DEC2 * DEG_TO_RAD
      DELAPH = RA1 - RA2
      SINDEL = DSIN(DEC1)*DCOS(DEC2)-DCOS(DEC1)*
     $DSIN(DEC2)*DCOS(DELAPH)
      DELDEL = DASIN(SINDEL)
      COSDEL = DCOS(DELDEL)
      SINARA = DSIN(DELAPH) * DCOS(DEC1)
      COSARA = DSIN(DEC1)*DSIN(DEC2)+DCOS(DEC1)*
     $DCOS(DEC2)*DCOS(DELAPH)
      ARAARA = DATAN2(SINARA,COSARA)
      ARAARA = ARAARA / DEG_TO_RAD
      DELDEL = DELDEL / DEG_TO_RAD
C       PRINT *,'ARAARA =',ARAARA,'DELDEL =',DELDEL
      RETURN
      END
*
      SUBROUTINE SETUP
C*****************************************************************
C
C
C
C*****************************************************************
*
      INCLUDE 'radec.inc'
      INCLUDE 'const.inc'
*
      REAL*8   AL,TO,T,ZETAO,Z,THETA
      REAL*8   ALAM
      REAL*8   CSL,OMEGA,ARG
      REAL*8   T2,T3
      NYRI = EPOCHFR
      NYRF = INT(EPOCHTO)
      MO = 1                             ! because NDA is the total number of day
      NDA = INT((EPOCHTO-DFLOAT(NYRF))*365.2422)
      AL = 30 * (MO - 1) + NDA
      TO = DFLOAT(NYRI - 1900)/100.0D0
      T = (DFLOAT(NYRF - NYRI) + AL / 365.2421988D0) / 100.0D0
*
      ZETAO=(2304.25D0+1.396D0*TO)*T+0.302D0*T**2+0.018D0*T**3
      Z=ZETAO+0.791D0*T**2
      THETA=(2004.682D0-0.853D0*TO)*T-0.426D0*T**2-0.042D0*T**3
      AM = (ZETAO+Z) * SEC_TO_RAD
      AN = THETA * SEC_TO_RAD
*
      T2=T*T
      T3=T*T2
      IF (ABS(EPOCHFR-1950D0).LT.ABS(EPOCHFR-2000D0)) THEN
*
* J1950 corrections from Almanach 1981  (p B18)
*
         ZETAO = 0.6402633D0*T + 0.0000839D0*T2 + 0.0000050D0*T3
         Z     = 0.6402633D0*T + 0.0003036D0*T2 + 0.0000050D0*T3
         THETA = 0.5567376D0*T - 0.0001183D0*T2 - 0.0000117D0*T3
      ELSE
*
* J2000 corrections from Almanach 1995  (p B18)
*
         ZETAO = 0.6406161D0*T + 0.0000839D0*T2 + 0.0000050D0*T3
         Z     = 0.6406161D0*T + 0.0003041D0*T2 + 0.0000051D0*T3
         THETA = 0.5567530D0*T - 0.0001185D0*T2 - 0.0000116D0*T3
      ENDIF
*
      AM = (ZETAO+Z) * DEG_TO_RAD
      AN = THETA * DEG_TO_RAD
*
      ALAM=(0.985647D0*AL+278.5D0)* DEG_TO_RAD
      SNL=DSIN(ALAM)
      CSL=DCOS(ALAM)
      OMEGA=259.183275D0-1934.142D0*(TO+T)
      ARG=OMEGA * DEG_TO_RAD
      DLONG=-8.3597D-5*DSIN(ARG)
      DOBLQ=4.4678D-5*DCOS(ARG)
      F3=0.91745051D0*CSL
      F4=0.39784993D0*CSL
      F5=0.39784993D0*DLONG
      RETURN
      END
*
      SUBROUTINE DANGLE(ALONG,ALAT,A)
c***********************************************************************
c version 1.0  mpifr cyber edition  22 may 1977.                       *
c see nod 2 manual page m 9.1                                          *
c c.g haslam       mar  1972.                                          *
c this routine recovers the longitude and latitude angles, along       *
c (0 - 360) and alat ( +.- 180) in deg_to_radrees, which correspond to the    *
c vector a.                                                            *
c***********************************************************************
*
      REAL*8   A(3)
      REAL*8   ALONG,ALAT
      REAL*8   AA
*
      INCLUDE 'const.inc'
*
      AA     = A(1)*A(1)+A(2)*A(2)
      AA     = DSQRT(AA)
      ALAT   = 90.0D0
      IF (AA.GE.0.000001D0) THEN
         AA = A(3)/AA
         ALAT = DATAN(AA)/DEG_TO_RAD     ! latitude in deg_to_radrees
      ENDIF
*
      IF (A(2).NE.0.D0 .OR. A(1).NE.0.D0) THEN
         ALONG=DATAN2(A(2),A(1))/DEG_TO_RAD    ! longitude in deg_to_radrees
      ELSE
         ALONG = 0.D0
      ENDIF
*
      IF (ALONG.LT.0.0) ALONG=ALONG+360.0D0
      RETURN
      END
*
*
      SUBROUTINE DCOSIN(ALONG,ALAT,A)    ! direction cosinus
c***********************************************************************
c version 1.0  mpifr cyber edition  22 may 1977.                       *
c see nod 2 manual page m 9.1                                          *
c c.g haslam       mar  1972.                                          *
c this routine computes the vector a for the angles, along and alat,   *
c which should be specified in degrees.                                *
c***********************************************************************
*
      REAL*8   A(3)
      REAL*8   ALONG,ALAT
      REAL*8   X,Y,Z
*
      INCLUDE 'const.inc'
*
      X    = DEG_TO_RAD*ALONG            ! x in radians
      Y    = DEG_TO_RAD*ALAT             ! y in radians
      Z    = DCOS(Y)
      A(1) = DCOS(X)*Z
      A(2) = DSIN(X)*Z
      A(3) = DSIN(Y)
      RETURN
      END
*
      SUBROUTINE MATMUL(MAT,A,B,N)
c***********************************************************************
c version 1.0  mpifr cyber edition  22 may 1977.  G.Haslam
c version 2.0  iram
c
c this routine provides the transformation of vector a to vector b     *
c using the n dimensional direction cosine array, mat.                *
c if n<0, multiply by transpose matrix
c***********************************************************************
*
      REAL*8    A(1),B(1),MAT(3,3)
      INTEGER   I,J,N
*
      IF (N.GT.0) THEN
         DO I=1,3
            B(I) = 0.D0
            DO J=1,3
               B(I) = B(I)+MAT(I,J)*A(J)
            ENDDO
         ENDDO
      ELSE
         DO I=1,3
            B(I) = 0.D0
            DO J=1,3
               B(I) = B(I)+MAT(J,I)*A(J)
            ENDDO
         ENDDO
      ENDIF
      RETURN
      END