      PROGRAM SMALLDIG
C     Digitize graphed curves with a Summagraphics MM 1201.
C     Does the following:
C        1.  Reads coordinates of points on curve.
C        2.  Gives area under curve by trapezoidal rule.
C        3.  Gives area enclosed by curve.
C     Written August 1992 by R. Thrun, Code R14, NSWC
C     Language: Microsoft Fortran 5.0
C
C     Uses the following routines:
C        COMMSUB library from the Forlib-Plus package to read the
C           comm port. This library works with MS Fortran 3.30, does 
C           not work with 3.31, but seems to work with MS Fortran 5.0.
C           Forlib-Plus is published by:
C               Alpha Computer Service
C               P.O. Box 2517
C               Cypress, Calif. 90630
C        KEYBRD -- by R. Thrun, checks for keyboard input
C        TTONE --- by R. Thrun, makes beep of specified freq and duration
C
      CHARACTER*30 OFILE
      CHARACTER*60 TITLE
      INTEGER*2 I, IB(2000), IBUT, N
      INTEGER*2 CERR, CSTAT, SSTAT
      LOGICAL*2 BREAK
      REAL*4 X(2000), Y(2000)
C
      WRITE (*,*) '   SMALLDIG IV Digitizing Program'
      WRITE (*,*) '  This will give coordinates of points on a curve'
      WRITE (*,*) '  and the area under or enclosed by the curve.'
      WRITE (*,*) 'What is the name of the output file? '
      READ (*, 1) OFILE
    1 FORMAT (A30)
C     STATUS='NEW' erases existing file due to bug in MS Fortran 3.30
C     Unspecified status gives an error on write (MS Fortran 3.30 bug)
C     For other compilers, 'NEW' will check for existing file.
      OPEN (2, FILE=OFILE, STATUS='NEW')
      WRITE (*,*) 'Give a title for the data: '
      READ (*, 2) TITLE
    2 FORMAT (A60)
      WRITE (2,2) TITLE
C
C     Initialize the serial port with Forlib+ routines.
C     Requires INTEGER*2 arguments, but INTEGER*4 constants work.
C     Arguments, in order, set 9600 baud, odd parity, 1 stop bit, 
C     8 bits, comm 1. These are factory settings on the tablet.
      CALL STBAUD (7, 1, 1, 8, 0, CSTAT, CERR)
      IF (CERR .EQ. 0) THEN
         WRITE (*,*) ' Error in initializing comm. port'
         STOP
      ENDIF
C     Now go to interrupt-driven input.
      CALL COMON
C
C     Initialize the Summagraphics MM 1201 tablet. Point mode, 1000 lpi.
C     We do not care about SSTAT, the word of status bits.
      CALL SNDBYT ('B', 0, SSTAT)
      CALL SNDBYT ('j', 0, SSTAT)
C
      WRITE (*,*)'Digitize an origin point with known coordinates.'
   10 CALL GETPNT (XORG, YORG, IBUT, BREAK)
      IF (BREAK) GO TO 100
      WRITE (*,*)'What are the (X,Y) coordinates of the point? '
      READ (*,*) XZERO, YZERO
C
   20 WRITE (*, 25)
   25 FORMAT (' Digitize another point a known distance away' /
     +        ' on the same horizontal line as the origin.')
      CALL GETPNT (X2, Y2, IBUT, BREAK)
      IF (BREAK) GO TO 100
      DX = X2 - XORG
      DY = Y2 - YORG
      S = SQRT(DX**2+DY**2)
      IF (S .GT. 0.25) GO TO 30 
      WRITE (*,*) 'The points should be further apart.'
      GO TO 20
C     Cosine and sine of rotation angle in tablet coordinates.
   30 CA = DX/S
      SA = DY/S
      WRITE (*,*) 'What is the horizontal distance between points? '
      READ (*,*) XSEP
      XSCALE = XSEP/S
C
      WRITE(*,*) 'Digitize two points a known vertical distance apart'
      WRITE(*,*) '   If the points are the same or the distance between'
      WRITE(*,*) '   is given as zero, then the horizontal and vertical'
      WRITE(*,*) '   scales will be the same.'
      WRITE(*,*) '   Get point 1'
      CALL GETPNT (X1, Y1, IBUT, BREAK)
      IF (BREAK) GO TO 100
      WRITE (*,*) '   Get point 2'
      CALL GETPNT (X2, Y2, IBUT, BREAK)
      IF (BREAK) GO TO 100
C     The vertical points are not necessarily on the same vertical line
C     Get the vertical difference, then the vertical scale.
      YR1 = -SA*(X1-XORG) + CA*(Y1-YORG)
      YR2 = -SA*(X2-XORG) + CA*(Y2-YORG)
      DY = YR2 - YR1
      IF (ABS(DY) .LT. 0.05) GO TO 40
      WRITE (*,*)'What is the vertical distance between points? '
      READ (*,*) YSEP
      IF (YSEP .EQ. 0.0) GO TO 40
      YSCALE = YSEP/DY
      GO TO 45
   40 YSCALE = XSCALE
C
      WRITE (*,*) 'Digitize points along a curve.'
   45 N = 0
   50 CALL GETPNT (XD, YD, IBUT, BREAK)
      IF (BREAK) GO TO 100
C     Cursor not in proximity to tablet.
      IF (IBUT .EQ. 0) GO TO 50
C     Template at upper right corner of tablet.
      IF (XD .LT. 11.) GO TO 60
C     Stop box
      IF (YD .GT. 10.7) GO TO 160
C     Delete last point
      IF (YD .GT. 9.7) GO TO 140
C     Normal additional point
   60 N = N + 1
      DX = XD - XORG
      DY = YD - YORG
      X(N) = ( CA*DX + SA*DY)*XSCALE + XZERO
      Y(N) = (-SA*DX + CA*DY)*YSCALE + YZERO
      IB(N) = IBUT
      WRITE (*, 65) X(N), Y(N), IB(N)
   65 FORMAT (1X, 1P2G13.5, I3)
      GO TO 50
C
  100 WRITE (*,*)'User break from keyboard.'
      GO TO 160
C
C     Delete last point.
  140 IF (N .GT. 0) N = N - 1
      GO TO 50
C
C     Write points to file.
  160 WRITE (*,*) 'Writing output file.'
      DO 170 I = 1,N
      WRITE (2, 165) X(I), Y(I), IB(I)
  165 FORMAT (1P2G13.5, I3)
  170 CONTINUE
C
C     Calculate area by trapezoidal rule.
      A = 0.0
      DO 180 I = 2,N
      A = A + (X(I)-X(I-1))*(Y(I)+Y(I-1))
  180 CONTINUE
      A = A/2.0
      WRITE (*, 185) A
      WRITE (2, 185) A
  185 FORMAT (5X, 'Area under curve =       ', 1PG15.5)
C     Since the curve will never close, because we will not hit exactly
C     on the starting point, we force a closure.
      A = A + (X(1)-X(N))*(Y(1)+Y(N))/2.
      WRITE (*, 190) A
      WRITE (2, 190) A
  190 FORMAT (5X, 'Area enclosed by curve = ', 1PG15.5)
C
C     Reset comm port, close file, and stop.
  200 CALL OFFINT
      STOP
      END
      SUBROUTINE GETPNT (X, Y, IBUT, BREAK)
C     Read a point from MM 1201 tablet using Forlib+ routines.
C     With the factory default jumper settings, the tablet sends
C     each point as a series of 5 binary-encoded bytes.
C
C     Output:
C        X,Y ---- coordinates in inches from lower left corner of tablet
C        IBUT --- number of button pressed
C        BREAK -- is .TRUE. if user pressed ESC or ^C
C
      INTEGER*2 IX, IY, IBUT
      CHARACTER*1 B(5)
      INTEGER*2 AVAIL, GETC, I, IB(5), IC, IER, NBYTES
      LOGICAL*2 BREAK, KEYHIT
C
      BREAK = .FALSE.
      NBYTES = 0
C     Alternately test for bytes in serial input buffer 
C     and check keyboard for ESC or ^C, indicating user break.
   10 CALL BFSTAT (AVAIL)
      IF (AVAIL .EQ. 0) GO TO 20
      NBYTES = NBYTES + 1
      CALL FCHCHR (B(NBYTES), IER)
      IF (NBYTES .EQ. 5) GO TO 50
C     Keyboard check. Swallow any bytes except ESC or ^C.
   20 IF (.NOT. KEYHIT(0)) GO TO 10
      IC = MOD(GETC(0),256)
      BREAK = (IC .EQ. 3) .OR. (IC .EQ. 27)
      IF (BREAK) RETURN
      GO TO 10
C
C     Convert the 5 binary-mode bytes sent by the tablet.
   50 DO 55 I = 1,5
      IB(I) = ICHAR(B(I))
   55 CONTINUE
C     First byte has proximity indicator and button number.
      IF (MOD(IB(1),128) .GE. 64) THEN
         WRITE (*,*)' Cursor was not in proximity to tablet'
         CALL TTONE (100, 9)
         X = 0.0
         Y = 0.0
         IBUT = 0
         RETURN
      ENDIF
      IBUT = MOD(IB(1),8)
      IX = IB(2) + IB(3)*128
      IY = IB(4) + IB(5)*128
      X = FLOAT(IX)/1000.
      Y = FLOAT(IY)/1000.
      CALL TTONE (900, 1)

      RETURN
      END
