program imstgs
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      integer*4 kode, modblk(8)
c
      external close_maf1  !$pragma C( close_maf1)
c
c-----------------------------------------------------------------------
c     *** initialization data collection parameters ***
c-----------------------------------------------------------------------
c
      write (12, '(''inside IMSTGS'')')
      call initialize (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 (*,'('' end of file'')')
      if (iret .eq. 1) write (*,'('' no more subperiods'')')
c
c---------------------------------------------------------------------
c     *** normal exit ***
c---------------------------------------------------------------------
c
  500 continue
      call close_maf1
      stop
      end

      subroutine initialize (kode, modblk)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      character cyorn
c
      integer*2 intp
c
      integer*4 ieof, modblk(8)
c
      common    /interpolate/intp
c
      external open_maf1  !$pragma C( open_maf1 )
c
c-----------------------------------------------------------------------
c     *** open the input data file ***
c-----------------------------------------------------------------------
c
      call open_maf1 (ieof)
c
c-----------------------------------------------------------------------
c     *** get start and stop time ***
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
  200 continue
      write (*,'('' Enter start and stop phase angles (deg): '', $)')
      read (*,*) iang1, iang2
      if (iabs(iang1) .gt. 180) then
         write (*, '('' ERROR: start angle out of range (+,- 180)'',)')
         go to 200
      end if
      if (iabs(iang2) .gt. 180) then
         write (*, '('' ERROR: start angle out of range (+,- 180)'',)')
         go to 200
      end if
      if (iang2 .lt. iang1) then
         write (*, '('' ERROR: stop angle less than start angle'')')
         go to 200
      end if
      modblk(7) = iang1
      modblk(8) = iang2
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 intp
c
      integer*4 cnt(32,2), iicx, lkuims(2,32), modblk(8)
      integer*4 ftms1(4), ieof, ifirst(2), iorboft, iorbtime, iorbyd,
     *          iwindow, istrip_times(180,2), itime(2,13), irec_yd,
     *          irec_ms
      integer*4 doy_date, immddyy, iyyddd
c
      real      cts(32,2), ims_data(2,180,32), orbs(4,32)
c
      real*8    ftms(4), stms(4)
c
      common    /comlim/stms,
     *          /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
      external open_maf1  !$pragma C( open_maf1 )
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     *** 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)
c         if (kode .eq. -10) go to 3000
         if (kode .eq. -10) then
            call open_maf1 (ieof)
            if (ieof .eq. 1) go to 3000
         end if
         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)
c         if (kode .eq. -10) go to 3000
         if (kode .eq. -10) then
            call open_maf1 (ieof)
            if (ieof .eq. 1) go to 3000
         end if
         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.dat', access = 'sequential',
     *      form = 'formatted', 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
      iyyddd = ftms(1)
      immddyy = doy_date (iyyddd)
      immdd = immddyy / 100
      imm = immdd / 100
      idd = immdd - (imm * 100)
      write (2,*) (modblk(i),i=1,8),iicx,(ftms1(i),i=1,4)
      write (2,*) imm, idd
      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 iat(8), iaz(8), ical(8), icde(512,2), 
     *          idata(2812), ifxs(8), iit(8), iiz(8), imdf(8), imsf(8), 
     *          irpac(8), irpf(8), itr(8), 
     *          iw7(16,8), izms(512), jcr(512,2)
c
      integer*4 cnt(32,2), imsh(32,2,8), irpa(32,2,8), 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
      call deffgs
      call getmd1
      call is_mass_scan (isms)
      if (isms .eq. 0) 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