;+
; NAME:
;     PNTCHEK
; PURPOSE:
;     Main routine to analyze antenna pointing data.  Responds to scans
;     of type POINT (scan code 16) or TEST DATA (scan code 1).
; CATEGORY:
;     OVRO APC INSTRUMENT CALIBRATION
; CALLING SEQUENCE:
;     pntchek[,filename][,hrec | ,after=after],output[,debug][,/cmdfile]
; INPUTS:
;     filename   the name of the file containing POINT data.  If
;                  omitted, the file DAILY.ARC in directory
;                  !DEFAULTS.WORKDIR is assumed.
;     hrec       the record number of the POINT scan header.  Data
;                  will be processed up to the next EOS segment.
;                  If this argument is given, then AFTER keyword
;                  is ignored.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     after      an optional time string of the standard form
;                  [yyyy.ddd ]hh:mm[:ss] after which to start
;                  looking for a valid POINT or TEST DATA scan.
;     cmdfile    a switch that indicates that the routine was called
;                  from a command file, so issue no modal messages, and
;                  print any error messages to <WORKDIR>\POINT.MSG.
;     debug      a switch that activates some debugging statements
;                  for getting an indication of intermediate steps.
;                  This switch is not compatible with /CMDFILE switch.
; ROUTINES CALLED:
;     openarc, getdata, lasthrec, tl_decode, decode, pntfit, nint
; OUTPUTS:
;     output     a 4-element string array containing formatted lines of
;                  output of the following form.  The two header lines are
;                  not included.  These lines are also appended to
;                  the file boresite.pnt in the directory !DEFAULTS.DBDIR.
;
;      Source Name     Date      Time   Rec#  HA(deg)  DEC(deg) Ant Axis Offset
;      ------------ ---------- -------- ----- -------- -------- --- ---- ------
;      3C454.3      1999/12/27 00:06:00    23   -5.236   16.147  1    HO   0.110
;      3C454.3      1999/12/27 00:06:00    23   -5.236   16.147  1    DO  -0.008
;      3C454.3      1999/12/27 00:15:36    23   -2.830   16.147  2    HO  -0.004
;      3C454.3      1999/12/27 00:15:36    23   -2.830   16.147  2    DO  -0.028
;
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 13-Jun-1999 by Dale E. Gary
;     20-Jul-1999  DG
;       Convert to work for either POINT or ENCODER ZERO data (NCODRZERO).
;       Now senses trajectory offsets and frequency used, through changes
;       to arguments to PNTFIT routine
;     30-Dec-1999  DG
;       Extensive changes to output format, to enable combining with
;       MOUNTCAL and WPOINT routines
;     05-Feb-2000  DG
;       Removed now unnecessary sin and cos lookup tables in GAINCOR call.
;     23-May-2000  DG
;       Changed to work with 6 or 2 elements.  The standard mode should
;       be two elements from now on, since that gives the best S/N.
;     10-Mar-2000  DG
;       Added check for short scans, so that routine can exit gracefully
;       without crashing.
;     13-Oct-2006  DG
;       Changed polarization channel from 2 (I) to 0 (R), due to new RL
;       polarization scheme.  Analyzing old data will now no longer work,
;       but there seems no reason to do that anyway.
;-
pro pntchek,filename,hrec,output,after=after,debug=debug,cmdfile=cmdfile

   ; If the /cmdfile switch is set, open a "new" POINT.MSG before doing
   ; anything.  If it exists it will be cleared.
   msglun = 0   ; Indicates POINT.MSG does not exist
   IF (keyword_set(cmdfile)) THEN $
                  openw,msglun,/get_lun,!defaults.workdir+'POINT.MSG'

   ; If no filename, assume standard DAILY.ARC file
   IF (n_elements(filename) EQ 0) THEN $
                  filename=!defaults.workdir+'DAILY.ARC'

   ; Open the file
   lun = openarc(filename,a,nrec)

   IF (n_elements(hrec) EQ 0) THEN BEGIN
      IF (keyword_set(after)) THEN hrec = lasthrec(a,nrec,tstr=after) $
                              ELSE hrec = lasthrec(a,nrec)
   ENDIF

   ; Read header info, etc.
   rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm,geometry,pRecent,pOld,cstruct)-1

   ; Verify that the specified record has either a POINT or TEST DATA
   ; scan code
   sc = header.tls.scancode
   IF (sc NE !SCAN.TEST AND sc NE !SCAN.POINT AND sc NE !SCAN.NCODRZERO) THEN BEGIN
      errstr = 'PNTCHEK: Header found is of wrong type.'
      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                                ELSE ans = widget_message(errstr,/ERROR)
      return
   ENDIF

   nant = header.nant
   bigdata = (*cstruct.pbigdata)
   pol = (*cstruct.ppol)
   ; The POINT scan has NBLK blocks of data, each of duration NOBS seconds.
   ; These parameters are set by the TRAJECTORY.  Currently this hardcodes
   ; the values for 81 TRAJECTORY, which are 18 blocks, 60 s per block.
   nblk = 18
   nobs = 60
   npersec = 1000./cfg.sampintms/obseq.nx    ; Number of cycles per second (may be <1)
   npts = nint(nobs*nblk*100./cfg.sampintms) ; Number of points expected

   ; Declare the arrays, a little larger than expected to be safe
   sn12 = fltarr(2,npts+100)
   t12  = lonarr(npts+100)
   npt = -1

   ; Get first cycle of data
   cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle)

   n = n_elements(cycle)

   id = progmeter(/init,Label='Working on scan at record: '+string(hrec,format='(I5)'))

   while(n_elements(cycle) eq n) do begin

      res = progmeter(id,npt*1.0/npts)
      npt = npt + 1
      ; Cycle is full, so save data
      bigdata(0) = cycle    ; Convert CYCLE to channelized floating array

      flagool,bigdata,header,obseq   ; Flag data for out of lock

      ; Correct for complex gain and attenuation
      smldata = gaincor(bigdata,pol,header,cfg,obseq,tcycle,GParm,Geometry)
;      if (npt lt 1080) then sn12[*,npt] = smldata[7:8,0,2] $
;                       else sn12[*,npt] = smldata[11:12,0,2]
      sn12[*,npt] = smldata[nant:nant+1,0,0]
      t12[npt] = tcycle
      cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle)
      while (n_elements(cycle) eq 1) do begin
         if (cycle eq !SEGM.HEADER or cycle eq !SEGM.EOS) then goto,break
         cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle)
      endwhile
   endwhile

break:
   free_lun,lun

   res = progmeter(id,/destroy)

   ; Require that at least half of the last block be present:
   if (npt lt npts-nobs/2) then begin
      print,'PNTCHEK: Scan at record ',hrec,' is too short.'
      output = ['Nofit','Nofit','Nofit','Nofit']
      return
   endif

   ; Truncate to actual size of data
   t12 = t12[0:npt]
   sn12 = sn12[*,0:npt]


   ; We allow NAQC seconds to acquire the source, and allow NTOL seconds
   ; on either end for "sloppy timing."
   nacq = 12
   ntol = 2

   ; Determine the start and stop indexes into the SN12 array where
   ; good data are expected.
   n = indgen(nblk)
   k1 =  nint((n*nobs+nacq+ntol)*npersec) - 1    ; Array of end indexes of bad data
   k2 = nint([0,((n+1)*nobs-ntol)]*npersec)<npt      ; Array of start indexes of bad data

   ; Set bad data parts to zero
   for i = 0, nblk-1 do begin
      sn12(*,k2(i):k1(i)) = !values.f_nan
   endfor

   k1 = nint((n*nobs+nacq+ntol)*npersec)          ; Array of start indexes of good data
   k2 = nint(((n+1)*nobs -ntol-1)*npersec)<(npt-1)  ; Array of end indexes of good data

   ; Create array to hold the averaged data
   avdat = fltarr(nblk)
   resid = fltarr(nblk)

   ; Loop over the 24 blocks (pointings) of data
   for i = 0, nblk-1 do begin

      sout = moment(sn12(0,k1[i]:k2[i]),/nan)
      cout = moment(sn12(1,k1[i]:k2[i]),/nan)
      n = n_elements(sout)
      avdat[i] = sqrt(sout[0]^2+cout[0]^2)
      resid[i] = (sqrt(sout[1]*sout[0]^2+cout[1]*cout[0]^2)/avdat[i])/sqrt(n)  ; To get standard error in the mean

   endfor

   ; Perform fit to data contained in AVDAT, returning the peak level
   ; PK, the pointing offsets in degrees PO(2), and rms error of the
   ; fit RMSERR(2)
   pntfit,avdat,resid,po,rmserr,obseq,traj,/debug

   if (sc EQ !SCAN.POINT OR sc EQ !SCAN.TEST) then begin
      olun = 0
      openu,olun,/get_lun,!defaults.dbdir+'boresite.pnt',/append

      if ((FSTAT(olun)).CUR_PTR eq 0) then begin
         out = 'Source Name     Date      Time   Rec#  HA(deg)  DEC(deg) Ant Axis Offset'
         printf,olun,out,format='(a)'
         out = '------------ ---------- -------- ----- -------- -------- --- ---- ------'
         printf,olun,out,format='(a)'
      endif

      on_ioerror,NULL

      line = string(traj.srcname,format='(A12)')+' '+header.tls.date
      ; Get best offset values
      print,traj.srcname
      print,"Pointing offsets for this location are:"
      print,po[0,*],format='("HA: ",F6.3," DEC: ",F6.3)'
      print,po[1,*],format='("HA: ",F6.3," DEC: ",F6.3)'

      ; Duration of measurement is 19 min, 12 sec, so offset time in the
      ; header by 1/4 of this in [msec], to get time of center of observation
      ; for antenna 1, and 3/4 for antenna 2
      dtmsec1 = ((19*60. + 12)*1000.)/4
      dtmsec2 = 3*((19*60. + 12)*1000.)/4
      tmsec1 = header.tls.msec+dtmsec1
      tmsec2 = header.tls.msec+dtmsec2

      ; Get Right Ascension and Declination at this time for antenna 1
      tarray = [1D,tmsec1,tmsec1^2]
      ra = total(geometry.a_ra*tarray)     ; [msec]
      dec = total(geometry.a_dec*tarray)   ; [masec]

      ; Calculate HA
      ha = lst(header.tls.year,header.tls.day,tmsec1) - ra   ; [msec]
      ha = ha*15./3600.D3
      if (ha lt -180.) then ha = 360.+ha else if (ha gt 180.) then ha = ha-360.
      dec = dec/3600.D3

      out = line+strmid(msec2str(tmsec1),0,9)+string(hrec,format='(I6)')$
                    +string(ha,dec,format='(2F9.3)')
      out2 = out+'  1    HO  '+string(po[0,0],format='(F6.3)')
      printf,olun,out2
      output = out2
      out2 = out+'  1    DO  '+string(po[0,1],format='(F6.3)')
      printf,olun,out2
      output = [output,out2]

      ; Get Right Ascension and Declination at this time for antenna 2
      tarray = [1D,tmsec2,tmsec2^2]
      ra = total(geometry.a_ra*tarray)     ; [msec]
      dec = total(geometry.a_dec*tarray)   ; [masec]

      ; Calculate HA
      ha = lst(header.tls.year,header.tls.day,tmsec2) - ra   ; [msec]
      ha = ha*15./3600.D3
      if (ha lt -180.) then ha = 360.+ha else if (ha gt 180.) then ha = ha-360.
      dec = dec/3600.D3

      out = line+strmid(msec2str(tmsec2),0,9)+string(hrec,format='(I6)')$
                    +string(ha,dec,format='(2F9.3)')
      out2 = out+'  2    HO  '+string(po[1,0],format='(F6.3)')
      printf,olun,out2
      output = [output,out2]
      out2 = out+'  2    DO  '+string(po[1,1],format='(F6.3)')
      printf,olun,out2
      output = [output,out2]
      free_lun,olun
   endif else if (sc EQ !SCAN.NCODRZERO) then begin
      print,"Encoder zero adjustments for optimum pointing are:"
      print,nint(po[0,*]*(2.^17/360.)),format='("HA: ",I4," DEC: ",I4)'
      print,nint(po[1,*]*(2.^17/360.)),format='("HA: ",I4," DEC: ",I4)'
      print,"Encoder zeroes after correction are:"
      print,nint(po[0,*]*(2.^17/360.))+cfg.zencoder[*,0],format='("HA: ",I4," DEC: ",I4)'
      print,nint(po[1,*]*(2.^17/360.))+cfg.zencoder[*,1],format='("HA: ",I4," DEC: ",I4)'
   endif

return
END