ACOMT.FTN
FTN77 
$FILES 0,1
      PROGRAM ACOMT(3,65), COMET SPECTRAL LINE PROGRAM <861128.0825>
********************************************************************
*     256 channel spectral line observing program.
*     Programs ACFIN, ACFPM, ACEPH, ALARM must be compiled and loaded.
*     Subroutines other than fileread are in &ACSUB 
*     Modified for binary write to tape file
* 
C     Spectral line source file parameters for ACOMT: 
C 
C     line 1 (free format): 
C     --------------------- 
C     name     character*20  max - source name, in quotes 
C     filename character*20  max - name of comet ephemeris parameter file 
C     nacc     integer*2     number of correlator accumulations 
C                            in each position 
C     nobs     integer*2     number of observations of this position pair 
C     kount    integer*2     octal correlator thumbwheel count, 
C                            written as  @37777 for example 
C     bw       real*4        correlator bandwidth in MHz eg 0.32
C     synofs   real*4        synthesizer offset in Hz applied to center 
C                            the line in passband.
C                            If zero, 0 km/s Vlsr is center of spectrum.
C     synrf    real*4        synthesizer offset in Hz applied to reference
C                            position ONLY, in addition to synofs, for
C                            frequency switching, and must NOT be zero. 
C     rf       real*4        line rest frequency in MHz 
C 
* 
      integer   idata(570), itime(5),  ipram(5),  buffer(10), 
     *          acfin(3),   aceph(3),  acfpm(3),  alarm(3), 
     *          isynth(4),  idvmsl(4), ndsel(4),
     *          inbuf(10),  ioutbuf(20) 
* 
      character filename*10, name*20, ktype*2, cbuffer*20,
     *          acephfile*20
* 
      logical   aok, fault, ifbrk, error, file_end, overflow, 
     *          lochange, newsource 
* 
      real      data(285), diff(256), ref(256), avrg(256),
     *          wave(4), rfmax(4), habm(4), decbm(4), xmc(4), dnoise(4) 
* 
      real*8    relcor, vlight
* 
*     data / idata arrays hold the spectrum and housekeeping info.
*     ioutbuf passes the comet element file name to ACEPH 
*     inbuf   retrieves the comet ra, dec and velocity from ACEPH 
* 
      equivalence          (data(1),idata(1)),   (name,idata(517)), 
     * (iy,idata(527)),    (id,idata(528)),      (ih,idata(529)), 
     * (im,idata(530)),    (is,idata(531)), 
     * (iw,idata(532)),    (indx1,idata(533)),   (ipoly,idata(534)),
     * (rf,idata(535)),    (fof,idata(537)),
     * (syon,idata(539)),  (syof,idata(541)),    (bw,idata(543)), 
     * (tint,idata(545)),  (tsys,idata(547)),    (ra50st,idata(549)), 
     * (dec50st,idata(551)),(rast,idata(553)),  (decst,idata(555)), 
     * (rxmc,idata(557)),  (yscale,idata(559)),  (yofset,idata(561)), 
     * (adav,idata(563)),  (n1,idata(565)),      (n2,idata(566)), 
     * (buffer,cbuffer),   (ioutbuf,acephfile), 
     * (inbuf(1),ra1950),  (inbuf(3),dec1950),
     * (inbuf(5),ranow),   (inbuf(7),decnow),    (inbuf(9),velocity)
* 
*     initialise parameters 
*     --------------------- 
* 
      data  n/256/, lu/7/, mtlu/8/, lu8660/34/, lu3336/33/, ludvm/35/,
     *      acfin/ 'ACFIN ' /, aceph/ 'ACEPH ' /, 
     *      acfpm/ 'ACFPM ' /, alarm/ 'ALARM ' /, 
     *      dtr  /0.017453293/, vlight /299792.5d0/, nchan /4/, 
* 
*      ifreq :       1        2        3        4 
     * rfmax   / 10000.,   6000.,   4000.,   2000./,
     * wave    /    3.6,     6.0,    13.0,    18.0/,
     * habm    / +0.053,  +0.150,  +0.070,  +0.075/,
     * decbm   / -0.212,  -0.566,  +0.450,  +1.750/,
     * xmc     /   384.,     64.,     96.,     72./,
     * ndsel   /      1,       2,       4,       5/,
     * dnoise  /   10.0,    32.0,    3.96,    11.2/,
     * idvmsl  /    77b,     77b,     77b,     77b/,
     * isynth  /   3330,       0,    3330,    3330/ 
* 
*      nextop   0 at start (no plot of previous spectrum),
*               1 to continue on present source 
*               2 to observe next source in file (bit 14 set) 
*                 or select new source file / non-sequential source (bit 11)
*               3 to end. 
*      iplot    2 for individual, 3 for average spectrum plot 
*      ndif     number of spectra plotted on current page so far. 
*      adav     number of current source spectra averaged so far. 
*      intvl    plotter vertical offset between spectra 
*      intzro   plotter zero level for first plot 
*      ipage    plotter page size 
*      ipoly    order of polynomial fitted to spectrum baseline.
*      iw       0 default to rectangular smoothing window only. 
* 
      call rmpar(ipram) 
      if (ipram(1) .eq. 1) lu = ipram(1)
   10 intvl = 1100
      intzro = 550
      ipage = 1 
      ipoly = 1 
      index = 0 
      iw = 0
      n1 = 1
      n2 = n
      fault = .false. 
      nfilerrors = 0
      newsource = .false. 
      call exec(11,itime,iy)
      if (lu .eq. 7) write (lu,'(a)') char(14b) 
      write (lu,1000) iy,itime(5),itime(4),itime(3) 
 1000 format(//'SPECTRA',i6,i4,2i3/23'='/)
   20 write (lu,'(''Program control options are :''/
     * ''Set S register bit 15 on to end''/ 
     * 19x,''14 to start the next source''/ 
     * 19x,''13 or type *BR,ACOMT to suspend.''/
     * 19x,''11 to start a new source file or another source''/ 
     * ''This program shuts down when the observations are complete''/
     *                     24x,''or the source sets.'')') 
* 
*     set up magtape at the end of previous data
*     ------------------------------------------
      call mtpos (lu,mtlu,mtn,mtfile) 
      write (lu,1020) mtfile
 1020 format (6x,'Starting file ',i4) 
* 
*     read parameters for disc source file restart here if bit 11 set 
*     --------------------------------------------------------------- 
   30 filename = ' '
      write (lu,'(6x,''Name of the source file ? _'')')
      read  (lu,'(a10)',err=30) filename
   40 write (lu,'(6x,a10,'' : start at source number ? _'')') filename 
      read  (lu,*,err=40) nsource 
          if (nsource .lt. 1)  go to 30 
      if (index .eq. 0) then
   50     write (lu,'(6x,''Number of the first observation ? _'')')
          read  (lu,*,err=50) ifirst
          if (ifirst .lt. 1 .or. ifirst .gt. 9999) then 
              go to 50
          else
              index = ifirst
          endif 
      end if
* 
*     read from source file : restart here if bit 14 set or next source due 
*     --------------------------------------------------------------------- 
  100 nextop = 0
      iplot  = 2
      ndif   = 0
      iadav  = 0
      adav   = 0.0
* 
      call filerd(lu,filename,nsource,error,file_end, 
     *    name,acephfile,nacc,nobs,kount,bw,synofs,synrf,rf)
      if (error) then 
          nfilerrors = nfilerrors + 1 
          if (nfilerrors .lt. 10) then
              write (lu,'(''trying next source... '')') 
              nsource = nsource + 1 
              go to 100 
          else
              fault = .true.
              go to 990 
          end if
      endif 
      if (file_end) then
          write (lu,'(''end of source file'')') 
          go to 990 
      endif 
* 
*     error check on count and bandwidth
      if (kount .gt. 37777b) kount = 37777b 
      count = kount + 1 
      if (bw .gt. 5.12) bw = 5.12 
*     reset plot position count if more than four obs 
      if (nobs .gt. 4) ndif = 0 
* 
*     find frequency index
      do j = 1,nchan
          if (rf .lt. rfmax(j)) ifreq = j 
      end do
      write (lu,1500) nsource, name, bw, kount, rf, synofs, synrf,
     *                dnoise(ifreq), wave(ifreq)
 1500 format(//'SOURCE',i3,' : ',a20/32'='/ 
     *         'Set the correlator bandwidth to',f5.2,' MHz ' 
     *         'and the thumbwheels to ',@6/
     *         'Line rest frequency =',f13.6,' MHz, '/
     *         'synthesizer offset applied =',f8.0,' Hz.'/
     *         'synthesizer offset applied to reference =',f8.0,' Hz.'/ 
     *         'noise calibration =',f7.3,' K at',f5.1,' cm.'/) 
* 
*     call ACEPH for current coordinates of comet - observability check 
      if (index .eq. ifirst) then 
          iprint = 1
      else
          iprint = 0
      end if
      call exec (9,aceph,lu,iprint,0,0,0,ioutbuf,20)
      call exec (14,1,inbuf,10) 
      onra  = ranow + habm(ifreq) / cos (decnow * dtr)
      ondec = decnow + decbm(ifreq) 
* 
*     find command hour angle, limit to -180 < comha < +180 
      comha = time(2)/240.0 - onra
      if (comha .gt. +180.0) comha = comha - 360.0
      if (comha .lt. -180.0) comha = comha + 360.0
      write (lu,1800) comha, ondec
1800  format(/6x,'Command angles are HA',f8.3,' DEC',f8.3)
* 
*     on-source integration time
      tint = 52.5 * 5.12/bw * count/16384.0 * nacc
*     observation time in seconds, approximately
      obstime = tint * 2 + 180
* 
*     check western horizon limit 
      if (comha .gt. halim(+90.0,ondec) - obstime/240.0) then 
          write (lu,'(''the source has set ! '')') 
          go to 990 
      end if
* 
*     check eastern horizon limit 
      if (comha .lt. halim(-90.0,ondec)) then 
*         source below east horizon , wait until it has risen 
          iswait =  (comha - halim(-90.0,ondec)) * 240
          write (lu,'(''Source has not yet risen ! ''/ 
     *                ''wait '',i5,'' minutes.'')') iswait / 60
          call exec (12,0,2,0,iswait) 
      endif 
* 
      write (lu,'(6x,''On-source integration time = '', 
     *            f12.3,'' seconds'')')  tint 
      write (lu,1850) nobs, obstime / 60
 1850 format (i8,' observations each of',f5.1,' minutes') 
* 
*     start of off-source observing loop
*     ----------------------------------
*     clear the switch register 
  200 call issr(0)
* 
*     store the start time of the previous observation
      is = itime(2) 
      im = itime(3) 
      ih = itime(4) 
      id = itime(5) 
*     store the coordinates of the previous observation 
      ra50st  = ra50
      dec50st = dec50 
      rast    = aonra 
      decst   = aondec
* 
*     obtain the start time of the current observation
      call exec(11,itime,iy)
      write (lu,2100) index,iy,(itime(j), j = 5,2,-1) 
 2100 format('observation ',i4,' starts at ',2i4,3i3) 
* 
*     run ACEPH for the comet position and velocity 
      call exec(9,aceph,lu,0,0,0,0,ioutbuf,20)
      call exec (14,1,inbuf,10) 
*     aonra, aondec are actual current coordinates, 
*     onra, ondec allow for the beam offsets
      ra50   = ra1950 
      dec50  = dec1950
      aonra  = ranow
      aondec = decnow 
      onra   = aonra + habm(ifreq) / cos (aondec * dtr) 
      ondec  = aondec + decbm(ifreq)
* 
*     relativistic doppler shift correction 
*     restart calculation here if lo multiplication factor has changed
  220 relcor = velocity / vlight
      freq = rf * dsqrt ((1d0 - relcor) / (1d0 + relcor)) 
      rxmc = xmc(ifreq) 
*     synthesizer frequency to center the line in the passband  (MHz) 
      fsyn = (freq - 60.0 + 0.5 * bw) / rxmc
*     synthesizer frequency offset for on & ref (MHz) 
      fof = synofs * 1.e-6 * rxmc 
* 
*     store previous synthesizer settings, find new ones
*     synof, synon are new synthesizer settings in MHz
      syof  = synof 
      synof = fsyn + (synofs + synrf) * 1.0e-6
      syon  = synon 
      synon = fsyn + synofs * 1.0e-6
      write (lu,2400) synof 
 2400 format (6x,'L.O. Synthesizer setting:',f11.6,' MHz ') 
* 
*     select dvm input signals
      input = idvmsl(ifreq) 
      call exec(2,113b,input,1,5,5) 
* 
*     set up fluke 8840 dvm: remote, reset, medium sample rate
      call rmote (ludvm)
      write (ludvm,'(''*S1'')') 
* 
*     tests before startup :
      if (index .eq. ifirst) then 
          write (lu,'(6x,''System check :'')') 
*         check that lo is computer controllable
          call chklo (lu,lu8660,lu3336,ludvm,isynth,
     *                xmc,ifreq,synof,ndsel,dnoise,lochange)
          if (lochange) then
*             lo multiplier has changed - recalculate synth frequency 
              go to 220 
          endif 
      endif 
* 
*     start of off-source track 
*     ------------------------- 
      call setlo (isynth,ifreq,synof,lu8660,lu3336) 
*     track the current (precessed) coordinates 
      call track (lu,4,onra,ondec)
*     start the correlator program
      call exec (10,acfin,n,nacc,kount,lu,0)
*     if continuing on the present source,
*     then store the previous spectrum on tape, and plot it 
      if (nextop .eq. 1) go to 500
* 
*     collect reference spectrum and store it.
  250 call exec(23,acfin) 
      call exec(14,1,idata,512) 
      do i = 1,n
          ref(i) = data(i)
      end do
* 
*     start of on source run
*     ----------------------
*     print synthesizer setting, track on-source coordinates. 
      write (lu,2400) synon 
      call setlo (isynth,ifreq,synon,lu8660,lu3336) 
*     restart track if integration time > 30 minutes
      if (tint .gt. 1800.0) call track(lu,4,onra,ondec) 
*     start correlator program, get spectrum on completion
      call exec(10,acfin,n,nacc,kount,lu,0) 
      call exec(23,acfin) 
      call exec(14,1,idata,512) 
      do i = 1,n
          diff(i) = data(i) 
      end do
* 
*     measure the system temperature
      call exec (2,113b,input,1,5,5)
      call exec (12,0,2,0,-1) 
      call calt (lu,ludvm,ifreq,ndsel,dnoise,isynth,synon,lu8660,lu3336,
     *           10,tson,fault) 
      if (fault) go to 960
*     obtain difference spectrum and normalise to kelvins 
      do i = 1,n
          diff(i) = (diff(i)/ref(i) - 1.0) * tson 
      end do
* 
*     check termination options at end of off and on source runs
*     ----------------------------------------------------------
*     increment run number, assume continuing with next observation 
      index = index + 1 
      if (index .eq. 10000) then
          index = 1 
          ifirst = 0
      endif 
      nextop = 1
* 
*     read antenna ha and dec via dio 
  360 call exec(1,11,buffer,6,1,3)
*     mask out tach sign and magnitude bits for ha, then dec
      buffer(1) = buffer(1) .and. 171777b 
      buffer(4) = buffer(4) .and. 171777b 
      read (cbuffer,'(2f6.3)',err=360) antha,antdc
      if (antdc .gt. 180.0) antdc = antdc - 360.0 
      write (lu,3600) antha,antdc 
 3600 format (6x,'Antenna HA ',f7.3,'  DEC ',f7.3)
* 
*     check bit 14 and if number of obs on current source done
      if ( issw(14) .lt. 0 .or. nobs .eq. (iadav + 1) ) nextop = 2
*     check bit 11 for new file / source option 
      if (issw(11) .lt. 0) then 
          nextop = 2
          newsource = .true.
      end if
* 
*     omit west horizon check if unnecessary
      if (antha .gt.  0.0 .and. antha .lt. 180.0) then
*         check if horizon reached on next run
          scnang = tint / 240.0 + 0.3 
          if ( (antha+2*scnang) .gt. halim(+90.0,antdc) ) then
              nextop = 3
              write (lu,'(''Next run on horizon'')')
          endif 
      endif 
* 
*     check bit 15 for end, over-rides previous options 
      if ( issw(15) .lt. 0) nextop = 3
* 
*     check bit 13 or break for program suspension
      if (issw(13) .lt. 0 .or. ifbrk(0)) then 
*         clear the switch register 
          call issr(0)
          pause 'ACOMT suspending - type GO,ACOMT to restart' 
* 
*         restart - form feed and chart advance, check cal & lo control 
          if (lu .eq. 7) write (lu,'(a)') char(14b) 
          write (lu,'(''ACOMT restarted''/14''='')')
*         force a chart advance before next plot
          ndif = 0
          call exec (2,113b,input,1,5,5)
          call setlo (isynth,ifreq,synon,lu8660,lu3336) 
          call calt (lu,ludvm,ifreq,ndsel,dnoise,isynth,synon,
     *               lu8660,lu3336,10,tson,fault) 
          if (fault) go to 960
      endif 
* 
*     clear s reg, continue if nextop not reset 
      call issr(0)
      if (nextop .eq. 1) go to 200
* 
*     else store time, synth freqs for last spectrum plot 
      is = itime(2) 
      im = itime(3) 
      ih = itime(4) 
      id = itime(5) 
      syon = synon
      syof = synof
      ra50st  = ra50
      dec50st = dec50 
      rast    = aonra 
      decst   = aondec
* 
*     store spectrum on magtape then plot it
*     --------------------------------------
*     store the spectrum for transfer to acfpm
  500 do i = 1,n
          data(i) = diff(i) 
      end do
*     store housekeeping for write to tape
      indx1 = index - 1 
      tsys = tson 
      call WRTTP (lu,data)
* 
*     check plotter parameters
      ndif = ndif + 1 
      iadvnc = 0
*     reset plot count if page full 
      if (ndif .eq. 9) ndif = 1 
*     chart advance before plot if first plot on page 
      if (ndif .eq. 1) iadvnc = 1 
*     set individual plot type, y axis offset and autoscale 
      iplot = 2 
      yofset = float( intvl * ndif + intzro ) 
      yscale = 0.0
* 
*     set up plotter control, call plotter program
      iout = iplot * 1000 + ipage * 10 + iadvnc 
      call exec(24,acfpm,iout,0,0,0,lu,idata,570) 
* 
*     include spectrum in the average and plot it 
*     ------------------------------------------- 
      do i = 1,n
          avrg(i) = (avrg(i) * adav + diff(i)) / (adav + 1.0) 
      end do
      avts = (avts * adav + tson) / (adav + 1.0)
      adav = adav + 1.0 
      iadav = ifix(adav)
      tsys = avts 
*     remove slope on average spectrum
      call rmv (n,avrg) 
*     index for average plot type 
      iplot = 3 
*     store average array for data transfer 
      do i = 1,n
          data(i) = avrg(i) 
      end do
*     concatenate plotter controls and plot 
      iout = iplot * 1000 + ipage * 10
      call exec(24,acfpm,iout,0,0,0,lu,idata,570) 
      if (nextop .eq. 1) go to 250
* 
*     start new source or terminate (nextop = 2 or 3) 
*     ----------------------------------------------- 
*     advance magtape past first eof
  900 call exec (3,1300b+mtlu)
      if (nextop .eq. 2) then 
*         update parameters for new source
          if (lu .eq. 7) then 
              write (lu,'(a)') char(14b)
          else
             write (lu,'(///)') 
          endif 
          nsource = nsource + 1 
          mtfile = mtfile + 1 
          call exec (11,itime,iy) 
          write (lu,1000) iy,itime(5),itime(4),itime(3) 
          write (lu,'(''mag tape '',i6,'' :_'')') mtn 
          write (lu,1020) mtfile
          if (newsource) then 
*             start new file / source 
              write (lu,'(''Now at source '',i4,'' in file '',a)')
     *            nsource, filename 
              newsource = .false. 
              go to 30
          else
*             restart on next source
              go to 100 
          end if
      endif 
* 
*     else nextop = 3 so park telescope and end 
*     ----------------------------------------- 
  960 call track (lu,2,0.0,334.114) 
  990 if (fault) call exec (10,alarm) 
      stop 'bye from ACOMT !'
      end 
* 
      subroutine filerd (lu,filename,nsource,error,file_end,
     *    name,acephfile,nacc,nobs,kount,bw,synofs,synrf,rf)
****************************************************************
*     read source data from disc file 
      character filename*10, name*20, acephfile*20, 
     *          cbuffer*80
      logical   error, file_end 
* 
      error = .false. 
      file_end   = .false.
      open (99,file=filename,status='old',iostat=ios,err=100) 
  100 if (ios .ne. 0) then
          write (lu,'(''Error '',i4,'' opening file '',a10)') 
     *        ios, filename 
          stop 'bye !'
      endif 
      rewind (99) 
      do i = 1,nsource
          cbuffer(1:80) = ' ' 
          read (99,'(a)',iostat=ios,err=800,end=900) cbuffer(1:80)
      end do
      close (99)
      read (cbuffer,*,iostat=ios,err=800,end=800) 
     *    name,acephfile,nacc,nobs,kount,bw,synofs,synrf,rf 
      return
* 
*         error returns 
  800     write (lu,'(''Error'',i5,'' at source'',i3//)',err=810) 
     *        ios, nsource
  810     error = .true.
          return
* 
  900     write (lu,'(''End of file at source'',i4)') nsource - 1 
          file_end = .true. 
          return
      end 
* 
      INCLUDE &ACSUB
      INCLUDE &HALIM
*********************************************************************** 
