c------------------------------------------------------------------------------
c
c    TRACE - Extraction of spectral traces from bitmaps of scanned spectra
c
c
C   Version 17.I.2009                             ----- Zbigniew KISIEL -----
C
C                          __________________________________________________
C                         | Institute of Physics, Polish Academy of Sciences |
C                         | Al.Lotnikow 32/46, Warszawa, POLAND              |
C                         |                             kisiel@ifpan.edu.pl  |
C                         |     http://info.ifpan.edu.pl/~kisiel/prospe.htm  |
C_________________________/--------------------------------------------------
C
c  PRELIMINARIES:
c
c  1. Scan the spectrum to lossless bitmap: recommended mode is 
c     300 dpi LZW compressed TIFF 
c
c  2. Convert to 8-bit (ie. indexed 256 colour) BMP standard. Convenient conversion
c     is possible with the batch convert mode of Irfanview.
c
c  3. Establish the RGB colours and its range for the two traces of interest.  
c     A useful tool for this is Gimp.
c
c  4. Write the colour values and their tolerances to file TRACE.INP
c
c     Then run TRACE
c
c
c  IMPORTANT:
c
c  1. Compilation with CVF compiler requires the use of /assume:byterecl option,
c     otherwise the RECL=32 is interpreted as 32 long words. Hence compile with:
c
c     df -static -fast -assume:byterecl trace.for
c
c------------------------------------------------------------------------------
C
C  Modification history:
C
c  11.11.08:  adapted from EXTRNEW as part of a project of digitising Nizhnii
c             Novgorod RAD spectra of H2O...HF
c  19.11.08:  batch operation, smoothing, differentiation
c  17.01.09:  generalisation: declarable 1 or 2 traces and Y-axis limits
c
c------------------------------------------------------------------------------
c  Files in use:
c
c    input: 
c
c     'filnam.bmp'         256 colour spectrum                               2
c     'trace.inp'          configuration file for TRACE                      1
c
c    output: 
c
c     'filnam_cal.xy'      averaged xy points for the calibration spectrum   3
c     'filnam_calall.xy'   all xy points the calibration spectrum (SO2)      4
c     'filnam_mol.xy'      xy points for the molecule spectrum               3
c     'filnam_molall.xy'   all xy points for the molecule spectrum           4
c
c------------------------------------------------------------------------------
c
      USE DFLIB
c      
      implicit integer*4 (i-n)
      real*8 ylow,yhigh,diffsc
      character*100 filin,line
      common /colors/icalr,icalg,icalb,icalrt,icalgt,icalbt,
     *               imolr,imolg,imolb,imolrt,imolgt,imolbt
      common /contrl/ioutb
c
      logical*4 fsys
c
      WRITE(*,3344)
3344  FORMAT(1X//'  ',76(1H_)/' |',T79,'|'/
     * ' |         T R A C E  -  Extraction of spectral traces',
     * ' from',T79,'|'/
     * ' |                       bitmaps of scanned spectra',
     * T79,'|'/
     * ' |',76(1H_),'|'/'  version 17.I.2009',T64,'Zbigniew KISIEL'/)
c
c
c...read the steering file TRACE.INP
c
      ioutb=0
      open(1,file='trace.inp',err=50,status='OLD')
25      read(1,'(a)')line
        if(line(1:1).eq.'!')goto 25
        read(line,26)iauto
        read(1,'(a)')line
        read(line,26)icalr,icalg,icalb
        read(1,'(a)')line
        read(line,26)icalrt,icalgt,icalbt
        read(1,'(a)')line
        read(line,26)imolr,imolg,imolb
        if(imolr.ne.0.or.imolg.ne.0.or.imolb.ne.0)ioutb=1
        read(1,'(a)')line
        read(line,26)imolrt,imolgt,imolbt
        read(1,'(a)')line
        read(line,26)nsmo,ntimes
        read(1,'(a)')line
        read(line,26)ndiff,i
        if(ndiff.ne.0.and.i.eq.0)then
          write(*,'(1x//'' ***** ERROR:  Zero scaling parameter '',
     *      ''not allowed for differentiation''//)')
          stop
        else
          diffsc=i
        endif
        read(1,'(a)')line
        read(line,26)iphasa,iphasb
        read(1,'(a)')line
        read(line,26)iminy,imaxy
        ylow=iminy
        yhigh=imaxy
26      format(43x,3i6)
      close(1)
c
c   
c...Single file input
C 
      if(iauto.ne.1)then
1       write(*,2)'         Picture file for analysis:   '
2       format(1x/1x,a,$)      
        read(*,'(a)',err=1)filin
c       
        call SCANSP(filin,nsmo,ntimes,ndiff,iphasa,iphasb,ylow,yhigh,
     *                                                        diffsc)   <-----
      endif
c
c
c...Automatic input: list and go through all .BMP files in the current
c   directory
c
c   Note that the DIR response is variable depending on installation:
c
c008/11/18  11:32 PM         8,956,678 11a.bmp
c008-11-18  23:32         8,956,678 11a.bmp
c
c   This option will only work if the first character of info for a file is 
c   '2', this assuming that the date begins with a full year. 
c   File names containing spaces are also not allowed.
c
      if(iauto.eq.1)then
        fsys=systemqq('dir *.bmp/Od>bmp_spec.lst')
        if(fsys.neqv..TRUE.)then
          write(*,'(1x/a//)')' ***** ERROR: Cannot do dir >bmp_spec.lst'
          stop
        endif        
c
        open(6,file='bmp_spec.lst',status='old')        
c
3       read(6,'(a)',end=4)line
        if(line(1:1).ne.'2')goto 3
        i=len_trim(line)-1
5       if(line(i:i).ne.' ')then
          i=i-1
          goto 5
        endif
        filin=line(i+1:len_trim(line))
        call SCANSP(filin,nsmo,ntimes,ndiff,iphasa,iphasb,ylow,yhigh,
     *                                                        diffsc)   <-----
c       write(*,*)filin
        goto 3
c
4       close(6)
        fsys=systemqq('del bmp_spec.lst')      
      endif
c
c      
      stop
c
50    write(*,'(1x//'' ***** ERROR: cannot read TRACE.INP''//)')
      stop
      end
c
c-----------------------------------------------------------------------------
c
      integer*2 function ihex(j)
      integer*1 j
c
c...elimination of interpretation of bytes greater than #7F as negative
c
      ihex=j
      if(j.lt.0)ihex=int2(j)+256
c
      return
      end
c
c-----------------------------------------------------------------------------
c
      subroutine SCANSP(filin,nsmo,ntimes,ndiff,iphasa,iphasb,
     *                                                ylow,yhigh,diffsc)
c
c    ROUTINE FOR SCANNING a .BMP FILE and identifying traces
c    Can be used in batch mode
c
C   FILIN - name of the file to be processed, part before the extension is
C           is used as generic name for all output files
C    NSMO - length of smoothing window, no smoothing if set to zero
C  NTIMES - the number of repeats for smoothing
C   NDIFF - length of differentiation window, no differentiation if set to zero
C  IPHASA - phase multiplier for trace A (+-1) to produce upward pointing
C           central peak
C  IPHASB - phase multiplier for trace B (+-1)
C
C       Y - vector with y coordinates of acceptable pixels in the current column
c           of the bitmap
c
      implicit integer*4 (i-n)
      real*8 yzero,diffsc,ylow,yhigh
c
      parameter (maxx=10000,maxy=5000,yzero=-11111.d0)
c
      integer*2 i,ihex
      integer*1 b,b1,B2,row(maxx),header(1500),backgr,curcol,
     *          buf(32)
      integer*2 bmap(maxx,maxy),pal(256,4)
      character*100 filin,fcal,fcala,fmol,fmola,fgle,fcals,fmols,fgles,
     *              fcalsd,fmolsd,fglesd
      real*8 y(maxy),ya(maxx),yb(maxx),ysum,yymed,ywork(maxx),ymin,ymax,
     *       aphase,bphase
c
      common /colors/icalr,icalg,icalb,icalrt,icalgt,icalbt,
     *               imolr,imolg,imolb,imolrt,imolgt,imolbt
      common /contrl/ioutb
c
      inperr=0      
      minpxa=100000
      maxpxa=-100000
      minpxb=100000
      maxpxb=-100000
      aphase=iphasa
      bphase=iphasb
c
c...Set up output files
c
      do 100 i=1,len_trim(filin)
        if(filin(i:i).eq.'.')then
          ii=i-1
          goto 101
        endif
100   continue
101   fcal=filin(1:ii)//'_A.xy'
      fcala=filin(1:ii)//'_all_A.xy'
      fmol=filin(1:ii)//'_B.xy'
      fmola=filin(1:ii)//'_all_B.xy'
      fgle=filin(1:ii)//'.gle'
c
      fcals=filin(1:ii)//'_sm_A.xy'
      fmols=filin(1:ii)//'_sm_B.xy'
      fgles=filin(1:ii)//'_sm.gle'
c
      fcalsd=filin(1:ii)//'_dif_A.xy'
      fmolsd=filin(1:ii)//'_dif_B.xy'
      fglesd=filin(1:ii)//'_dif.gle'
c
c_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\
c
c...Open and check the input BMP file
c
      open(2,file=filin,access='direct',RECL=32,err=1,status='OLD')
c
c
c......................BMP f o r m a t..................................
c
C...Structure of a Windows RGB .BMP file:
C
c   __byte number
c  |
c  |  byte length
c  |   | 
c  |   |     BITMAPFILEHEADER (14 bytes)
c
c  1   2     BM                    
c  3   4     file size             
c  7   4     zeroed and not used   
c 11   4     offset of bitmap      
c
c            BITMAPINFOHEADER (40 bytes)
c
c 15   4     size of this block 28hex=40dec
c 19   4     picture width
c 23   4     picture height
c 27   2     no of planes (=1)
c 29   2     bits per pixel (=8 for 256 cols)
c 31   4     compression type (=0 for no compression)     
c 35   4     size of image (may be 0 if not compressed)
c 39   4     zeroed
c 43   4     zeroed
c 47   4     no of colours in palette (100hex for 256 cols.)
c 51   4     no of important colours in palette (100hex for 256 cols.)
c
c            RGBQUAD (4x256 bytes for 256 colours)
c
c 55   4     quad of bytes for red, green, blue, fill which is repeated for
c            no of colours in palette
c
c            BITMAP
c
c  variable  for 256 colours contains byte picture elements in rows, with
c            possible filling of each row so that it occupies a multiple of
c            4 bytes (32 bits).  The rows are written with picture pixels
c            going from left to right, and rows from BOTTOM to TOP.
c_________
C         
c            NOTE that in graphics programs indices for x,y pixel coordinates
c            and index values in the bitmap run from zero.  Thus the 256
c            bitmap indices are 0 to 255
c      
c......................BMP f o r m a t..................................
c
c
c...read the whole header up to the bitmap, assuming it is 1078 bytes long 
c   (ie for BMP with 256 colours)
c
      nr=0
      do 30 i=0,1024,32
        nr=nr+1
        read(2,REC=nr)(header(i+j),j=1,32)
30    continue
      nr=nr+1
      read(2,REC=nr)(header(j),j=1057,1078),(buf(j),j=1,10)
c
c...Check whether Windows .BMP
c
      b=header(1)
      b1=header(2)
      b2=header(15)
      if(b.eq.'B'.and.b1.eq.'M'.and.b2.eq.#28)then
        write(*,4)filin
4       format(1x/1x,78(1h-)/
     * ' Windows .BMP format confirmed for:  ',a/)
      else
        write(*,5)CHAR(7),filin,b,b1,B2
5       format(1x/1x,2a,' is not Windows .BMP'//
     * 10x,'----> first two bytes are: ',2a/
     * 10x,'      and last byte of biSize is ',i4,' and not 40dec)'//)
        close(2)
        inperr=1
        return
      endif
c
c...bitmap offset
c
      b1=header(11)
      b2=header(12)
      i=int2(b2)*256+ihex(b1)                                           <-----
      write(*,7)'         bitmap offset = ',i,b1,b2
      ioffs=i
C
C...BITS PER PIXEL
C
      i=int2(header(30))*256+ihex(header(29))                           <-----
      write(*,7)'        bits per pixel = ',i,header(29),header(30)
      write(*,'(1x)')
7     format(1x,a,i4,10x,z2,',',z2)
      if(i.ne.8)then
        inperr=1
        return
      endif
c
c...UNCOMPRESSED?
c
      if(header(31).ne.0)then
        write(*,'(1x/'' THIS PICTURE IS COMPRESSED --- ''/
     *     '' biCompression bytes (hex) are '',4z2//)')
     *      (header(j),j=31,34)
        inperr=1
        return
      endif
c
c...picture width and height
c
      i=int2(header(20))*256+ihex(header(19))                           <-----
      write(*,7)'         picture width = ',i,header(19),header(20)
      iwidth=i
      iw32=(iwidth/4)*4
      if(iw32.ne.i)iw32=iw32+4
      write(*,7)'         rounded width = ',iw32
c
      i=int2(header(24))*256+ihex(header(23))                           <-----
      write(*,7)'        picture height = ',i,header(23),header(24)
      ihight=i
c
      if(iwidth.gt.maxx.or.ihight.gt.maxy)then
        write(*,'(1x/'' Bitmap limits exceeded: current settings: '',
     *  2i7//)')maxx,maxy
        stop
      endif
           
c
c...set up the pallette: byte order is GRB so reverse it for RGB
c
      nn=54
      do 10 n=1,256
        nn=nn+1
        pal(n,3)=ihex(header(nn))                                       <-----
        nn=nn+1
        pal(n,2)=ihex(header(nn))                                       <-----
        nn=nn+1
        pal(n,1)=ihex(header(nn))                                       <-----
        nn=nn+1
        pal(n,4)=ihex(header(nn))                                       <-----
10    continue      
c
c...dump pallette
c
      open(7,file='h.pal',status='unknown')
      do 40 n=1,256
        write(7,*)n,pal(n,1),pal(n,2),pal(n,3),pal(n,4)
40    continue
      close(7)
c
c
c......................SCAN THE PICTURE..........................
c
c...loop over all rows; background colour is set to colour of first element in
c   first row of picture)
c
      write(*,'(1x/'' READING ROW:'')')    
      ifirst=10
      backgr=buf(1)
      do 20 ii=1,ihight
c
c...read row 
c
        if(mod(ii,10).eq.0)write(*,'(i5,$)')ii
c
        do 221 j=1,ifirst
          row(j)=buf(j)
221     continue          
        ilast=iw32-ifirst-((iw32-ifirst)/32)*32
        do 21 jj=ifirst,iw32-ilast-32,32
          nr=nr+1
          read(2,rec=nr)(row(jj+j),j=1,32)
21      continue
        nr=nr+1
        ifirst=32-ilast
        if(ii.ne.ihight)then
          read(2,rec=nr,err=200)
     *      (row(iw32-j),j=ilast-1,0,-1),(buf(j),j=1,ifirst)
        else
          read(2,rec=nr,err=200)
     *      (row(iw32-j),j=ilast-1,0,-1)
        endif
c
        do 22 n=1,iwidth
          bmap(n,ii)=ihex(row(n))                                       <-----
22      continue
c
20    continue
c
      close(2)
c
c...dump selected column
c
c     open(7,file='hcol.dmp',status='unknown')
c     ncol=
c     do 40 n=ihight,1,-1
c       write(7,*)n,pal(n,1),pal(n,2),pal(n,3),pal(n,4)
c40    continue
c      close(7)
c
c_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\_/\
c
c...Extract the curves from picture in BMAP containing with coordinates
C
c          1,ihight ... iwidth,ihight
c          1,1 ........ iwidth,1
c
c   which correspond to the standard orientation
c   Bitmaps are scaned columnwise
c 
c
c
c...Trace A spectrum (usually the calibrant)....................................
c
      do 51 n=1,maxx
        ya(n)=yzero
51    continue        
c      
      open(3,file=fcal,err=50,status='UNKNOWN')
      open(4,file=fcala,err=50,status='UNKNOWN')
      write(4,35)
35    format('!'/'!    X     Y          Xpix  Ypix       R   G   B'/'!')

c
      minr=icalr-icalrt
      if(minr.lt.0)minr=0
      maxr=icalr+icalrt
      if(maxr.gt.255)maxr=255
      ming=icalg-icalgt
      if(ming.lt.0)ming=0
      maxg=icalg+icalgt
      if(maxg.gt.255)maxg=255
      minb=icalb-icalbt
      if(minb.lt.0)minb=0
      maxb=icalb+icalbt
      if(maxb.gt.255)maxb=255
c
      write(*,
     * '(1x//'' Extracting trace A spectrum, '',
     * ''central colour and limits:'')')
      write(*,29)'   red:  ',icalr,minr,maxr
      write(*,29)' green:  ',icalg,ming,maxg
      write(*,29)'  blue:  ',icalb,minb,maxb
29    format(1x,a,i4,4x,2i4)
c
      do 24 i=1,iwidth
        ysum=0
        npts=0
        do 23 j=1,ihight
          n=bmap(i,j)+1
          ire=pal(n,1)
          igr=pal(n,2)
          ibl=pal(n,3)
          if(ire.lt.minr.or.ire.gt.maxr)goto 23
          if(igr.lt.ming.or.igr.gt.maxg)goto 23
          if(ibl.lt.minb.or.ibl.gt.maxb)goto 23
c
c...write all identified pixels as found
c
          write(4,'(2i6,8x,2i6,4x3i4)')i,j,i-1,ihight-j,ire,igr,ibl
c
          npts=npts+1
          y(npts)=j
          ysum=ysum+j
23      continue 
c
        if(npts.eq.0)goto 24
c
        if(npts.gt.3)then
36        n=npts/2
          if(n*2.eq.npts)then
            yymed=y(n+1)-y(n)
          else
            yymed=0.5d0*(y(n+2)-y(n))
          endif
c
          if(y(2)-y(1).gt.3.d0*yymed)then
            ysum=0
            do 37 nn=1,npts-1
              y(nn)=y(nn+1)
              ysum=ysum+y(nn)
37          continue
            npts=npts-1
            if(npts.gt.3)then
              goto 36
            else
              goto 38
            endif
          endif              
c
          if(y(npts)-y(npts-1).gt.3.d0*yymed)then
            ysum=ysum-y(npts)
            npts=npts-1
            if(npts.gt.3)goto 36
          endif              
c                    
        endif
c
c...write one averaged pixel for a column
c
38      yj=ysum/npts
        j=yj
        write(3,'(i6,f9.2,8x,2i6,4x3i4)')i,yj,i-1,ihight-j
        ya(i)=yj
        if(i.lt.minpxa)minpxa=i
        if(i.gt.maxpxa)maxpxa=i
c        
24    continue
c
      close(3)
      close(4)
c 
c...Trace B spectrum, usually the studied molecule..............................
c
      if(ioutb.ne.1)goto 105
c
      do 52 n=1,maxx
        yb(n)=yzero
52    continue        
c
      open(3,file=fmol,err=50,status='UNKNOWN')
      open(4,file=fmola,err=50,status='UNKNOWN')
      write(4,35)
c
      minr=imolr-imolrt
      if(minr.lt.0)minr=0
      maxr=imolr+imolrt
      if(maxr.gt.255)maxr=255
      ming=imolg-imolgt
      if(ming.lt.0)ming=0
      maxg=imolg+imolgt
      if(maxg.gt.255)maxg=255
      minb=imolb-imolbt
      if(minb.lt.0)minb=0
      maxb=imolb+imolbt
      if(maxb.gt.255)maxb=255
c
      write(*,'(1x//'' Extracting trace B spectrum, '',
     * ''central colour and limits:'')')
      write(*,29)'   red:  ',imolr,minr,maxr
      write(*,29)' green:  ',imolg,ming,maxg
      write(*,29)'  blue:  ',imolb,minb,maxb
c
      do 28 i=1,iwidth
        ysum=0
        npts=0
        do 27 j=1,ihight
          n=bmap(i,j)+1
          ire=pal(n,1)
          igr=pal(n,2)
          ibl=pal(n,3)
          if(ire.lt.minr.or.ire.gt.maxr)goto 27
          if(igr.lt.ming.or.igr.gt.maxg)goto 27
          if(ibl.lt.minb.or.ibl.gt.maxb)goto 27
c
          write(4,'(2i6,8x,2i6,4x,3i4)')i,j,i-1,ihight-j,ire,igr,ibl
c
          npts=npts+1
          y(npts)=j
          ysum=ysum+j
27      continue         
c
        if(npts.eq.0)goto 28
c
        if(npts.gt.3)then
39        n=npts/2
          if(n*2.eq.npts)then
            yymed=y(n+1)-y(n)
          else
            yymed=0.5d0*(y(n+2)-y(n))
          endif
c
          if(y(2)-y(1).gt.3.d0*yymed)then
            ysum=0
            do 42 nn=1,npts-1
              y(nn)=y(nn+1)
              ysum=ysum+y(nn)
42          continue
            npts=npts-1
            if(npts.gt.3)then
              goto 39
            else
              goto 41
            endif
          endif              
c
          if(y(npts)-y(npts-1).gt.3.d0*yymed)then
            ysum=ysum-y(npts)
            npts=npts-1
            if(npts.gt.3)goto 39
          endif              
c                    
        endif
c
c...write one averaged pixel for a column
c
41      yj=ysum/npts
        j=yj
        write(3,'(i6,f9.2,8x,2i6,4x,3i4)')i,yj,i-1,ihight-j
        yb(i)=yj
        if(i.lt.minpxb)minpxb=i
        if(i.gt.maxpxb)maxpxb=i
c
28    continue
c
      close(3)
      close(4)
c
c.......................................................................
c
c   From here on deal with both traces simultaneously
c
c
c...write gle script to allow previewing traces as extracted
c
105   minpx=minpxa
      maxpx=maxpxa
      if(ioutb.eq.1)then
        if(minpxb.lt.minpx)minpx=minpxb
        if(maxpxb.gt.maxpx)maxpx=maxpxb
      endif
      call wrigle(fgle,fcal,fmol,ylow,yhigh,minpx,maxpx)                <-----
c
c
c...List missing x values
c
      write(*,'(1x/'' Trace A missing pixels: '',$)')
      do 56 n=minpxa,maxpxa
        if(ya(n).eq.yzero)write(*,'(i5,$)')n
56    continue
c
      if(ioutb.eq.1)then
        write(*,'(1x//'' Trace B missing pixels: '',$)')
        do 57 n=minpxb,maxpxb
          if(yb(n).eq.yzero)write(*,'(i5,$)')n
57      continue
      endif
c
c
c...Interpolate missing x values for plot A
c      
      last=0
      nzeros=0
      do 53 n=1,maxpxa
c
        if(ya(n).ne.yzero)then
          if(nzeros.eq.0)then
            last=n
            goto 53
          endif
c
          if(nzeros.gt.0)then
            if(last.eq.0)then
              do 55 i=1,n-1
                ya(i)=0.d0
55            continue
              last=n
              goto 53
            endif                              
c
            do 54 nn=last+1,n-1
              ya(nn)=ya(last)+
     *              real(nn-last)*(ya(n)-ya(last))/real(n-last)
54          continue
            nzeros=0
            last=n
            goto 53
          endif
        endif
c
        if(ya(n).eq.yzero)then
          nzeros=nzeros+1
        endif
c
53    continue
c
c.......................................................................
c
c...Interpolate missing x values for plot B
c      
      if(ioutb.ne.1)goto 106
c      
      last=0
      nzeros=0
      do 63 n=1,maxpxb
c
        if(yb(n).ne.yzero)then
          if(nzeros.eq.0)then
            last=n
            goto 63
          endif
c
          if(nzeros.gt.0)then
            if(last.eq.0)then
              do 65 i=1,n-1
                yb(i)=0.d0
65            continue
              last=n
              goto 63
            endif                              
c
            do 64 nn=last+1,n-1
              yb(nn)=yb(last)+
     *              real(nn-last)*(yb(n)-yb(last))/real(n-last)
64          continue
            nzeros=0
            last=n
            goto 63
          endif
        endif
c
        if(yb(n).eq.yzero)then
          nzeros=nzeros+1
        endif
c
63    continue
c
c.......................................................................
c
c...Check whether no missing x values left
c
106   write(*,'(1x//'' Trace A missing pixels: '',$)')
      do 66 n=minpxa,maxpxa
        if(ya(n).eq.yzero)write(*,'(i5,$)')n
66    continue
c
      if(ioutb.ne.1)goto 107
      write(*,'(1x/'' Trace B missing pixels: '',$)')
      do 67 n=minpxb,maxpxb
        if(yb(n).eq.yzero)write(*,'(i5,$)')n
67    continue
c
c.......................................................................
c
c...smooth NTIMES times with smoothing length NSMO (from TRACE.inp)
c
107   if(nsmo.eq.0)goto 108
c
      write(*,73)nsmo,ntimes
73    format(1x//i4,' point smoothing, repeated',i3,' times'/)
c
      n=0
      do 69 nn=minpxa,maxpxa
        n=n+1
        ywork(n)=ya(nn)
69    continue
      call smooth(ywork,n,nsmo,ntimes)                                  <-----
      i=minpxa-1
      do 70 nn=1,n
        i=i+1
        ya(i)=ywork(nn)
70    continue
c
      if(ioutb.eq.1)then
        n=0
        do 71 nn=minpxb,maxpxb
          n=n+1
          ywork(n)=yb(nn)
71      continue
        call smooth(ywork,n,nsmo,ntimes)                                <-----
        i=minpxb-1
        do 72 nn=1,n
          i=i+1
          yb(i)=ywork(nn)
72      continue
      endif
c
c...write data
c
68    open(3,file=fcals,status='unknown')
      do 58 n=minpxa,maxpxa
        write(3,'(i6,f9.2,8x,2i6)')n,ya(n),n-1,ihight-int(ya(n))
58    continue        
      close(3)
c
      if(ioutb.eq.1)then
        open(4,file=fmols,status='unknown')
        do 158 n=minpxb,maxpxb
          write(4,'(i6,f9.2,8x,2i6)')n,yb(n),n-1,ihight-int(yb(n))
158     continue
        close(4)
      endif
c
c...write gle script for previewing smoothed traces
c
      call wrigle(fgles,fcals,fmols,ylow,yhigh,minpx,maxpx)             <-----
c
c.......................................................................
c
c...differentiate with window length NDIFF (from TRACE.inp)
c
c   Note that the scaling parameter DIFFSC is somewhat arbitrary. The value 
c   of 2550 rescales (13pt) derivatives of H2OHF spectra into the pixel
c   range of 300 dpi bitmaps
c
108   if(ndiff.eq.0)return
c
      write(*,74)ndiff
74    format(i4,' point differentiation'/)
c
      n=0
      do 75 nn=minpxa,maxpxa
        n=n+1
        ywork(n)=ya(nn)
75    continue
      call differ(ywork,n,ndiff,1)                                      <-----
      i=minpxa-1      
      ymax=-1.d10
      ymin= 1.d10
      do 76 nn=1,n
        i=i+1
        ya(i)=aphase*ywork(nn)/diffsc
        if(ya(i).gt.ymax)ymax=ya(i)
        if(ya(i).lt.ymin)ymin=ya(i)
76    continue
c
      if(ioutb.eq.1)then
        n=0
        do 77 nn=minpxb,maxpxb
          n=n+1
          ywork(n)=yb(nn)
77      continue
        call differ(ywork,n,ndiff,1)                                    <-----
        i=minpxb-1
c       ymax=-1.d10
c       ymin= 1.d10
        do 78 nn=1,n
          i=i+1
          yb(i)=bphase*ywork(nn)/diffsc
          if(yb(i).gt.ymax)ymax=yb(i)
          if(yb(i).lt.ymin)ymin=yb(i)
78      continue
      endif
c
c...write data
c
      open(3,file=fcalsd,status='unknown')
      do 79 n=minpxa,maxpxa
        write(3,'(i6,f13.3)')n,ya(n)
79    continue        
      close(3)
c
      if(ioutb.eq.1)then
        open(4,file=fmolsd,status='unknown')
        do 179 n=minpxb,maxpxb
          write(4,'(i6,f13.3)')n,yb(n)
179     continue        
        close(4)
      endif
c
c...write gle script for previewing differentiated traces
c
      call wrigle(fglesd,fcalsd,fmolsd,ymin,ymax,minpx,maxpx)           <-----
c
      return
c
c.......................................................................
c
c...error conditions
c
200   write(*,'(1x//'' ***** ERROR: Problem reading the bitmap:''//
     *     7x,a//)')filin(1:len_trim(filin))
      return
c
50    write(*,'(1x//'' ***** ERROR: Cannot open output file''//)')
      return
c
1     write(*,'(1x//'' ***** ERROR: Cannot open input file:''//7x,a//)')
     * filin(1:len_trim(filin))
c     
      return
      end
c
c-----------------------------------------------------------------------------
c
      subroutine wrigle(fgle,fila,filb,ymin,ymax,minpx,maxpx)
c
      implicit integer*4 (i-n)
      character*100 fgle,fila,filb
      character ysh*12
      real*8 ymin,ymax
      common /contrl/ioutb
c
      bord=0.05d0*(ymax-ymin)
      yshift=3.0*bord
      write(ysh,'(f12.1)')yshift
1     if(ysh(1:1).eq.' ')then
        ysh=ysh(2:len_trim(ysh))
        goto 1
      endif      
      open(1,file=fgle(1:len_trim(fgle)),status='unknown')
c
      write(1,'(a)')'size 29.5 21.0'
      write(1,'(a)')' '
      write(1,'(a)')'lwd=0.03'
      write(1,'(a)')'set lwidth lwd'
      write(1,'(a)')'set join round'
      write(1,'(2a)')'yshift=',ysh
      write(1,'(a)')' '
      write(1,'(a)')'amove -1.5  1.5'
      write(1,'(a)')'begin graph'
      write(1,'(a)')'      nobox'
      write(1,'(a)')'      size 35 21'
      write(1,'(2a)')'      xtitle "\it x\rm\ pixel" hei 0.8 dist 0.4 ',
     *               'font texcmr'
      write(1,'(a,i6,a,i6,a)')'      xaxis min ',minpx-50,
     *          '  max ',maxpx+50,'  dticks  1000 dsubticks   100'
      write(1,'(a)')'      xticks length 0.3'
      write(1,'(a)')'      xsubticks length 0.1'
      write(1,'(a)')'      xlabels hei 0.6 '
      write(1,'(a)')'! '
      write(1,'(a,f14.2,a,f14.2)')'      yaxis min ',ymin-bord,
     *                                      '  max ',ymax+bord
      write(1,'(a)')'      yticks length 0.3'
      write(1,'(a)')'      ysubticks length 0.1'
      write(1,'(a)')'      ylabels hei 0.6'
      write(1,'(2a)')'      ytitle "\it y\rm\ pixel" hei 0.8 ',
     *              'dist .4 font texcmr'
      write(1,'(a)')'! '
      write(1,'(3a)')'      data ',fila(1:len_trim(fila)),' d1=c1,c2'
      write(1,'(a)')'      d1 lstyle 1 lwidth lwd color black'
      write(1,'(a)')'  '
c
      if(ioutb.eq.1)then
        write(1,'(3a)')'      data ',filb(1:len_trim(filb)),' d2=c1,c2'   
c
        n=len_trim(fgle)
        if(fgle(n-4:n-4).eq.'d')then       
          write(1,'(a)')'      let d3=d2+yshift'
          write(1,'(a)')'      d3 lstyle 1 lwidth lwd color red'      
        else
          write(1,'(a)')'      d2 lstyle 1 lwidth lwd color red'
        endif
      endif
c      
      write(1,'(a)')'  '
      write(1,'(a)')'end graph'
c
      close(1)
c
      return
      end
c
c--------------------------------------------------------------------------
c
      subroutine SMOOTH(ydata,npts,nspt,ntimes)
c
c   NSPT - number of points in smoothing polynomial
c  YDATA - buffer containing the data, from 1 to NPTS
c NTIMES - the number of times that the smooth is to be repeated
c
      PARAMETER (maxx=10000,maxsmo=1001)
c
      implicit real*8 (a-h,o-z)
      real*8 ydata(maxx),temp(maxx),spol(maxsmo)
c
C...Least squares smoothing using the Savitsky,Golay method
c
C    For smoothing interval of length 2m+1 the elements of the smoothing
c    (cubic) polynomial are given by:
C
C           3(3m**2 + 3m -1 - 5s**2)
C   c(s) =  ------------------------
C           (2m+1) (2m-1) (2m+3)
C
C   where s runs from -m to +m (T.H.Edwards and P.D.Wilson, Applied
C                               Spectroscopy 28,541-545(1974))
C
C
c...set up coefficients in smoothing polynomial
c
      M=NSPT/2
      T1=3.D0/((2*M+1)*(2.D0*M-1.D0)*(2*M+3))
      T2=3*M*M+3.D0*M-1.D0
      DO 1103 j=1,NSPT
        IS=j-M-1
        SPOL(j)=T1*(T2-5*IS*IS)
1103  CONTINUE
C
c...Smooth NTIMES times a TEMP copy of data in YDATA into YDATA.
c   For long smoothing polynomials only smooth once.
c
      ISTRT=M+1
      IFIN=npts-M
      DO 1104 k=1,ntimes
        do 1105 j=1,npts
          temp(j)=ydata(j)
1105    continue
        DO 543 I=1,npts
          SUM=0.
          DO 544 J=1,NSPT
            IS=J-M-1
            II=I+IS
            IF(II.LT.1)II=-II+1
            IF(II.GT.npts)II=npts-(II-npts-1)
              SUM=SUM+temp(II)*SPOL(J)
544       CONTINUE
          ydata(I)=sum
543     CONTINUE
1104  continue
C
c
      return
      end
c
c--------------------------------------------------------------------------
c
      subroutine DIFFER(ydata,npts,nspt,ntimes)
c
C...Least squares 1st derivative (with smoothing) using cubic-quartic 
c   coefficients from the Savitsky,Golay paper (A.Savitsky, M.J.E.Golay,
c   Analytical Chemistry 36,1627-1639(1964))
c
c   NSPT - number of points in differentiation polynomial
c  YDATA - buffer containing the data, from 1 to NPTS
c NTIMES - the number of times that the differentiation is to be repeated
c
      PARAMETER (maxx=10000,maxsmo=25)
c
      implicit real*8 (a-h,o-z)
      real*8 ydata(maxx),temp(maxx),spol(maxsmo)
      integer*4 pol5(5),pol7(7),pol9(9),pol11(11),pol13(13),pol15(15),
     *          pol17(17),pol19(19),pol21(21),pol23(23),pol25(25)
c
c  01 coefficient for 9 points corrected from 129 to 126 on the basis of all
c  even coefficients for 01
c
c  03 coefficient for 21 points changed from 79504 to 79564 by assuming 
c  typo 6 to 0 (but this is a rather shaky basis)
c
      data pol5/1,-8,0,8,1/
      data pol7/22,-67,-58,0,58,67,-22/
      data pol9/86,-142,-193,-126,0,
     *         126, 193, 142, -86/
      data pol11/300,-294,-532,-503,-296,0,
     *           296, 503, 532, 294,-300/
      data pol13/1133,-660,-1578,-1796,-1489, -832,0,
     *            832,1489, 1796, 1578,  660,-1133/
      data pol15/12922,-4121,-14150,-18334,-17842,-13843, -7506,0,
     *            7506,13843, 17842, 18334, 14150,  4121,-12922/
      data pol17/748,-98,-643,-930,-1002,-902,-673,-358,0,
     *           358,673, 902,1002,  930, 643,  98,-748/
      data pol19/6936,  68,-4648,-7481,-8700,-8574,-8179,-5363,-2816,0,
     *           2816,5363, 8179, 8574, 8700, 7481, 4648,  -68,-6936/
      data pol21/ 84075, 10032,-43284,-78176,-96947,-101900,-95338,
     *           -79564,-56881,-29592,0,
     *            29592, 56881, 79564, 95338,101900,  96947, 78176,
     *            43284,-10032, 84075/
      data pol23/ 3938,  815,-1518,-3140,-4130,-4567,-4530,
     *           -4098,-3350,-2365,-1222,0, 
     *            1222, 2365, 3350, 4098, 4530, 4567, 4130, 
     *            3140, 1518, -815,-3938/
      data pol25/ 30866,  8602, -8525,-20982,-29236,-33754,-35003,
     *           -33450,-29562,-23806,-16649, -8558,0,
     *             8558, 16649, 23806, 29562, 33450, 35003, 33754,
     *            29236, 20982,  8525, -8602,-30866/
C
C...legality checks
C      
      if(nspt.lt.5.or.nspt.gt.25.or.((nspt/2)*2.eq.nspt))then
        write(*,1)nspt
1       format(1x//' SORRY!  Differentiation window width has to be an',
     *   ' ODD number'/
     *             '         between 5 and 25, while you have ',
     *   'specified',i10//)
        stop
      endif
C
C
c...set up coefficients in smoothing polynomial
c
c
      M=NSPT/2
      do 1103 j=1,nspt
        if(nspt.eq. 5)spol(j)=pol5(j)
        if(nspt.eq. 7)spol(j)=pol7(j)
        if(nspt.eq. 9)spol(j)=pol9(j)
        if(nspt.eq.11)spol(j)=pol11(j)
        if(nspt.eq.13)spol(j)=pol13(j)
        if(nspt.eq.15)spol(j)=pol15(j)
        if(nspt.eq.17)spol(j)=pol17(j)
        if(nspt.eq.19)spol(j)=pol19(j)
        if(nspt.eq.21)spol(j)=pol21(j)
        if(nspt.eq.23)spol(j)=pol23(j)
        if(nspt.eq.25)spol(j)=pol25(j)
1103  continue        
C
c...Smooth NTIMES times a TEMP copy of data in YDATA into YDATA.
c   For long smoothing polynomials only smooth once.
c
      ISTRT=M+1
      IFIN=npts-M
      DO 1104 k=1,ntimes
        do 1105 j=1,npts
          temp(j)=ydata(j)
1105    continue
        DO 543 I=1,npts
          SUM=0.
          DO 544 J=1,NSPT
            IS=J-M-1
            II=I+IS
            IF(II.LT.1)II=-II+1
            IF(II.GT.npts)II=npts-(II-npts-1)
              SUM=SUM+temp(II)*SPOL(J)
544       CONTINUE
          ydata(I)=sum
543     CONTINUE
1104  continue
C
c
      return
      end
c
c-----------------------------------------------------------------------------
c-----------------------------------------------------------------------------