* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* nic_xycoord.f
*
*
      SUBROUTINE INIT_XYCOORD(SRA,SDEC,AZI,ELE,XCOORD,YCOORD,
     $ST,SITE_LATITUDE,NPOINTS)
*
c**************************************************************************
c   This subroutine converts az and el first into az, el offsets from     *
c   field centre and then into nasmyth coords.
c   (cf subroutine dazel (Bonn MPIfR))
c
c   It should take refraction, pointing and homology into account (?)
c**************************************************************************
*
      INCLUDE 'const.inc'
*
      REAL*4 AZI(1),ELE(1),XCOORD(1),YCOORD(1),ST(1)
      REAL*8 SITE_LATITUDE
      INTEGER NPOINTS
*
      REAL*8 BMAT(9),AVEC(3),BVEC(3),SRA,SDEC,SAZ,SEL,SHA,LATITUDE
*     REAL*4 COL,PNUL,P14         ! variable used in previous version
*     REAL*4 DELAY
*
*
      REAL*4 CELE,SELE,TELE,CAZI,SAZI,DAZI
      REAL*8 DIFF_AZI,DIFF_ELE,DAZ,DEL
*
      INTEGER I
*
*  Initialize matrix bmat
*
      LATITUDE = SITE_LATITUDE/DEG_TO_RAD! latitude in degrees
      CALL AMSET(LATITUDE,BMAT)
*
      DO I=1,NPOINTS
*
*  Compute SAZ and SEL (Azimuth and Elevation from the field centre)
*
         SHA = DBLE (ST(I)) / 240.D0 - SRA     ! hour angle in degrees
*
         CALL DCOSIN(SHA,SDEC,AVEC)
         CALL MATMUL(BMAT,AVEC,BVEC,-3)
         CALL DANGLE(SAZ,SEL,BVEC)
*
         CELE=COS(DEG_TO_RAD*ELE(I))
         SELE=SIN(DEG_TO_RAD*ELE(I))
         TELE=SELE/CELE
         CAZI=COS(DEG_TO_RAD*AZI(I))
         SAZI=SIN(DEG_TO_RAD*AZI(I))
         DAZI=SIN(2.*DEG_TO_RAD*AZI(I))
*
*
c         daz      = -1877.33 - 2.79/cele + 5.94*tele -
c     $        4.42*tele*cazi -
c     $        9.88*tele*sazi +
c     $        2.*dazi/cele +
c     $        (col - 2000.*0.786/10501000./sec_to_rad+0.5*0.01)/cele
*
         DAZ = 0
         DIFF_AZI = AZI(I)-SAZ - DAZ/3600.D0
*
* --- map from -180.d0 to 180.d0
*
         DO WHILE (DIFF_AZI.GT.180.D0)
            DIFF_AZI = DIFF_AZI - 360.D0
         ENDDO
         DO WHILE (DIFF_AZI.LT. -180.D0)
            DIFF_AZI = DIFF_AZI + 360.D0
         ENDDO
*
c         del     = 4.42*sazi - 9.88*cazi + 344.83 - 82.65*cele -
c     $        52.23*sele + p14/tele - 0.002/tele**3 + pnul +
c     $        (30180.*cele+1680.)*0.786/10501000./sec_to_rad +
c     $        (-1.36*cele- 0.25)*0.02
*
         DEL = 0
         DIFF_ELE = ELE(I)-SEL- DEL/3600.D0
*
* --- map from -180.d0 to 180.d0
*
         DO WHILE (DIFF_ELE.GT.180.D0)
            DIFF_ELE = DIFF_ELE - 360.D0
         ENDDO
         DO WHILE (DIFF_ELE.LT. -180.D0)
            DIFF_ELE = DIFF_ELE + 360.D0
         ENDDO
*
         XCOORD(I)=(DIFF_AZI)*36000.D0*CELE
         YCOORD(I)=(DIFF_ELE)*36000.D0
*	type *,I,XCOORD(I),YCOORD(I)
*
      ENDDO
      RETURN
      END
*
*
      SUBROUTINE FIND_SUBSCAN_LIMIT (LST,SUBSCAN_TIME,
     $SUBSCAN_IFIRST,SUBSCAN_ILAST,XCOORD,YCOORD,
     $OAZM,SINT,NSUBSCAN,NRECORD)
C**************************************************************
C
C     FIND_SUBSCAN_LIMIT() :
C        Identify first record and last record for each subscan
C
C**************************************************************
*
      REAL*4    LST(1)
      REAL*4    SUBSCAN_TIME(1)
      INTEGER*4 SUBSCAN_IFIRST(1),SUBSCAN_ILAST(1)
      REAL*4    XCOORD(1),YCOORD(1)
      REAL*4    OAZM,SINT
      INTEGER   NSUBSCAN,NRECORD
*
      INTEGER   II,ISUBSCAN
      INTEGER   FIRST_RECORD,LAST_RECORD
      REAL      SUBSCAN_START
      REAL      OFFSET_MIDNIGHT
      REAL      XLIMIT,DELTA_Y,DELTA_INF
*
      II = 1
      OFFSET_MIDNIGHT = 0.
      XLIMIT = ABS(OAZM)*3600.*10 + SINT/10.   ! in deciarcseconds
      DELTA_INF = (SINT*10)/4            ! delta_inf= sint/4 (in deciarcseconds)
*
      FIRST_RECORD=1
      DO WHILE ( ABS(XCOORD(II)).GT.XLIMIT .AND. II.LE.NRECORD)
         II = II+1
      ENDDO
      DO ISUBSCAN=1,NSUBSCAN
         FIRST_RECORD = 0
         LAST_RECORD  = 0
*
         SUBSCAN_START = SUBSCAN_TIME(ISUBSCAN)- 0.01 + OFFSET_MIDNIGHT
         IF (OFFSET_MIDNIGHT.EQ.0 .AND. LST(II).GT.86400.)
     $   SUBSCAN_START = SUBSCAN_START+86400.
*
         DO WHILE ((FIRST_RECORD.EQ.0 .OR. LAST_RECORD.EQ.0)
     $      .AND. II.LE.NRECORD)
            IF (II.GT.1) THEN
               IF ((LST(II).LT.LST(II-1)).AND.(LST(II).LT.100)) THEN
                  IF (FIRST_RECORD.EQ.0) THEN
                     SUBSCAN_START = SUBSCAN_START + 86400.
                  ENDIF
                  OFFSET_MIDNIGHT = 86400.
               ENDIF
            ENDIF
*
            IF (LST(II)+OFFSET_MIDNIGHT .LT. SUBSCAN_START) THEN
               II = II+1
            ELSE
               IF (FIRST_RECORD.EQ.0) THEN
                  DELTA_Y = 0
                  FIRST_RECORD = II
               ELSE
                  DELTA_Y = YCOORD(II)-YCOORD(II-1)
               ENDIF
               IF ( ABS(XCOORD(II)).LE.XLIMIT .AND.
     $         ABS(DELTA_Y).LT.DELTA_INF) THEN
                  II = II+1
               ELSE
                  LAST_RECORD = II-1
               ENDIF
            ENDIF
         ENDDO
         IF (II.GT.NRECORD) LAST_RECORD = NRECORD
*
         IF (FIRST_RECORD.EQ.0 .OR.
     $   LAST_RECORD - FIRST_RECORD.LE.0) THEN
            WRITE(6,*) 'no data in the subscan ',ISUBSCAN
            FIRST_RECORD = 0
            LAST_RECORD = 0
         ENDIF
         SUBSCAN_IFIRST(ISUBSCAN) = FIRST_RECORD
         SUBSCAN_ILAST(ISUBSCAN)  = LAST_RECORD
         DO WHILE ( ABS(XCOORD(II)).GT.XLIMIT .AND. II.LE.NRECORD)
            II = II+1
         ENDDO
      ENDDO
      RETURN
      END
*
*
      SUBROUTINE FIND_SUBSCAN_RAWLIMIT (SUBSCAN_IFIRST,SUBSCAN_ILAST,
     $YCOORD,SINT,NSUBSCAN,NRECORD)
C**************************************************************
C
C     FIND_SUBSCAN_RAWLIMIT() :
C        Identify first record and last record for each subscan
C
C        This routine finds CONTINUOUS subscans (no hole between subscans)
C                           ----------
C
C**************************************************************
*
      INTEGER*4 SUBSCAN_IFIRST(1),SUBSCAN_ILAST(1)
      REAL*4    YCOORD(1)
      REAL*4    SINT
      INTEGER   NSUBSCAN,NRECORD
*
      INTEGER   II,ISUBSCAN
      INTEGER   FIRST_RECORD,LAST_RECORD
      REAL      DELTA_Y,DELTA_INF,MEAN_Y
      INTEGER   N
      CHARACTER STRING*30
*
      II = 1
      DELTA_INF = (SINT*10)/2*1.3        ! delta_inf= sint/2 (in deciarcseconds)
*
      DO ISUBSCAN=1,NSUBSCAN
         FIRST_RECORD = 0
         LAST_RECORD  = 0
*
         DO WHILE ((FIRST_RECORD.EQ.0 .OR. LAST_RECORD.EQ.0)
     $      .AND. II.LE.NRECORD)
*
            IF (FIRST_RECORD.EQ.0) THEN
               FIRST_RECORD = II
               DELTA_Y = 0
               MEAN_Y = YCOORD(II)
               N      = 1
            ELSE
               DELTA_Y = YCOORD(II)-MEAN_Y
               MEAN_Y = (N*MEAN_Y + YCOORD(II))/(N+1)
               N      = N+1
            ENDIF
            IF (ABS(DELTA_Y).LT.DELTA_INF) THEN
               II = II+1
            ELSE
               LAST_RECORD = II-1
            ENDIF
         ENDDO
         IF (II.GT.NRECORD) LAST_RECORD = NRECORD
*
         IF (FIRST_RECORD.EQ.0 .OR.
     $   LAST_RECORD - FIRST_RECORD.LE.0) THEN
            WRITE (STRING,'(a,i4)') 'No data in subscan ', ISUBSCAN
            CALL MESSAGE (2,1,'BASE',STRING)
            FIRST_RECORD = 0
            LAST_RECORD = 0
         ENDIF
         SUBSCAN_IFIRST(ISUBSCAN) = FIRST_RECORD
         SUBSCAN_ILAST(ISUBSCAN)  = LAST_RECORD
      ENDDO
      SUBSCAN_ILAST(NSUBSCAN) = NRECORD  ! all records must be considered
      RETURN
      END
