program spntgs
c
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c
      byte      icx
c
      integer*4 kode, modblk(8)
c
      external close_maf1  !$pragma C( close_maf1 )
c
c-----------------------------------------------------------------------
c     *** initialization data colection parameters ***
c-----------------------------------------------------------------------
c
      call initialize (kode, modblk, icx)
      if (kode .eq. -1) go to 500
c
c-----------------------------------------------------------------------
c     *** process the data ***
c-----------------------------------------------------------------------
c
      call proces (modblk, iret, kflag, icx)
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, icx)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      byte      icx, text(4)
c
      integer*4 ieof, modblk(8)
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 times ***
c-----------------------------------------------------------------------
c
      call lmssetr (kode)
      if (kode .eq. -1) go to 99
c
c-----------------------------------------------------------------------
c     *** get duration of each frame ***
c-----------------------------------------------------------------------
c
      call frasetr (kode)
c
c-----------------------------------------------------------------------
c     *** define which head ***
c-----------------------------------------------------------------------
c
      call hedset (modblk(1))
      if (modblk(1) .eq. 4) go to 10
c
c-----------------------------------------------------------------------
c     *** define mass setting range ***
c-----------------------------------------------------------------------
c
      call masset (jm1, jm2)
c
c-----------------------------------------------------------------------
c     * single mass det define ilohi lo/hi flag ***
c-----------------------------------------------------------------------
c
      call mastxt (jm1, jm2, modblk(2), text, text)
      modblk(3)=jm1
      modblk(4)=jm2
      call mascnv (modblk(3), modblk(4))
   10 continue
c
c-----------------------------------------------------------------------
c     *** define rpa range ***
c-----------------------------------------------------------------------
c
  100 continue
      write (*,'('' Enter minimum and maximum RPA settings: '', $)')
      read (*,*) irpa1, irpa2
      if (irpa1 .lt. 0  .or.  irpa1 .gt. 1000) then
         write (*, '('' ERROR: minimum rpa out of range (0->1000)'')')
         go to 100
      end if
      if (irpa2 .lt. 0  .or.  irpa2 .gt. 1000) then
         write (*, '('' ERROR: maximum rpa out of range (0->1000)'')')
         go to 100
      end if
      if (irpa2 .lt. irpa1) then
         write (*, '('' ERROR: maximum rpa less than minimum'')')
         go to 100
      end if
      modblk(5) = irpa1
      modblk(6) = irpa2
c
c-----------------------------------------------------------------------
c     *** define aperture bias ***
c-----------------------------------------------------------------------
c
      call apset (icx)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
   99 return
      end

      subroutine proces (modblk,kode,lflag,icx)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      logical   ifirst1
c
      byte      icx
c
      integer*4 ctn(32), iicx, modblk(8)
c
      integer*4 ftms1(4),ieof,ifirst(2),iorboft,iorbtime,iorbyd,
     *          iwindow,istrip_times(180,2),itime(2,13),
     *          pitch_angle(180,2,3)
      integer*4 doy_date, immddyy, iyyddd
c
      real     cts(32),orbs(4,13),spn_data(180,32)
c
      real*8   stms(4),ftms(4)
c
      common    /comlim/stms,
     *          /fratms/ftms,ndiv,
     *          /orbdat/itime,orbs,
     *          /orbit/iorboft,iwindow,iorbyd,iorbtime,inum,ifirst
c
      data      ifirst1/.true./
c
      external open_maf1  !$pragma C( open_maf1 )
c
c-----------------------------------------------------------------------
c     *** initialize ***
c-----------------------------------------------------------------------
c
      lflag=0
      mflag=0                  ! flag for pitch angle ident
      do i=1,180
         do j=1,32
            spn_data(i,j) = -1.0
         end do
         istrip_times(i,1) = 0
         istrip_times(i,2) = 0
      end do
c
c-----------------------------------------------------------------------
c     *** initialise the orbit parameters and  ***
c     ***    initialize save information       ***
c-----------------------------------------------------------------------
c
      call orboftr
      do i=1,13
         itime(1,i)=-1
         itime(2,i)=-1
         do j=1,4
            orbs(j,i)=0.0
         end do
      end do
c
c-----------------------------------------------------------------------
c     *** if first time in routine read first record ***
c-----------------------------------------------------------------------
c
      if (ifirst1) then
         call getmf1 (kode)
c         if (kode .eq. -10) go to 3000            ! end of file
         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     *** loop through the number of divisions for this time span ***
c-----------------------------------------------------------------------
c
      indiv=0
      do 2000 idiv=1,ndiv
         indiv=indiv+1
c
c-----------------------------------------------------------------------
c        *** initialize working arrays ***
c-----------------------------------------------------------------------
c
         minang=-1000
         nang=32
         jflag=0
         do i=1,32
            cts(i)=0.0
            ctn(i)=0
         end do
c
c-----------------------------------------------------------------------
c        *** record within time limits ? ***
c-----------------------------------------------------------------------
c
  200    call lmstssr (kode,i1,i2)
         if (kode) 400,600,1800     ! before,during,after
c
c-----------------------------------------------------------------------
c        *** get next record ***
c-----------------------------------------------------------------------
c
  400    call getmf1 (kode)
c         if (kode .eq. -10) go to 3000            ! end of file
         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 in current 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        *** test if required aperture bias ***
c-----------------------------------------------------------------------
c
         call aptst (icx,kode)
         if (kode .eq. -1) go to 800     ! no
c
c-----------------------------------------------------------------------
c        *** proper aperture bias - continue processing ***
c-----------------------------------------------------------------------
c
         call deffgs                        ! define flags
         call getmd1                        ! and instrument mode
c
c-----------------------------------------------------------------------
c        *** collect the data ***
c-----------------------------------------------------------------------
c
         if (minang .eq. -1000) call pitpha (minang)
         call srtang (modblk,cts,ctn,nang,i1,i2,kflag)
         jflag=max0(jflag,kflag)
         lflag=max0(jflag,lflag)
  800    continue
c
c-----------------------------------------------------------------------
c        *** see if need more data for this division ***
c-----------------------------------------------------------------------
c
         if (i2 .eq. 512) go to 400
c
c-----------------------------------------------------------------------
c        *** average data ***
c-----------------------------------------------------------------------
c
         do i=1,32
            if (ctn(i) .le. 0) then
               cts(i)=-1.0
            else
               cts(i)=cts(i)/float(ctn(i))
            end if
         end do
c
c-----------------------------------------------------------------------
c        *** put data in storage arrays ***
c-----------------------------------------------------------------------
c
 1000    continue
         iflag=1
         istrip_times(idiv,1)=stms(2)/1000.0d0
         istrip_times(idiv,2)=stms(2)/1000.0d0
         do ispn=1,nang
            spn_data(idiv,ispn)=cts(ispn)
         end do
c
c-----------------------------------------------------------------------
c        *** store the pitch angle (radial head only) ***
c-----------------------------------------------------------------------
c
         if (modblk(1) .ne. 1) go to 1200
         if (minang .eq. -1000) go to 1200
         pitch_angle(idiv,1,1)=minang
         pitch_angle(idiv,1,2)=stms(2)/1000.0d0
         pitch_angle(idiv,1,3)=stms(4)/1000.0d0
         mflag=1
         maxang=minang+180
         if (maxang .gt. 180) maxang=maxang-360
         pitch_angle(idiv,2,1)=maxang
         pitch_angle(idiv,2,2)=stms(2)/1000.0d0
         pitch_angle(idiv,2,3)=stms(4)/1000.0d0
 1200    continue
c
c-----------------------------------------------------------------------
c        *** get time of next division ***
c-----------------------------------------------------------------------
c
 1800    call lmsnxtr (kode)
         if (kode .ne. 0) go to 3000
 2000 continue
c
c-----------------------------------------------------------------------
c     *** process the last division ***
c-----------------------------------------------------------------------
c
 3000 continue
      do i=1,32
         if (ctn(i) .le. 0) then
            cts(i)=-1.0
         else
            cts(i)=cts(i)/float(ctn(i))
         end if
      end do
      istrip_times(idiv,1)=stms(2)/1000.0d0
      istrip_times(idiv,2)=stms(2)/1000.0d0
      do ispn=1,nang
         spn_data(idiv,ispn)=cts(ispn)
      end do
      if (modblk(1) .ne. 1) go to 3100
      if (minang .eq. -1000) go to 3100
      pitch_angle(idiv,1,1)=minang
      pitch_angle(idiv,1,2)=stms(2)/1000.0d0
      pitch_angle(idiv,1,3)=stms(4)/1000.0d0
      mflag=1
      maxang=minang+180
      if (maxang .gt. 180) maxang=maxang-360
      pitch_angle(idiv,2,1)=maxang
      pitch_angle(idiv,2,2)=stms(2)/1000.0d0
      pitch_angle(idiv,2,3)=stms(4)/1000.0d0
 3100 continue
c
c-----------------------------------------------------------------------
c     *** output data array to output file ***
c-----------------------------------------------------------------------
c
c      open (unit = 2, file = 'spntgs.dat', access = 'sequential',
c     *      form = 'formatted', status = 'new', err = 3200)
      open (unit = 2, file = 'spntgs.dat', access = 'sequential',
     *      form = 'formatted', err = 3200)
      go to 3300
 3200 continue
      type '('' error opening output file'')'
      return
 3300 continue
      ftms1(1)=ftms(1)
      ftms1(2)=ftms(2)/1000.0d0
      ftms1(3)=ftms(3)
      ftms1(4)=ftms(4)/1000.0d0
      if (icx .eq. '0') then
         iicx = 0
      else if (icx .eq. '2') then
         iicx = 2
      else if (icx .eq. '4') then
         iicx = 4
      else if (icx .eq. '8') then
         iicx = 8
      else
         iicx = 10
      end if
      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,*) iorboft,ifirst(1),ifirst(2)
      write (2,*) inum
      do i=1,inum
         write (2,*) (itime(j,i),j=1,2),(orbs(j,i),j=1,4)
      end do
      write (2,*) ndiv
      write (2,*) nang
      do i=1,ndiv
c*         write (2,*) istrip_times(i,1),istrip_times(i,2),
c*     *               (spn_data(i,j),j=1,nang)
         write (2,*) (spn_data(i,j),j=1,nang)
      end do
      if (modblk(1) .eq. 1) then
         do i=1,ndiv
            write (2,*) ((pitch_angle(i,j,k),j=1,2),k=1,3)
         end do
      end if
      close (2)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end

      subroutine srtang(modblk,cts,ictn,nang,i1,i2,jflag)
c
c----------------------------------------------------------------------|
c                                                                      |
c     this subroutine sorts data for a given head and given mass in    |
c     terms of instrument spin phase angle.                            |
c                                                                      |
c----------------------------------------------------------------------|
c
      integer*2 icde(512,2),idat(2812),jcr(512,2),kcel(512)
c      integer*2 icde(512,2),ictn(nang),idat(2812),jcr(512,2),jmsh(512)
c     *          jrpa(512),kcel(512),modblk(8)
c
      integer*4 modblk(8), jrpa(512), jmsh(512), ictn(nang)
      integer*4 idcodc, ic
c
      real      cts(nang)
c
      common    /maf1/idat,
     *          /rpamsh/jrpa,jmsh
c
      equivalence (icde(1,1),idat(1789)),(jcr(1,1),idat(253)),
     *            (kcel(1),idat(1277))
c
c-----------------------------------------------------------------------
c     *** initialize ***
c-----------------------------------------------------------------------
c
      jflag=0
c
c-----------------------------------------------------------------------
c     *** first call refang, to define phase angle reference ***
c-----------------------------------------------------------------------
c
      call refang(degsam,ramang)
c
c-----------------------------------------------------------------------
c     *** see which head ***
c-----------------------------------------------------------------------
c
      if(modblk(1).eq.4)go to 4000     ! electrometer
      if(modblk(1).eq.1)go to 401     ! radial
c
c-----------------------------------------------------------------------
c     *** z head processing ***
c-----------------------------------------------------------------------
c
      idet=(modblk(1)-2)*2+modblk(2)            ! get detector
c
      do 430 i=i1,i2
         if (jrpa(i) .ge. modblk(5)) then     ! check energy
            if (jrpa(i) .le. modblk(6)) then
               if (jmsh(i) .ge. modblk(3)) then     ! check mass
                  if (jmsh(i) .le. modblk(4)) then
                     izms=igtzms(idet,i)     ! have data ?
                     if (izms .le. 0) go to 430
                     jflag=1
                     ang=ramang+(i-1)*degsam     ! get angle
                     ang=ang+900.0
                     ang=amod(ang,360.0)
                     iang=1+ang*nang/360.0
                     temp = float(idcodc(icde(i,izms),ic))
                     cts(iang)=cts(iang)+float(idcodc(icde(i,izms),ic))
                     ictn(iang)=ictn(iang)+ic
                  end if
               end if
            end if
         end if
  430 continue
      go to 9000
c
c-----------------------------------------------------------------------
c     *** radial head processing ***
c-----------------------------------------------------------------------
c
  401 continue
      do 100 i=i1,i2
         if (jrpa(i) .ge. modblk(5)) then     ! check energy
            if (jrpa(i) .le. modblk(6)) then
               if (jmsh(i) .ge. modblk(3)) then     ! check mass
                  if (jmsh(i) .le. modblk(4)) then
                     jflag=1
                     ang=ramang+(i-1)*degsam     ! get angle
                     ang=ang+900.0
                     ang=amod(ang,360.0)
                     iang=1+ang*nang/360.0
                     cts(iang)=cts(iang)+
     *                                float(idcodc(jcr(i,modblk(2)),ic))
                     ictn(iang)=ictn(iang)+ic
                  end if
               end if
            end if
         end if
  100 continue
       go to 9000
c
c-----------------------------------------------------------------------
c    *** electrometer channel ***
c-----------------------------------------------------------------------
c
 4000 continue
      do 5000 i=i1,i2
         if (jrpa(i) .ge. modblk(5)) then     ! check energy
            if (jrpa(i) .le. modblk(6)) then
               jflag=1
               ang=ramang+(i-1)*degsam     ! get angle
               ang=ang+900.0
               ang=amod(ang,360.0)
               iang=1+ang*nang/360.0
               cts(iang)=cts(iang)+float(idcodc(kcel(i),ic))
               ictn(iang)=ictn(iang)+ic
            end if
         end if
 5000 continue
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
 9000 continue
      return
      end