program create_spin_rpa
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
c
      integer*4 header(24), ihave_data, ii, kode, maf1_err, modblk(8),
     *          nang, irpa(32), ndif
c
      real      oa_data(6), re
c
      common    /hedsav/header
      common    /blktim/iyd, it
c
      external close_maf1  !$pragma C( close_maf1 )
c
c-----------------------------------------------------------------------
c     *** initialize the data collection parameters ***
c-----------------------------------------------------------------------
c
      call initialize (modblk, nang, kode)
c
c-----------------------------------------------------------------------
c     *** search for start time ***
c-----------------------------------------------------------------------
c
      kode = -1
  100 continue
      do while (kode .eq. -1)
         call getmf1 (maf1_err)
         if (maf1_err .lt. 0) go to 9000
c*         call mshmsm (it, ih, im, is, ims)
c*         write (*,*) 'record time = ', iyd, ih, im, is, ims
         call lmstst (kode)
      end do
c
c-----------------------------------------------------------------------
c     *** if record within processing time start processing ***
c-----------------------------------------------------------------------
c
  200 continue
      if (kode .eq. 0) then
c
c-----------------------------------------------------------------------
c        *** collect the data ***
c-----------------------------------------------------------------------
c
         call cnts_eng_spn (modblk, kode, maf1_err, ihave_data, irpa,
     *                      ndif, nang)
c
c-----------------------------------------------------------------------
c        *** if have data process it and save it ***
c-----------------------------------------------------------------------
c
         if (ihave_data .eq. 1) then
c
c-----------------------------------------------------------------------
c           *** save the OA data ***
c-----------------------------------------------------------------------
c
            re = 0.0
            do ii = 1, 3
               re = re +
     *              float(header(ii+8)) * float(header(ii+8)) * 1e-6
            end do
            re = sqrt(re)
            oa_data(1) = re
            vel = 0.0
            do ii = 1, 3
               vel = vel +
     *               float(header(ii+11)) * float(header(ii+11)) * 1.e-4
            end do
            vel = sqrt(vel)
            oa_data(2) = vel
            oa_data(3) = float(header(19)) * 0.01
            oa_data(4) = float(header(17)) * 0.01
            oa_data(5) = float(header(15)) * 0.01
            oa_data(6) = float(header(16)) * 0.01
c
c-----------------------------------------------------------------------
c           *** output the data ***
c-----------------------------------------------------------------------
c
            call output_data (modblk, oa_data, irpa, ndif, nang)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
         end if     ! ihave_data
      end if     ! kode
c
c-----------------------------------------------------------------------
c     *** see if more data to process ***
c-----------------------------------------------------------------------
c
      if (maf1_err .lt. 0) go to 9000
c
c-----------------------------------------------------------------------
c     *** see if more subperiods to collect ***
c-----------------------------------------------------------------------
c
  300 continue
      call lmsnxt (kode)
      if (kode .eq. 0) then
         call lmstst (kode)
         if (kode .eq. -1) then
            go to 100
         else if (kode .eq. 0) then
            go to 200
         else
            go to 300
         end if
      else
         write (*, '('' No more subperiods to process: END OF JOB'')')
         stop
      end if
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
 9000 continue
      if (maf1_err .eq. -1) then
         write (*, '('' MAF1 read error: END OF JOB'')')
         stop
      end if
      if (maf1_err .eq. -10) then
         write (*, '('' MAF1 end of file: END OF JOB'')')
         stop
      end if
c
      call close_maf1
      close (10)
c
      end



      subroutine cnts_eng_spn (modblk, kode, maf1_err, ihave_data, irpa,
     *                         ndif, nang)
c
c-----------------------------------------------------------------------
c
c this routine does all the data accumulation for the three heads
c
c-----------------------------------------------------------------------
c
      integer*4 ictnm_rpa(360), ictnp_rpa(360), ictnr_rpa(360),
     *          ictnr_spn(360), ie, ipangs(360,2),
     *          maf1_err, imsh(512),
     *          irpa(32), is, jflag, jrpa(512), kode, ihave_data,
     *          modblk(8),
     *          nang, ndif
c
      real      ctsm_rpa(360), ctsm_rpa_sqr(360), ctsp_rpa(360), 
     *          ctsp_rpa_sqr(360), ctsr_rpa(360), ctsr_rpa_sqr(360),
     *          ctsr_spn(360), ctsr_spn_sqr(360)
c
      common    /head_data/ctsm_rpa, ctsm_rpa_sqr, ictnm_rpa,
     *                     ctsp_rpa, ctsp_rpa_sqr, ictnp_rpa,
     *                     ctsr_rpa, ctsr_rpa_sqr, ictnr_rpa,
     *                     ctsr_spn, ctsr_spn_sqr, ictnr_spn,
     *                     ipangs,
     *          /rpamsh/jrpa, imsh
      common    /blktim/iyd, it
c
c-----------------------------------------------------------------------
c     *** initialize data values ***
c-----------------------------------------------------------------------
c
      ihave_data = 0
      call zero_data
c
c-----------------------------------------------------------------------
c     *** see if record is within time limits ***
c-----------------------------------------------------------------------
c
  100 call lmstss (kode, is, ie)
      if (kode) 200,300,1000
c
c-----------------------------------------------------------------------
c     *** get next record ***
c-----------------------------------------------------------------------
c
  200 call getmf1 (maf1_err)
      if (maf1_err .eq. -10) go to 1000           ! end of file
      if (maf1_err .eq. -1) go to 1000
c*      call mshmsm (it, ihr, imn, isc, ims)
c*      write (*,*) 'record time = ', iyd, ihr, imn, isc, ims
      go to 100
c
c-----------------------------------------------------------------------
c     *** this record in current time limits so process ***
c-----------------------------------------------------------------------
c
  300 continue
      if (ihave_data .eq. 0) call orbsav
      call deffgs                        ! define flags
      call getmd1                        ! and instrument mode
c
c-----------------------------------------------------------------------
c     ***  determine the unique RPA sequence  ***
c     *** assumed present in first 32 samples ***
c-----------------------------------------------------------------------
c
      call fndval (jrpa, 32, irpa, ndif)         ! get unique values
      if (ndif .le. 0) go to 200                 ! check we have some
      call sort (irpa, ndif)                     ! and sort them
c
c-----------------------------------------------------------------------
c     *** collect the RADIAL head spin data ***
c-----------------------------------------------------------------------
c
      modblk(1) = 1
      call srtang (modblk, ctsr_spn, ctsr_spn_sqr, ictnr_spn, nang, is,
     *             ie, jflag)
      call parset (ipangs, nang)
      ihave_data = max0(jflag, ihave_data)
c
c-----------------------------------------------------------------------
c     *** collect the RADIAL head RPA data ***
c-----------------------------------------------------------------------
c
      modblk(1) = 1
      call srtrp2 (modblk, ctsr_rpa, ctsr_rpa_sqr, ictnr_rpa, is, ie,
     *             jflag, irpa, ndif)
      ihave_data = max0(jflag, ihave_data)
c
c-----------------------------------------------------------------------
c     *** collect the +Z head RPA data ***
c-----------------------------------------------------------------------
c
      modblk(1) = 2
      call srtrp2 (modblk, ctsp_rpa, ctsp_rpa_sqr, ictnp_rpa, is, ie,
     *             jflag, irpa, ndif)
      ihave_data = max0(jflag, ihave_data)
c
c-----------------------------------------------------------------------
c     *** collect the -Z head RPA data ***
c-----------------------------------------------------------------------
c
      modblk(1) = 3
      call srtrp2 (modblk, ctsm_rpa, ctsm_rpa_sqr, ictnm_rpa, is, ie,
     *             jflag, irpa, ndif)
      ihave_data = max0(jflag, ihave_data)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      go to 200                       ! get next record
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
 1000 return
      end



      subroutine zero_data
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      integer*4 ictnm_rpa(360), ictnp_rpa(360), ictnr_rpa(360),
     *          ictnr_spn(360), ii, ipangs(360,2)
c
      real      ctsm_rpa(360), ctsm_rpa_sqr(360), ctsp_rpa(360), 
     *          ctsp_rpa_sqr(360), ctsr_rpa(360), ctsr_rpa_sqr(360),
     *          ctsr_spn(360), ctsr_spn_sqr(360)
c
      common    /head_data/ctsm_rpa, ctsm_rpa_sqr, ictnm_rpa,
     *                     ctsp_rpa, ctsp_rpa_sqr, ictnp_rpa,
     *                     ctsr_rpa, ctsr_rpa_sqr, ictnr_rpa,
     *                     ctsr_spn, ctsr_spn_sqr, ictnr_spn,
     *                     ipangs
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      do ii = 1, 360
         ictnm_rpa(ii) = 0
         ctsm_rpa(ii) = 0.0
         ctsm_rpa_sqr(ii) = 0.0
         ictnp_rpa(ii) = 0
         ctsp_rpa(ii) = 0.0
         ctsp_rpa_sqr(ii) = 0.0
         ictnr_rpa(ii) = 0
         ctsr_rpa(ii) = 0.0
         ctsr_rpa_sqr(ii) = 0.0
         ictnr_spn(ii) = 0
         ctsr_spn(ii) = 0.0
         ctsr_spn_sqr(ii) = 0.0
      end do
      call parini (ipangs, nang)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end



      subroutine initialize (modblk, nang, kode)
c
c-----------------------------------------------------------------------
c     This routine initializes all the collection values.
c-----------------------------------------------------------------------
c
      byte      text(4)
c
      integer*4 iang1, iang2, irpa1, irpa2, jm1, jm2, kode, modblk(8),
     *          nang
c
      external open_maf1  !$pragma C( open_maf1 )
c
c-----------------------------------------------------------------------
c     *** open the input data file***
c-----------------------------------------------------------------------
c
      call open_maf1
c
c-----------------------------------------------------------------------
c     *** get the various time limits ***
c-----------------------------------------------------------------------
c
      call lmsset (kode)
c
c-----------------------------------------------------------------------
c     *** get which head ***
c-----------------------------------------------------------------------
c
      modblk(1) = 0
c
c-----------------------------------------------------------------------
c     *** get the mass setting range and channel ***
c-----------------------------------------------------------------------
c
      call masset (jm1, jm2)
      call mastxt (jm1, jm2, modblk(2), text, text)
      modblk(3) = jm1
      modblk(4) = jm2
      call mascnv (modblk(3), modblk(4))
c
c-----------------------------------------------------------------------
c     *** get the 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     *** get the 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-----------------------------------------------------------------------
c
  300 continue
      write (*, '('' Enter number of angle bins (0 = 72): '',$)')
      read (*, *) nang
      if (nang .lt. 0  .or.  nang .gt. 360) then
         write (*, '('' Invalid number of angle bins'')')
         go to 300
      end if
      if (nang .eq. 0) nang = 72
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end



      subroutine output_data (modblk, oa_data, irpa, ndif, nang)
c
c-----------------------------------------------------------------------
c
c     This routine calculates the RPA values, counts/sample and standard
c     deviation for the Z head RPA curves and the center angle, counts/
c     sample and standard deviation for the RADIAL head spin curve.
c
c-----------------------------------------------------------------------
c
      byte      ilohi(2), txt(4,2)
c
      character buf*80
c
      integer*4 icollect_times(4), ictnm_rpa(360), ictnp_rpa(360),
     *          ictnr_rpa(360), ictnr_spn(360), idur, ii, inc,
     *          ipangs(360,2), irpa(32),
     *          irun_times(4), modblk(8), nang, ndif
c
      real      ctsm_rpa(360), ctsm_rpa_sqr(360), ctsp_rpa(360), 
     *          ctsp_rpa_sqr(360), ctsr_rpa(360), ctsr_rpa_sqr(360),
     *          ctsr_spn(360), ctsr_spn_sqr(360)
      real      oa_data(6), del_ang, del_ang2, out(12), vrpa(360)
c
      common    /com_times/irun_times, icollect_times, idur, inc,
     *          /head_data/ctsm_rpa, ctsm_rpa_sqr, ictnm_rpa,
     *                     ctsp_rpa, ctsp_rpa_sqr, ictnp_rpa,
     *                     ctsr_rpa, ctsr_rpa_sqr, ictnr_rpa,
     *                     ctsr_spn, ctsr_spn_sqr, ictnr_spn,
     *                     ipangs
c
      data      ilohi/'L', 'H'/
c
c-----------------------------------------------------------------------
c     *** output first line of text ***
c-----------------------------------------------------------------------
c
      iyr = (icollect_times(1) / 1000)
      idn = mod(icollect_times(1), 1000)
      call mshmsm (icollect_times(2), ihr1, imn1, isc1, ms1)
      call mshmsm (icollect_times(4), ihr2, imn2, isc2, ms2)
      write (buf, 100) iyr, idn, ihr1, imn1, isc1, ihr2, imn2, isc2
  100 format ('DE RIMS  ', i2, '/', i3.3, ' ', i2.2, 2(':', i2.2),
     *         ' - ', i2.2, 2(':',i2.2))
      write (10, *) (buf(i:i), i = 1, 35)
c
c-----------------------------------------------------------------------
c     *** output second line of text ***
c-----------------------------------------------------------------------
c
      jm1 = modblk(3) * (((modblk(2) - 1) * 5) - 4)
      jm2 = modblk(4) * (((modblk(2) - 1) * 5) - 4)
      call mastxt (jm1, jm2, idum, txt, txt(1,2))
      write (buf, 200) ilohi(modblk(2)), (txt(i, modblk(2)), i = 1,4),
     *                 modblk(3), modblk(4)
  200 format (a1, '/', 4a1, '   mass range = ', i4, ' -> ', i4)
      write (10, *) (buf(i:i), i = 1, 34)
c
      ichan = modblk(2)
      if (ichan .eq. 1) then
         if (txt(2,ichan)  .eq. '+') then
            ion_num = 1   ! L/H+
         else if (txt(4,ichan) .eq. '+') then
            ion_num = 5   ! L/He++
         else
            ion_num = 3   ! L/He+
         end if
      else
         if (txt(1,ichan) .eq. 'H') then
            ion_num = 2   ! H/He+
         else if (txt(1,ichan) .eq. 'N') then
            ion_num = 7   ! H/N+
         else if (txt(1,ichan) .eq. 'M') then
            ion_num = 8   ! H/Mo+
         else if (txt(3,ichan) .eq. '+') then
            ion_num = 6   ! H/O++
         else
            ion_num = 4   ! H/O+
         end if
      end if
c
c-----------------------------------------------------------------------
c     *** output third and forth line of text ***
c-----------------------------------------------------------------------
c
      rpa1 = modblk(5) * 0.05
      rpa2 = modblk(6) * 0.05
      write (buf, 300) modblk(5), rpa1, modblk(6), rpa2
  300 format ('rpa range = ', i4, ' (', f5.2, 'ev)', ' -> ', i4, ' (',
     *        f5.2, 'ev)')
      write (10, *) (buf(i:i), i = 1, 44)
      write (buf, 400) modblk(7), modblk(8)
c
  400 format ('angle range = ', i4, ' -> ', i4, ' deg')
      write (10, *) (buf(i:i), i = 1, 30)
c
c-----------------------------------------------------------------------
c     *** output fifth line of text ***
c-----------------------------------------------------------------------
c
      write (buf, 500) oa_data
  500 format ('re = ', f4.2, '  vel = ', f4.2, '  mlt = ', f5.2,
     *        '  mlat = ', f6.2, '  L = ', f6.2, '  invlat = ', f6.2)
      write (10, *) (buf(i:i), i = 1, 78)
c
c-----------------------------------------------------------------------
c     *** output number of angle bins and number of rpa values ***
c-----------------------------------------------------------------------
c
      write (10, '(1x, i3, 1x, i2, 1x, i2)') nang, ndif, ion_num
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      do ii = 1, 360
         vrpa(ii) = -1.0
      end do
      do ii = 1, ndif
         vrpa(ii) = float(irpa(ii)) * 0.05
      end do
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      del_ang = 360.0 / float(nang)
      del_ang2 = del_ang / 2.0
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      do ii = 1, nang
c
         out(1) = del_ang2 + (del_ang * float(ii -1)) - 180.0
c
         if (ictnr_spn(ii) .gt. 0) then
            out(2) = ctsr_spn(ii) / float(ictnr_spn(ii))
            avg_sqr = ctsr_spn_sqr(ii) / float(ictnr_spn(ii))
            dif = avg_sqr - out(2)**2
            if (ictnr_spn(ii) .gt. 1) then
               factor = ictnr_spn(ii) / (ictnr_spn(ii) - 1)
            else
               factor = 1.0
            end if
            out(3) = sqrt(factor * dif)
         else
            out(2) = -1.0
            out(3) = -1.0
         end if
c
         out(4) = vrpa(ii)
c
         if (ictnr_rpa(ii) .gt. 0) then
            out(5) = ctsr_rpa(ii) / float(ictnr_rpa(ii))
            avg_sqr = ctsr_rpa_sqr(ii) / float(ictnr_rpa(ii))
            dif = avg_sqr - out(5)**2
            if (ictnr_rpa(ii) .gt. 1) then
               factor = ictnr_rpa(ii) / (ictnr_rpa(ii) - 1)
            else
               factor = 1.0
            end if
            out(6) = sqrt(factor * dif)
         else
            out(5) = -1.0
            out(6) = -1.0
         end if
c
         if (ictnm_rpa(ii) .gt. 0) then
            out(7) = ctsm_rpa(ii) / float(ictnm_rpa(ii))
            avg_sqr = ctsm_rpa_sqr(ii) / float(ictnm_rpa(ii))
            dif = avg_sqr - out(7)**2
            if (ictnm_rpa(ii) .gt. 1) then
               factor = ictnm_rpa(ii) / (ictnm_rpa(ii) - 1)
            else
               factor = 1.0
            end if
            out(8) = sqrt(factor * dif)
         else
            out(7) = -1.0
            out(8) = -1.0
         end if
c
         if (ictnp_rpa(ii) .gt. 0) then
            out(9) = ctsp_rpa(ii) / float(ictnp_rpa(ii))
            avg_sqr = ctsp_rpa_sqr(ii) / float(ictnp_rpa(ii))
            dif = avg_sqr - out(9)**2
            if (ictnp_rpa(ii) .gt. 1) then
               factor = ictnp_rpa(ii) / (ictnp_rpa(ii) - 1)
            else
               factor = 1.0
            end if
            out(10) = sqrt(factor * dif)
         else
            out(9) = -1.0
            out(10) = -1.0
         end if
c
         out(11) = ipangs(ii,1)
         out(12) = ipangs(ii,2)
c
         write (10,
     *          '(1x, f7.2, f10.2, f8.2, f6.2, 3(f10.2, f8.2)), 2i5')
     *          out
c
      end do
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      return
      end