Program AUTO
#
#
# Build as :
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~ RAT AUTO.RAT >AUTO.FTN   ~
#~ FOR/F4P AUTO/LIST        ~
#~ LINK/NOCHEC/F4P/LAB AUTO ~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# AUTO interactively screens data from an auto-analyser instrument
# and using peak scanning routines (PEAK), an analysis of elemental
# concentrations is made.
#
# AUTO allows for multiple background readings to be taken throughout
# the sample run, without the need of re-measuring standards. Results
# displayed on the VDU are calculated using the first measured background
# readings, but those displayed on the printer use a linear, least squares
# evaluation of all backgrounds entered during the course of the
# run.
#
# AUTO allows parallel input from the VDU and the the Analog recorder. This
# strategy allows commands to be entered while a sample run is in progress.
# Unfortunately FORTRAN does not normally allow this, so a MACRO subroutine
# called FORIO was written to queue a VDU read with interrupts set. 
# This procedure prevents normal FORTRAN I/O to the VDU. 
# In order to use fortran read or writes to the VDU, the FORIO read
# must be killed (call forio('K',rbuf,nchar,rflg).
#
# Two start options can be given either via the parameter starting
# mechanism or via keyboard prompt.
# They are :
#
# /rej  Reject peaks below deignated value
# /sim  For simulated input (64 word blocks) of raw A/D data
# /isw  Switch on diagnostics (only when /sim)
# /rec  For crash recovery. In this case, the file 'AUTO.DAT'
#       is opened , the header information read and positoned
#       to the rear of file. Standards do not have to be re-read
#       as the previous calibration factor is used.
# /dmp  Dump data blocks into file auto.dmp
#
# A recovery file 'auto.dat' is created 'on the fly'. Its format is :
#
# Header record 1
# ---------------
#   5A2 Unit of measure
# F10.5 Calibration factor
#  F7.0 First background reading
#    I2 Replication of samples
# F10.0 Start time (secs)
#   8I1 Channel map
#  F6.4 Corr. Coeff.
#    I8 Peak height top standard
#    I8 Peak height bot standard
#    I8 Peak height of reject sample
#
# Data records
# ------------
# A2    Data type ('BA'- background, 'SA'- sample, 'EN'- end record)
#                 ('RA'- Out of Range warning, 'ST'- Standard )
#                 ('RS'- Start repelication Sequence.)
# F6.0  Position of peak
# F7.0  Total height of peak
#
#
#
logical start,stds,end,pause,eof,backf,bgflg,sim,rec,dmp,rflg,isw
logical lfst,rej,seq,rejflg
byte fn(20),idat(9),icr
integer in(64),units(5),nostds(10),chn(8),mcr(2),pcmd(7),rbuf(40)
real valstd(10),back(10),backx(10),peastd(50),stdx(50),val1(10),val2(10)
real pst(50)
common /abuf/ibuf(1040)
data mcr/3RMCR,3R.../,icr/13/
data pcmd/'pr','i ','= ','au','to','.p','rt'/
data backf/.false./,ibel/7/,sim/.false./,rec/.false./,dmp/.false./
data start/.false./,stds/.true./,end/.false./,eof/.false./,pause/.false./
data rflg/.false./,rej/.false./,rejflg/.false./
#
write (5,('0            ** START OF PROGRAM AUTO **'//))
#
call getops(sim,rec,dmp,isw,rej)
if (.not.(isw.and.sim)) isw=.false.
if (rej) rejflg=.true.
#
stime=secnds(0.)
#
if (dmp) open (unit=3,name='auto.dmp',form='UNFORMATTED')
#
#
if (sim)    # Simulating input ?
{
   write(5,('$Enter the simulation filename ->'))
   read(5,(20a1))fn
   fn(20)=0
   open(unit=1,name=fn,type='OLD',form='UNFORMATTED',readonly)
   write(5,('$Enter the sampling time ->'))
   read(5,(f7.4))time
   ithrsh=.33/time
}
#
#
if (rec)    # Crash recovery ?
{
   open (unit=2,name='auto.dat',type='OLD')
   read (2,(1x,5a2,f10.5,f7.0,i2,f10.0,8i1,f6.4,3i8))_
        units,fac,back(1),nosams,stime,chn,corr,itop,ibot,irejpk
   close (unit=2)
   call labini (1,chn,ist)
   call labwt (idum,time,ibn)
   if (time.gt..03)
   {
      write (5,('0Sampling time TOO SLOW'))
      write (5,(' Type "LAB SET .03" to achieve at least 30sps'))
      write (5,('0          !!! AUTO ABORT !!!',a1,/))ibel
      call exit (133)
   }
   write (5,('0Recovery done - waiting on next sample',a1,/))ibel
   stds=.false.
   ihl=0
   itop=0
   ibot=5000
   irejpk=0
   ibacn=1
   ithrsh=.33/time
}
else
{
   #
   # Question Time
   #
   write(5,('0Start of the AUTO-ANALYSER program'//))
   write(5,('$Enter the unit of measure (e.g. %,mg,ppm) ->'))
   read(5,(5a2))units
   write(5,('0Enter the replication and value of standards'))
   write(5,(' Examples:'))
   write(5,(' 3,10      Implies 3 standards at value 10'))
   write(5,(' 3,22.5    Implies 3 standards at value 22.5'))
   write(5,(' 0,0       Finished entering standards'/))
   nost=0
   do i=1,10
   {  
      write(5,('$->'))
      read(5,*)nostds(i),valstd(i)
      if(nostds(i).eq.0)break
   }
   #
   # Split standards into individual arrays
   #
   for ({i=1;j=1};valstd(i).ne.0.;i=i+1)
   {
      do k=1,nostds(i)
      {
         stdx(j)=valstd(i)
         j=j+1
      }
   }
   nstds=j-1
   no1=nstds
   write(5,('$Enter the replication of samples ->'))
   read(5,*)nosams
   write(5,('$Enter computer channel number ->'))
   read(5,*)chn(1)
   ichn=chn(1)
   #
   #
   # Start reading background
   #
   write(5,('0When the chart is reading background'))
   repeat
   {
      write(5,('      enter the BA command',/))
      write (5,('$AUTO>'))
      read(5,(a2))icmd
      call cnvtuc(2,icmd)
      if(icmd.ne.'BA') write(5,('      ** Invalid Command',a1))ibel
   }
   until (icmd.eq.'BA')
   start=.false.
   end=.false.
   pause=.false.
   ihl=0
   itop=0
   ibot=5000
   ibacn=1
   #
   # Initialilize the LAB unit
   # and wait for background
   #
   if (sim)
   {
      ifst=input(in)
      ba=0.
      do i=1,32 ; ba=ba+in(i)
      back(1)=ba/32.
   }
   else
   {
      call labini(1,chn,ist)
      if(ist.ne.0)
      {
         write(5,('0Error ',i4,' from labini - ABORT !'))ist
         call exit(255)
      }
      write (5,('0Please Wait'))
      call labwt(idum,time,ibn)
      if (time.gt..03)
      {
         write (5,('0Sampling time TOO SLOW'))
         write (5,(' Type "LAB SET .03" to achieve at least 30sps'))
         write (5,('0          !!! AUTO ABORT !!!',a1,/))ibel
         call exit (133)
      }
      ithrsh=.33/time
      call labwt(idum,time,ibn)      #Make sure of complete buffer
      #
      call backg(ichn,ibn,bck)
      x=secnds(stime)
      x=x-16.*time
      backx(1)=x
      back(1)=bck
      write (5,('0Background has been recorded',a1))ibel
   }
   #
   if(rej) write(5,('0Waiting for REJECT peak',A1,/))ibel
   else    write (5,('0Starting Data Collection'/))
   #
}
#
# Wait for data or peak
#
if (.not.sim) call forio('R',rbuf,'AUTO>',5,nchar,rflg)    # Post read on TI:
#
repeat
{
   if(.not.sim)
   {
      do i=1,64 ;in(i)=0
      while (pause)
      {
         start=.false.
         if (end) break 2
         if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)
         write(5,('0!!paused!!',a1,/))ibel
         call forio('R',rbuf,'AUTO>',5,nchar,rflg)
         while (rflg)
         {
            do i=1,10
            {
               if(.not.rflg) break 2
               call labwt(idum,time,idn)
            }
            if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)            
            write(5,('+!! paused!!',a1))ibel
            call forio('R',rbuf,'AUTO>',5,nchar,rflg)
         }
         call cmd(rbuf,pause,end,seq,bgflg,rflg)
      }
      ict=0
      do i=1,32
      {
         call labwt(idum,time,ibn)
         if (.not.rflg)
         {
            call cmd(rbuf,pause,end,seq,bgflg,rflg)
            if (.not.(pause.or.backf.or.bgflg))
            call forio('R',rbuf,'AUTO>',5,nchar,rflg)
         }
         if (seq)
         {
            open (unit=2,name='auto.dat',type='OLD',access='APPEND')
            write(2,(' RS',f6.0,f7.0))xdum,ydum
            close (unit=2)
            seq=.false.
            break
         }
         if(end.or.pause) break
         if (backf)
         {
            call backg(ichn,ibn,bck)
            x=secnds(stime)
            x=x-16.*time
            backx(ibacn)=x
            back(ibacn)=bck
            if (ibacn.gt.1.and..not.stds)
            {
               #
               # Record backgrounds on file
               #
               open (unit=2,name='auto.dat',type='OLD',access='APPEND')
               write (2,(' BA',f6.0,f7.0))x,bck
               close (unit=2)
            }
            if (rflg) call forio('K',rbuf,,,,rflg)
            write (5,('0Background has been recorded',a1,/))ibel
            call forio('R',rbuf,'AUTO>',5,nchar,rflg)
            backf=.false.
          }
         if (bgflg)
         {
            # 
            # Background requested - clean up current buffer
            # and return later
            #
            ibacn=ibacn+1
            if (ibacn.gt.10)
            {
               write (5,('0** Limit of 10 background readings exceeded'))
               write (5,('0          Value ignored !',a1))ibel
            }
            else
            {
               if (rflg) call forio('K',rbuf,,,,rflg)
               write (5,('0Start of background - please wait'))
                backf=.true.
                bgflg=.false.
            }
         bgflg=.false.
         }
          #
          # Average 16 samples and
          # Pack the buffer
          #
          j=ibn*260+5+ichn
          do ii=1,2
          {
             xac=0.
             do jj=1,16
             {
                xac=xac+ibuf(j)
                j=j+8
             }
             ict=ict+1
             in(ict)=xac/16.
          }
      }
      if(dmp) write(3)in
   }
   else
   {
      if(start)
      {
         inpt=input(in)
         if (inpt.eq.-10) break
         if (inpt.ne.0)
         {
            write(5,('0File error ',i4,'  !! Abort !!',a1))inpt,ibel
            stop 133
         }
      }
   }
   repeat
   {
      kct=ict/3
      kct=kct*3
      call peak(isw,in,kct,ithrsh,.not.start,ipos,ihight,istat,izps)
      start=.true.
      is=istat.and.3
      if (is.gt.0)
      {
         #
         # Peak found
         #
         if (ihight.gt.ihl)
         {
            #
            # Collect maximum height and position
            #

            ihl=ihight
            posx=secnds(stime)-64.*time+ipos*time
         }
         if (is.eq.1)
         {
            if (rej)
            {
               irejpk=ihl
               if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)  #kill read
               write (5,('0Reject sample read - waiting for Standards'))
               if (.not.(rflg.or.sim)) call forio('R',rbuf,'AUTO>',5,nchar,rflg) 
               rej=.false.
               start=.false.
               ihl=0
               next
            }
            if (ihl.lt.irejpk)
            {
               if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)  #kill read
               write (5,(' Peak below threshold - rejected',a1))ibel
               if (.not.(rflg.or.sim)) call forio('R',rbuf,'AUTO>',5,nchar,rflg) 
               next
            }
            if (stds)
            {
               #
               # Standards come here
               #
               if(ihl.lt.ibot) ibot=ihl
               if(ihl.gt.itop) itop=ihl
               nstds=nstds-1
               pst(no1-nstds)=ihl
               peastd(no1-nstds)=ihl-bck    # standard peak height
               if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)  #kill read
               write (5,(' standard ',f8.3,' has been read'))stdx(no1-nstds)
               if (.not.(rflg.or.sim)) call forio('R',rbuf,'AUTO>',5,nchar,rflg) 
               if (izps.ne.0.and..not.rejflg)
               {
                  if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)  #kill read
                  write (5,('0FATAL - last standard too low for detection'_
                             ,a1))ibel
                  call exit(113)
               }
               if (nstds.le.0)
               {
                  #
                  # Calibrate the standards
                  #
                  open (unit=3,name='AUTO.PRT')
                  call fit(0,no1,peastd,stdx,yint,fac)
                  if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)  #kill read
                  write (5,('0std-val     std-peak     est-val'))
                  write (5,(' -------     --------     -------'))
                  write (3,('1std-val     std-peak     est-val'))
                  write (3,(' -------     --------     -------'))
                  sx=0. ; sx2=0. ; sy2=0.
                  do j=1,no1
                  {
                     y=stdx(j)
                     yest=peastd(j)*fac
                     sx=sx+y*yest
                     sx2=sx2+(y*y)
                     sy2=sy2+(yest*yest)
                     write (5,(1x,f7.3,3x,f10.3,3x,f9.4))stdx(j),peastd(j),yest
                     write (3,(1x,f7.3,3x,f10.3,3x,f9.4))stdx(j),peastd(j),yest
                  }
                  write (5,('0Equation is :  value = ',f8.4,' X peak'))fac
                  corr=sx/sqrt(sx2*sy2)
                  write (5,('0Correlation Coeff = ',f6.4,/))corr
                  write (5,('0Standards completed - waiting for samples.'))
                  write (3,('0Equation is :  value = ',f8.4,' X peak'))fac
                  corr=sx/sqrt(sx2*sy2)
                  write (3,('0Correlation Coeff = ',f6.4,/))corr
                  close(unit=3)
                  if (.not.(rflg.or.sim)) call forio('R',rbuf,'AUTO>',5,nchar,rflg)  
                  #
                  # Open data file and write header
                  #
                  open (unit=2,name='auto.dat')
                  write (2,(1x,5a2,f10.5,f7.0,i2,f10.0,8i1,f6.4,3i8))_
                        units,fac,back(1),nosams,stime,chn,corr,itop,_
                        ibot,irejpk
                  ibw=1
                  while (ibw.lt.ibacn)
                  {
                      ibw=ibw+1
                      write (2,(' BA',f6.0,f7.0))backx(ibw),back(ibw)
                  }
                  istw=1
                  while (istw.le.no1)
                  {
                     write(2,(' ST',f7.2,f6.0))stdx(istw),pst(istw)
                     istw=istw+1
                  }
                  close (unit=2)
                  stds=.false.
               }
            }
            else
            {
               #
               # Samples come here
               #
               if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)  #kill read
               if (sim) lfst=.false.
               else lfst=.true.
               open (unit=2,name='auto.dat',type='OLD',access='APPEND')
               if (izps.ne.0.and..not.rejflg)
                  #
                  # If zero height peaks
                  #
                  do  i=1,izps
                  {
                      isamp=isamp+1
                      val=0.
                      if (lfst) write (5,('+',a1'RA   ',i4,3X,f11.4,a1))icr,_
                                isamp,val,icr
                      else write (5,(1x,'RA   ',i4,3X,f11.4))isamp,val
                      lfst=.false.
                      write(2,(' RA',f6.0,f7.0))posx,val
                  }
               isamp=isamp+1
               val=(ihl-bck)*fac               # Calibrate sample
               iflg='SA'
               if (ihl.gt.itop.or.ihl.lt.ibot) iflg='RA'
               xhl=ihl
               write(2,(1x,a2,f6.0,f7.0))iflg,posx,xhl
               if (iflg.eq.'SA') iflg='  '
               if (lfst)write (5,('+',a1,a2,'   ',i4,3X,f11.4,a1))_
                  icr,iflg,isamp,val,icr
               else write (5,(1x,a1,a2,'   ',i4,3X,f11.4))icr,iflg,isamp,val   
               close (unit=2)
               if (.not.(rflg.or.sim)) call forio('R',rbuf,'AUTO>',5,nchar,rflg)  
            }
            ihl=0
         }
      }
      is=istat.and.64
   }
   until (is.eq.0)
}
until (end)        
#
if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)  #kill read
open (unit=2,name='auto.dat',type='OLD',access='APPEND')
x=0.
write (2,(' EN',f6.0,f7.0))x,x
write (5,('0Processing the report - please wait',a1))ibel
close (unit=1)
close (unit=2)
open (unit=1,name='auto.dat',type='OLD')
open (unit=2,name='auto.prt',type='OLD',access='APPEND')
#
# Count the number of backgrounds
#
ibak=1
read (1,(1x,5a2,f10.5,f7.0,i2,f10.0,8i1,f6.4,3i8))_
     units,fac,back(1),nosams,stime,chn,corr,itop,ibot,irejpk
repeat
{
   read (1,(1x,a2,f6.0,f7.0))iflg,xpos,yhight
   if (iflg.eq.'EN') break
   if (iflg.eq.'BA')
   {
      ibak=ibak+1
      back(ibak)=yhight
      backx(ibak)=xpos
   }
}
if(ibak.gt.1) call fit(1,ibak,backx,back,yint,bfac)
ict=1
ilin=0
call date(idat(1))
call header(idat,units,corr)
close (unit=1)
open (unit=1,name='auto.dat',type='OLD')
read (1,(1x,a2))i
repeat
{
   v1=0.
   v2=0.
   i=0
   while (i.lt.nosams)
   {
      i=i+1
      read(1,(1x,a2,f6.0,f7.0))ityp,posx,yhl
      if (ityp.eq.'RS')
      {
         i=i-1
         break
      }
      if (ityp.eq.'EN') break 2
      if (ityp.eq.'BA')read(1,(1x,a2,f6.0,f7.0))ityp,posx,yhl
      if (ityp.eq.'ST') next 2
      val1(i)=(yhl-back(1))*fac
      if (val1(i).lt.0.)val1(i)=0.
      v1=v1+val1(i)
      val2(i)=(yhl-(yint+posx*bfac))*fac
      if (val2(i).lt.0.)val2(i)=0.
      v2=v2+val2(i)
   }
   if (i.gt.0)
   {
      inum=i
      xi=i
      avv1=v1/xi
      avv2=v2/xi
      iflg='  '
      if (ityp.eq.'RA')iflg='RA'
      if (ibak.gt.1) write (2,(1x,a2,1x,i4,3x,f11.4,9x,9f11.4))_
                     iflg,ict,avv2,(val2(i),i=1,inum)
      else write (2,(1x,a2,1x,i4,3x,f11.4,9x,9f11.4))_
                 iflg,ict,avv1,(val1(i),i=1,inum)
      ict=ict+1
      ilin=ilin+1
      if (ilin.gt.55)
      {
         ilin=0
         call header(idat,units,corr)
      }
   }
}
#
close (unit=1)
close (unit=2)
call spawn (mcr,,,1,,,,pcmd,14)
write (5,('0File ',4a2,' sent to printer'))(pcmd(j),j=4,7)
write (5,('0To cleanup type "DELETE AUTO.*" then "LOGOUT"',/))
write (5,('0All finished - bye.',a1))ibel
stop
end
#
#
#
#
Subroutine backg(ichn,ibn,back)
#~~~~~~~~~~~~~~~~~~~~~~~~~~
# Calculates the average background over a full block of data
#
# ichn - Supplied A/D channel number
# bck  - back ground value
#
common /abuf/ibuf(1040)
ba=0.
inc=ichn
ib=260*ibn+5
do i=1,32
{
   ba=ba+ibuf(ib+inc)
   inc=inc+8
}
back=ba/32.
return
end
#
#
#
subroutine fit(itype,inum,x,y,yint,fact)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# itype = Input (0=linear thru orig., 1=linear)
# inum  = Input = number of points
# x     = Input = x array
# y     = Input = y array
# yint  = Returned y intercept
# fact  = Returned gradient
#
#
real x(50),y(50)
#
f1=0.;f2=0.;f3=0.;f4=0.;f5=inum
do i=1,inum
{
   f1=f1+x(i)
   f2=f2+x(i)*x(i)
   f3=f3+y(i)
   f4=f4+x(i)*y(i)
}
if(itype.eq.0)
{
   yint=0.
   fact=f4/f2
}
else
{
   yint=(f3*f2-f1*f4)/(f5*f2-f1*f1)
   fact=(f5*f4-f1*f3)/(f5*f2-f1*f1)
}
return
end
#
#
#
function input(in)
#~~~~~~~~~~~~~~~~~
#
integer in(64)
#
call errsns
input=0
read(1,end=10,err=20)in
return
10    input=-10
return
20    call errsns(input)
return
end
#
#
#
Subroutine header(idat,units,corr)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~
byte idat(9)
integer units(5)
#
write (2,('1      Auto-Analyser Report For ',9a1,/))idat
write (2,('         The unit of measure  = ',5a2))units
write (2,('         The correlation coef = ',f6.4,/))corr
write (2,('    Sample  Average-result       Individual-results'))
write (2,('    ------  --------------       ------------------'))
return
end
#
#
#
subroutine getops (sim,rec,dmp,isw,rej)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
logical sim,rec,dmp,isw,rej
byte cmd(80),opts(20)
data ncmds/5/,opts/'/','S','I','M','/','R','E','C','/','D','M','P'_
,'/','I','S','W','/','R','E','J'/
call getmcr(cmd,ilen)
if (ilen.eq.-80)
{
   write(5,('$Options ->'))
   read(5,(Q,80a1))ilen,cmd
}
call cnvtuc(ilen,cmd)
ic=0
for (i=1 ; i.le.4*ncmds ; i=i+4)
{
   ic=ic+1
   call stfind(opts(i),4,cmd,ilen,ipos,jstat)
   if (jstat.eq.0)
   {
      if (ic.eq.1) sim=.TRUE.
      if (ic.eq.2) rec=.TRUE.
      if (ic.eq.3) dmp=.TRUE.
      if (ic.eq.4) isw=.TRUE.
      if (ic.eq.5) rej=.TRUE.
   }
}
return
end
#
#
#
subroutine cmd(rbuf,pause,end,seq,bgflg,rflg)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
#
logical pause,end,bgflg,rflg,seq
integer rbuf(40),cmdtab(7)
data cmdtab/'BA','PA','CO','EN','RS','  ',0/,ibel/7/
#
call cnvtuc(2,rbuf)
call stfind(rbuf,2,cmdtab,14,ipos,jst)
if (jst.eq.0)
{
   if(ipos.eq.1)bgflg=.true.
   if(ipos.eq.3)pause=.true.
   if(ipos.eq.5)pause=.false.
   if(ipos.eq.7)end=.true.
   if(ipos.eq.9)seq=.true.
}
else
{
   if (rflg) call forio('K',rbuf,'AUTO>',5,nchar,rflg)
   write (5,('0Command Error',a1,/))ibel
}
rbuf(1)='  '
return
end