      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 use 
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.0 - 1989 Feb 23')
      PARAMETER (OUTC=6)      ! terminal output unit
      PARAMETER (INC=5)      ! terminal input unit
      PARAMETER (INSTA=90)      ! unit for station catalog
      PARAMETER (INSORC=90)      ! unit for source catalog
      PARAMETER (PR=8)      ! printer output
      PARAMETER (NPARS=45)      ! 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 (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(MAXS),ISTERM,LIMITS,TICKS,DOTS
C
      CHARACTER*1   ANSWER
      CHARACTER*3   LABEL(MAXB)
      CHARACTER*8   FMT,TIMEC
      CHARACTER*8   SNAME(MAXSRC),XSNAME,STNAME(MAXS),XSTNAME
      CHARACTER*165 LEGND1
      CHARACTER*30  LEGND2
      CHARACTER*64  UNITS
      CHARACTER*128 PLTFIL,LSTFIL,POSFIL
C
      INTEGER*4 MNAME,DNAME,IS1(MAXB),IS2(MAXB),PGBEGIN
      INTEGER*2 LRA(16),LDEC(16)
C
      REAL*4    TCTR,PCTR
      REAL*4    HL(MAXS),HAL(MAXS),ZL(MAXS),CLCD(MAXS),CLSD(MAXS)
      REAL*4    SLCD(MAXS),XPOS(MAXS),YPOS(MAXS)
      REAL*4    SLSD(MAXS),H(MAXB,240),HA(MAXS),ZA(MAXS),AZ(MAXS)
      REAL*4    T(MAXS),P(MAXS),BZ(MAXB),BE(MAXB),BH(MAXB)
      REAL*4    CDEC(MAXSRC),SDEC(MAXSRC)
C
      REAL*8    DATE,GMST,GMSTP,RR,TT,PP,FREQ,WL,XCTR,YCTR,ZCTR
      REAL*8    RA(MAXSRC),DEC(MAXSRC)
      REAL*8    CLAT(MAXS),SLAT(MAXS),X(MAXS),Y(MAXS),Z(MAXS)
      REAL*8    PARS(NPARS),VALS(NPARS),ENDMRK,BLANK
C
      DATA   PARS/ 'PRINT', 'NOTICKS', 'NOINVERT', 'NOBLID',
     1                 'TITLE',3*' ','NUP','FREQ', 'UVMAX', 'DATE',
     2                 10*'SOU????', 'RA', 'DEC', 'DOTS',
     3                 15*'STA?????', 'DT', 2*'GST', 'LINEW???',
     4                   'FONT' /
      DATA   VALS/ -1D0, -1D0, -1D0, -1D0,
     1                 'HAZI',3*' ',2D0, 0D0, 0D0, 0D0,
     2                 10*'  ', 0D0, 0D0, -1D0,
     3                 2*'***UNSET',13*' ', 0.1D0,2*-1D0, 1D0,
     4                   1D0 /
      DATA   ENDMRK/ '/' /,BLANK/' '/,MINUS/'-'/
C
C-----------------------------------------------------------------------
C            Read file names and parameters
C-----------------------------------------------------------------------
C
      WRITE (OUTC,500) VERSN
    5      WRITE (OUTC,501) '$Graphics device (default: no plot): '
      READ (INC,501) PLTFIL
      PLOT = PLTFIL.NE.' '
      IF (PLOT) THEN
            CALL STR$UPCASE(PLTFIL,PLTFIL)
            IER = PGBEGIN(9,PLTFIL,1,1)
            IPLOT = 1
            IF (IER.NE.1) GOTO 5
      END IF
C
   10      WRITE (OUTC,501) '$File name for listing (default: terminal): '
      READ (INC,501) LSTFIL
      IF (LSTFIL.EQ.' ') LSTFIL = 'SYS$OUTPUT'
      OPEN (UNIT=PR, NAME=LSTFIL, STATUS='NEW',FORM='FORMATTED',ERR=10)
C
      WRITE (OUTC,501) ' Control parameters (end with /)'
      WRITE (OUTC,501) ' eg SOURCES (or RA, DEC), STATIONS, FREQ:'
      MODE = 1
      IF (ISTERM('SYS$INPUT')) 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(25).EQ.0D0
      INVERT = VALS(3).EQ.-1D0
      BLID   = VALS(4).EQ.-1D0
      MINUP  = VALS(9)
      FREQ   = VALS(10)      ! frequency/MHz
      IF (FREQ.LE.0D0) THEN
          FREQ = 0D0
          WL = 1000            ! 1 km/m
          UNITS = '(kilometers)'
      ELSE
          WL = VELC/FREQ       ! 1 million wavelengths /m
          WRITE (UNITS,6021) NINT(FREQ)
      END IF
      UVMAX  = VALS(11)
      DT     = VALS(41)  ! time between points (hours)
      IF(DT.LT.0.1) DT=0.1
      GST1   = VALS(42)
      GST2   = VALS(43)
      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(44)
        NFONT  = VALS(45)
      CALL PGSLW(LINEWD)
      CALL PGSCF(NFONT)
C
      WRITE (PR,601) VERSN
      IF (FREQ.GT.0.0) WRITE (PR,6022) FREQ,WL
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-----------------------------------------------------------------------
      MDATE = VALS(12)
      IY = MOD(MDATE/3600,100)
      IM = MOD(MDATE,3600)/60
      ID = MOD(MDATE,60)
      IF (MDATE.LE.0) CALL IDATE(IM,ID,IY)
      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(12+ISRC).NE.BLANK) THEN
            NSRC = NSRC+1
            WRITE (SNAME(NSRC),501) VALS(12+ISRC)
            CALL STR$UPCASE(SNAME(ISRC),SNAME(ISRC))
            FOUND(NSRC) = .FALSE.
          END IF
 1040      CONTINUE
C
      IF(VALS(23).NE.0D0 .OR. VALS(24).NE.0D0) THEN
            FOUND(1) = .TRUE.
            RA(1) = VALS(23)/3600.0
            DEC(1) = VALS(24)/3600.0
            NSRC = MAX(NSRC,1)
      END IF
      IF(NSRC.EQ.1 .AND. FOUND(1)) GOTO 1200
C
      OPEN(UNIT=INSORC, NAME='SOURCES', READONLY,
     1           FORM='FORMATTED', TYPE='OLD', ERR=1900)
 1100      READ(INSORC,511,END=1800,ERR=1900) XSNAME,IHR,IMR,SSR,
     1           ISIGN,IDD,IMD,SSD
      DO 1150 ISRC=1,NSRC
          IF(SNAME(ISRC).EQ.XSNAME) THEN
            RA(ISRC) = IHR +(IMR/60.0 + SSR/3600.0)
            DEC(ISRC) = IDD +(IMD/60.0 + SSD/3600.0)
            IF(ISIGN.EQ.MINUS) DEC(ISRC)=-DEC(ISRC)
            FOUND(ISRC) = .TRUE.
          END IF
 1150      CONTINUE
      GOTO 1100
 1800      CLOSE (UNIT=INSORC)
C
 1200      OK = .TRUE.
      DO 1250 ISRC=1,NSRC
      IF(FOUND(ISRC)) THEN
            CALL RFORM8(RA(ISRC)*PI/12D0,LRA,16,3)
            CALL SFORM8(DEC(ISRC)*PI/180D0,LDEC,16,3)
            WRITE(PR,640) SNAME(ISRC),LRA,LDEC
            CDEC(ISRC) = DCOS(DEC(ISRC)*PI/180D0)
            SDEC(ISRC) = DSIN(DEC(ISRC)*PI/180D0)
      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
      WRITE(OUTC,501) '$Read telescope positions from a file? (N)'
      READ(INC,501) ANSWER
      IF(ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y') THEN
          WRITE(OUTC,501) '$File with telescope positions: '
          READ(INC,501) POSFIL
          OPEN (UNIT=INSTA, NAME=POSFIL, READONLY,
     1           FORM='FORMATTED', STATUS='OLD', ERR=2900)
          WRITE (PR,620)
C
 2200     READ (INSTA,521,END=2400,ERR=2900) XSTNAME,RR,TD,TM,TS,
     1           PD,PM,PS, XZL,XHL
          CALL STR$UPCASE(XSTNAME,XSTNAME)
      ELSE
          NS=0
          WRITE(OUTC,501) '$Site name, latitude (d, m, s), 
     1           longitude (d, m, s), zenith angle lim (d), ',/,
     2           '$hour angle lim (h) : '
          READ(INC,501) STNAME(1), TD, TM, TS, PD, PM, PS, ZL, HL
          WRITE(OUTC,501) '$Telescope positions by R (R(meters), 
     1           theta(deg CCW from north))',/,
     2           '$or X (X(meters east), Y(meters north))? (default: X) '
          READ(INC,501) ANSWER
          IF (ANSWER.EQ.'R'.OR.ANSWER.EQ.'r') THEN
2250          WRITE(OUTC,501) '$R, theta (R>10000=>stop): '
              NS=NS+1
              READ(INC,*) RADIUS, THETA
              XPOS(NS)=-RADIUS*SIN(THETA*RPDEG)
              YPOS(NS)=RADIUS*COS(THETA*RPDEG)
              IF(RADIUS.LT.10000.) GOTO 2250
          ELSE
              WRITE(OUTC,501) 'X, Y ((ABS(X) or ABS(Y) > 10000 => stop): '
              NS=NS+1
              READ(INC,*) XPOS(NS),YPOS(NS)
              IF(ABS(XPOS(NS)).LT.10000..AND.ABS(YPOS(NS)).LT.10000.) GOTO 2260
2280      ENDIF
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 (PR,621) XSTNAME,RR,TD,TM,TS,PD,PM,PS
      TT = TD + TM/60D0 + TS/3600D0      ! latitude (deg)
      PP = PD + PM/60D0 + PS/3600D0      ! w longitude (deg)
      IF (RR.LT.6D6) CALL GEOID(RR,TT)
      SLAT = DSIN(TT*RPDEG)
      CLAT = DCOS(TT*RPDEG)
      TCTR = TT
      PCTR = PP/15D0
      XCTR = RR*CLAT*DCOS(PP*RPDEG)
      YCTR = RR*CLAT*DSIN(PP*RPDEG)
      ZCTR = RR*SLAT
      WRITE (PR,623) RR,TT,PP,XCTR,YCTR,ZCTR,HL,ZL
C
 2400      CLOSE(UNIT=INSTA)
C
C-----------------------------------------------------------------------
C            Coordinates of telescopes
C-----------------------------------------------------------------------
      DO 2500 I=1,NS
          X(M)=RR*DCOS(PP*RPDEG)*DCOS(TT*RPDEG)
     1           +XPOS(NS)*DSIN(PP*RPDEG)
     2           -YPOS(NS)*DCOS(PP*RPDEG)*DSIN(TT*RPDEG)
          Y(M)=RR*DSIN(PP*RPDEG)*DCOS(TT*RPDEG)
     1           -XPOS(NS)*DCOS(PP*RPDEG)
     2           -YPOS(NS)*DSIN(PP*RPDEG)*DSIN(TT*RPDEG)
          Z(M)=RR*DSIN(TT*RPDEG)
     1           +YPOS(NS)DCOS(TT*RPDEG)
2500  CONTINUE
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(M1)-X(M2))/WL
            BY = (Y(M1)-Y(M2))/WL
            BZ(J) = (Z(M1)-Z(M2))/WL
            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)
            FMT = '(I2,1X)'
            IF (J.GT.9) FMT='(I3)'
            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 at each station
C
      DO 4200 IS=1,NS
          CLCD(IS) = CLAT(IS)*CDEC(ISRC)
          CLSD(IS) = CLAT(IS)*SDEC(ISRC)
          SLCD(IS) = SLAT(IS)*CDEC(ISRC)
          SLSD(IS) = SLAT(IS)*SDEC(ISRC)
          COSZAL   = DCOS(ZL(IS)*RPDEG)
          COSHAL   = (COSZAL-SLSD(IS))/CLCD(IS)
          IF (COSHAL.GE.1.0) THEN
C            ! Source is never visible
            HAL(IS) = -0.01
            WRITE (PR,501) ' ++WARNING++ Source '//SNAME(ISRC)//
     1                  ' cannot be seen from station '//STNAME(IS)
            WRITE (OUTC,501) ' ++WARNING++ Source '//SNAME(ISRC)//
     1                  ' cannot be seen from station '//STNAME(IS)
          ELSE IF (COSHAL.LE.-1.0 .AND. HL(IS).GE.12.0) THEN
C            ! Source is circumpolar
            HAL(IS) = 12.0
          ELSE IF (COSHAL.LE.-1.0) THEN
C            ! Circumpolar, but HA limit is effective
            HAL(IS) = HL(IS)
          ELSE
C            ! ZA limit is effective; HA limit may be
            HAL(IS) = ACOS(COSHAL)/RPDEG/15.0
            HAL(IS) = MIN(HAL(IS),HL(IS))
          END IF
 4200      CONTINUE
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 for each station
C
          DO 5200 IS=1,NS
            HA(IS) = AMOD(GHA-P(IS),24.0)
            IF (HA(IS).LT.-12.0) HA(IS) = HA(IS)+24.0
            IF (HA(IS).GT.12.0)  HA(IS) = HA(IS)-24.0
            CHA = COS(HA(IS)*RPHR)
            IF (ABS(HA(IS)).LE.HAL(IS)) THEN
                ZA(IS) = ACOS(SLSD(IS)+CLCD(IS)*CHA)/RPDEG
                DEN = CDEC(ISRC)*DSIN(HA(IS)*RPHR)
                AZX = ATAN2(CLSD(IS)-SLCD(IS)*CHA, DEN)/RPDEG
                AZ(IS) = AMOD(AZX-90.0,360.0)
                IF (AZ(IS).LT.0.0) AZ(IS)=AZ(IS)+360.0
                UP(IS) = .TRUE.
                NUP = NUP+1
            ELSE
                HA(IS) = IMPOSS
                ZA(IS) = IMPOSS
                AZ(IS) = IMPOSS
                UP(IS) = .FALSE.
            END IF
 5200          CONTINUE
C
C            Compute IHA for each baseline
C
          DO 5300 IB=1,LJ
            IF (NUP.GE.MINUP .AND. UP(IS1(IB)).AND.UP(IS2(IB))) THEN
                H(IB,K) = AMOD(GHA-BH(IB)-6.0, 24.0)
                IF (H(IB,K).LT.-12.0) H(IB,K) = H(IB,K)+24.0
                IF (H(IB,K).GT.12.0)  H(IB,K) = H(IB,K)-24.0
            ELSE
                H(IB,K) = IMPOSS
            END IF
 5300          CONTINUE
C
C            Print section 1: HA, ZA, AZ for each station
C
          IF (PRINT1 .AND. NUP.GE.MINUP) 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) (STNAME(M),M=M1,M2)
                WRITE (PR,652)
            END IF
            IF (MOD(NL,10).EQ.1) WRITE(PR,90)
            WRITE (PR,655,IOSTAT=JUNK) UT,GST,
     1                  (HA(M),ZA(M),AZ(M),M=M1,M2)
          END IF
 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 6100 IB=1,LJ
                  IF (H(IB,K).NE.IMPOSS) GOTO 6200
 6100                CONTINUE
                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(IB,K),IB=J1,J2)
 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
      LEGND1 = SNAME(ISRC)
      DO 7104 IS=1,NS
          I = LEN1(LEGND1)+2
          LEGND1(I:) = STNAME(IS)
 7104      CONTINUE
      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(J,K)
            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)
                CALL GRMARK(IPLOT,.TRUE.,'03'X,.FALSE.,1,UTIK,VTIK)
                IF(INVERT) CALL GRMARK(IPLOT,.TRUE.,'03'X,
     1                  .FALSE.,1, -UTIK, -VTIK)
                CALL GRVECT(IPLOT,2,.FALSE.,1,U,V)
                END IF
 7300      CONTINUE
      IF (.NOT.INVERT) GOTO 7200
      DO 7310 K=1,LK
            HAI = H(J,K)
            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-----------------------------------------------------------------------
 1900 CALL ERROR('Cannot read SOURCES catalog')
 2900 CALL ERROR('Cannot read STATIONS catalog')
C-----------------------------------------------------------------------
C            Format statements
C-----------------------------------------------------------------------
  500 FORMAT(' HAZI computes u,v-tracks for VLBI experiments'/
     1       ' (Version ',A,')'/)
  501 FORMAT(A)
  601 FORMAT(' Program HAZI'/' (Version ',A,')'/)
 6021 FORMAT('(million wavelengths,',I6,' MHz)')
 6022 FORMAT(/5X,'Frequency',F11.3,' MHz'/
     1       5X,'Wavelength',-4PF10.3,' cm')
  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,A8/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,T121,
     1       2F6.1)
  630 FORMAT(//' Baselines (millions of wavelengths):'
     1       //5X,'#',15X,'Stations',20X,
     2       'Bx',13X,'By',13X,'Bz',14X,'Beq',11X,'GHA (h)')
  631 FORMAT(3X,I3,8X,A8,' -- ',A8,7X,3F15.6,2X,2F14.6)
  651 FORMAT(//3X,'UT',5X,'GST',3X,<NSP>(9X,A8,5X))
  652 FORMAT(16X,        