C @(#)plot.for	19.1 (ES0-DMD) 02/25/03 13:28:47
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C  @(#)plot.for	5.3  (ESO-IPG)  3/28/93  12:53:39
      SUBROUTINE PLOT(NPTS,X,Y,SYMBOL)
c
c	Copyright (C) Andrew T. Young, 1990
C
C  If you have any comments, or problems with this program,
C  please report them to:
C                         A.T.YOUNG                        28 DEC. 1985
C
C
C  Plots on printer or crt screen, in area ncols x nrows.
C   (Default is 80 x 24 for terminal screen.)
C
C   (Default o/p unit is 6; can be reset -- see below.)
C
C  ARGUMENTS: NPTS is no.of points in X,Y arrays; SYMBOL plots at each point.
C   (NPTS is negative if more data sets follow, but 1st set fixes limits.)
C
C    PLOT determines X- and Y-limits unless forced -- see below.
C
C  If NPTS = 0 and...
C
C    SYMBOL = 'I'   Inverts X if X(1)=1., Inverts Y if Y(1)=1.
C             'L'   sets LIMITS in X-Y plane to (X(1),X(2)), (Y(1),Y(2)).
C             'P'   sets PAGE SIZE to  NCOLS=X(1), NROWS=Y(1).
C             'S'   all points plot  SAME  SYMBOL(1).
C             'D'   all points plot  DIFFERENT  SYMBOLs(NP).
C             'O'   sets treatment for points Outside plot area.
C             'U'   sets ftn.output unit to INT(X(1)).  (DEFAULT = 6.)
C
C
      IMPLICIT NONE
C
      INTEGER NPTS, NCOLS, NROWS, NUNIT, NP, N, IX, IY, NROW, LAST, NCOL
      REAL X, Y, XMIN, XMAX, YMIN, YMAX, CMX, RMX, XX, DELX, UNITX,DELY, 
     1         UNITY, DY, YROW, DX, ABSXX
C
      SAVE
      EXTERNAL ROUND
C
      LOGICAL NOTYET,LIMITS,READY,INVX,INVY,SAME,XOUT,YOUT
      DIMENSION X(1000),Y(1000),SYMBOL(1000)
      INTEGER MXROWS,MXCOLS
      PARAMETER (MXROWS=32, MXCOLS=80)
      CHARACTER ROW(MXROWS)*80,SYMB,SYMBOL,WORD*5,W7*7
C
      DATA LIMITS/.FALSE./,READY/.FALSE./,INVX/.FALSE./,INVY/.FALSE./,
     1 SAME/.TRUE./, XOUT/.FALSE./,YOUT/.FALSE./
      DATA NCOLS/80/, NROWS/24/, NUNIT/6/
C
C
      NP=NPTS

      IF(NP)1,2,15
    1 NOTYET=.TRUE.
C    Neg.NPTS flags more data.
      NP=-NP
      GO TO 16
C
C  SPECIAL if NPTS=0.
C
    2 NOTYET=.TRUE.
      IF(SYMBOL(1).NE.'L') GO TO 3
C  L = LIMITS in X-Y plane.
      XMIN=X(1)
      XMAX=X(2)
      YMIN=Y(1)
      YMAX=Y(2)
      LIMITS=.TRUE.
      RETURN
C
C  P = PAGE SIZE -- set NROWS, NCOLS.
    3 IF(SYMBOL(1).NE.'P') GO TO 5
      IF(X(1).LT.7. .OR. Y(1).LT.5.) GO TO 90
      CMX=MXCOLS
      IF(X(1).GT.CMX)THEN
          IF(NUNIT.EQ.6)THEN
              WRITE(ROW(1),4)'NCOLS = ',X(1),CMX
              CALL TV(ROW(1))
          ELSE
              WRITE(NUNIT,4)'NCOLS = ',X(1),CMX
          END IF
      END IF
    4 FORMAT(/' NOT ENOUGH ROOM FOR ',A8,F5.0//' Reduced to ',F4.0/)
      RMX=MXROWS
      IF(Y(1).GT.RMX)THEN
          IF(NUNIT.EQ.6)THEN
              WRITE(ROW(1),4)'NROWS = ',Y(1),RMX
              CALL TV(ROW(1))
          ELSE
              WRITE(NUNIT,4)'NROWS = ',Y(1),RMX
          END IF
      END IF
      NCOLS=MIN(X(1),CMX)
      NROWS=MIN(Y(1),RMX)
      RETURN
C
C  I = INVERTED plots.
    5 IF(SYMBOL(1).NE.'I') GO TO 7
      INVX=.FALSE.
      IF(X(1).NE.1.) GO TO 6
      INVX=.TRUE.
    6 INVY=.FALSE.
      IF(Y(1).NE.1.)RETURN
      INVY=.TRUE.
      RETURN
C
C  S = SAME SYMBOL for all points.
    7 IF(SYMBOL(1).NE.'S') GO TO 8
      SAME=.TRUE.
      RETURN
C
C  D = DIFFERENT SYMBOL for each point.
    8 IF(SYMBOL(1).NE.'D') GO TO 9
      SAME=.FALSE.
      RETURN
C
C  O = treatment for points OUTSIDE plot area:
C      flagged by ^ or v if Y(1)=1.  (Not flagged if 0.)
C      flagged by > or < if X(1)=1.  (Not flagged if 0.)
    9 IF(SYMBOL(1).NE.'O')GO TO 10
      XOUT=.FALSE.
      YOUT=.FALSE.
      IF (X(1).GT.0.) XOUT=.TRUE.
      IF (Y(1).GT.0.) YOUT=.TRUE.
      RETURN
C
C  U = UNIT number.
   10 IF(SYMBOL(1).NE.'U')GO TO 91
      NUNIT=X(1)
      RETURN
C
C
C   NORMAL ENTRY.
C
   15 NOTYET=.FALSE.
   16 IF(READY) GO TO 34
      IF(LIMITS) GO TO 18
C
C  SET LIMITS on first entry, if not set before.
C
      XMIN=X(1)
      XMAX=X(1)
      YMIN=Y(1)
      YMAX=Y(1)
      IF(NP.EQ.1) GO TO 92
      DO 17 N=2,NP
      XMIN=MIN(XMIN,X(N))
      XMAX=MAX(XMAX,X(N))
      YMIN=MIN(YMIN,Y(N))
   17 YMAX=MAX(YMAX,Y(N))
C    INVERT limits:
      IF(INVX)THEN
      XX=XMIN
      XMIN=XMAX
      XMAX=XX
      END IF
      IF(INVY)THEN
      XX=YMIN
      YMIN=YMAX
      YMAX=XX
      END IF
C    Expand limits to next ROUND values.
   18 IF(XMIN.EQ.XMAX .OR. YMIN.EQ.YMAX) GO TO 95
      CALL ROUND(XMIN,XMAX,NCOLS-5,DELX,IX,UNITX)
      CALL ROUND(YMIN,YMAX,NROWS-3,DELY,IY,UNITY)
C
C  CLEAR IMAGE.
C
      DO 19 NROW=1,NROWS-2
   19 ROW(NROW)(:NCOLS)='     I'
C
C  LABEL AXES, vert.first.
C
      DY=(DELY/(IY))*UNITY
      YROW=YMIN/UNITY
C   Ticks and values.
      DO 20 NROW=NROWS-3,1,-IY
      WRITE(ROW(NROW),'(1X,F4.0,''+'')') YROW
      LAST=NROW
   20 YROW=YROW+DELY
      IF(UNITY-1.)21,32,24
C   Scale label.
   21 WRITE(WORD,'(F4.2)')UNITY
      IF(UNITY.GT..009)GO TO 25
C    Exponent.
   22 WRITE(W7,'(1PE7.1)')UNITY
      WORD(2:5)=W7(4:7)
      GO TO 25
   24 WRITE(WORD,'(F5.0)')UNITY
      IF(UNITY.GT.150.)GO TO 22
   25 READ(WORD,'(1X,A4)')WORD
      IF(LAST.GT.2) LAST=0
      LAST=LAST+1
      ROW(LAST+1)(2:5)=WORD(:4)
      IF(UNITY.LT.150. .AND. UNITY.GT..009)ROW(LAST)(3:3)='X'
C
C   X-AXIS...
C
   32 ROW(NROWS-3)(7:NCOLS)='-------------------------------------------
     1------------------------------------------------------------------
     2-----'
      DX=(DELX/(IX))*UNITX
      XX=XMIN/UNITX
C   Ticks and values.
      DO 33 NCOL=6,NCOLS-2,IX
C       ABSXX=ABS(XX)+.4
C       Find longest label to print:
        ABSXX=MAX(XMIN, XMAX, -10.*XMAX, -10.*XMIN)/UNITX+.4
        IF (ABSXX.LT.99.5) THEN
          WRITE(ROW(NROWS-2)(NCOL-2:NCOL+1),'(F4.0)') XX
        ELSE IF(IX.GE.4 .AND. ABSXX.LT.999.5) THEN
          WRITE(ROW(NROWS-2)(NCOL-2:NCOL+2),'(F5.0)') XX
        ELSE IF(IX.GT.4 .AND. ABSXX.LT.9999.5) THEN
          WRITE(ROW(NROWS-2)(NCOL-3:NCOL+2),'(F6.0)') XX
        ELSE
C         Label every other tick.
          IF(MOD(ABS(XX)+.4,ABS(DELX+DELX)).LT.0.5)
     1       WRITE(ROW(NROWS-2)(NCOL-4:NCOL+2),'(F7.0)') XX
        END IF
        ROW(NROWS-3)(NCOL:NCOL)='+'
   33 XX=XX+DELX
      IF(UNITX.EQ.1.) GO TO 34
C    Exponent.
      WRITE(W7,'(1PE7.1)')UNITX
      ROW(NROWS-2)(NCOLS-3:NCOLS)=W7(4:7)
   34 READY=.TRUE.
      IF(NP.EQ.0) RETURN
C
C  PUT POINTS in IMAGE.
C
      DO 35 N=1,NP
        IF (X(N).EQ.3.E33 .OR. Y(N).EQ.3.E33) GO TO 35
        NROW=(YMAX-Y(N))/DY + 1.5
        NCOL=(X(N)-XMIN)/DX + 6.5
C    Treat pts.Outside area:
        IF(NROW.LT.1) THEN
C         Y is above area.
          IF (YOUT) THEN
C           Put edge marker in plot.
            IF(NCOL.LT.6) THEN
C             Both X & Y out.  Mark upper left corner.
              IF (XOUT) ROW(1)(6:6)=CHAR(92)
            ELSE IF(NCOL.GT.NCOLS) THEN
C             Both X & Y out.  Mark upper right corner.
              IF (XOUT) ROW(1)(NCOLS:NCOLS)='/'
            ELSE
C             X is inside.  Mark correct column.
              ROW(1)(NCOL:NCOL)='^'
            END IF
          ELSE
C           Ignore.
          END IF
          GO TO 35
        ELSE IF(NROW.GT.NROWS-3) THEN
C         Y is below area.
          IF (YOUT) THEN
C           Put edge marker in plot.
            IF(NCOL.LT.6) THEN
C             Both X & Y out.  Mark lower left corner.
              IF (XOUT) ROW(NROWS-3)(6:6)='/'
            ELSE IF(NCOL.GT.NCOLS) THEN
C             Both X & Y out.  Mark lower right corner.
              IF (XOUT) ROW(NROWS-3)(NCOLS:NCOLS)=CHAR(92)
            ELSE
C             X is inside.  Mark correct column.
              ROW(NROWS-3)(NCOL:NCOL)='v'
            END IF
          ELSE
C           Ignore.
          END IF
          GO TO 35
        ELSE
C         Y is inside area, so NROW is OK.  Check X:
          IF(NCOL.LT.6) THEN
C           X is off left edge.
            IF (XOUT) ROW(NROW)(6:6)='<'
            GO TO 35
          ELSE IF(NCOL.GT.NCOLS) THEN
C           X is off right edge.
            IF (XOUT) ROW(NROW)(NCOLS:NCOLS)='>'
            GO TO 35
          ELSE
C           Both are inside.  Plot.
          END IF
        END IF
C
        SYMB=ROW(NROW)(NCOL:NCOL)
        IF ((XOUT.OR.YOUT).AND.(SYMB.EQ.'^'.OR.SYMB.EQ.'v'.OR.
     1           SYMB.EQ.'>'.OR.SYMB.EQ.'<'.OR.SYMB.EQ.'/'.OR.
     2           SYMB.EQ.CHAR(92)))THEN
        ELSE IF(SYMB.EQ.' '.OR.SYMB.EQ.'-'.OR.
     1     SYMB.EQ.' '.OR.SYMB.EQ.'-'.OR.SYMB.EQ.'I'.OR.SYMB.EQ.'+')THEN
          IF(SAME)THEN
            ROW(NROW)(NCOL:NCOL)=SYMBOL(1)
          ELSE
            ROW(NROW)(NCOL:NCOL)=SYMBOL(N)
          END IF
        ELSE
          ROW(NROW)(NCOL:NCOL)='$'
        END IF
   35 CONTINUE
      IF(NOTYET) RETURN
C
C  PRINT PLOT.
C
      IF(NUNIT.EQ.6)THEN
          DO 40 N=1,NROWS-2
              CALL TVN(ROW(N)(:NCOLS))
   40     CONTINUE
      ELSE
          WRITE(NUNIT,'(A)') (ROW(N)(:NCOLS),N=1,NROWS-2)
      END IF
      LIMITS=.FALSE.
      READY=.FALSE.
      RETURN
C
C  ERROR RETURNS.
C
   90 IF(NUNIT.EQ.6)THEN
          WRITE(ROW(1),'(2(A,G10.5))')'NCOLS = ',X(1),', NROWS = ',Y(1)
          CALL TVN(ROW(1))
          WRITE(ROW(1),'(A)')'PAGE TOO SMALL TO PLOT IN.'
          CALL TVN(ROW(1))
      ELSE
          WRITE(NUNIT,'(2(A,G10.5))')'NCOLS = ',X(1),', NROWS = ',Y(1)
          WRITE(NUNIT,'(A)')'PAGE TOO SMALL TO PLOT IN.'
C          WRITE(NUNIT,*)'NCOLS = ',X(1),', NROWS = ',Y(1)
C          WRITE(NUNIT,*)'PAGE TOO SMALL TO PLOT IN.'
      END IF
      GO TO 99
C
   91 IF(NUNIT.EQ.6)THEN
      WRITE(ROW(1),'(A)')SYMBOL(1)//
     1       ' IS NOT A RECOGNISED CODE FOR NPTS = 0.'
       CALL TVN(ROW(1))
       WRITE(ROW(1),'(A)')
     1    'POSSIBLE ATTEMPT TO PLOT A ZERO-LENGTH ARRAY.'
       CALL TVN(ROW(1))
      ELSE
       WRITE(NUNIT,'(A)')SYMBOL(1),' IS NOT A RECOGNISED CODE'//
     1     ' FOR NPTS = 0.'
       WRITE(NUNIT,'(A)')'POSSIBLE ATTEMPT TO PLOT A ZERO-LENGTH ARRAY.'
C       WRITE(NUNIT,*)SYMBOL(1),' IS NOT A RECOGNISED CODE FOR NPTS = 0.'
C       WRITE(NUNIT,*)'POSSIBLE ATTEMPT TO PLOT A ZERO-LENGTH ARRAY.'
      END IF
      GO TO 99
C
   92 IF(NUNIT.EQ.6)THEN
        WRITE(ROW(1),'(2(A,G10.5),A)')'ONE POINT (',X(1),',',Y(1),
     1     ') CANNOT SET LIMITS.'
        CALL TVN(ROW(1))
      ELSE
        WRITE(NUNIT,'(2(A,G10.5),A)')'ONE POINT (',X(1),',',Y(1),
     1     ') CANNOT SET LIMITS.'
C        WRITE(NUNIT,*)'ONE POINT (',X(1),',',Y(1),') CANNOT SET LIMITS.'
      END IF
      GO TO 99
C
   95 IF(XMAX.EQ.XMIN) SYMB='X'
      IF(YMAX.EQ.YMIN) SYMB='Y'
      IF(NUNIT.EQ.6)THEN
          WRITE(ROW(1),'(A)')'ALL '//SYMB//' VALUES EQUAL.'
          CALL TVN(ROW(1))
      ELSE
          WRITE(NUNIT,'(A)')'ALL '//SYMB//' VALUES EQUAL.'
      END IF
      LIMITS=.FALSE.
C
   99 IF(NUNIT.EQ.6)THEN
          CALL TV('NO PLOT PRODUCED.')
          CALL TVN(' ')
      ELSE
          WRITE(NUNIT,'(A)') ' '
          WRITE(NUNIT,'(A)')'NO PLOT PRODUCED.'
          WRITE(NUNIT,'(A)') ' '
C          WRITE(NUNIT,*)
C          WRITE(NUNIT,*)'NO PLOT PRODUCED.'
C          WRITE(NUNIT,*) 
      END IF
      RETURN
C
      END
      SUBROUTINE XAXIS(XLIMS)
C
C   Adds dashed X-axis to plot.  Limits MUST be established first!
C
      IMPLICIT NONE
C
      REAL XLIMS(2), XVAL, STEP
C
      INTEGER I
C
C
      STEP=(XLIMS(2)-XLIMS(1))*0.05
      XVAL=XLIMS(1)+STEP/2.
      DO 10 I=1,20
        CALL PLOT(-1, XVAL, 0., '-')
        XVAL=XVAL+STEP
   10 CONTINUE
C
C
      RETURN
      END
