
      PROGRAM LAZI
C-----------------------------------------------------------------------
C LAZI calculates (u,v) tracks for an interferometer with all telescopes
C located at one site.  This is a modification of T. J. Pearson's program
C HAZI, version 2.1, which is loosely based on Jon Romney's program HAZI.
C
C For HAZI:
C T.J.Pearson,  VAX-11 Fortran
C
C Version 1.0:  1979 May 10 - new program (T.J. Pearson)
C Version 1.1:  1979 Oct 15
C Version 1.2:  1979 Nov  6 - more than 1 source.
C Version 1.3:  1980 Mar 14 - ask questions.
C Version 1.4:  1980 Mar 25 - options NOTICKS, DOTS, GST.
C Version 1.5:  1980 Jun 30 - increase to 15 stations.
C Version 1.6:  1980 Nov 23 - cure bug in tick marks.
C Version 1.7:  1981 Jun 11 - Fortran-77 standard.
C Version 1.8:  1982 Jun  9 - new graphics devices.
C Version 2.0:  1982 Dec 23 - correct bug for invisible sources
C                             or sources circumpolar about S pole;
C                             add parameter NUP; compare names in upper
C                             case; scale in km if frequency not given.
C Version 2.1:  1984 Sep 25 - add FONT and LINEWIDTH parameters
C                             to match UVPLOT (TJP).
C
C For LAZI:
C J. T. Armstrong
C
C Version 1.0:  1989 Feb 23 - Converted to single-site optical use 
C         1.1:  1989 June 1 - Ability to read telescope position file
C-----------------------------------------------------------------------
      CHARACTER*(*) VERSN
      INTEGER       OUTC,INC,INSTA,INSORC,PR,NPARS,MAXSRC,MAXS,MAXB
      REAL*4        IMPOSS
      REAL*8        VELC,SIDR,PI,RPDEG,RPHR
C
      PARAMETER (VERSN='1.1 - 1989 May 12')
      PARAMETER (OUTC=6)      ! terminal output unit
      PARAMETER (INC=5)      ! terminal input unit
      PARAMETER (INSTA=11)      ! unit for station catalog
      PARAMETER (INSORC=12)      ! unit for source catalog
      PARAMETER (PR=8)      ! printer output
      PARAMETER (NPARS=94)      ! number of control parameters
      PARAMETER (MAXSRC=10)      ! maximum number of sources
      PARAMETER (MAXS=15)      ! maximum no of stations
C                        ! NB dimension of FOUND should exceed
C                        ! MAXS and MAXSRC
      PARAMETER (MAXB=MAXS*(MAXS-1)/2) ! maximum no of baselines
      PARAMETER (MAXCOL=4) ! maximum number of colors
      PARAMETER (IMPOSS=1E10)      ! flag for source-not-up
      PARAMETER (VELC=2.99792458D8)      ! velocity of light m/s
      PARAMETER (SIDR=1.0027379093D0) ! sidereal-solar ratio
      PARAMETER (PI=3.14159265359D0)
      PARAMETER (RPDEG=PI/180D0)
      PARAMETER (RPHR=PI/12D0)
C
      LOGICAL PRINT1,PRINT2,PLOT,INVERT,BLID,OK,FOUND(20)
      LOGICAL UP,ISTERM,LIMITS,TICKS,DOTS
C
      CHARACTER*1   LRA(16), LDEC(16)
      CHARACTER*1   ANSWER
      CHARACTER*3   LABEL(MAXB)
      CHARACTER*8   FMT,TIMEC
      CHARACTER*10  SNAME(MAXSRC),XSNAME
      CHARACTER*10  XSTNAME,STNAME(MAXS)
      CHARACTER*165 LEGND1
      CHARACTER*30  LEGND2
      CHARACTER*64  UNITS
      CHARACTER*128 PLTFILE,LSTFILE,POSFILE,SRCFILE
C
      INTEGER   ICOORD(MAXS)
      INTEGER*4 MNAME,DNAME,IS1(MAXB),IS2(MAXB),PGBEGIN
C
      REAL*8    HL,HAL,ZL,CLCD,CLSD,SLCD,SLSD
      REAL*8    CDEC(MAXSRC),SDEC(MAXSRC)
      REAL*8    H(240,MAXS),HA,ZA,AZ
      REAL*8    T(MAXS),P,BZ(MAXB),BE(MAXB),BH(MAXB)
      REAL*8    TCTR,PCTR,RADIUS
      REAL*8    DATE,GMST,GMSTP,GHA
      REAL*8    RR,TT,PP,FREQ,WL6,XCTR,YCTR,ZCTR
      REAL*8    RA(MAXSRC),DEC(MAXSRC),XRA,XDEC
      REAL*8    CHA,DEN,AZX
      REAL*8    CLAT,SLAT,X(MAXS),Y(MAXS),Z(MAXS)
      REAL*8    PARS(NPARS),VALS(NPARS),ENDMRK,BLANK
      REAL*8    EPOS(MAXS),NPOS(MAXS),ELEV(MAXS),XSTELEV
C
      DATA   PARS/ 'PRINT', 'NOTICKS', 'NOINVERT', 'NOBLID',
     1             'TITLE', 3*' ', 'SITE', 9*' ', 'LAT', 'LONG',
     2             'NUP', 'ZENlim', 'HAlim', 'WAVEL', 'UVMAX', 
     3             'DAY', 'MONTH', 'YEAR',
     4             10*'SOU????', 'RA', 'DEC', 'DOTS',
     5             'DT', 2*'GST', 'LINEW???',
     6             'FONT','PLotfile',11*' ','LSTfile',11*' ',
     7             'POSfile',11*' ','SRCfile',11*' ' /
C     
      DATA   VALS/ 1D0, -1D0, -1D0, -1D0,
     1             'LAZI', 3*' ', '***UNSET',9*' ', 0D0, 0D0,
     2             2D0, 50, 5, 500D0, 0D0, 
     3             0D0, 0D0, 0D0, 
     4             10*'  ', 0D0, 0D0, -1D0,
     5             0.1D0, 2*-1D0, 1D0,
     6             1D0, 12*' ',  12*' ',
     7             12*' ',12*' ' /
C
      DATA   ENDMRK/ '/' /,BLANK/' '/,MINUS/'-'/
C
C-----------------------------------------------------------------------
C            Read file names and parameters
C-----------------------------------------------------------------------
C
      WRITE (OUTC,501) ' Control parameters (end with /)'
      WRITE (OUTC,501) ' eg SOURCES (or RA, DEC), WAVEL:'
      MODE=2
      CALL KEYIN(PARS,VALS,NPARS,ENDMRK,MODE,INC,OUTC)
C
      PRINT1 = VALS(1).EQ.0D0 .OR. VALS(1).EQ.1D0
      PRINT2 = VALS(1).EQ.0D0 .OR. VALS(1).EQ.2D0
      TICKS  = VALS(2).EQ.-1D0
      DOTS   = VALS(41).EQ.0D0
      INVERT = VALS(3).EQ.-1D0
      BLID   = VALS(4).EQ.-1D0
      DO 8 J=9,18
          XSTNAME(J-8:)=CHAR(VALS(J))
    8 CONTINUE
C      WRITE (OUTC,'(1X,10A1)') XSTNAME
      TT     = VALS(19)
      PP     = VALS(20)
      MINUP  = VALS(21)
      ZL     = VALS(22)
      HL     = VALS(23)
      WL6    = VALS(24)*1D-3      ! one million wavelengths in m
C
      IF (WL6.LE.0D0) WL6=5D-1  ! 10E6*(default lambda [= 500 nm])=0.5 m
      FREQ=(VELC/WL6)*1D-6       ! frequency in THz
C
      WRITE (OUTC,500) VERSN
      WRITE (PLTFILE, '(12A8)') (VALS(J),J=47,58)
      PLOT = PLTFILE.NE.' '
      IF (PLOT) THEN
            CALL STR$UPCASE(PLTFILE,PLTFILE)
            IER = PGBEGIN(9,PLTFILE,1,1)
            IPLOT = 1
            IF (IER.NE.1) GOTO 900
      END IF
C
      WRITE (LSTFILE, '(12A8)') (VALS(J),J=59,70)
    9 IF (LSTFILE.NE.' ') THEN
            OPEN (UNIT=PR, FILE=LSTFILE, STATUS='NEW',
     1           FORM='FORMATTED',ERR=10)
      ELSE
            OPEN (UNIT=PR, FILE='LAZI.OUT', STATUS='OLD',
     1           FORM='FORMATTED',ERR=10)
      END IF
      GOTO 11
   10 LSTFILE='LAZI.OUT'
      GOTO 9
   11 CONTINUE
C
      WRITE (POSFILE,'(12A8)') (VALS(J), J=71,82)
      WRITE (SRCFILE,'(12A8)') (VALS(J), J=83,94)
C
      UVMAX  = VALS(25)
      DT     = VALS(42)  ! time between points (hours)
      IF(DT.LT.0.1) DT=0.1
      GST1   = VALS(43)
      GST2   = VALS(44)
      LIMITS = GST1.GT.0.0 .AND. GST2.GT.0.0
      IF (GST2.LT.GST1) GST2 = GST2+24.0
      IF (GST2-GST1.GT.24.0) LIMITS = .FALSE.
      LINEWD = VALS(45)
      NFONT  = VALS(46)
      CALL PGSLW(LINEWD)
      CALL PGSCF(NFONT)
C
      WRITE (PR,601) VERSN
      WRITE (PR,6022) FREQ,VALS(24)
C-----------------------------------------------------------------------
C      U.T. - G.S.T. conversion
C      If date not supplied, use today's date
C      DATE = Julian Date
C      GMST = Greenwich Mean S.T. at 0 U.T. (sec)
C-----------------------------------------------------------------------
      IY = VALS(28)
      IM = VALS(27)
      ID = VALS(26)
      IF (VALS(28).LE.0.) CALL IDATE(IM,ID,IY)
      IF(IY.LT.100) IY = 1900 + IY
      CALL TDATE(IY,IM,ID,JD,MNAME,DNAME)
      CALL JULDA(IY,IM,ID,DATE,GMST)
      GMSTP = DMOD(GMST,1D0)
      WRITE(PR,610) IY,MNAME,ID,DNAME,DATE,TIMEC(IDINT(GMST)),GMSTP
      GMST = GMST/3600D0
C-----------------------------------------------------------------------
C      Sources:
C      If RA, DEC not supplied, look up source in catalog
C-----------------------------------------------------------------------
      NSRC = 0
      SNAME(1) = ' '
      DO 1040 ISRC=1,MAXSRC
          IF (VALS(28+ISRC).NE.BLANK) THEN
            NSRC = NSRC+1
            WRITE (SNAME(NSRC),501) VALS(28+ISRC)
            CALL STR$UPCASE(SNAME(ISRC),SNAME(ISRC))
            FOUND(NSRC) = .FALSE.
          END IF
 1040      CONTINUE
C
      IF(VALS(39).NE.0D0 .OR. VALS(40).NE.0D0) THEN
            FOUND(1) = .TRUE.
            RA(1) = VALS(39)
            DEC(1) = VALS(40)
            NSRC = MAX(NSRC,1)
      END IF
      IF(NSRC.EQ.1 .AND. FOUND(1)) GOTO 1200
C
      OPEN (UNIT=INSORC, FILE=SRCFILE,   
     1           FORM='FORMATTED', STATUS='OLD', ERR=1900)
 1100 DO 1150, JCHK=1,NSRC
          CALL GTSRC(INSORC,XSNAME,XRA,XDEC) 
          DO 1150 ISRC=1,NSRC
              IF(SNAME(ISRC).EQ.XSNAME) THEN
                RA(ISRC) = XRA
                DEC(ISRC) = XDEC
                FOUND(ISRC) = .TRUE.
              END IF
 1150     CONTINUE
C 1100      READ(INSORC,511,END=1190,ERR=1900) XSNAME,IHR,IMR,SSR,
C     1           ISIGN,IDD,IMD,SSD
C      DO 1150 ISRC=1,NSRC
C          IF(SNAME(ISRC).EQ.XSNAME) THEN
C            RA(ISRC) = IHR +(IMR/60.0 + SSR/3600.0)
C            DEC(ISRC) = IDD +(IMD/60.0 + SSD/3600.0)
C            IF(ISIGN.EQ.MINUS) DEC(ISRC)=-DEC(ISRC)
C            FOUND(ISRC) = .TRUE.
C          END IF
C 1150      CONTINUE
C      GOTO 1100
 1190      CLOSE (UNIT=INSORC)
C
 1200      OK = .TRUE.
      DO 1250 ISRC=1,NSRC
      IF(FOUND(ISRC)) THEN
            CALL RFORM8(RA(ISRC)*RPHR,LRA,16,3)
            CALL SFORM8(DEC(ISRC)*RPDEG,LDEC,16,3)
            WRITE(PR,640) SNAME(ISRC),LRA,LDEC
            CDEC(ISRC) = DCOS(DEC(ISRC)*RPDEG)
            SDEC(ISRC) = DSIN(DEC(ISRC)*RPDEG)
      ELSE
            WRITE (OUTC,501) ' +++ERROR+++ Source not found'//
     1                  ' in catalog: '//SNAME(ISRC)
            OK = .FALSE.
      END IF
 1250      CONTINUE
C
C-----------------------------------------------------------------------
C            Station coordinates
C-----------------------------------------------------------------------
C
      IF (POSFILE.NE.' ') THEN
          OPEN (UNIT=INSTA, FILE=POSFILE, 
     1           FORM='FORMATTED', STATUS='OLD', ERR=2245)
          WRITE (PR,620)
C
          CALL GTSTA(INSTA,NS,XSTNAME,TT,PP,XSTELEV,0)
C                First time: no. of stations, site name, 
C                latitude, longitude, site elevation
C          WRITE (OUTC,2201) NS,XSTNAME,TT,PP,XSTELEV
 2201     FORMAT(1X,'NS:',I3,,2X,'SITE:',A10,2X,
     1             'LAT:',F6.2,2X,'LONG:',F6.2,2X,'ELEV:',F6.2)
          DO 2200, I=1,NS
C               WRITE (OUTC,*) I
               CALL GTSTA(INSTA,ICOORD(I),STNAME(I),
     1                     EPOS(I),NPOS(I),ELEV(I),0)
C                    Second and subsequent times: polar-or-rectangular
C                    flag, station name, position, position, elevation
C               WRITE (OUTC,2202) ICOORD(I),STNAME(I),
C     1                     EPOS(I),NPOS(I),ELEV(I)
 2202     FORMAT(1X,'ICOORD:',I3,2X,'STA:',A10,2X,'EAST:',F5.0,2X,
     1              'NORTH:',F6.2,2X,'ELEV:',F6.2)
               IF (ICOORD(I).EQ.1) THEN    ! convert polar to x-y
                    RADIUS=EPOS(I)
                    EPOS(I)=-RADIUS*SIN(NPOS(I)*RPDEG)
                    NPOS(I)=RADIUS*COS(NPOS(I)*RPDEG)
               END IF
 2200     CONTINUE
          CALL STR$UPCASE(XSTNAME,XSTNAME)
 2240     CLOSE(UNIT=INSTA)
C
      ELSE
C
 2245     WRITE(OUTC,501) 
     1       '$Polar [P] or X-Y [X] telescope positions? (default: X): '
          READ(INC,501) ANSWER
          IF (ANSWER.EQ.'P'.OR.ANSWER.EQ.'p') THEN
              WRITE(OUTC,502) '$R in meters, theta CCW from north.'
 2250         WRITE(OUTC,502) '$Name, R, theta (^Z to stop): '
              READ(INC,*,ERR=2250,END=2255) 
     1            (STNAME(I),EPOS(I),NPOS(I),I=1,MAXS)
 2255         NS=I-1
              DO 2259 I=1,NS
                  RADIUS=EPOS(I)
                  EPOS(I)=-RADIUS*SIN(NPOS(I)*RPDEG)
                  NPOS(I)=RADIUS*COS(NPOS(I)*RPDEG)
 2259         CONTINUE
          ELSE
              WRITE(OUTC,502) 'X, Y in meters east and north.'
 2260         WRITE(OUTC,502) '$Name, X, Y (^Z to stop): '
              READ(INC,*,ERR=2260,END=2265) 
     1            (STNAME(I),EPOS(I),NPOS(I),I=1,MAXS)
 2265         NS=I-1
          END IF
      END IF
C
C-----------------------------------------------------------------------------
C             coordinates of array center
C-----------------------------------------------------------------------------
      IF(HL.EQ.0.0) HL=12.0
      HL = ABS(HL)
      IF (ZL.EQ.0.0) ZL=90.0
      WRITE (OUTC,621) XSTNAME,RR,TT,PP
      WRITE (PR,621) XSTNAME,RR,TT,PP
      IF (RR.LT.6D6) CALL GEOID(RR,TT)
      SLAT = DSIN(TT*RPDEG)
      CLAT = DCOS(TT*RPDEG)
      TCTR = TT
      PCTR = PP/15D0                          ! Longitude in HOURS
      XCTR = (RR+XSTELEV)*CLAT*DCOS(PP*RPDEG)
      YCTR = (RR+XSTELEV)*CLAT*DSIN(PP*RPDEG)
      ZCTR = (RR+XSTELEV)*SLAT
      WRITE (PR,623) (RR+XSTELEV),TT,PP,XCTR,YCTR,ZCTR,HL,ZL
      WRITE (OUTC,623) (RR+XSTELEV),TT,PP,XCTR,YCTR,ZCTR,HL,ZL
C
C
C-----------------------------------------------------------------------
C            Coordinates of telescopes
C
C            Coordinate system has X pointing outward from center of
C            Earth to 0 deg longitude, 0 deg latitude; Y pointing to
C            90 deg west longitude, 0 deg latitude; Z pointing outward 
C            from North Pole.  A left-handed coordinate system!
C-----------------------------------------------------------------------
      WRITE(PR,2499)
      WRITE(OUTC,2499)
 2499 FORMAT(' Telescope positions:')
      DO 2500 I=1,NS
          X(I)=(RR+ELEV(I))*DCOS(PP*RPDEG)*DCOS(TT*RPDEG)
     1           +EPOS(I)*DSIN(PP*RPDEG)
     2           -NPOS(I)*DCOS(PP*RPDEG)*DSIN(TT*RPDEG)
          Y(I)=(RR+ELEV(I))*DSIN(PP*RPDEG)*DCOS(TT*RPDEG)
     1           -EPOS(I)*DCOS(PP*RPDEG)
     2           -NPOS(I)*DSIN(PP*RPDEG)*DSIN(TT*RPDEG)
          Z(I)=(RR+ELEV(I))*DSIN(TT*RPDEG)+NPOS(I)*DCOS(TT*RPDEG)
          WRITE(OUTC,2501) I,X(I),Y(I),Z(I)
          WRITE(PR,2501) I,X(I),Y(I),Z(I)
2501  FORMAT(1X,I2,1X,'X= ',F15.6,1X,'Y= ',F15.6,1X,'Z= ',F15.6) 
2500  CONTINUE               
      WRITE (OUTC,501) ' Telescope positions calculated'
C-----------------------------------------------------------------------
C            Baselines
C-----------------------------------------------------------------------
      J = 0
      WRITE(PR,630)
      DO 3100 M2=2,NS
      NS1 = M2-1
      DO 3100 M1=1,NS1
            J = J+1
            IS1(J) = M1
            IS2(J) = M2
            BX = (X(M2)-X(M1))/WL6
            BY = (Y(M2)-Y(M1))/WL6
            BZ(J) = (Z(M2)-Z(M1))/WL6
            BE(J) = SQRT(BX**2+BY**2)
            BH(J) = ATAN2(BY,BX)/RPHR
            WRITE(PR,631) J,STNAME(M1),STNAME(M2),
     1             BX,BY,BZ(J),BE(J),BH(J)
            WRITE(OUTC,631) J,STNAME(M1),STNAME(M2),
     1             BX,BY,BZ(J),BE(J),BH(J)
            FMT = '(I2,1X)'
            IF (J.GT.9) FMT='(I3)'
            WRITE(LABEL(J),FMT) J
C           ENCODE(3,FMT,LABEL(J)) J
 3100      CONTINUE
      LJ=J
C
C-----------------------------------------------------------------------
C            Begin loop for sources
C-----------------------------------------------------------------------
C
      DO 8000 ISRC=1,NSRC
C
C            Determine hour-angle limits
C
      CLCD = CLAT*CDEC(ISRC)
      CLSD = CLAT*SDEC(ISRC)
      SLCD = SLAT*CDEC(ISRC)
      SLSD = SLAT*SDEC(ISRC)
      COSZAL   = DCOS(ZL*RPDEG)
      COSHAL   = (COSZAL-SLSD)/CLCD
      IF (COSHAL.GE.1.0) THEN
C        ! Source is never visible
         HAL = -0.01
         WRITE (PR,501) ' ++WARNING++ Source '//SNAME(ISRC)//
     1                  ' cannot be seen from station '//XSTNAME
         WRITE (OUTC,501) ' ++WARNING++ Source '//SNAME(ISRC)//
     1                  ' cannot be seen from station '//XSTNAME
         ELSE IF (COSHAL.LE.-1.0 .AND. HL.GE.12.0) THEN
C        ! Source is circumpolar
             HAL = 12.0
         ELSE IF (COSHAL.LE.-1.0) THEN
C        ! Circumpolar, but HA limit is effective
             HAL = HL
         ELSE
C        ! ZA limit is effective; HA limit may be
             HAL = ACOS(COSHAL)/RPDEG/15.0
             HAL = MIN(HAL,HL)
         END IF
C
      GST = DMOD(12.0+RA(ISRC)-DT,24D0)
      IF (GST.LT.0.0) GST = GST+24.0
      IF (LIMITS) GST = GST1      ! supplied start time
      GST0 = IFIX(GST/DT)*DT
C
C-----------------------------------------------------------------------
C            Print Section 1:  HA, ZA, AZ
C-----------------------------------------------------------------------
C
      LK = 24.0/DT
      IF (LIMITS) LK = (GST2-GST1)/DT      ! supplied time range
      NPAGE = NS/5
      IF(NS.GT.NPAGE*5) NPAGE=NPAGE+1
      IF(.NOT.PRINT1) NPAGE = 1
      DO 5100 IPAGE=1,NPAGE
      M1 = (IPAGE-1)*5 + 1
      M2 = MIN0(M1+4,NS)
      NSP = M2 - M1 + 1
      NL = 0
      DO 5100 K=1,LK
          GST = GST0 + K*DT
          IF(GST.GE.24.0) GST=GST-24.0
          GHA = GST-RA(ISRC)      ! Greenwich Hour Angle of source
          NUP = 0
C
C            Compute HA, ZA, AZ 
C
        HA = DMOD(GHA-PCTR,24D0)
        IF (HA.LT.-12.0) HA = HA+24.0
        IF (HA.GT.12.0)  HA = HA-24.0
        CHA = COS(HA*RPHR)
        IF (ABS(HA).LE.HAL) THEN
            ZA = ACOS(SLSD+CLCD*CHA)/RPDEG
            DEN = CDEC(ISRC)*DSIN(HA*RPHR)
            AZX = ATAN2(CLSD-SLCD*CHA, DEN)/RPDEG
            AZ = DMOD(AZX-90.0,360D0)
            IF (AZ.LT.0.0) AZ=AZ+360.0
            IF(ABS(ZA).LT.ZL) THEN
                UP = .TRUE.
            ELSE
                UP = .FALSE.
            END IF
        ELSE
            HA = IMPOSS
            ZA = IMPOSS
            AZ = IMPOSS
            UP = .FALSE.
        END IF
C
C            Compute IHA for each baseline
C
        DO 5101 IB=1,LJ
            IF (UP) THEN
                H(K,IB) = DMOD(GHA-BH(IB)-6.0, 24D0)
                IF (H(K,IB).LT.-12.0) H(K,IB) = H(K,IB)+24.0
                IF (H(K,IB).GT.12.0)  H(K,IB) = H(K,IB)-24.0
            ELSE
                H(K,IB) = IMPOSS
            END IF
 5101   CONTINUE
C
C            Print section 1: HA, ZA, AZ 
C
        IF (PRINT1 .AND. UP) THEN
          UT = DMOD((GST-GMST)/SIDR,24D0)
          IF (UT.LT.0.0) UT = UT+24.0
          NL = NL+1
          IF (MOD(NL,50).EQ.1) THEN
              WRITE (PR,91)  SNAME(ISRC)
              WRITE (PR,651) XSTNAME
              WRITE (PR,652)
          END IF
          IF (MOD(NL,10).EQ.1) WRITE(PR,90)
          WRITE (PR,655,IOSTAT=JUNK) UT,GST,
     1                  HA,ZA,AZ
        END IF
C
 5100   CONTINUE
C
C-----------------------------------------------------------------------
C            Print Section 2: IHA
C-----------------------------------------------------------------------
      IF (PRINT2) THEN
          NPAGE = LJ/15
          IF (LJ.GT.NPAGE*15) NPAGE = NPAGE + 1
          DO 6500 IPAGE = 1,NPAGE
            J1 = (IPAGE-1)*15 + 1
            J2 = MIN(J1+14,LJ)
            LJP = J2 - J1 + 1
            NL = 0
            DO 6400 K=1,LK
              DO 6400 IB=1,LJ
                IF (H(K,IB).NE.IMPOSS) GOTO 6200
                GOTO 6400
C
 6200           GST = GST0 + K*DT
                IF (GST.GT.24.0) GST = GST-24.0
                NL = NL+1
                IF (MOD(NL,50).EQ.1) THEN
                  WRITE (PR,91) SNAME(ISRC)
                  WRITE (PR,851) (IB,IB=J1,J2)
                END IF
                IF (MOD(NL,5).EQ.1) WRITE(PR,90)
                WRITE (PR,855,ERR=6400) GST,H(K,IB)
 6400            CONTINUE
 6500          CONTINUE
      END IF
C-----------------------------------------------------------------------
C            Plot section
C-----------------------------------------------------------------------
C
      IF (.NOT.PLOT) GOTO 8000
      IF (UVMAX.NE.0.0) THEN
          UMAX = UVMAX
          VMAX = UVMAX
      ELSE
          UMAX = 0.0
          VMAX = 0.0
          DO 7100 IB=1,LJ
            UMAX = AMAX1(UMAX,BE(IB))
            VMAX = AMAX1(VMAX, ABS(BZ(IB)*CDEC(ISRC)) +
     1                         ABS(BE(IB)*SDEC(ISRC)) )
 7100          CONTINUE
          UMAX = 1.2*UMAX
          VMAX = 1.2*VMAX
      END IF
C
C WRITE STATEMENTS TO CHECK TO SEE IF THE PROPER VALUES EXIST AT THIS
C POINT:
C
      WRITE(OUTC,501) UNITS 
C
      LEGND1(1:) = XSTNAME
      LEGND1(15:) = SNAME(ISRC)
      WRITE (LEGND2,702) DEC(ISRC)
      IF (LIMITS) WRITE (LEGND2(15:),703) GST1,GST2
      CALL PGSCI(5)
      CALL PGENV(UMAX,-UMAX,-VMAX,VMAX,1,1)
      CALL PGLABEL(UNITS,' ',LEGND1)
      DO 7200 J=1,LJ
            U1 = BE(J)
            V0 = BZ(J)*CDEC(ISRC)
            V1 = BE(J)*SDEC(ISRC)
            IPP = 2
            CALL PGSCI(3)
      DO 7300 K=1,LK
            HAI = H(K,J)
            IF(HAI.EQ.IMPOSS) THEN
                  IPP = 2
                  GOTO 7300
            END IF
            HAIR = HAI*RPHR
            U = U1*COS(HAIR)
            V = V0 + V1*SIN(HAIR)
            IF(IPP.EQ.2 .AND. BLID) CALL GRCHAR(IPLOT,.TRUE.,0.,
     1                  .FALSE.,U,V,LABEL(J))
            IF (DOTS) IPP=3
            CALL GRVECT(IPLOT,IPP,.FALSE.,1,U,V)
            IPP = 1
            IF(TICKS .AND. ABS(HAI-NINT(HAI)).LT.DT/2.0) THEN
                HAIR = FLOAT(NINT(HAI))*RPHR
                UTIK = U1*COS(HAIR)
                VTIK = V0 + V1*SIN(HAIR)
C                CALL GRMARK(IPLOT,.TRUE.,'03'X,.FALSE.,1,UTIK,VTIK)
C                IF(INVERT) CALL GRMARK(IPLOT,.TRUE.,'03'X,
C     1                  .FALSE.,1, -UTIK, -VTIK)
                CALL GRMARK(IPLOT,.TRUE.,'03',.FALSE.,1,UTIK,VTIK)
                IF(INVERT) CALL GRMARK(IPLOT,.TRUE.,'03',
     1                  .FALSE.,1, -UTIK, -VTIK)
                CALL GRVECT(IPLOT,2,.FALSE.,1,U,V)
                ENDIF
 7300      CONTINUE
      IF (.NOT.INVERT) GOTO 7200
      DO 7310 K=1,LK
            HAI = H(K,J)
            IF(HAI.EQ.IMPOSS) GOTO 7310
            HAIR = HAI*RPHR
            U = U1*COS(HAIR)
            V = V0 + V1*SIN(HAIR)
            CALL GRVECT(IPLOT,3,.FALSE.,1,-U,-V)
 7310    CONTINUE
 7200 CONTINUE
      CALL GRTERM
C
 8000 CONTINUE
C
      IF(PLOT) CALL PGEND
      CLOSE (UNIT=PR)
      CALL EXIT
C-----------------------------------------------------------------------
C            Error exits
C-----------------------------------------------------------------------
  900 CALL ERROR('Cannot open PLTFILE')
 1900 CALL ERROR('Cannot read SOURCES catalog')
 2900 CALL ERROR('Cannot read STATIONS catalog')
C-----------------------------------------------------------------------
C            Format statements
C-----------------------------------------------------------------------
  500 FORMAT(' LAZI computes u-v tracks for a single-site 
     1       interferometer'/' (Version ',A,')'/)
  501 FORMAT(A)
  502 FORMAT((A)/)
  601 FORMAT(' Program LAZI'/' (Version ',A,')'/)
 6021 FORMAT('(million wavelengths,',I6,' MHz)')
 6022 FORMAT(/5X,'Frequency',F11.3,' THz'/
     1       5X,'Wavelength',F10.3,' nm')
  610 FORMAT(//5X,'Date: ',I4,1X,A3,I3,' (',A3,')'/
     1       5X,'Julian date: ',F12.1/
     2       5X,'Greenwich Mean S.T. at 0h U.T.: ',A8,F4.3)
  511 FORMAT(2X,A8,5X,2I2,F6.3,A1,2I2,F6.3)
  640 FORMAT(/5X,'Source: ',A8,5X,'RA:',16A1,5X,'Dec:',16A1)
  521 FORMAT(A8,2X,F10.0,2(F16.0,F3.0,F6.3),2F5.0)
  620 FORMAT(//' Stations:'//4X,'Name',9X,'Rad./elev.',
     1       6X,'Latitude',17X,'Longitude',19X,'X',12X,'Y',12X,'Z',
     2       7X,'HA/ZA limits')
  621 FORMAT(/2X,A10/3X,'Input',10X,F9.1,2(F14.6,F4.0,F6.2))
  623 FORMAT(3X,'Converted',6X,F9.1,F14.6,10X,F14.6,/,13X,3F13.1,
     1       T121,2F6.1)
  630 FORMAT(//' Baselines (millions of wavelengths):'
     1       //3X,'#',8X,'Stations',11X,
     2       'Bx',9X,'By',9X,'Bz',9X,'Beq',5X,'GHA (h)')
  631 FORMAT(1X,I3,2X,A8,'-',A8,2X,3F11.6,2X,2F10.6)
  651 FORMAT(//3X,'UT',5X,'GST',12X,A8)
  652 FORMAT(22X,'HA',4X,'ZA',5X,'AZ')
  655 FORMAT(1X,F5.1,2X,F5.1,7X,F5.1,2X,F4.1,2X,F5.1)
  851 FORMAT('   GST',10X,'Interferometer Hour Angles'
     1       //10X,/(I4,3X))
  855 FORMAT(1X,F5.1,4X,/(F5.1,2X))
   90 FORMAT(1X)
   91 FORMAT(1H1,10X,A8)
  701 FORMAT(A8)
  702 FORMAT('DEC =',F6.1)
  703 FORMAT('GST =',F5.1,',',F5.1)
C-----------------------------------------------------------------------
      END
C
C=======================================================================
C
      SUBROUTINE STR$UPCASE ( A, B )
      CHARACTER *(*) A, B
      INTEGER *4 I, K
      K = LEN(A)
      DO 100 I = 1, K
          IF ( ( A(I:I) .LE. 'z' ) .AND. ( A(I:I) .GE. 'a' ) ) THEN
              B(I:I) = CHAR( ICHAR ( A(I:I) ) - 32 )
          ELSE
              B(I:I) = A(I:I)
          END IF
  100 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE GTSTA(INSTA,ITYPE,STANAME,P1,P2,P3,OUTSTA)
C-----------------------------------------------------------------------
C This is a modification of subroutine GTMOD in VISPLOT.
C
C GTMOD:
C T.J.Pearson   1979 FEB 18
C Modified for VAX version  1979 JULY 5
C Modified to allow comments 1981 February 9
C Modified to copy input comments to an output file 1982 September 22
C
C GTSTA:
C Adapted to read a list of optical interferometer siderostat positions
C for LAZI:
C J.T.Armstrong  1989 May 13
C
C Read one model component (card) from unit 'INSTA'.
C Each card carries a character string and up to 4 numbers.
C The first card carries the number of stations (ITYPE), the
C site name (STANAME), the site elevation (P1), latitude (P2),
C and longitude (P3).  The subsequent cards carry
C the coordinate system type ITYPE, the station name STANAME,  
C and 3 floating-point parameters P1-P3 representing the station 
C coordinates.
C      ITYPE = 1 => P1, P2 are distance and position angle from
C                   array center (PA = CCW from north)
C      ITYPE = 2 => P1, P2 are distance east and north from array
C                   center
C      P3 = elevation (default = 0)
C Copy comments following a '!' character to unit 'OUTSTA'
C Ignore zero-flux components (or blank lines)
C-----------------------------------------------------------------------
      REAL*8       PARS(3),P1,P2,P3
      INTEGER      INSTA,OUTSTA,ITYPE
      INTEGER*2    CARD(80),BLANK,EXCL
      CHARACTER*1  ANSWER
      CHARACTER*10 STANAME
      CHARACTER*64 FNAME
C
      DATA BLANK/'  '/
      DATA EXCL/'!'/
C
C Trim off trailing blanks:
   5  READ(INSTA, '(80A1)', END=50) CARD
      IF(CARD(1).EQ.EXCL) GOTO 5
      L = 80
      DO WHILE (CARD(L).EQ.BLANK)
            L = L-1
            IF (L.EQ.0) GOTO 5
      END DO
C 
      I = 1
C Read type or number of stations:
      CALL SKIPBL(CARD,L,I)
      CALL CTOI2(CARD,L,I,ITYPE)
C Read site or station name: 
      J = 1
      CALL SKIPBL(CARD,L,I)
      CALL CTOSTR(CARD,L,I,STANAME)
C      DO WHILE (CARD(I).NE.BLANK.AND.CARD(I).NE.EXCL)
C           STANAME(J:) = CHAR(CARD(I))
C           I = I+1
C           J = J+1
C      END DO
      IF(CARD(I).EQ.EXCL) GOTO 50
C Read parameters: 
      P3 = 0.
      DO 30 J=1,3
          CALL SKIPBL(CARD,L,I)
          CALL CTOR2(CARD,L,I,PARS(J))
   30 CONTINUE
      IF (I.GT.L) GOTO 40
      IF (CARD(I).EQ.BLANK) GOTO 40
      IF (CARD(I).EQ.EXCL) THEN
            IF (OUTSTA.GT.0) WRITE (OUTSTA,'(80A1)') (CARD(J),J=I,L)
            GOTO 40
      END IF
      INQUIRE(UNIT=INSTA,NAME=FNAME)
      CALL ERROR('Invalid model format: file '//FNAME)
   40 CONTINUE
      IF (ITYPE.EQ.0) GOTO 5
      P1 = PARS(1)
      P2 = PARS(2)
      P3 = PARS(3)
   50 RETURN 
      END
C
C==========================================================================

        SUBROUTINE GTSRC(INSORC,SNAME,RA,DEC)
C
C                   J.T.ARMSTRONG JUNE 1 1989  NDP FORTRAN
C
C                   READ SOURCE NAME, RA, DEC FROM A FILE
C
        REAL*8       RA,DEC
        CHARACTER*10 SNAME
        INTEGER      INSORC
        INTEGER*2    CARD(80)
C
        DATA BLANK/' '/
        DATA EXCL/'!'/
C Trim off trailing blanks
   5  READ(INSORC, '(80A1)', END=50) CARD
      IF(CARD(I).EQ.EXCL) GOTO 5
      L=80
      DO WHILE (CARD(L).EQ.BLANK)
            L = L-1
            IF (L.EQ.0) GOTO 5
      END DO
C
      I=1
C Read source name
      CALL SKIPBL(CARD,L,I)
      CALL CTOSTR(CARD,L,I,SNAME)
      WRITE(6,1000) SNAME
 1000 FORMAT(' GTSRC: SNAME: ',A8)
      IF (I.GT.L) GOTO 50
C Read RA
      CALL SKIPBL(CARD,L,I)
      CALL CTOR2(CARD,L,I,RA)
      WRITE(6,1001) I,RA
 1001 FORMAT(' GTSRC: I, RA = 'I3,2X,F20.10)
      IF (I.GT.L) GOTO 50
C Read DEC
      CALL SKIPBL(CARD,L,I)
      CALL CTOR2(CARD,L,I,DEC)
      WRITE(6,1002) I,DEC
 1002 FORMAT(' GTSRC: I,DEC = 'I3,2X,F20.10)
   50 CONTINUE
      RETURN
      END
C=======================================================================
        SUBROUTINE CTOR2(S,N,I,R)
C
C               T.J.PEARSON  1978 MAY 11
C               MODIFIED 1979 FEB 6  VAX/VMS FORTRAN 4
C               MODIFIED 1989 MAY 26  CONVERT FR. FUNCTION TO SUBROUTINE
C
C               EVALUATE NUMBER AT S(I); INCREMENT I
C               RETURN AS REAL*8
C
        REAL*8    R,RR,FRACTN
        INTEGER   IEXPON,SIGN
        INTEGER   N,I,IEXPON,ICOLON
        INTEGER*2 S(N),COLON,PLUS,MINUS,POINT,E,ELC
        CHARACTER*2 CCOLON,CPLUS,CMINUS,CPOINT,CE,CELC
        EQUIVALENCE (CCOLON,COLON), (CPLUS,PLUS), (CMINUS,MINUS), 
     +              (CPOINT,POINT), (CE,E), (CELC,ELC)
        DATA      CCOLON/':'/,CPLUS/'+'/,CMINUS/'-'/,CPOINT/'.'/,
     1            CE/'E'/,CELC/'e'/
        DATA      TEN/10.0/
C
        R=0D0
        IF(I.GT.N) GOTO 270
        SIGN=+1
        IF(S(I).NE.MINUS) GOTO 130
            I=I+1
            SIGN=-1
            CALL SKIPBL(S,N,I)
            GOTO 140
  130   IF(S(I).NE.PLUS) GOTO 140
            I=I+1
            CALL SKIPBL(S,N,I)
  140   ISTART=I
        ICOLON=0
  150   CALL CTOD2(S,N,I,RR)
        R = RR*SIGN/(60D0**ICOLON) + R
        IF(I.GT.N) GOTO 270
        IF(S(I).NE.COLON) GOTO 160
            I=I+1
            ICOLON=ICOLON+1
            GOTO 150
C
C               LOOK FOR FRACTIONAL PART, IF ANY
C
  160   IF(S(I).NE.POINT) GOTO 170
        J=I+1
        CALL CTOD2(S,N,J,RR)
        FRACTN=RR*SIGN
        R=R + FRACTN*10D0**(I-J+1)
        I=J
        IF(I.GT.N) GOTO 270
C
C               Look for exponent, if any; an exponent may not stand
C               alone
C
  170   IF(.NOT.(S(I).EQ.E.OR.S(I).EQ.ELC) .OR. I.EQ.ISTART) GOTO 270
        I=I+1
        IF(I.GT.N) GOTO 270
        SIGN=+1
        IF(S(I).NE.MINUS) GOTO 180
            I=I+1
            SIGN=-1
            GOTO 190
  180   IF(S(I).NE.PLUS) GOTO 190
            I=I+1
  190   CALL CTOI2(S,N,I,IEXPON)
        R=R * 10D0**(SIGN*IEXPON)
  270   CONTINUE
        RETURN
        END
C=======================================================================
        SUBROUTINE CTOD2(S,N,I,R)
C
C               T.J. PEARSON 1979 JUN 4  VAX/VMS FORTRAN 4
C               MODIFIED 1989 MAY 26  CONVERT FR. FUNC. TO SUBROUTINE
C               EVALUATE INTEGER AT S(I); INCREMENT I
C               RETURN AS REAL*8
C
        INTEGER   N,I
        REAL*8    R
        INTEGER*2 S(N),DIGITS(10)
        CHARACTER*2 CDIGIT(10)
        EQUIVALENCE (CDIGIT(1),DIGITS(1))
        DATA      CDIGIT/'0','1','2','3','4','5','6','7','8','9'/
C
        R=0D0
   10   IF(I.GT.N) RETURN
        DO 20 K=1,10
        IF(S(I).EQ.DIGITS(K)) GOTO 30
   20   CONTINUE
        RETURN
   30   R=R*10D0 + (K-1)
        I=I+1
        GOTO 10
        END
C=======================================================================

        SUBROUTINE CTOI2(S,N,I,J)
C               TJP 1977 DECEMBER 4  CIT 370/158
C               MODIFIED 1979 FEB 5  VAX/VMS FORTRAN 4
C               MODIFIED 1989 MAY 26 CONV. FR. FUNCT. TO SUBROUTINE
C
C               EVALUATE INTEGER AT S(I); INCREMENT I
C
        INTEGER   J,N,I
        INTEGER*2 S(N),DIGITS(10)
        CHARACTER*2 CDIGIT(10)
        EQUIVALENCE (CDIGIT(1),DIGITS(1))
        DATA      CDIGIT/'0','1','2','3','4','5','6','7','8','9'/
C
        J=0
   10   IF(I.GT.N) RETURN
        DO 20 K=1,10
        IF(S(I).EQ.DIGITS(K)) GOTO 30
   20   CONTINUE
        RETURN
   30   J=J*10+K-1
        I=I+1
        GOTO 10
        END
C========================================================================

        SUBROUTINE CTOSTR(S,N,I,CSTR)
C                 J.T.ARMSTRONG JUNE 1 1989  NDP FORTRAN
C
C                 READ CHARACTERS INTO A 10-CHARACTER STRING, 
C                 INCREMENT I
C
        INTEGER       N,I
        INTEGER*2     S(N),BLANK
        CHARACTER*10  CSTR
C
        DATA BLANK/'  '/
C
        STRING='          '
        IF(I.GT.N) RETURN
        J=1
        DO WHILE (S(I).NE.BLANK)
             CSTR(J:) = CHAR(S(I))
             I=I+1
             J=J+1
        END DO
        WRITE(6,1000) CSTR
 1000   FORMAT(' CTOSTR: CSTR: ',A10)
        RETURN
        END



