program counts_time
c
c-----------------------------------------------------------------------
c counts vs time for any selected head(s) for all species
c uses rims mf1
c
c v1.0 written by     rlwest     bcss     19-apr-1988
c
c rcts,ictn,kflag: (8,3)  1=L/H+
c                         2=H/HE+
c                         3=L/HE+
c                         4=H/O+
c                         5=L/HE++
c                         6=H/O++
c                         7=H/N+
c                         8=L/Molecular Ions
c
c                         1=radial head
c                         2=+z head
c                         3=-z head
c-----------------------------------------------------------------------
c
      integer*4 idata, iheads, modblk(8), iavemax, kflag(8,3)
c
      external close_maf1  !$pragma C( close_maf1 )
      external close_ctm_out  !$pragma C( close_ctm_out )
c
c-----------------------------------------------------------------------
c     *** initialization ***
c-----------------------------------------------------------------------
c
   49 continue
      call input (kode, modblk, iavemax, iheads)
      if (kode .eq. -1) go to 99
c
c-----------------------------------------------------------------------
c     *** processing ***
c-----------------------------------------------------------------------
c
   40 continue
c
c-----------------------------------------------------------------------
c     *** accumulate data (kflag set to 1 if any collected) ***
c-----------------------------------------------------------------------
c
      call proces (modblk, iret, kflag, iavemax, iheads)
      idata = 0
      do i = 1,8
         do j = 1,3
            idata = max(idata, kflag(i,j))
         end do
      end do
      if (idata .eq. 0) go to 42
c
c-----------------------------------------------------------------------
c     *** write header information and data ***
c-----------------------------------------------------------------------
c
      call wrtout (modblk, kflag, iheads, iavemax)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
   42 if (iret .eq. 0) go to 40
c
c-----------------------------------------------------------------------
c     *** end of file logic ***
c-----------------------------------------------------------------------
c
  510 if (iret .eq. -10) write (*, 520)
  520 format (' end of file')
      if (iret .eq. 1) write (*,521)
  521 format (' no more subperiods')
      go to 49
c
c-----------------------------------------------------------------------
c     *** normal exit ***
c-----------------------------------------------------------------------
c
   99 continue
      call close_ctm_out
      call close_maf1
      call exit(0)
      stop
      end



      subroutine input (kode, modblk, iavemax, iheads)
c
c-----------------------------------------------------------------------
c input subroutine for ctstim
c
c the modblk array stores a given data selection, as follows
c   element   functions            values
c      1      head           1=radial,2=+z,3=-z
c      2      ilohi          1=lo mass,2=hi mass
c      3      im1            lower mass setting
c      4      im2            upper mass setting
c      5      irpa1          lower rpa setting
c      6      irpa2          upper rpa setting
c      7      iang1          lower phase angle (+/- 180)
c      8      iang2          upper phase angle (+/- 180)
c-----------------------------------------------------------------------
c
      character ians*1
c
      integer*4 iavemax, iheads, kode, modblk(8)
c
      external open_maf1  !$pragma C( open_maf1 )
      external open_ctm_out  !$pragma C( open_ctm_out )
c
c-----------------------------------------------------------------------
c     *** open the input data file ***
c-----------------------------------------------------------------------
c
      call open_maf1 (ieof)
      call open_ctm_out (ieof)
c
c-----------------------------------------------------------------------
c     *** get start and stop times ***
c-----------------------------------------------------------------------
c
      call lmssetr (kode)
      if (kode .eq. -1) go to 99
c
c-----------------------------------------------------------------------
c     *** get duration of each frame and resolution ***
c-----------------------------------------------------------------------
c
      call frasetr (kode)
c
c-----------------------------------------------------------------------
c     *** define which head ***
c-----------------------------------------------------------------------
c
      modblk(1) = 0
      iheads = 0
   25 continue
      type '('' which head do you wish to process: radial=1, -z=2,''
     *       '' +z=4'')'
      type '('' (if more than one head, sum the numbers)? '',$)'
      accept *, iheads
      if (iheads .eq. 0) then
         go to 25
      else if (iheads .eq. 1) then
         type '('' processing radial head only, correct (y/n) ? '',$)'
      else if (iheads .eq. 2) then
         type '('' processing -z head only, correct (y/n) ? '',$)'
      else if (iheads .eq. 3) then
         type'('' processing radial and -z heads, correct (y/n) ? '',$)'
      else if (iheads .eq. 4) then
         type '('' processing +z head only, correct (y/n) ? '',$)'
      else if (iheads .eq. 5) then
         type'('' processing radial and +z heads, correct (y/n) ? '',$)'
      else if (iheads .eq. 6) then
         type '('' processing -z and +z heads, correct (y/n) ? '',$)'
      else if (iheads .eq. 7) then
         type '('' processing all three heads, correct (y/n) ? '',$)'
      end if
      accept '(a)', ians
      if (ians .ne. 'Y'  .and.  ians .ne. 'y') go to 25
      modblk(2)=0
      modblk(3)=0
      modblk(4)=0
c
c-----------------------------------------------------------------------
c     *** define rpa range ***
c-----------------------------------------------------------------------
c
      call rpaset (modblk(5), modblk(6))
c
c-----------------------------------------------------------------------
c     *** define angle range ***
c-----------------------------------------------------------------------
c
      call angset (modblk(7), modblk(8))
c
c-----------------------------------------------------------------------
c     *** define whether averages or maximums needed ***
c-----------------------------------------------------------------------
c
      call avemax (iavemax)
c
   99 return
      end



      subroutine proces (modblk, kode, kflag, iavemax, iheads)
c
c-----------------------------------------------------------------------
c this routine does all the accumulation for thes current time interval
c-----------------------------------------------------------------------
c
      logical   ifirst
      integer*2 idat(2812)
c
      integer*4 iheads, jchan(8), jmassh(8), jmassl(8), kflag(8,3),
     *          modblk(8), iavemax, jflag, kode, ndiv
c
      integer*4 iydoy, imil
c
      real      rcte(600,8,3), ctn(8,3), cts(8,3), rcts(600,8,3)
c
      real*8    run_times(4), inc, dur, frame_times(4), strip_times(4),
     *          recyd, rect, del_time, div_time, save_time(600)
c
      common    /blktim/iydoy, imil,
     *          /comlim/strip_times,
     *          /comtms/run_times, inc, dur,
     *          /data/save_time, rcts, rcte,
     *          /fratms/frame_times, ndiv
      common    /maf1/idat
c
      data      ifirst/.true./,
     *          jchan/1, 2, 1, 2, 1, 2, 2, 1/,
     *          jmassl/2*3500, 2*810, 2*1690,  980, 410/,
     *          jmassh/2*3730, 2*935, 2*1870, 1050, 520/
c
c-----------------------------------------------------------------------
c    *** zero out arrays ***
c-----------------------------------------------------------------------
c
      do i = 1,600
         do j = 1,8
            do k = 1,3
               rcts(i,j,k) = 0.0
               rcte(i,j,k) = 0.0
            end do
         end do
      end do
      lflag=0
c
c-----------------------------------------------------------------------
c     *** calculate time per division ***
c-----------------------------------------------------------------------
c
      lflag = 0
      if (frame_times(4) .gt. frame_times(2)) then
         stop = frame_times(4)
      else
         stop = frame_times(4) + 86400000000.0d+00
      end if
      del_time = (stop - frame_times(2)) / float(ndiv)
c
c-----------------------------------------------------------------------
c     *** calculate center time of division ***
c-----------------------------------------------------------------------
c
       div_time = strip_times(2) + (del_time / 2.0d+00)
       if (div_time .gt. 86400000000.0d+00)
     *    div_time = div_time - 86400000000.0d+00
c
c-----------------------------------------------------------------------
c     *** define and output new frame limits ***
c-----------------------------------------------------------------------
c
      call deffrar (kode)
c
c-----------------------------------------------------------------------
c     *** get first block in ***
c-----------------------------------------------------------------------
c
      if (ifirst) then
         call getmf1 (kode)
         recyd = iydoy
         rect = float(imil) * 1000.0d0
         if (kode .eq. -10) go to 510
         ifirst = .false.
      end if
      nsamp = dur / 15625.0d0
      jdiv = ndiv / 12
      do 500 idiv = 1,ndiv
         itsamp = 0
c
c-----------------------------------------------------------------------
c     *** initialise working arrays ***
c-----------------------------------------------------------------------
c
         do i = 1,8
            do j = 1,3
              cts(i,j) = 0.0
              ctn(i,j) = 0.0
            end do
         end do
c
c-----------------------------------------------------------------------
c     *** see if within time limits ***
c-----------------------------------------------------------------------
c
   40    call lmstssr (kode, i1, i2)
         if (kode) 10, 200, 30
c
c-----------------------------------------------------------------------
c     *** get next record ***
c-----------------------------------------------------------------------
c
   10    call getmf1 (kode)
         recyd = iydoy
         rect = float(imil) * 1000.0d0
         if (kode .eq. -10) go to 510            ! end of file
         go to 40
c
  200    continue
c
c-----------------------------------------------------------------------
c     *** this block in current limits so process ***
c-----------------------------------------------------------------------
c
         call deffgs                        ! define flags
         call getmod                        ! and instrument mode
c
c-----------------------------------------------------------------------
c     *** process the data for counts panel ***
c-----------------------------------------------------------------------
c
         itsamp = itsamp + (i2 - i1) + 1
c                  radial head
         if (mod(iheads,2) .eq. 1) then
            modblk(1) = 1
            do i = 1,8
               modblk(2) = jchan(i)
               modblk(3) = jmassl(i)
               modblk(4) = jmassh(i)
               if (iavemax .eq. 1) then
                  call getcts (modblk,cts(i,1),ctn(i,1),i1,i2,jflag)
               else
                  call getmax (modblk,cts(i,1),ctn(i,1),i1,i2,jflag)
               end if
               kflag(i,1) = max0(jflag, kflag(i,1))
            end do
         end if
c                  +z head
         if (iheads .gt. 3) then
            modblk(1) = 2
            do i = 1,8
               modblk(2) = jchan(i)
               modblk(3) = jmassl(i)
               modblk(4) = jmassh(i)
               if (iavemax .eq. 1) then
                  call getcts (modblk,cts(i,2),ctn(i,2),i1,i2,jflag)
               else
                  call getmax (modblk,cts(i,2),ctn(i,2),i1,i2,jflag)
               end if
               kflag(i,2) = max0(jflag, kflag(i,2))
            end do
         end if
c                  -z head
         if (iheads.eq.2 .or. iheads.eq.3 .or. iheads.gt.5) then
            modblk(1) = 3
            do i = 1,8
               modblk(2) = jchan(i)
               modblk(3) = jmassl(i)
               modblk(4) = jmassh(i)
               if (iavemax .eq. 1) then
                  call getcts (modblk,cts(i,3),ctn(i,3),i1,i2,jflag)
               else
                  call getmax (modblk,cts(i,3),ctn(i,3),i1,i2,jflag)
               end if
               kflag(i,3) = max0(jflag, kflag(i,3))
            end do
         end if
c
c-----------------------------------------------------------------------
c     *** see if done with this subperiod (averaging interval) ***
c-----------------------------------------------------------------------
c
         if (i2 .eq. 512  .and.  itsamp .lt. nsamp) go to 10
c
c-----------------------------------------------------------------------
c     *** average the data ***
c-----------------------------------------------------------------------
c
         do j = 1,3
            do i = 1,8
               if (ctn(i,j) .gt. 0) then
                  rcts(idiv,i,j) = cts(i,j) / ctn(i,j)
                  rcte(idiv,i,j) = sqrt(cts(i,j)) / ctn(i,j)
               end if
            end do
         end do
         save_time(idiv) = div_time
c
c-----------------------------------------------------------------------
c     *** get next limits ***
c-----------------------------------------------------------------------
c
   30    call lmsnxtr (kode)
         if (kode .ne. 0) go to 510
c
c-----------------------------------------------------------------------
c     *** calculate center time of division ***
c-----------------------------------------------------------------------
c
       div_time = strip_times(2) + (del_time / 2.0d+00)
       if (div_time .gt. 86400000000.0d+00)
     *    div_time = div_time - 86400000000.0d+00
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
  500 continue
c
c-----------------------------------------------------------------------
c     *** done with this time interval so let's write it ***
c-----------------------------------------------------------------------
c
  510 continue
      return
      end



      subroutine wrtout (modblk,kflag,iheads,iavemax)
c
c-----------------------------------------------------------------------
c     *** does the writting ***
c-----------------------------------------------------------------------
c
      character buf1(39), buf2(34), buf3(77)
      character avemax(2)*5, line1*39, line2*34, line3*77
c
      integer*4 ihr, ii, imn, ims, isc, jhr, jmn, jms, jsc,
     *          modblk(8), ndiv
c
      integer*4 itmdm
c
      real*4    rcte(600,8,3), rcts(600,8,3)
c
      real*8    frame_times(4), save_time(600)
c
      data      avemax/'  ave','  max'/
c
      common    /fratms/frame_times, ndiv,
     *          /data/save_time,rcts,rcte
c
      external write_ctm_data  !$pragma C( write_ctm_data )
c
c-----------------------------------------------------------------------
c     *** create title line1 ***
c-----------------------------------------------------------------------
c
      encode (39,100,line1) avemax(iavemax)
  100 format('DE RIMS     counts vs time (v1.0) ',a5)

      do ii = 1,39
         buf1(ii) = line1(ii:ii)
      end do
c
c-----------------------------------------------------------------------
c     *** determine times and create line2 ***
c-----------------------------------------------------------------------
c
      itmdm=frame_times(2)/1000.0d0
      call mshmsm(itmdm,ihr,imn,isc,ims)
      itmdm=frame_times(4)/1000.0d0
      call mshmsm(itmdm,jhr,jmn,jsc,jms)
      iyr=frame_times(1)/1000.0d0
      idn=dmod(frame_times(1),1000.0d0)
      if(iyr.lt.10)iyr=iyr+80
      encode (34,200,line2) iyr,idn,ihr,imn,isc,ims,jhr,jmn,jsc,jms
  200 format (i2,'/',i3.3,i3.2,2('.',i2.2),'.',i3.3,' - ',i2.2,
     *        2('.',i2.2),'.',i3.3)

      do ii = 1,34
         buf2(ii) = line2(ii:ii)
      end do
c
c-----------------------------------------------------------------------
c     *** create line3 ***
c-----------------------------------------------------------------------
c
      e1=modblk(5)*0.050
      e2=modblk(6)*0.050
      encode (77,300,line3) modblk(5),e1,modblk(6),e2,modblk(7),
     *                      modblk(8)
  300 format ('rpa range = ',i4,' (',f4.1,'ev) - ',i4,' (',f4.1,'ev)',
     *        10x,'angle range = ',i4,' to ',i4)

      do ii = 1,77
         buf3(ii) = line3(ii:ii)
      end do
c
c-----------------------------------------------------------------------
c     *** do the writting ***
c-----------------------------------------------------------------------
c
      call write_ctm_data (line1, line2, line3, iheads, ndiv,
     *                     save_time, rcts, rcte)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end



      subroutine resltn (dummy,ndiv)
c
c-----------------------------------------------------------------------
c dummy routine resltn that returns ndiv=180 regardless of dummy value
c remember that resltn is a line plot program.
c-----------------------------------------------------------------------
c
      integer*4 ndiv
      integer*4 dummy
      data      itm/0/
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      if (itm .eq. 0) then
         type*,'note! keep data resolution in mind when ',
     *         'selecting ndiv'
         ndiv=600
         call iask (ndiv,'enter number of division')
         nsav=ndiv
         itm=1
      else
         ndiv=nsav
      end if
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
        return
      end



      subroutine getmax (modblk,cts,ctn,i1,i2,jflag)
c
c-----------------------------------------------------------------------
c this subroutine gets maximum data count for a given head and given 
c mass. cts and ctn are updated with all the samples i1 thru i2, that
c satisfy the modblk requirements.
c
c input arguments are:
c      i1    - start sample index of data to be included
c      i2    - stop sample index. i1, i2 are in range [1,512]
c output arguments, in addition to cts, ctn,
c      jflag - set to 0 if no data added to arrays, otherwise 1
c
c note *** before this routine is called, the information in the
c following common blocks must be defined
c      common      updated by calling
c      ------      ------------------
c      rpamsh      getmod.  contains instrument rpa, ims settings
c      i7flgs      deffgs.  contains instrument mode flags.
c
c v1.0 jfe johnson 13 aug 82
c v1.1 jfej        17 aug 82 added idcodc call
c v1.2 blg         27 jan 83 modified getcts into getmax
c v1.3 rlwest       8 jul 86 incorporated the new version if idcodc
c
c-----------------------------------------------------------------------
c
      integer*2 icde(512,2), idat(2812), jcr(512,2)
c
      integer*4 jmsh(512), jrpa(512), modblk(8)
c
      integer*4 idcodc, ictn
c
      equivalence (icde(1,1),idat(1789)), (jcr(1,1),idat(253))
c
      common    /rpamsh/jrpa,jmsh,
     *          /maf1/idat
c
      jflag=0
c
c-----------------------------------------------------------------------
c first call refang, to define phase angle reference
c-----------------------------------------------------------------------
c
      call refang (degsam,ramang)
c
      if (modblk(1) .eq. 1) go to 401
c
c-----------------------------------------------------------------------
c z head processing
c-----------------------------------------------------------------------
c
      idet=(modblk(1)-2)*2+modblk(2)
      do 430 i=i1,i2
c
c-----------------------------------------------------------------------
c                   check that rpa setting is in the required range
c-----------------------------------------------------------------------
c
         if(jrpa(i).lt.modblk(5).or.jrpa(i).gt.modblk(6))go to 430
c
c-----------------------------------------------------------------------
c                  check that this mass is in the range required
c-----------------------------------------------------------------------
c
         if(jmsh(i).lt.modblk(3).or.jmsh(i).gt.modblk(4))go to 430
c
c-----------------------------------------------------------------------
c                  check this sample is in the angle range required
c-----------------------------------------------------------------------
c
         iang=isangl(degsam,ramang,i)
         if(iang.lt.modblk(7).or.iang.gt.modblk(8))go to 430
c
c-----------------------------------------------------------------------
c                  check we have data for this z head detector
c-----------------------------------------------------------------------
c
         izms=igtzms(idet,i)
         if(izms.le.0)go to 430
c
c-----------------------------------------------------------------------
c                  have data
c-----------------------------------------------------------------------
c
         jflag=1
         if (cts .lt. float(idcodc(icde(i,izms),ictn)))
     *       cts=float(idcodc(icde(i,izms),ictn))
         ctn=1.0
c-----------------------------------------------------------------------
  430 continue
c
c-----------------------------------------------------------------------
c done z head procesing
c-----------------------------------------------------------------------
c
      go to 400
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
  401 continue
c
c-----------------------------------------------------------------------
c radial head processing
c-----------------------------------------------------------------------
c
      do 100 i=i1,i2
c
c-----------------------------------------------------------------------
c                  check that rpa setting is in the required range
c-----------------------------------------------------------------------
c
         if(jrpa(i).lt.modblk(5).or.jrpa(i).gt.modblk(6))go to 100
c
c-----------------------------------------------------------------------
c                  check that this mass is in the range required
c-----------------------------------------------------------------------
c
         if(jmsh(i).lt.modblk(3).or.jmsh(i).gt.modblk(4))go to 100
c
c-----------------------------------------------------------------------
c                  check this sample is in the angle range required
c-----------------------------------------------------------------------
c
         iang=isangl(degsam,ramang,i)
         if(iang.lt.modblk(7).or.iang.gt.modblk(8))go to 100
c
c-----------------------------------------------------------------------
c                  have data
c-----------------------------------------------------------------------
c
         jflag=1
         if (cts .lt. float(idcodc(jcr(i,modblk(2)),ictn)))
     *       cts=float(idcodc(jcr(i,modblk(2)),ictn))
         ctn=1.0
c-----------------------------------------------------------------------
  100 continue
c
c-----------------------------------------------------------------------
c done
c-----------------------------------------------------------------------
c
  400 continue
      return
      end



      subroutine avemax (iavemax)
c
c-----------------------------------------------------------------------
c  this subroutine requests whether average counts or a maximum
c  count is required.  iavemax = 1 for average count
c                             iavemax = 2 for maximum count
c-----------------------------------------------------------------------
c
      integer*4 iavemax
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      write(*,100)
  100 format($,' plot: averages=1, maximums=2 ? ')
      read(*,*,end=99)iavemax
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
   99 call exit(0)
      stop
      end