$DEBUG
C
      SUBROUTINE SQUINT(DELTA)
c
c   SQUEEZE of data points (if point number/frequency characteristic too 
c   shallow) and linearisation through quadratic interpolation
c
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*2 (I-N)
      PARAMETER (NMAXPT=32000, maxfre=27000)
      INTEGER*4 FREQS,MINF,MAXF,FTURN(10),INTERM
      INTEGER*2 NTURN(10)
C
      COMMON /SPEC/freqs(maxfre),ISPEC(NMAXPT),NPTS,ISMALL,ILARGE
     1       /PLOTC/IPLOT(NMAXPT)
C
      WRITE(*,10)
10    FORMAT(1X/' **** NOTE: this option should only be run once ',
     * 'all unlocks have been eliminated'/)
c
c...SQUEEZE closely spaced data points: such pairs of points are replaced
c   by a single point of mean frequency and intensity
c
      MINF=1000000000
      MAXF=0
      DO 1 N=1,NPTS
        IF(FREQS(N).LT.MINF)THEN
          MINF=FREQS(N)
          NMIN=N
        ENDIF
        IF(FREQS(N).GT.MAXF)THEN
          MAXF=FREQS(N)
          NMAX=N        
        endif
1     CONTINUE
C
C...first scan up, second down
      IF(MAXF.GT.FREQS(1).AND.MAXF.GT.FREQS(NPTS))THEN
        DELTA=DBLE(MAXF)-DBLE(FREQS(1))+DBLE(MAXF)-
     *      DBLE(FREQS(NPTS))
      ENDIF
C
C...first scan down, second up
      IF(MINF.LT.FREQS(1).AND.MINF.LT.FREQS(NPTS))THEN
        DELTA=DBLE(FREQS(1))-DBLE(MINF)+DBLE(FREQS(NPTS))-
     *      DBLE(MINF)
      ENDIF
C
      NITER=0
C
2     HALFAV=0.5D0*DELTA/NPTS
      ICHANG=0
      NITER=NITER+1
C
      I=1
      ii=1
3     IF(ABS(dble(FREQS(I))-dble(FREQS(I+1))).LT.HALFAV)THEN
        FREQS(ii)=0.5d0*(dble(freqs(i))+dble(freqs(i+1)))
        ISPEC(II)=0.5d0*(dble(ISPEC(I))+dble(ispec(i+1)))
        i=i+2
        ii=ii+1
        ICHANG=1         
      ELSE
        ISPEC(II)=ISPEC(I)
        FREQS(II)=FREQS(I)
        i=i+1
        ii=ii+1
      ENDIF
      if(i.lt.npts)goto 3
      IF(I.EQ.NPTS)THEN
        FREQS(II)=freqs(I)
        ispec(II)=ispec(I)
        ii=ii+1       
      ENDIF
      WRITE(*,'('' Iteration number'',I3,''  average step'',F8.4,
     * ''  squeezed pairs'',I4)')
     * NITER,(DELTA*0.001D0)/NPTS,NPTS-II+1
      NPTS=II-1
c
      IF(NITER.GT.30)GOTO 20
      IF(ICHANG.NE.0)GOTO 2
C
C...Prelude to linearisation - determine frequencies at turning points ie.
c   locate changes in sign of the gradient, to keep things simple the central
c   turnover point is discarded and its neighbours taken.  This account for 
c   the various possibilities of a jump in frequency at joint between two
c   spectral directions
C
20    I=1
      NTURN(I)=1
      FTURN(I)=FREQS(1)
      DELTA=DBLE(FREQS(2))-DBLE(FREQS(1))
C
      aversp=ispec(1)+ispec(npts)
      DO 11 N=2,NPTS-1
        DELTA1=DBLE(FREQS(N+1))-DBLE(FREQS(N))
        IF(DELTA1/DELTA.LT.0.D0)THEN
          I=I+1
          NTURN(I)=N-1
          FTURN(I)=FREQS(N-1)
          I=I+1
          NTURN(I)=N+1
          FTURN(I)=FREQS(N+1)
          IF(I.GE.9)GOTO 12
        ENDIF
        DELTA=DELTA1
        aversp=aversp+ispec(n)
11    CONTINUE
C
      I=I+1
      NTURN(I)=NPTS
      FTURN(I)=FREQS(NPTS)
12    WRITE(*,14)
14    FORMAT(1X/' Turning points:')
      DO 15 N=1,I
        WRITE(*,16)N,NTURN(N),DBLE(FTURN(N))*0.001D0
15    CONTINUE
16    FORMAT(I20,'.   point, frequency:  ',I5,F11.3)
C
C...Linearisation is only allowed if there are four entries in the
C   turning point table
C
      IF(I.EQ.4)THEN
23      WRITE(*,22)
22      FORMAT(1X/' DO YOU WANT TO LINEARIZE ?   ',$)
        READ(*,'(I1)',ERR=23)IFLAG
        IF(IFLAG.EQ.1)GOTO 24
        GOTO 17
      ELSE
        WRITE(*,21)
21      FORMAT(1X/' THIS NUMBER HAS TO EQUAL 4 FOR LINEARISATION')
        WRITE(*,'(1X/'' .... Press ENTER to continue .... '',$)')
        READ(*,'(I1)',ERR=17)II
        GOTO 17
      ENDIF
C
C...LINEARISATION proper of the two spectral sections between the turning
C   points
C
24    ISMALL=32000
      ILARGE=-32000
      aversp=aversp/npts
c
      iturn=1
      DELTA=(DBLE(FTURN(2))-dble(FTURN(1)))/
     *    dble(NTURN(2)-NTURN(1))
      fend=dble(fturn(2))-delta
      fstart=dble(fturn(1))+delta
      X=fstart-delta
      J=0
      I=2
      WRITE(*,'(1X/'' Linearisation now at point:      '',$)')
C
C----------------
C   INTERPOLATION LOOP
C
801   X=X+DELTA
      IF(DELTA.GT.0.0)THEN
        IF(X.ge.FEND)GOTO 43
      ELSE
        IF(X.le.FEND)GOTO 43
      ENDIF
      J=J+1
C
C...centre the three data points to be used for interpolation on the
C   interpolation point
c
30    dd=abs(x-freqs(i))
      if(dd.gt.abs(x-freqs(i+1)))then
        i=i+1
        goto 30
      endif
      IF(I.EQ.100*(I/100))WRITE(*,'(I6,$)')I
C
C...quadratic interpolation for unequal intervals, using three pairs of points
C   (x1,y1), (x2,y2), (x3,y3) using Newton's formula:
C              f(x) = y1 + (x-x1) D12 + (x-x1)(x-x2) D123
C               D12 = (y2-y1)/(x2-x1)
C               D23 = (y3-y2)/(x3-x2)
C              D123 = (D23-D12)/(x3-x1)
c
      x1=dble(freqs(i-1))
      x2=dble(freqs(i))
      x3=dble(freqs(i+1))
      y1=ispec(i-1)
      y2=ispec(i)
      y3=ispec(i+1)
      D12=(y2-y1)/(x2-x1)
      D23=(y3-y2)/(x3-x2)
      D123=(D23-D12)/(x3-x1)
      INTERM=y1+(x-x1)*D12+(x-x1)*(x-x2)*D123
      IF(INTERM.GT.32000)INTERM=32000
      IF(INTERM.LT.-32000)INTERM=-32000
      IPLOT(J)=INTERM
      if(abs(x1-x2).gt.5.*abs(delta).and.
     *      ((x-x2.lt.0.d0.and.delta.gt.0.d0).or.
     *       (x-x2.gt.0.d0.and.delta.lt.0.d0) ))iplot(j)=aversp
      if(abs(x3-x2).gt.5.*abs(delta).and.
     *       ((x-x2.gt.0.d0.and.delta.gt.0.d0).or.
     *        (x-x2.lt.0.d0.and.delta.lt.0.d0) ))iplot(j)=aversp
      IF(IPLOT(J).LT.ISMALL)ISMALL=IPLOT(J)
      IF(IPLOT(J).GT.ILARGE)ILARGE=IPLOT(J)
      GOTO 801
C----------------
c...bookkeeping and backtrack to linearize the return spectrum.  To
C   facilitate coaddition it is ensured that points on the return spectrum
c   coincide in frequency with those on the forward spectrum.
c
43    if(iturn.eq.1)then
        fturn(1)=fstart
        nturn(2)=j
        fturn(2)=x-delta
        delta=-delta
        i=nturn(3)+1
        nturn(3)=j+1       
        fturn(3)=fturn(3)+delta
        N=NINT((FTURN(3)-FTURN(2))/DELTA)
        FTURN(3)=FTURN(2)+N*DELTA
        x=fturn(3)-delta
        fend=fturn(4)-delta
        iturn=2
        goto 801
      else
        nturn(4)=j
        fturn(4)=x-delta
      endif  
      NPTS=J
c     
c...transfer linearized spectrum into appropriate arrays
c
      delta=-delta
      freq=fturn(1)-delta
      DO 42 N=1,nturn(2)
        freq=freq+delta
        ISPEC(N)=IPLOT(N)
        freqs(n)=freq
42    CONTINUE
      delta=-delta
      freq=fturn(3)-delta
      do 53 n=nturn(3),nturn(4)
        freq=freq+delta
        ISPEC(N)=IPLOT(N)
        freqs(n)=freq     
53    continue
      WRITE(*,'(1X)')
      DELTA=ABS(DELTA)
C
17    RETURN
      END
C
C------------------------------------------------------------------------
C
      SUBROUTINE READSP(DELTAL,FILNAM)
c      
      IMPLICIT INTEGER*2 (I-N)
      IMPLICIT REAL*4 (A-H,O-Z)
      PARAMETER (NMAXPT=32000, maxfre=27000, maxcal=200, 
     *  maxspl=500, maxvol=500)
c
      COMMON /CALVOL/VKM,IVKM,NVKM,PPS
      COMMON /LINFRE/FSTART,FEND,FINCR
      COMMON /CALFRE/XCALPT,YCALPT,NCALPT,COM1,LINEAR
      COMMON /SPEC/freqs(maxfre),ISPEC(NMAXPT),NPTS,ISMALL,ILARGE,
     1       /INFO/COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
     1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
     1 SCANSP,FRMOD,FRAMPL
      COMMON /OLD/IOLD(NMAXPT),nptold,nvkmol,ivkmol(maxvol),v1old,v2old,
     1   vkmold(maxvol),ncalol,fstold,fenold,finold,xcalol(maxcal),
     1   ycalol(maxcal)
      COMMON /PRETTY/BOLD,NORMAL
      COMMON /SPL/NSPLIN,X,Y,A,B,C,D,GRA,CA,GRB,CB
c
      INTEGER*2 IVKM(maxvol),NVKM,LINEAR
      INTEGER*4 freqs,MAXF,MINF
      REAL*4 VKM(maxvol),PPS
      REAL*8 X(maxspl),Y(maxspl),A(maxspl),B(maxspl),C(maxspl),
     * D(maxspl),GRA,CA,GRB,CB
      REAL*8 xcalol,ycalol,fstold,fenold,finold
      REAL*8 XCALPT(maxcal),YCALPT(maxcal),FSTART,FEND,FINCR
      REAL*8 ASCMIN,ASCMAX,XASC,YASC,SCASC,DELTAL,DELTA1,DELTA2
      CHARACTER FILNAM*25,CHARFL,BOLD*4,NORMAL*4
      CHARACTER*72 COMENT,COM1
      CHARACTER*6 LAMP,SCANSP
      CHARACTER*20 SAMPLE
C
C...Input of spectrum stored in the standard format
C
499   WRITE(*,'(1X//'' NAME OF DATA FILE ?  ''\)')
      READ(*,'(A)',ERR=499)FILNAM
      OPEN(2,FILE=FILNAM,FORM='BINARY',ERR=499)
      READ(2,ERR=494)COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
     1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
     1 SCANSP,PPS,FRMOD,FRAMPL
      READ(2,ERR=494)NPTS,ISMALL,ILARGE
      IF(IYEAR.LT.1980.OR.IMON.GT.12.OR.IDAY.GT.31.OR.IMON.LT.1.
     *    OR.IDAY.LT.1.OR.PPS.EQ.0.0)GOTO 494
C
      WRITE(*,155)CHAR(27),'[2J'
      WRITE(*,'(1H+,A1,''[7m'',2X,A72,2X,A1,''[0m''/)')CHAR(27),
     1 COMENT,CHAR(27)
      WRITE(*,'(''          Time and Date:  '',A4,I2,'':'',I2,'':'',
     1 I2,'', '',I2,''/'',I2,''/'',I4,A4)')BOLD,IHOUR,IMIN,ISEC,
     1 IDAY,IMON,IYEAR,NORMAL
      WRITE(*,'(1X/''                   Lamp:   '',A4,A,A4)')
     1 BOLD,LAMP,NORMAL
      WRITE(*,'(''         Vkm limits /mV:  '',A4,F6.0,''-'',
     1 F6.0,A4)')BOLD,VKMST,VKMEND,NORMAL
      WRITE(*,'(''                Grid /V:  '',A4,F6.1,A4)')
     1 BOLD,GRID,NORMAL
      WRITE(*,'(1X/''                 Sample:  '',A4,A,A4)')
     1 BOLD,SAMPLE,NORMAL
      WRITE(*,'('' Sample pressure /mTorr:  '',A4,F6.1,A4)')
     1 BOLD,SAMPRE,NORMAL
      WRITE(*,'(1X/''               Gain /mV:  '',A4,F6.1,A4)')
     1 BOLD,GAIN,NORMAL
      WRITE(*,'(''       Time constant /s:  '',A4,F6.3,A4)')
     1 BOLD,TIMEC,NORMAL
      WRITE(*,'(''             Phase /deg:  '',A4,F6.1,A4)')
     1 BOLD,PHASE,NORMAL
      WRITE(*,'(''      SIPLOW scan speed:    '',A4,A,A4)')
     1 BOLD,SCANSP,NORMAL
      WRITE(*,'(''      Points per second:  '',A4,F6.1,A4)')
     1 BOLD,PPS,NORMAL
      WRITE(*,'(1X/''  Modulation freq. /KHz:  '',A4,F8.3,A4)')
     1 BOLD,FRMOD,NORMAL
      WRITE(*,'(''   Modulation ampl. /mV:  '',A4,F6.1,A4)')
     1 BOLD,FRAMPL,NORMAL
C
      DO 120 N=1,NPTS
        READ(2,ERR=494)ISPEC(N)
        IOLD(N)=ISPEC(N)
120   CONTINUE
      NPTOLD=NPTS
C
C...read BWO voltages (present if PPS negative)
C
      IF(PPS.LT.0)THEN
        READ(2)NVKM
        WRITE(*,'(/1X,A4,I5,A4,'' measurements of BWO voltage'')')
     *      BOLD,NVKM,NORMAL         
        NVKMOL=NVKM
        V1OLD=VKMST
        V2OLD=VKMEND
        DO 4561 N=1,NVKM
          READ(2)IVKM(N),VKM(N)
          IVKMOL(N)=IVKM(N)
          VKMOLD(N)=VKM(N)
4561    CONTINUE
        OPEN(3,FILE='VOLTS',STATUS='UNKNOWN')
        DO 4560 N=1,NVKM
          WRITE(3,4562)IVKM(N),VKM(N)
4560    CONTINUE
4562    FORMAT(I5,F15.5)
        CLOSE(3)
        NSPLIN=NVKM
        DO 637 N=1,NVKM
          X(N)=IVKM(N)
          Y(N)=VKM(N)
637     CONTINUE
        CALL SPLINE
      ENDIF
C
C...Read frequency linearized file (present if PPS > 50 or comment is
C   as set by FREQLIN), this may contain frequency calibration points
c   which are a leftover from FREQLIN
C
      IF(COMENT.EQ.COM1.or.PPS.GT.50.)THEN
        READ(2)FSTART,FEND,FINCR
          WRITE(*,'(1x/''  Frequency from '',A4,F9.2,A4,'' to '',
     1  A4,F9.2,A4,'' MHz,  every:'',F9.6,'' MHz'')')
     2  BOLD,FSTART,NORMAL,BOLD,FEND,NORMAL,FINCR
          FSTOLD=FSTART
          FENOLD=FEND
          FINOLD=FINCR
          READ(2)NCALPT
          NCALOL=NCALPT
          DO 700 N=1,NCALPT
            READ(2)XCALPT(N),YCALPT(N)
            XCALOL(N)=XCALPT(N)
            YCALOL(N)=YCALPT(N)
700       CONTINUE
          LINEAR=1
      ENDIF
C
C...Read the frequency channel (present if PPS < -500), and decide whether
C   spectrum has not been linearized
C
      IF(PPS.LT.-500.)THEN
        MAXF=0
        MINF=1000000000
        NMINF=0
        NMAXF=0
        DO 710 N=1,NPTS
          READ(2)FREQS(N)
          IF(FREQS(N).GT.MAXF)THEN
            MAXF=FREQS(N)
            NMAXF=N
          ENDIF
          IF(FREQS(N).LT.MINF)THEN
            MINF=FREQS(N)
            NMINF=N
          ENDIF
710     CONTINUE
        NLIM=MIN(NMINF,NMAXF)
        IF(NLIM.EQ.1)NLIM=MAX(NMINF,NMAXF)
c
c...single gradient spectrum
        if(nlim.eq.npts)then
          DELTA1=ABS(DBLE(FREQS(NLIM)-FREQS(1))/(NLIM-1))
          n=npts/2
          delta2=abs(dble(freqs(n)-freqs(1))/(n-1))
          IF(MAX(DELTA1,DELTA2)/MIN(DELTA1,DELTA2).GT.1.001D0)GOTO 20
          deltal=delta1
          goto 20
        endif          
c
c...forward/backward spectrum
        DELTA1=ABS(DBLE(FREQS(NLIM-1)-FREQS(1))/(NLIM-2))
        DELTA2=ABS(DBLE(FREQS(NPTS)-FREQS(NLIM+1))/(NPTS-NLIM-1))       
        IF(MAX(DELTA1,DELTA2)/MIN(DELTA1,DELTA2).GT.1.001D0)GOTO 20
        DELTAL=DELTA1
      ENDIF
C
20    IF(PPS.GT.-500.)THEN
         WRITE(*,'(1X,A4,I5,A4,'' points,  ranging from''
     1    ,A4,I6,A4,'' to'',A4,I6,A4)')BOLD,NPTS,NORMAL,BOLD,
     1    ISMALL,NORMAL,BOLD,ILARGE,NORMAL
      ELSE
         WRITE(*,'(1X,A4,I5,A4,'' points,  ranging from''
     1    ,A4,I6,A4,'' to'',A4,I6,A4,'' and '',A,F8.1,''-'',F8.1,
     1    A,'' MHz'')')BOLD,NPTS,NORMAL,BOLD,ISMALL,NORMAL,BOLD,
     1    ILARGE,NORMAL,BOLD,DBLE(MINF)*0.001D0,DBLE(MAXF)*0.001D0,
     1    NORMAL
      ENDIF
C        
      GOTO 493
C    
C...Input from an ASCII file which is executed once binary read fails
C   (data is treated as linearized file and only the frequency values
C    of the first and last points matter)
C
494   WRITE(*,492)char(7)
492   FORMAT(1X/'---- This is not a file in ZK binary format: will now',
     * ' try reading it'/
     *          '     as two column ASCII, assuming that first line '
     * ' contains the'/
     *          '     number of points and second is the comment',A)
      CLOSE(2)
      OPEN(2,FILE=FILNAM,ERR=489)
      READ(2,*,ERR=489)NPTS
      READ(2,1300,ERR=489)COMENT
1300  FORMAT(A)
      ASCMIN=1.E20
      ASCMAX=-1.E20
      DO 488 N=1,NPTS
        READ(2,*,ERR=489)XASC,YASC
        IF(N.EQ.1)FSTART=XASC
        IF(N.EQ.NPTS)FEND=XASC
        IF(YASC.GT.ASCMAX)ASCMAX=YASC
        IF(YASC.LT.ASCMIN)ASCMIN=YASC
488   CONTINUE
      FINCR=(FEND-FSTART)/(NPTS-1)
      FSTOLD=FSTART
      FENOLD=FEND
      FINOLD=FINCR
      NCALPT=0
      NCALOL=0
      LINEAR=1
      REWIND(2)
      READ(2,1300)CHARFL
      READ(2,1300)CHARFL      
      SCASC=20000.D0/(ASCMAX-ASCMIN)
      DO 487 N=1,NPTS
        READ(2,*)XASC,YASC
        ISPEC(N)=(YASC-ASCMIN)*SCASC-10000.D0
        IOLD(N)=ISPEC(N)
487   CONTINUE
      NPTOLD=NPTS
      ISMALL=-10000
      ILARGE=10000
      PPS=51.
      IDAY=1
      IMON=1
      IYEAR=2000
      WRITE(SAMPLE,'(1P(2E10.3))')ASCMIN,ASCMAX
      SCANSP=' '
      LAMP=' '
      WRITE(*,155)CHAR(27),'[2J'
155   FORMAT(1H+,A1,A)
      WRITE(*,'(1H+,A1,''[7m'',2X,A72,2X,A1,''[0m''/)')CHAR(27),
     1 COMENT,CHAR(27)
      WRITE(*,'(1X/1X,A4,I5,A4,'' points,  ranging from''
     1 ,A4,1PE12.4,A4,'' to'',A4,E12.4,A4)')BOLD,NPTS,NORMAL,BOLD,
     1 ASCMIN,NORMAL,BOLD,ASCMAX,NORMAL
       WRITE(*,'(''  Frequency from '',A4,F9.2,A4,'' to '',A4,F9.2,
     1  A4,'' units'')')BOLD,FSTART,NORMAL,BOLD,FEND,NORMAL
      GOTO 493
C
489   WRITE(*,'(1X/'' ***** INPUT PROBLEMS!!!!''/)')
C
493   CLOSE(2)
      WRITE(*,491)
491   FORMAT(1X/' .... Press ENTER to continue.....  ',$)
      READ(*,1300,ERR=490)CHARFL
490   WRITE(*,'(1X)')
C
      return
      end
C
C------------------------------------------------------------------------
C
      SUBROUTINE APPEND(PPS)
c      
C  - if spectrum with frequency channel ie. PPS<500 then another spectrum
c    of the same type is appended
c  - if linearized spectrum ie. PPS>50 then another similar spectrum of the
c    same type is merged (in overlapped sections points from the second
c    spectrum are used)
c
      IMPLICIT INTEGER*2 (I-N)
      IMPLICIT REAL*4 (A-H,O-Z)
      PARAMETER (NMAXPT=32000, maxfre=27000, maxcal=200, 
     *  maxspl=500, maxvol=500)
c
      COMMON /LINFRE/FSTART,FEND,FINCR
      COMMON /SPEC/freqs(maxfre),ISPEC(NMAXPT),NPTS,ISMALL,ILARGE,
     1       /INFO/COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
     1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
     1 SCANSP,FRMOD,FRAMPL
c
      INTEGER*2 NVKMA,IVKMA
      INTEGER*4 freqs
      REAL*4 PPS
      REAL*8 FSTART,FEND,FINCR,FSNEW,FENEW,FINEW
      CHARACTER FILNAM*25,CDUM*6,CDUM1
      CHARACTER COMENT*72,LAMP*6,SCANSP*6,SAMPLE*20
C
C...Input of spectrum stored in the standard format, skip the header
c   information, the second spectrum is read into ISPEC directly
c   following the first spectrum
C
499   WRITE(*,'(1X//'' NAME OF DATA FILE ?  ''\)')
      READ(*,'(A)',ERR=499)FILNAM
      OPEN(2,FILE=FILNAM,FORM='BINARY',ERR=499)
C     READ(2,ERR=493)COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
C    1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
C    1 SCANSP,PPS,FRMOD,FRAMPL
      READ(2,ERR=493)CDUM,CDUM,CDUM,CDUM,CDUM,CDUM,CDUM,CDUM,CDUM,
     1 CDUM,CDUM,CDUM,IDUM,IDUM,IDUM,IDUM,IDUM,IDUM,
     1 CDUM,RDUM,RDUM,RDUM,CDUM,CDUM,CDUM,CDUM1,CDUM1,RDUM,RDUM,RDUM,
     1 RDUM,CDUM,RDUM,RDUM,RDUM
C
      READ(2,ERR=493)NPTSA,ISMALA,ILARGA
      IF( (NPTSA+NPTS.GT.MAXFRE.AND.PPS.LT.-500.0) .OR.
     1    (NPTSA+NPTS.GT.NMAXPT.AND.PPS.GT.-500.0)     )THEN
        WRITE(*,'(1X/'' **** spectrum for appending too long'',A/)')
     1   CHAR(7)
        GOTO 493
      ENDIF
      IF(ISMALA.LT.ISMALL)ISMALL=ISMALA
      IF(ILARGA.GT.ILARGE)ILARGE=ILARGA
C
      DO 120 N=NPTS+1,NPTS+NPTSA
        READ(2,ERR=493)ISPEC(N)
120   CONTINUE
C
C...Action for merging a linearized spectrum, taken when:
c    - spectrum 2 follows spectrum 1
c    - both spectra have the same point spacing which is contiguous from one
c      to another
C
      IF(PPS.GT.50.)THEN
        READ(2,ERR=493)FSNEW,FENEW,FINEW
c
c...tests for same increment and contiguity 
c
        if(ABS(1.D0-finew/FINCR).GT.1.D-8)goto 200       
        rsteps=(fsnew-fend)/finew
        nsteps=nint(rsteps)
        if(float(iabs(nsteps))-abs(rsteps).gt.1.d-5)goto 200
        if(fsnew-fend.gt.fincr)goto 200
        if(fenew.le.fstart)goto 200
        goto 201
c
200     write(*,'(1x/'' Spectra will not merge - check bounds and'',
     *  '' frequency increments!'',a1/)')char(7)
        goto 493
c
c...if both spectra have a bridging point
c
201     if(nsteps.eq.0)then
          do 202 n=npts+1,npts+nptsa
            ispec(n-1)=ispec(n)
202       continue
          npts=npts+nptsa-1
          fend=fenew
          WRITE(*,203)NPTSA-1,filnam
203       FORMAT(1X/' ....',I5,' Points added from:  ',a)
          goto 493
        endif
c
c...if spectrum 2 only replaces some of spectrum 1
c
        if(fsnew.ge.fstart.and.fenew.le.fend)then
          nn=npts
          nnn=nint((fsnew-fstart)/fincr)
          do 121 n=1,nptsa
            nnn=nnn+1
            nn=nn+1
            ispec(nnn)=ispec(nn)
121       continue
          WRITE(*,122)NPTSA,filnam
122       FORMAT(1X/' ....',I5,' Points replaced from:  ',a)
        endif
        GOTO 493
      ENDIF
C
C...Action for appending spectrum with frequency channel: skip BWO voltages
C
      IF(PPS.LT.0)THEN
        READ(2)NVKMA
        DO 4561 N=1,NVKMA
          READ(2)IVKMA,VKMA
4561    CONTINUE
      ENDIF
C
C...Read the frequency channel (if present)
C
      IF(PPS.LT.-500.)THEN
        DO 710 N=NPTS+1,NPTS+NPTSA
          READ(2,ERR=493)FREQS(N)
710     CONTINUE
      ENDIF
      NPTS=NPTS+NPTSA
c
      WRITE(*,491)NPTSA,NPTS
491   FORMAT(1X/' ....',I5,' Points appended to a total of',I5,
     *  ' points')
C
493   CLOSE(2)
C
      return
      end
C
C------------------------------------------------------------------------
C
      SUBROUTINE SAVFRE(DELTA)
C
      IMPLICIT INTEGER*2 (I-N)
      IMPLICIT REAL*4 (A-H,O-Z)
      PARAMETER (NMAXPT=32000, maxfre=27000, maxvol=500, maxcal=200)
      COMMON /LINFRE/FSTART,FEND,FINCR
      COMMON /CALFRE/XCALPT,YCALPT,NCALPT,COM1,LINEAR
      COMMON /CALVOL/VKM(maxvol),IVKM(maxvol),NVKM,PPS
      COMMON /SPEC/freqs(maxfre),ISPEC(NMAXPT),NPTS,ISMALL,ILARGE
      INTEGER*4 FREQS
      CHARACTER*72 COM1
      CHARACTER FILNAM*25
      REAL*8 XCALPT(maxcal),YCALPT(maxcal),FSTART,FEND,FINCR
      REAL*8 DELTA
C
      IF(DELTA.EQ.0.D0)THEN
        WRITE(*,'(1X/'' **** Run this option AFTER linearisation'')')
        RETURN
      ENDIF
C
507   WRITE(*,577)1,npts
577   FORMAT(1X/'     Spectral point limits are:  ',2i6/
     * '     IXMIN,IXMAX (range of points to be output):  '\)
      READ(*,*,ERR=507)IXMIN,IXMAX
      IF(IXMAX.LE.IXMIN)GOTO 579
      IF(IXMIN.LT.0)GOTO 579
      IF(IXMAX.GT.NPTS)GOTO 579
      IF( DBLE(FREQS(IXMIN+1)-FREQS(IXMIN))/
     *    DBLE(FREQS(IXMAX)-FREQS(IXMAX-1)). LT.0.D0)THEN
        WRITE(*,'(1X/'' **** This is not a spectrum with a single'',
     *    '' gradient'')')
        RETURN
      ENDIF
      GOTO 5578
579   WRITE(*,600)CHAR(7)
600   FORMAT(' **** ERROR in bounds',A/)
      GOTO 507
C
5578  LINEAR=1
      FSTART=DBLE(FREQS(IXMIN))*0.001D0
      FEND=DBLE(FREQS(IXMAX))*0.001D0
      FINCR=(FEND-FSTART)/(IXMAX-IXMIN)
      NCALPT=0
      PPS=ABS(PPS)
      J=0
      DO 1 N=IXMIN,IXMAX
        J=J+1
        ISPEC(J)=ISPEC(N)
1     CONTINUE
C
1302  WRITE(*,'(1X/'' OUTPUT FILE NAME :  ''\)')
      READ(*,'(A)',ERR=1302)FILNAM
      CALL SAVESP(FILNAM)
C
      RETURN
      END
C
C------------------------------------------------------------------------
C
      SUBROUTINE SAVESP(FILNAM)
C
C   This routine saves the recorded data and various descriptive
C   information into a file written in binary format.
C
      IMPLICIT INTEGER*2 (I-N)
      IMPLICIT REAL*4 (A-H,O-Z)
      PARAMETER (NMAXPT=32000, maxfre=27000, maxvol=500, maxcal=200)
      COMMON /CALVOL/VKM,IVKM,NVKM,PPS
      COMMON /LINFRE/FSTART,FEND,FINCR
      COMMON /CALFRE/XCALPT,YCALPT,NCALPT,COM1,LINEAR
      COMMON /SPEC/freqs(maxfre),ISPEC(NMAXPT),NPTS,ISMALL,ILARGE,
     1       /INFO/COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
     1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
     1 SCANSP,FRMOD,FRAMPL
      INTEGER*2 IVKM(maxvol),NVKM,LINEAR
      INTEGER*4 freqs
      REAL*4 VKM(maxvol)
      REAL*8 XCALPT(maxcal),YCALPT(maxcal),FSTART,FEND,FINCR
      CHARACTER FILNAM*25,CHARFL,STAT*7
      CHARACTER*72 COMENT,COM1
      CHARACTER*6 LAMP,SCANSP
      CHARACTER*20 SAMPLE
C
22    STAT='NEW'
26    OPEN(3,FILE=FILNAM,STATUS=STAT,FORM='BINARY',ERR=21)
      WRITE(3)COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
     1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
     1 SCANSP,PPS,FRMOD,FRAMPL
      WRITE(3)NPTS,ISMALL,ILARGE
      DO 20 N=1,NPTS
        WRITE(3)ISPEC(N)
20    CONTINUE
C
      IF(COMENT.EQ.COM1.or.PPS.GT.50.)THEN
        WRITE(3)FSTART,FEND,FINCR
        WRITE(3)NCALPT
        WRITE(3)(XCALPT(I),YCALPT(I),I=1,NCALPT)
      ENDIF
C
      IF(NVKM.NE.0)THEN
        WRITE(3)NVKM
        DO 25 N=1,NVKM
          WRITE(3)IVKM(N),VKM(N)
25      CONTINUE
      ENDIF
C
      IF(PPS.LT.-500.)THEN
        DO 30 N=1,NPTS
          WRITE(3)FREQS(N)
30      CONTINUE
      ENDIF          
C
      CLOSE(3)
      GOTO 23
21    WRITE(*,'(1X,a1/'' THIS FILE ALREADY EXISTS, OVERWRITE (Y/N)?  ''
     * ,$)')CHAR(7)
      READ(*,1300,ERR=21)CHARFL
      IF(CHARFL.EQ.'Y'.OR.CHARFL.EQ.'y')THEN
        STAT='UNKNOWN'
        GOTO 26
      ENDIF
      IF(CHARFL.EQ.'N'.OR.CHARFL.EQ.'n')THEN
27      WRITE(*,'('' New file name:  '',$)')
        READ(*,1300,ERR=27)FILNAM
1300    FORMAT(A)
        GOTO 22
      ENDIF
C
23    RETURN
      END
C
C------------------------------------------------------------------------
C
      SUBROUTINE SPLINE
C
C   DETERMINATION OF CUBIC SPLINE FUNCTION THROUGH A GIVEN SET OF POINTS
C
C     X,Y - vectors containing on input the X and Y coordinates of points
C           in order of increasing X
C A,B,C,D - vectors containing on output coefficients of the cubic spline
C           components
C       N - number of points
C
C  Cubic spline component i is only defined in the region X.i to X.i+1,
C  value of the spline function in this region should be calculated using
C  the reduced coordinate X-X.i
C
      IMPLICIT INTEGER*2 (I-N)
      parameter (maxspl=500)
      REAL*8 X(maxspl),Y(maxspl),A(maxspl),B(maxspl),C(maxspl),
     * D(maxspl),GRA,CA,GRB,CB
      COMMON /SPL/N,X,Y,A,B,C,D,GRA,CA,GRB,CB
C
C   Sort of data points in order of increasing X
C
192   DO 201 I=1,N-1
        J=I
206     J=J+1
        IF(X(J)-X(I))203,204,204
203     XX=X(I)
        X(I)=X(J)
        X(J)=XX
        XX=Y(I)
        Y(I)=Y(J)
        Y(J)=XX
204     IF(J.EQ.N)GOTO 201
        GOTO 206
201   CONTINUE
C
C   Spline function
C
      DO 6 I=1,N
        A(I)=0.
        B(I)=0.
        C(I)=0.0
6     D(I)=0.0
      DO 1 I=1,N-1
1     A(I)=X(I+1)-X(I)
      DO 2 I=2,N-1
        SU=A(I-1)+A(I)
        B(I)=A(I-1)/SU
        C(I)=A(I)/SU
2     D(I)=((Y(I+1)-Y(I))/A(I)-(Y(I)-Y(I-1))/A(I-1))/SU
      B(1)=2.
      DO 3 I=2,N-1
        K=I+1
        B(I)=2.-B(K)*C(I)/B(I-1)
3     D(K)=D(K)-B(K)*D(I)/B(I-1)
      C(N-1)=3.*D(N-1)/B(N-2)
      DO 4 I=1,N-2
        K=N-I
4     C(K)=(3.*D(K)-C(K)*C(K+1))/B(K-1)
      C(1)=0.
      C(N)=0.
      DO 5 I=1,N-1
        K=I+1
        B(I)=(Y(K)-Y(I))/A(I)-A(I)*(C(K)+2.*C(I))/3.
        D(I)=(C(K)-C(I))/3./A(I)
5     A(I)=Y(I)
C
      GRA=B(1)
      CA=Y(1)-GRA*X(1)
      T=X(N)-X(N-1)
      GRB=B(N-1)+2.*C(N-1)*T+3.*D(N-1)*T*T
      CB=Y(N)-GRB*X(N)
C
      RETURN
      END
C
C------------------------------------------------------------------------
C
      SUBROUTINE VALY(YY,XX,NFUNCT)
C
C   VALY calculates the value of cubic spline function for a given ordinate.
C
C     YY - contains on output the value of the spline function at XX
C     XX - value of X for evaluation of the spline function
C NFUNCT - contains on output the number of spline component used for the cal-
C          culation
C
C in addition:
C
C     A,B,C,D should contain coefficients of the cubic spline components
C           X should contain the X values at the beginnings of the cubic
C             components
C  GRA,CA and GRB,CB should contain the coefficients of the lower and
C             upper linear end components resp.
C
      IMPLICIT INTEGER*2 (I-N)
      PARAMETER (maxspl=500)
c
      REAL*8 X(maxspl),Y(maxspl),A(maxspl),B(maxspl),C(maxspl),
     * D(maxspl),GRA,CA,GRB,CB,XX,YY
      COMMON /SPL/NPTS,X,Y,A,B,C,D,GRA,CA,GRB,CB
C
      IF(XX.GT.X(1).AND.XX.LT.X(NPTS))GOTO 55
      IF(XX.GE.X(NPTS))GOTO 56
      YY=XX*GRA+CA
      NFUNCT=0
      GOTO 57
56    YY=XX*GRB+CB
      NFUNCT=NPTS
      GOTO 57
C
55    DO 27 I=NPTS-1,1,-1
      IF(XX.LT.X(I))GOTO 27
      NFUNCT=I
      GOTO 28
27    CONTINUE
C
28    T=XX-X(NFUNCT)
      YY=A(NFUNCT)+B(NFUNCT)*T+C(NFUNCT)*T*T+D(NFUNCT)*T*T*T
C
57    RETURN
      END
C
C------------------------------------------------------------------------
C
      SUBROUTINE PEPSON(FILNAM)
C      
C...Graphics printer output (for EPSON compatible printer)
C
      IMPLICIT INTEGER*2 (I-N)
      IMPLICIT REAL*4 (A-H,O-Z)
      PARAMETER (NMAXPT=32000, maxfre=27000, maxvol=500, maxcal=200)
C
      COMMON /SPEC/freqs(maxfre),ISPEC(NMAXPT),NPTS,ISMALL,ILARGE,
     *       /INFO/COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
     * LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
     * SCANSP,FRMOD,FRAMPL
     *       /PLOTC/IPLOT(NMAXPT)
     *       /PRETTY/BOLD,NORMAL
      COMMON /LINFRE/FSTART,FEND,FINCR
      COMMON /CALFRE/XCALPT,YCALPT,NCALPT,COM1,LINEAR
      COMMON /CALVOL/VKM(maxvol),IVKM(maxvol),NVKM,PPS
C
      REAL*8 D1,D2,FSTART,FEND,FINCR,XCALPT(maxcal),YCALPT(maxcal),FMARK
      INTEGER*4 INTERM,freqs
      CHARACTER FILNAM*25,BOLD*4,NORMAL*4
      CHARACTER*72 COMENT,COM1
      CHARACTER*6 LAMP,SCANSP
      CHARACTER*20 SAMPLE
C
C...X-axis limits
C
1283  WRITE(*,1420)1,NPTS
1420  FORMAT(1X/'       Overall points range in the spectrum :   ',
     *  I5,',',I5/' X-axis limits in the form IXMIN,IXMAX <CR> :   '\)
      READ(*,*,ERR=1283)NS,NE
      IF(NS.LT.1.OR.NS.GE.NE)GOTO 1283
      IF(NE.GT.NPTS)GOTO 1283
      IF(NS.EQ.1)NS=2
      IF(NE.GT.NPTS-24)NE=NPTS-24
      LPLOT=NE-NS
      RLPLOT=LPLOT/2400.
      IF(RLPLOT.LT.1.2)RLPLOT=LPLOT/2150.
      N=LPLOT
      IF(LINEAR.EQ.1)THEN
        D1=FSTART+(NS-1)*FINCR
        D2=FSTART+(NE-1)*FINCR
      ENDIF
      IK=1
1422  WRITE(*,1290)BOLD,NS,NE,N,RLPLOT,NORMAL
1290  FORMAT(1X/1X,A4,' Plot will be between data points ',I5,
     1' and ',I5,'  and will be'/15X,I5,' printer points and ',F5.2,
     1 ' A4 pages long.',A4/14X,'-1 - change endpoints'/
     1 15X,'0 - OK'/13X,'1-9 - various compression ratios',10X\)
      READ(*,1,ERR=1422)IFLAG
1     FORMAT(I5)
      IF(IFLAG.EQ.-1)GOTO 1283
      IF(IFLAG.EQ.0.AND.N.EQ.LPLOT)GOTO 1421
      IF(IFLAG.EQ.0.AND.N.LT.LPLOT)GOTO 1426
      IF(IFLAG.LT.1.OR.IFLAG.GT.9)GOTO 1422
      IK=IFLAG
      N=LPLOT/IK
      RLPLOT=N/2400.
      IF(RLPLOT.LT.1.2)RLPLOT=N/2150.
      GOTO 1422
C
C...Implementation of compression ratio defined by IK
C
1426  NP=NS-1
      II=0
      INTERM=0
      DO 56 N=NS,NE
        II=II+1
        INTERM=INTERM+ISPEC(N)
        IF(II.LT.IK)GOTO 56
        NP=NP+1
        IPLOT(NP)=INTERM/IK
        II=0
        INTERM=0
56    CONTINUE
      NE=NP
      LPLOT=NE-NS
      GOTO 1423
C
1421  DO 1425 N=NS,NE
        IPLOT(N)=ISPEC(N)
1425  CONTINUE
C
C...Marker spacing
C
1423  IF(LINEAR.EQ.1)THEN
1474    WRITE(*,'(1X/'' Marker spacing /MHz =  '',$)')
        READ(*,'(F10.4)',ERR=1474)FMARK
        WRITE(*,'('' .... marker will be printed every'',F5.0,
     * '' printer points, OK ?  '',$)')FMARK/(FINCR*IK)
        READ(*,1,ERR=1474)IFLAG
        IF(IFLAG.NE.1)GOTO 1474
      ENDIF
C
C...Y-axis scaling
C
1473  WRITE(*,1401)
1401  FORMAT(1X/' DO YOU WANT TO SELECT THE Y-AXIS RANGE ',
     1 'MANUALLY ?   '\)
      READ(*,1,ERR=1473)IFLAG
      IF(IFLAG.NE.1)GOTO 1405
1404  WRITE(*,1402)
1402  FORMAT(1X/' Y-axis limits in the form IYMIN,IYMAX <CR> :   ',$)
      READ(*,*,ERR=1404)MIN,MAX
      IF(MIN.GT.MAX)GOTO 1404
      GOTO 1406
C
C...Scale spectrum into range 1-800 for printing
C
1405  MIN=32000
      MAX=-32000
      DO 1287 N=NS,NS+LPLOT
        IF(IPLOT(N).GT.MAX)MAX=IPLOT(N)
        IF(IPLOT(N).LT.MIN)MIN=IPLOT(N)
1287  CONTINUE
1406  SCALE=799./(MAX-MIN)
      IYSHFT=1
      DO 1286 N=NS,NS+LPLOT
        INTERM=(IPLOT(N)-MIN)*SCALE
        IF(INTERM.LT.0)INTERM=0
        IF(INTERM.GT.799)INTERM=799
        IPLOT(N)=INTERM+IYSHFT
1286  CONTINUE
C
C...table of specifications (in condensed draft mode)
C
      OPEN(4,FILE='PRN',STATUS='NEW')
      WRITE(4,'(1X,A)')CHAR(15)
      WRITE(4,1650)COMENT,SAMPLE,SAMPRE,FILNAM,GAIN,PHASE,IHOUR,IMIN,
     * ISEC,IDAY,IMON,IYEAR,TIMEC,LAMP,VKMST,VKMEND,GRID
      WRITE(4,1651)NPTS,ISMALL,
     * ILARGE,FRMOD,FRAMPL,PPS,SCANSP,NS,NE,IK,MIN,MAX
      IF(LINEAR.EQ.1)WRITE(4,1652)D1,D2,FMARK
      WRITE(4,1653)
1650  FORMAT(2X,72(1H_)/1X,'I',T61,'IF/MHz  N3,NK I'
     * /1X,'I',A,T75,'I'/' I',T75,'I'
     * /1X,'I   Sample: ',A,
     *    T40,'    pressure/mTorr: ',F5.1,T75,'I'
     * /1X,'I Filename: ',A,
     *    T40,'gain/mV, phase/deg: ',F5.2,',',F5.0,T75,'I'
     * /1X,'I Recorded: ',I2,':',I2,':',I2,I4,'/',I2,'/',I4,
     *    T40,'   time constant/s: ',F5.3,T75,'I'
     * /1X,'I     Lamp: ',A,',',F8.2,'-Vkm/mV-',F7.2,
     *    T45,  '       grid/V: ',F5.0,T75,'I')
1651  FORMAT(
     *  1X,'INo of pts: ',I5,I7,'<y<',I6,
     *    T40,'        modulation: ',F4.1,'kHz,',F4.1,'mV',T75,'I'
     * /1X,'I     scan: ',F6.1,' points/s and speed of ',A,T75,'I'
     * /1X,'I',T75,'I'
     * /1X,'I PRINTOUT: x range (',I6,',',I6,')',
     *    T40,'       X-compression: ',I2,T75,'I'
     * /1X,'I           y range (',I6,',',I6,')',T75,'I')
1652  FORMAT(1X,
     *     'I     frequency/MHz (',F9.2,',',F9.2,')',T46,
     *             'Marker spacing:',F7.2,' MHz',T75,'I')
1653  FORMAT(1X,'I',72(1H_),'I'///)
c
c...cancel condensed printing if not a linearized spectrum
c
      IF(LINEAR.NE.1)WRITE(4,'(1X,A)')CHAR(18)
      CLOSE(4)
C
C...Print spectrum
C
      IF(LINEAR.EQ.1)THEN
        CALL PRNSPM(NS,LPLOT,D1,D2,FMARK,IK)
      ELSE
        CALL PRNSPE(NS,LPLOT)
      ENDIF
C
      RETURN
      END     
C
C_____________________________________________________________________________
C
      SUBROUTINE PRNSPE(NSTART,LPLOT)
C
C     PRINTOUT OF SPECTRUM IN PRINTER GRAPHICS MODE WITH FREQUENCY AXIS
C     ALONG THE LENGTH OF THE PAPER
C
C   Vertical axis resolution (across the paper) is 120 characters per inch -
C      double density graphics
C   Horizontal resolution (along the paper) is 216 characters  per inch - 
C      generated by the routine by interleaving three passes of the printhead
C   Precautions are taken against sending the EOF character (code 26) which
C      causes the printing to crash - this can take place on sending either N1
C      or the bit pattern
C
C      NSTART - first spectral point in IPLOT for plotting
C       LPLOT - length of spectrum to be plotted 
C        LINE - contains graphics vectors for three interleaved passes of
C               the printhead
C        NMIN - records the lowest byte in LINE that has been altered
C        NMAX - records the highest byte in LINE that has been altered
C
      IMPLICIT INTEGER*2 (I-N)
      PARAMETER (NMAXPT=32000)
      COMMON /PLOTC/IPLOT(NMAXPT)
      INTEGER*2 LINE(800,3),NMIN(3),NMAX(3),NSTART,LPLOT
      CHARACTER ESC,CONTR
c
c...L is width of graphics window for spectrum in horizontal printer dots,
c   it is assumed that values in ISPEC have been prescaled to range 1,L
C
      L=800
C
      ESC=CHAR(27)
      CONTR=CHAR(1)
      OPEN(4,FILE='PRN',STATUS='NEW')
c
c...set line spacing to 1/216 inch
C
      WRITE(4,'(1X,3A1)')ESC,'3',CHAR(1)
C
C...zero output vectors
C
      DO 5 NLINE=1,3
        NMIN(NLINE)=L
        NMAX(NLINE)=1
        DO 5 N=1,L
          LINE(N,NLINE)=0
5     CONTINUE
C
C...Loop through all points in groups of 24 at a time
C
      DO 8 NPT=NSTART+1,NSTART+LPLOT-24,24
        NSHIFT=1
        NB=IPLOT(NPT+24)
        NMIN(1)=NB
        NMAX(1)=NB
        NLINE=1
C
C...Set up the three output vectors in LINE by splitting each group of
C   24 points among three staggered passses of eight pins of the printhead
C   so that:   first pass points - 1  4  7  10  13  16  19  22 
C             second pass points -  2  5  8  11  14  17  20  23
C              third pass points -   3  6  9  12  15  18  21  24
C   The loop is in descending order from the point above the last plotted
C   spectral point to the point preceding the first plotted spectral
C   point where:
C
C     NB - value currently considered N'th spectral point
C     NA - value of the next higher (N+1)'st spectral point
C     NM - midpoint of NA and NB used for apportioning segments of NA-NB line
C          among columns of pixels
C  I1,I2 - pixel limits of pixels associated with line NA-NB to be set in
C          the current pixel column
C  I3,I4 - pixel limits of pixels associated with line NA-NB to be set in
C          the next higher pixel column
C NSHIFT - takes values 1,2,4,8,... and on addition to an element of
C          LINE sets a given bit 
C
        DO 1 N=NPT+23,NPT-1,-1
          NA=NB
          NB=IPLOT(N)
          NM=(NA+NB)/2
          I1=NA
          I2=NM
          I3=NM
          I4=NB
          IF(NB.GT.NA)GOTO 40
          I1=NM
          I2=NA
          I3=NB
          I4=NM
40        IF(N.EQ.NPT+23)GOTO 22
C
C... part of line NA-NB at the higher point (NA) - this is skipped for
C    the first value of N ie. NPT+23
C
          IF(I1.LT.NMIN(NLINE))NMIN(NLINE)=I1
          IF(I2.GT.NMAX(NLINE))NMAX(NLINE)=I2
          DO 2 NN=I1,I2
            IF(LINE(NN,NLINE).LT.NSHIFT)
     *        LINE(NN,NLINE)=NSHIFT+LINE(NN,NLINE)
2         CONTINUE
          IF(N.LT.NPT)GOTO 25
          NLINE=NLINE+1
          IF(NLINE.LE.3)GOTO 22
          NLINE=1
          NSHIFT=NSHIFT*2
C
C...part of line NA-NB at the current point (NB) - this is skipped for the
C   last value of N ie. NPT-1
C
22        IF(I3.LT.NMIN(NLINE))NMIN(NLINE)=I3
          IF(I4.GT.NMAX(NLINE))NMAX(NLINE)=I4
          DO 3 NN=I3,I4
            LINE(NN,NLINE)=NSHIFT+LINE(NN,NLINE)
3         CONTINUE
1       CONTINUE
C
C...Print the output vector for the group of 24 points - in three passses
C   If N1 is 26 then increase the length of the vector by one (zero) byte
C
25      DO 30 NLINE=3,1,-1
12        N2=NMAX(NLINE)/256
          N1=NMAX(NLINE)-256*N2
          IF(N1.NE.26)GOTO 37
          NMAX(NLINE)=NMAX(NLINE)+1
          GOTO 12
C
C...set printer to double density graphics mode (120pts/inch) to print
C   graphics line N1+N2*256 dots long
C
37        WRITE(4,'(1X/1X,5A1\)',ERR=15)ESC,'*',CONTR,CHAR(N1),CHAR(N2)
C
C...send the component bytes for each pass - if any byte is 26 (11010), it
C   is reset to 10 (1010 - leading bit reset to zero)
C
          NRESET=0
75        DO 14 N=1,NMAX(NLINE)
            IF(LINE(N,NLINE).EQ.26)LINE(N,NLINE)=10
            WRITE(4,'(A1\)',ERR=50)CHAR(LINE(N,NLINE))
14        CONTINUE
          goto 30
c
c...rescue loop in case of problems
50        IF(NRESET.GT.0)GOTO 53
          write(*,51)'in sending graphics line in loop 14',char(7)
51        format(1x/' ***** PRINTER PROBLEM ',2a,':'/
     *     17x,'1) switch printer off and on to reset'/
     *     17x,'2) press ENTER to resume            ..... ',$) 
55        read(*,'(i1)',err=52)nn
52        nreset=nreset+1
          WRITE(4,'(1X,3A1)',err=53)ESC,'3',CHAR(1)
          WRITE(4,'(1X/1X,5A1\)',ERR=53)ESC,'*',CONTR,CHAR(N1),
     *         CHAR(N2)
          goto 75
c
53        if(nreset.ge.3)then
            write(*,'(1x/'' ???? Too bad, I give up (check '',
     *        '' ERROR.DMP)''/)')
            goto 15
          endif
          WRITE(*,54)
54        format(1x/' ***** Still something wrong - try resetting',
     *     ' again, then ENTER .... ',$)
          goto 55
C
30      CONTINUE
C
C...Advance the paper by 21/216 inch for next three passes and reset parts of
C   vectors set in last three passes
C
        WRITE(4,'(3A1\)',err=56)ESC,'J',CHAR(21)
        GOTO 57
c
c...rescue code
56      nreset=0
        write(*,51)'on advancing by 21/216 inch',char(7)
60      read(*,'(i1)',err=58)nn
58      nreset=nreset+1
        WRITE(4,'(1X,3A1)',err=59)ESC,'3',CHAR(1)
        WRITE(4,'(3A1\)',err=59)ESC,'J',CHAR(21)
        goto 57
59      if(nreset.ge.3)then
          write(*,'(1x/'' ???? This is not getting anywhere - I '',
     *     ''give up !!!!''/)')
          close(4)
          return
        endif
        write(*,54)
        goto 60        
c        
57      DO 36 NLINE=1,3
          DO 6 N=NMIN(NLINE),NMAX(NLINE)
            LINE(N,NLINE)=0
6         CONTINUE
          NMIN(NLINE)=L
          NMAX(NLINE)=1
36      CONTINUE
8     CONTINUE
C
C...Reset line spacing to the default 36/216 = 1/6 inch
C
      WRITE(4,'(1X,3A1)')ESC,'3',CHAR(36)
      CLOSE(4)
      GOTO 16
C
C...Dump in case of problems
C
15    OPEN(3,FILE='ERROR.DMP',STATUS='NEW')
      WRITE(*,70)
70    FORMAT(1X/' **** File ERROR.DMP has been produced on problem'/)
      WRITE(3,71)
71    FORMAT(' The next line contains:'/' NPT, N, NLINE, NMAX(NLINE)',
     * ', N2, N1, (IPLOT(I),I=NPT+8,NPT-1,-1)'/)
      WRITE(3,17)NPT,N,NMAX(NLINE),N2,N1,(IPLOT(I),I=NPT+8,NPT-1,-1)
17    FORMAT(6I4,2X,10I6)
      WRITE(3,72)
72    FORMAT(1X/' The following are the contents of LINE(I,1),'
     * 'LINE(I,2),LINE(I,3)'/' where I is from 1 to NMAX(NLINE):')
      DO 31 NLINE=1,3
        WRITE(*,'(1X)')
        WRITE(3,'(1X,20I4)')(LINE(I,NLINE),I=1,NMAX(NLINE))
31    CONTINUE
      CLOSE(3)
C
16    RETURN
      END
C
C_____________________________________________________________________________
C
C
        SUBROUTINE PRNSPM(NSTART,LPLOT,F1,F2,FMARK,NCOMPR)
C
C     PRINTOUT OF SPECTRUM and FREQUENCY SCALE IN PRINTER GRAPHICS MODE
C     WITH FREQUENCY AXIS ALONG THE LENGTH OF THE PAPER
C
C   Vertical axis resolution (across the paper) is 120 characters per inch -
C      double density graphics
C   Horizontal resolution (along the paper) is 216 characters  per inch - 
C      generated by the routine by interleaving three passes of the printhead
C   Precautions are taken against sending the EOF character (code 26) which
C      causes the printing to crash - this can take place on sending either
C      NN1 or the bit pattern
C
C
      IMPLICIT INTEGER*2 (I-N)
      PARAMETER (NMAXPT=32000, maxcal=200)
      COMMON /PLOTC/IPLOT(NMAXPT)
      COMMON /LINFRE/FSTART,FEND,FINCR
      COMMON /CALFRE/X,Y,NCALPT,COM1,LINEAR
      REAL*8 FSTART,FEND,FINCR,X(maxcal),Y(maxcal)
      INTEGER*2 LINE(840,3),NMIN(3),NMAX(3)
      INTEGER*4 INTERM
      REAL*8 FMARK,F1,F2,FREQ,DELTAF,HDELTA,FMINV,FR,
     1 RN1,RN2
      CHARACTER ESC,CONTR,COM1*72
      L=840
C
      DELTAF=FINCR*NCOMPR
      HDELTA=DELTAF*0.5D0
      FMINV=1.D0/FMARK
C
      ESC=CHAR(27)
      CONTR=CHAR(1)
      OPEN(4,FILE='PRN',STATUS='NEW')
c
c...set linespacing to 1/216 inch
c
      WRITE(4,'(1X,3A1)')ESC,'3',CHAR(1)
C
      DO 5 NLINE=1,3
        NMIN(NLINE)=0
        NMAX(NLINE)=0
        LINE(1,NLINE)=255
        DO 5 N=2,L
          LINE(N,NLINE)=0
5     CONTINUE
C
C...M A I N    L O O P   the spectrum is considered 24 points at a time and 
C       split up into three interleaved passes of the printhead, each of which
C       prints 8 points
C
      FREQ=F1-24*DELTAF
      DO 8 NPT=NSTART+1,NSTART+LPLOT-24,24
        FREQ=FREQ+24*DELTAF
        NSHIFT=1
        mten=0
        NB=IPLOT(NPT+24)
        NMIN(1)=NB
        NMAX(1)=NB
        NLINE=1
C
C   Setting up the output vector LINE for three passes of the printhead
C
        DO 1 N=NPT+23,NPT-1,-1
          NA=NB
          NB=IPLOT(N)
          NM=(NA+NB)/2
          I1=NA
          I2=NM
          I3=NM
          I4=NB
          IF(NB.GT.NA)GOTO 40
          I1=NM
          I2=NA
          I3=NB
          I4=NM
40        IF(N.EQ.NPT+23)GOTO 22
C
C  .... part of line at the first point (NA) - this is skipped for the first
C       value of N ie. NPT+23
C
          IF(I1.LT.NMIN(NLINE))NMIN(NLINE)=I1
          IF(I2.GT.NMAX(NLINE))NMAX(NLINE)=I2
          DO 2 NN=I1,I2
            IF(LINE(NN,NLINE).LT.NSHIFT)
     *           LINE(NN,NLINE)=NSHIFT+LINE(NN,NLINE)
2         CONTINUE
C
C  .... markers
C
          FR=FREQ+(N+1-NPT)*DELTAF
          F1=FR+HDELTA
          F2=FR-HDELTA
          RN1=DINT(F1*FMINV)
          RN2=DINT(F2*FMINV)
          IF(RN1.EQ.RN2)GOTO 61
          MARLEN=4
          IF(RN1.EQ.DINT(RN1/10.)*10.)then
            MARLEN=9
            mten=nline
            ften=fr
          endif
          DO 71 NN=2,2+MARLEN
            LINE(NN,NLINE)=LINE(NN,NLINE)+NSHIFT
71        CONTINUE
C
61        IF(N.LT.NPT)GOTO 25
          NLINE=NLINE+1
          IF(NLINE.LE.3)GOTO 22
          NLINE=1
          NSHIFT=NSHIFT*2
C
C  .... part of line for the second point (NB) - this is skipped for the last
C       value of N ie. NPT-1
22        IF(NMAX(NLINE).EQ.0)GOTO 43
          IF(I3.LT.NMIN(NLINE))NMIN(NLINE)=I3
          IF(I4.GT.NMAX(NLINE))NMAX(NLINE)=I4
          GOTO 44
43        NMIN(NLINE)=I3
          NMAX(NLINE)=I4
44        DO 3 NN=I3,I4
            LINE(NN,NLINE)=NSHIFT+LINE(NN,NLINE)
3         CONTINUE
1       CONTINUE
C
C...Print the output vector, three lines of graphics at 1/216 inch
c   linespacing
C
25      DO 30 NLINE=3,1,-1
12        NN2=(NMAX(NLINE))/256
          NN1=NMAX(NLINE)-256*NN2
          IF(NN1.NE.26)GOTO 37
          NMAX(NLINE)=NMAX(NLINE)+1
          GOTO 12
c
c...Select double density graphics, NN1 and NN2 define the number of bytes of
c   data to follow = n1+n2*256
c
37        WRITE(4,'(1X/1X,5A1\)',ERR=15)ESC,'*',CONTR,CHAR(NN1),
     *                 CHAR(NN2)
          DO 14 N=1,NMAX(NLINE)
            IF(LINE(N,NLINE).EQ.26)LINE(N,NLINE)=10
            WRITE(4,'(A1\)',ERR=15)CHAR(LINE(N,NLINE))
14        CONTINUE
c
c...Write the first six digits of the marker frequency, first breaking the
c   carriage hold of graphics mode so that printer advances 1/216 inch, then
c   backspacing the printer by 1/216 inch and writing the number in condensed
c   draft mode
c
          if(nline.eq.mten)then
            INTERM=NINT(FTEN)
            write(4,'(1X/1X,3a,''  '',I6\)')ESC,'j',char(1),INTERM
          endif
30      CONTINUE
c
c...Advance paper by 21/216 inch to clear the region just printed
c
        WRITE(4,'(3A1\)')ESC,'J',CHAR(21)
c
c...Clear used parts of the output vector
c
        DO 36 NLINE=1,3
          DO 6 N=NMIN(NLINE),NMAX(NLINE)
            LINE(N,NLINE)=0
6         CONTINUE
          line(1,nline)=255
          DO 70 N=2,36
            LINE(N,NLINE)=0
70        CONTINUE
          NMIN(NLINE)=0
          NMAX(NLINE)=0
36      CONTINUE
8     CONTINUE
C
      WRITE(4,'(1X,3A1)')ESC,'3',CHAR(36)
      WRITE(4,'(1X,A)')CHAR(18)
      CLOSE(4)
      GOTO 16
C
15    OPEN(3,FILE='ERROR.DMP',STATUS='NEW')
      WRITE(3,17)NPT,N,NMAX(NLINE),RN2,RN1,(IPLOT(I),I=NPT+8,NPT-1,-1)
17    FORMAT(3I4,2F10.0,2X,10I6)
      DO 31 NLINE=1,3
        WRITE(*,'(1X)')
        WRITE(3,'(1X,20I4)')(LINE(I,NLINE),I=1,NMAX(NLINE))
31    CONTINUE
      CLOSE(3)
C
16    RETURN
      END
C_____________________________________________________________________________
C
C
      SUBROUTINE PROLAN
C      
C   Roland DXY-880A output (the plotter is assumed to be connected
C   to the default printer port)
C
      IMPLICIT INTEGER*2 (I-N)
      IMPLICIT REAL*4 (A-H,O-Z)
      PARAMETER (NMAXPT=32000, maxfre=27000)
      COMMON /SPEC/freqs(maxfre),ISPEC(NMAXPT),NPTS,ISMALL,ILARGE,
     *       /INFO/COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC,
     * LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE,
     * SCANSP,FRMOD,FRAMPL,
     *       /PLOTC/IPLOT(NMAXPT),
     *       /PRETTY/BOLD,NORMAL
      CHARACTER BOLD*4,NORMAL*4
      INTEGER*4 INTERM,freqs
      CHARACTER*72 COMENT
      CHARACTER*6 LAMP,SCANSP
      CHARACTER*20 SAMPLE
C
C...X-axis limits on data set
C
1284  ISTP=10
      WRITE(*,1420)1,NPTS
1420  FORMAT(1X/'       Overall points range in the spectrum :   ',
     *  I5,',',I5/' X-axis limits in the form IXMIN,IXMAX <CR> :   '\)
      READ(*,*,ERR=1284)NS,NE
      IF(NS.LT.1.OR.NS.GE.NE)GOTO 1284
      IF(NE.GT.NPTS)GOTO 1284
      IF(NS.EQ.1)NS=2
      LPLOT=NE-NS
      RLPLOT=LPLOT*0.01
      N=LPLOT
1922  WRITE(*,1990)BOLD,NS,NE,N,RLPLOT,NORMAL
1990  FORMAT(1X/1X,A4,' Plot will be between data points ',I5,
     1' and ',I5,'  and will be'/15X,I5,' printer points and ',F5.2,
     1 ' cm long.',A4/14X,'-1 - change endpoints'/
     1 15X,'0 - OK'/13X,'1-9 - various compression ratios',10X\)
      READ(*,1,ERR=1922)IFLAG
1     FORMAT(I5)
      IF(IFLAG.EQ.-1)GOTO 1284
      IF(IFLAG.EQ.0.AND.N.EQ.LPLOT)GOTO 1921
      IF(IFLAG.EQ.0.AND.N.LT.LPLOT)GOTO 1926
      IF(IFLAG.LT.1.OR.IFLAG.GT.9)GOTO 1922
C
C...X-axis compression ratio
C
      IK=IFLAG
      N=LPLOT/IK
      RLPLOT=N*0.01
      GOTO 1922
C
C...compress spectrum and transfer to plot buffer
C
1926  NP=NS-1
      II=0
      INTERM=0
      DO 1956 N=NS,NE
        II=II+1
        INTERM=INTERM+ISPEC(N)
        IF(II.LT.IK)GOTO 1956
        NP=NP+1
        IPLOT(NP)=INTERM/IK
        II=0
        INTERM=0
1956  CONTINUE
      NE=NP
      LPLOT=NE-NS
      GOTO 1923
C
C...transfer points to plot buffer for uncompressed spectrum
C
1921  DO 1925 N=NS,NE
        IPLOT(N)=ISPEC(N)
1925  CONTINUE
C
C...Y-axis limits on data set
C
1923  WRITE(*,1401)
1401  FORMAT(1X/' DO YOU WANT TO SELECT THE Y-AXIS RANGE ',
     1 'MANUALLY ?   '\)
      READ(*,1,ERR=1923)IFLAG
      IF(IFLAG.NE.1)GOTO 1905
1904  WRITE(*,1402)
1402  FORMAT(1X/' Y-axis limits in the form IYMIN,IYMAX <CR> :   ',$)
      READ(*,*,ERR=1904)MIN,MAX
      IF(MIN.GT.MAX)GOTO 1904
      GOTO 1906
C
C...find Y-axis limits on data for autoscaling
C
1905  MIN=32000
      MAX=-32000
      DO 1987 N=NS,NS+LPLOT
        IF(IPLOT(N).GT.MAX)MAX=IPLOT(N)
        IF(IPLOT(N).LT.MIN)MIN=IPLOT(N)
1987  CONTINUE
C
C...scale spectrum
C
1906  WRITE(*,1950)
1950  FORMAT(1X/' Height of spectrum and baseline offset (cm):  '\)
      READ(*,*,ERR=1906)HCM,YSHFT
      IYSHFT=YSHFT*100.
      SCALE=(HCM*100.)/(MAX-MIN)
      IHT=HCM*100.
      DO 1986 N=NS,NS+LPLOT
        INTERM=(IPLOT(N)-MIN)*SCALE
        IF(INTERM.LT.0)INTERM=0
        IF(INTERM.GT.IHT)INTERM=IHT
        IPLOT(N)=INTERM+IYSHFT
1986  CONTINUE
C
C...plot spectrum
C
      open(2,file='PRN',status='new')
c
      WRITE(2,*)'H'
      WRITE(2,1953)' M',0,IPLOT(NS)
      WRITE(*,'(1X/'' Ready to plot, press ENTER  '',A\)')CHAR(7)
      READ(*,'(I1)')I
      I=0
      IPLX=0
      IPLY=IPLOT(NS)
      DO 1952 N=NS,NS+LPLOT
        I=I+1
        CALL SLWPEN(I,IPLOT(N))
1952  CONTINUE
1953  FORMAT(a,i5,',',i5)
      WRITE(2,1953)' M',I,0
        WRITE(2,*)'H'
      CLOSE(2)
C
      RETURN
      END
C
C_____________________________________________________________________________
C
C
      SUBROUTINE SLWPEN(IX,IY)
C
C   Routine to draw a line on the Roland from the current pen position
C   to point IX,IY
C
C            -  It is assumed that the pen is lowered and positioned
C               at point IPLX,IPLY, after moving IPLX,IPLY are updated
C               so it is not necessary to do this externally
C            -  plotting is done at constant pen speed irrespective of
C               distance between the end points - additional plotter points
C               are interpolated between the end points at spacing ISTP along
C               the line
C
      IMPLICIT INTEGER*2 (I-N)
      IMPLICIT REAL*4 (A-H,O-Z)
      COMMON /PLOTER/IPLX,IPLY,ISTP
C
      DX=IX-IPLX
      DY=IY-IPLY
      IF(ABS(DX).LT.6.E-6)GOTO 3
C
      X=IPLX
      GR=DY/DX
      YINT=IPLY-GR*X
      R=SQRT(DX**2+DY**2)
      NT=R/ISTP
      IF(NT.EQ.0)NT=1
      STPX=DX/NT
      DO 1 N=1,NT
        X=X+STPX
        JY=YINT+GR*X
        JX=X
        WRITE(2,1953)' D',JX,JY
1     CONTINUE
1953  FORMAT(a,i5,',',i5)
      GOTO 5
C
3     NT=1.7*ABS(DY)/ISTP
      IF(NT.EQ.0)NT=1
      STPY=DY/NT
      Y=IPLY
      DO 6 N=1,NT
        Y=Y+STPY
        JY=Y
        WRITE(2,1953)' D',IPLX,JY
6     CONTINUE
C
5     IPLX=IX
      IPLY=IY
C
      RETURN
      END
C_____________________________________________________________________________
C_____________________________________________________________________________