program imstgs_pvwave
c
c-----------------------------------------------------------------------
c                                                                      |
c     this program creates an ascii file of the data to be plotted with|
c     a pv~wave program called 'something'.                            |
c                                                                      |
c----------------------------------------------------------------------|
c     v1.0                                                             |
c     written by     r.l.west     bscc     apr-11-1994                 |
c                    modification of imstgs                            |
c-----------------------------------------------------------------------
c
      integer*2 kode, modblk(8)
c
c-----------------------------------------------------------------------
c     *** initialization ***
c-----------------------------------------------------------------------
c
      lin=1
      lunti=5
      lunto=6
      call deluns (lin,lunti,lunto)
      write (lunto, '('' imstgs_pvwave v1.0'')')
c
c-----------------------------------------------------------------------
c     *** get the input data file ***
c-----------------------------------------------------------------------
c
      call deopen
c
c-----------------------------------------------------------------------
c     *** define how data is to be processed ***
c-----------------------------------------------------------------------
c
  100 continue
      call input (kode, modblk)
      if (kode .eq. -1) go to 500
c
c-----------------------------------------------------------------------
c     *** process the data ***
c-----------------------------------------------------------------------
c
      call process (modblk, iret, kflag)
c
c---------------------------------------------------------------------
c     *** end of file logic ***
c---------------------------------------------------------------------
c
      if (iret .eq. -10) write (lunto,'('' end of file'')')
      if (iret .eq. 1) write (lunto,'('' no more subperiods'')')
      call deffrar (kode)
      if (kode .ne. 0) then
         write (lunto,'('' no more frames'')')
         call rewded
         go to 100
      end if
      go to 100
c
c---------------------------------------------------------------------
c     *** normal exit ***
c---------------------------------------------------------------------
c
  500 continue
      stop
      end

      subroutine input (kode, modblk)
c
c----------------------------------------------------------------------|
c                                                                      |
c     input subroutine for rpatgs                                      |
c                                                                      |
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,4=electrometer    |
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----------------------------------------------------------------------|
c
      character cyorn
c
      integer*2 intp, modblk(8)
c
      common    /interpolate/intp
c
c-----------------------------------------------------------------------
c     *** get time limits ***
c-----------------------------------------------------------------------
c
      call lmssetr (kode)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      call frasetr (kode)
c
c-----------------------------------------------------------------------
c     *** define which head ***
c-----------------------------------------------------------------------
c
      call hedset (ihead)
      modblk(1)=ihead
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      do ii = 2, 6
         modblk(ii) = 0
      end do
c
c-----------------------------------------------------------------------
c     *** get angle range ***
c-----------------------------------------------------------------------
c
      call angset (modblk(7),modblk(8))
c
c-----------------------------------------------------------------------
c     *** see if wish to interpolate ***
c-----------------------------------------------------------------------
c
   25 continue
      type '('' do you wish to interpolate the data (y/n)? '',$)'
      accept '(a1)', cyorn
      if (cyorn .eq. 'y'  .or.  cyorn .eq. 'y') then
         intp = 1
      else if (cyorn .eq. 'n'  .or.  cyorn .eq. 'n') then
         intp = 0
      else
         go to 25
      end if
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end

      subroutine process (modblk, kode, lflag)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      logical   ifirst1, sclim2_flag
c
      integer*2 cnt(32,2), iicx, intp, lkuims(2,32), modblk(8)
c
      integer*4 ftms1(4), ifirst(2), iorboft, iorbtime, iorbyd, iwindow,
     *          istrip_times(180,2), itime(2,13),
     *          irec_yd, irec_ms
c
      real      cts(32,2), ims_data(2,180,32), orbs(4,32)
c
      real*8    ftms(4), stms(4)
c
      common    /comlim/stms,
     *          /deluns/lundat, lunti, lunto,
     *          /fratms/ftms, ndiv,
     *          /interpolate/intp,
     *          /orbdat/itime, orbs,
     *          /orbit/iorboft, iwindow, iorbyd, iorbtime, inum, ifirst,
     *          /scale/lkuims
      common    /blktim/irec_yd, irec_ms
c
      data      ifirst1/.true./
c
c-----------------------------------------------------------------------
c     *** initialize ***
c-----------------------------------------------------------------------
c
      lflag = 0
      do ii = 1, 180
         do jj = 1, 32
            ims_data(1,ii,jj) = -1.0
            ims_data(2,ii,jj) = -1.0
         end do
         istrip_times(ii,1) = 0
         istrip_times(ii,2) = 0
      end do
c
c-----------------------------------------------------------------------
c     *** output frame limits ***
c-----------------------------------------------------------------------
c
      call fraoptr (lunto)
c
c-----------------------------------------------------------------------
c     *** initialize the orbit parameters and ***
c     ***    initialize save information      ***
c-----------------------------------------------------------------------
c
      call orboftr
      do ii = 1, 13
         itime(1,ii) = -1
         itime(2,ii) = -1
         do jj = 1, 4
            orbs(jj,ii) = 0.0
         end do
      end do
c
c-----------------------------------------------------------------------
c     *** if first time in, read first record ***
c-----------------------------------------------------------------------
c
      if (ifirst1) then
         call getmf1 (kode)
         if (kode .eq. -10) go to 3000
         ifirst1 = .false.
      end if
c
c-----------------------------------------------------------------------
c     *** main loop to collect data ***
c-----------------------------------------------------------------------
c
      indiv = 0
      sclim2_flag = .true.
      do idiv = 1, ndiv
         indiv = indiv + 1
c
c-----------------------------------------------------------------------
c        *** initialise working arrays ***
c-----------------------------------------------------------------------
c
         jflag = 0
         do ii = 1, 32
            cts(ii,1) = 0.0
            cts(ii,2) = 0.0
            cnt(ii,1) = 0
            cnt(ii,2) = 0
         end do
c
c-----------------------------------------------------------------------
c        *** within time limits ? ***
c-----------------------------------------------------------------------
c
  200    call lmstssr (kode, i1, i2)
         if (kode) 400, 600, 1600     ! before, during, after
c
c-----------------------------------------------------------------------
c        *** get next record ***
c-----------------------------------------------------------------------
c
  400    call getmf1 (kode)
         if (kode .eq. -10) go to 3000
         go to 200
c
c-----------------------------------------------------------------------
c        *** this record within current time limits so process ***
c-----------------------------------------------------------------------
c
  600    continue
c
c-----------------------------------------------------------------------
c        *** see if time to save orbit parameters ***
c-----------------------------------------------------------------------
c
         call orb_save
c
c-----------------------------------------------------------------------
c        *** collect the data ***
c-----------------------------------------------------------------------
c
         call srtims (cts, cnt, i1, i2, modblk, sclim2_flag, jflag)
c
c-----------------------------------------------------------------------
c        *** see if need more data for this strip ***
c-----------------------------------------------------------------------
c
         if (i2 .eq. 512) go to 400
c
c-----------------------------------------------------------------------
c        *** done with this strip ***
c-----------------------------------------------------------------------
c        *** average the data ***
c-----------------------------------------------------------------------
c
         do ii = 1,32
            do jj = 1, 2
               if (cnt(ii,jj) .le. 0) then
                  cts(ii,jj) = -1.0
               else
                  cts(ii,jj) = cts(ii,jj) / float(cnt(ii,jj))
               end if
            end do
         end do
         if (intp .ne. 0) then
            call intrr1 (cts, 32)
            call intrr1 (cts(1,2), 32)
         end if
c
c-----------------------------------------------------------------------
c        *** put data in storage arrays ***
c-----------------------------------------------------------------------
 1000    continue
         istrip_times(idiv,1)=stms(2)/1000.0d0
         istrip_times(idiv,2)=stms(4)/1000.0d0
         do ii=1,32
            ims_data(1,idiv,ii)=cts(ii,1)
            ims_data(2,idiv,ii)=cts(ii,2)
         end do
c
c-----------------------------------------------------------------------
c      *** get time of next strip ***
c-----------------------------------------------------------------------
c
 1600    call lmsnxtr (kode)
         if (kode .ne. 0) go to 3000
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      end do     ! idiv = 1, ndiv
c
c-----------------------------------------------------------------------
c     *** process last division ***
c-----------------------------------------------------------------------
c
 3000 continue
      do ii = 1,32
         do jj = 1, 2
            if (cnt(ii,jj) .le. 0) then
               cts(ii,jj) = -1.0
            else
               cts(ii,jj) = cts(ii,jj) / float(cnt(ii,jj))
            end if
         end do
      end do
      if (intp .ne. 0) then
         call intrr1 (cts, 32)
         call intrr1 (cts(1,2), 32)
      end if
      istrip_times(idiv,1)=stms(2)/1000.0d0
      istrip_times(idiv,2)=stms(4)/1000.0d0
      do ii=1,32
         ims_data(1,idiv,ii)=cts(ii,1)
         ims_data(2,idiv,ii)=cts(ii,2)
      end do
c
c-----------------------------------------------------------------------
c     *** output data array to output file ***
c-----------------------------------------------------------------------
c
      open (unit = 2, file = 'imstgs_pvwave.dat', access = 'sequential',
     *      form = 'formatted', status = 'new', recordtype = 'variable',
     *      disp = 'keep', err = 3100)
      go to 3200
 3100 continue
      type '('' error opening output file'')'
      return
 3200 continue
      ftms1(1)=ftms(1)
      ftms1(2)=ftms(2)/1000.0d0
      ftms1(3)=ftms(3)
      ftms1(4)=ftms(4)/1000.0d0
      iicx = 10
      write (2,*) (modblk(i),i=1,8),iicx,(ftms1(i),i=1,4)
      write (2,*) inum,iorboft,ifirst(1),ifirst(2)
      do i=1,inum
         write (2,*) (itime(j,i),j=1,2),(orbs(j,i),j=1,4)
      end do
      ndif=32
      write (2,*) ndiv,ndif
      write (2,*) (lkuims(1,i),i=1,ndif)
      write (2,*) (lkuims(2,i),i=1,ndif)
      do i=1,ndiv
         write (2,*) (ims_data(1,i,j),j=1,ndif)
      end do
      do i=1,ndiv
         write (2,*) (ims_data(2,i,j),j=1,ndif)
      end do
      close (2)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end

      subroutine srtims (cts, cnt, i1, i2, modblk, sclim2_flag, jflag)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      logical   sclim2_flag
c
      integer*2 cnt(32,2), iat(8), iaz(8), ical(8), icde(512,2), 
     *          idata(2812), ifxs(8), iit(8), iiz(8), imdf(8), imsf(8), 
     *          imsh(32,2,8), irpa(32,2,8), irpac(8), irpf(8), itr(8), 
     *          iw7(16,8), izms(512), jcr(512,2), jmsh(512), jrpa(512), 
     *          modblk(8)
c
      real      cts(32,2)
c
      common    /maf1/idata
c
      equivalence (idata(125), iw7(1,1)),
     *            (idata(253), jcr(1)), 
     *            (idata(1789), icde(1)), 
     *            (imsh(1,1,1), jmsh(1)), 
     *            (irpa(1,1,1), jrpa(1))
c
c-----------------------------------------------------------------------
c     *** get the reference angle ***
c-----------------------------------------------------------------------
c
      call refang (degsam, ramang)
c
c-----------------------------------------------------------------------
c     *** get the word 7 values ***
c-----------------------------------------------------------------------
c
      call w7flgs (iw7, imdf, itr, iit, iiz, ical, ifxs, irpac, iat,
     *             iaz, imsf, irpf)
c
c-----------------------------------------------------------------------
c     *** return if not in mass scan mode ***
c-----------------------------------------------------------------------
c
      if (imsf(1) .eq. 1) then
         jflag = 0
         return
      end if
c
c-----------------------------------------------------------------------
c     *** this loop is executed once every second ***
c-----------------------------------------------------------------------
c
      do ir = 1, 8
c
c-----------------------------------------------------------------------
c        *** this loop every half second to define instrument mode ***
c-----------------------------------------------------------------------
c
         do it = 1, 2
            call expmd1 (imdf(ir), ical(ir), imsf(ir), irpf(ir),
     *                   irpac(ir), irpa(1,it,ir), imsh(1,it,ir))
c
c-----------------------------------------------------------------------
c           *** scale ims and rpa settings to bin numbers
c-----------------------------------------------------------------------
c
            if (sclim2_flag) call sclim2 (imsh(1,it,ir), ibad)
c
c-----------------------------------------------------------------------
c	    *** if ibad = 1, bad memory dump, get next record ***
c-----------------------------------------------------------------------
c
            if (ibad .eq. 1) then
               return 
            else
               sclim2_flag = .false.
            end if
c       
            call sclims (imsh(1,it,ir))
            call sclrpa (irpa(1,it,ir))
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
         end do     ! it = 1, 2
      end do     ! ir = 1, 8
c
c-----------------------------------------------------------------------
c     *** record is good so process it ***
c-----------------------------------------------------------------------
c
      if (modblk(1) .eq. 1) then
c
c-----------------------------------------------------------------------
c        *** radial head processing ***
c-----------------------------------------------------------------------
c
         do ii = i1, i2
            if (jrpa(ii) .gt. 0  .and.  jmsh(ii) .gt. 0) then
               imp = jmsh(ii)
               jflag = 1
               iang = isangl (degsam, ramang, ii)
               if (iang.ge.modblk(7) .and. iang.le.modblk(8)) then
                  do ilh = 1, 2
                     cts(imp, ilh) = cts(imp,ilh) +
     *                               idcodc(jcr(ii,ilh), ic)
                     cnt(imp, ilh) = cnt(imp, ilh) + ic
                  end do
               end if
            end if
         end do     ! ii = i1, i2
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      else
c
c-----------------------------------------------------------------------
c        *** z head processing ***
c-----------------------------------------------------------------------
c
         do ilh = 1, 2
            idet = (modblk(1)-2)*2 + ilh
            call defzms (iiz, iaz, itr, idet, izms)
            do ii = i1, i2
               if (jrpa(ii) .gt. 0  .and.  jmsh(ii) .gt. 0) then
                  imp = jmsh(ii)
                  if (izms(ii) .gt. 0) then
                     ich = izms(ii)
                     jflag = 1
                     iang = isangl (degsam, ramang, ii)
                     if (iang.ge.modblk(7) .and. iang.le.modblk(8)) then
                        cts(imp, ilh) = cts(imp,ilh) +
     *                                  idcodc(icde(ii,ich), ic)
                        cnt(imp, ilh) = cnt(imp, ilh) + ic
                     end if
                  end if
               end if
            end do
         end do     ! ilh = 1, 2
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      end if     ! if (modblk(1) .eq. 1
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end