PRO showak1,nf1,nf2
; This procedure makes a plots of boat wakes as seen in the current-
; velocity measurements from the boat-wake 
; setup at Tiburon.  A series of current-velocity files are read in one
; by one.  For each file, the appropriate file of ferry runs from the
; VTM system is read in, as well as files for observed tidal current
; and stage.  A graph of cvss1 is made for each recorded ferry run.
; (As of June 2001, the only Larkspur monohull transmitting is the
; Marin.)  The time of passage of the ferry from VTM data is
; indicated, and the tidal stage and current velocity are posted.  The
; graph is either displayed on the screen or written to
; a printer file.  The series of files is selected in the prodedure 
; read1.pro.  RWB, June 6, 2001.
  IF n_params() ne 2 THEN BEGIN
    print,'showak reads files of current velocity and files of VTM runs, '
    print,'and displays wakes made by logged ferries.'
    print,'Incorrect number of parameters.'
    print,'Calling sequence:  showak,nfile1,nfile2'
    print,'nfile1 and nfile2 are file numbers in cus-file sequences.'
    print,'To select the sequence, edit showak.pro.'
    print,'Example:  showak,0,21'
    RETURN
  ENDIF
  COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
  COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
     d0,d1,dw0,dw1
  nfile1=nf1;
  nfile2=nf2;
  init; Initialization.
;	Start reading in files of current velocities.
  FOR ifile=nfile1,nfile2 do begin
    readcv; reads one file of current-velocity data.
    readvtm; read in the dayfile of VTM data for this run.
;	Start of loop over ferry runs from VTM data.  We will exit the loop if
;	either there are no more ferry runs or if a ferry run takes place (at
; 	least partially) after the current file of cv data (requiring
;       reading in another file of cv data).
    irun=0; Run counter.
    nrun=n_elements(tpstcr)
    WHILE idone eq 0 DO BEGIN
      findwake; locate a wake from the VTM file in this cv file.
      IF idone eq 1 THEN BREAK; need to read another cv file.
      dspwake; display the wake, with annotations.
    ENDWHILE
    print,'end of printing of wakes for file number ',ifile
  ENDFOR
  endall; Close things out.
  SAVE
END; ****************************************************************

; ****** Subroutines ************

PRO init
; Initialization for showak1.  RWB. June 6 2001.
; nfile1, nfile2, ifile:  first, last and current .cur file number
; iprint:  = 1 if printout goes to a PostScript file.
; idone:  goes to 1 when there are no more wakes to display.
; imarksf:  = 1 if scheduled ferries are to be shown.
; dday0, dday1:  Time before and after VTM crossing to include in wake.
; itype:
; bday
; vtmfile:  file name of current dayfile of vtm data.
; irun, nrun:  current run number and number of runs in vtm dayfile.
; cvss1old, d1pstold: cv and time arrays for previous .cur file.
; cvss1cat,d1pstcat: concatenated from previous and current .cur file.
; tpstrcr:  time-of-crossing array for current VTM dayfile.
; idir: direction array for current VTM dayfile.
; d0, d1: first and last times from current .cur file.
; dwo, dw1: first and last times for wake from VTM dayfile.
  COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
  COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
     d0,d1,dw0,dw1
  COMMON getrtcw,rmcfile,rmwfile,vrmc,trmc,hrmw,trmw
  vtmfile=''
  rmcfile=''
  rmwfile=''
  jprint=1; printout level
  dday0=3./long(60.*24.); Show this much before time of ferry passage.
  dday1=8./long(60.*24.); Show this much after the ferry.
  dtj70=julday(1,1,1970,0,0,0); Add this to system time to get Julian day.
  imarksf=0; 1 -> mark scheduled ferries.
  itype=0
  bday=double(0.)
  iprint=0
  chr1=''
  read,'Output to printer?' + $
    '(enter 1 for printer, return for screen)',chr1
  iprint = fix(chr1)
  IF iprint eq 1 THEN BEGIN
    set_plot,'ps'
    device,/landscape
  ENDIF
  !p.multi=[0,2,2,0,0]; four plots per page
  !p.multi=0;	one plot per page
  RETURN
END; ***************************************************************

PRO readcv; ********************************************************
; This routine reads the next file of Tiburon current-velocity data.
; It is concatenated with the previous file, unless it is the first.
  COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
  COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
     d0,d1,dw0,dw1
   read1,ifile; read in one file of current-velocity data.
   restore
   infname=infile; put file name in common.
   idone=0; When we are finished with this file, idone will be set to 1.
   npts=n_elements(cvss1); number of points in arrays.
   xmax=float(npts)
   day1pst=day1-.3333333; change DOS time (GMT) to PST
; Fill concatenated plotting arrays.
   IF ifile eq nfile1 THEN BEGIN
      cvss1cat=cvss1
      d1pstcat=day1pst
   ENDIF ELSE BEGIN
      cvss1cat=[cvss1old,cvss1]
      d1pstcat=[d1pstold,day1pst]
   ENDELSE
   cvss1old=cvss1; Save data for next time.
   d1pstold=day1pst
   d0 = day1pst(0); start of current cv file.
   d1 = day1pst(npts-1); end
   IF jprint ge 1 THEN print,'d0, d1 = ',d0,d1
   RETURN
END; *****************************************************************

PRO readvtm; *********************************************************
; This routine opens the day-file of VTM data which corresponds to the
; run just opened.
  COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
  COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
     d0,d1,dw0,dw1
    timtx1, d0,iy0,im0,id0,tday0; starting year, month, day number.
    vtmf='/d4/Darmat/data/naut/vtm/'+string(iy0,format='(i4.4)')+ $
      '/'+string(im0,format='(i2.2)')+'/mar'+string(iy0 mod 10, $
      format='(i1.1)')+ $
      string(im0,format='(i2.2)')+string(id0,format='(i2.2)')+'.dat'
    IF jprint ge 1 THEN print,'vtmf = ',vtmf
    IF vtmf eq vtmfile THEN RETURN; The right data are already loaded.
    vtmfile=vtmf;
    list=findfile(vtmfile,count=nvtmfile); Look for the file.
    IF nvtmfile eq 0 THEN BEGIN
      print,'Cannot find vtm file ',vtmfile
      idone = 1; This should make us move on to the next file of
               ; current data.
    ENDIF ELSE BEGIN
      rdsk,tpstcr,vtmfile,18; These are the crossing times for all runs on
                          ; this date.
      rdsk,idir,vtmfile,8; Direction of runs; 1 means S->N, 0 means N->S.
      tpstcr=tpstcr/86400.+julday(im0,id0,iy0,0)-dtj70; convert to DOS days.
;      tpstcr=tpstcr/86400.+dosday(iy0,im0,id0); convert the array to Dos days.
      IF jprint ge 1 THEN print,tpstcr
   ENDELSE
   RETURN
END; *****************************************************************

PRO findwake; ********************************************************
; This routine looks in the array for VTM ferry crossings for one day
; to find a crossing contained within the cv data in memory.
;** Decision point.  We have to compare d0 and d1 with the various crossing
; times.  If any crossing time falls in this input file, plot it out.
; irun will be a pointer, which will be moved every time a run is plotted
; or a new vtm file is read.  
; One assumption will be made about the data - no ferry run overlapping two
; PST days will be considered.
;
;    If the end of the wake is before the start of this file, increment the
;    ferry run number.  Do this until all runs for this date have been
;    considered.
;    When out of ferry runs, a new .cur file needs to be read in.
;
;    If the end of the wake is after the end of this file, a new .cur file
;    needs to be read in.
;
;    Now we assume that the end of this wake is in this file.  We plot from
;    the special plotting arrays, with two files concatenated.
;
; Try and find a wake with the end of its window in this file.  We
; may inherit a value of irun from the last crossing.
  COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
  COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
     d0,d1,dw0,dw1
   dw0=tpstcr(irun)-dday0; start time for wake window
   dw1=dw0+dday0+dday1; end time for wake window
   IF jprint > 2 THEN print,'irun, dw0, dw1 = ',irun, dw0,dw1
; The next loop is to increment up to the first run which ends in the
; current data file (and thus requires the current data file).
   WHILE dw1 lt d0 and irun lt nrun-1 DO BEGIN
      irun=irun+1
      dw0=tpstcr(irun)-dday0; start time for wake window
      dw1=dw0+dday0+dday1; end time for wake window
      print,'irun, dw0, dw1 = ',irun, dw0,dw1
   ENDWHILE
   IF dw1 lt d0 THEN idone=1; no more wakes in this file.
   IF dw1 gt d1 THEN idone=1; need data from next file.
   IF jprint ge 1 THEN print,'irun,idone,dw0,dw1= ',irun,idone,dw0,dw1
   RETURN
END; ***************************************************************

PRO dspwake; ******************************************************
; This routine graphs the current velocity for a wake from a VTM
; run.
   COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
   COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
      d0,d1,dw0,dw1
   chr1=''; character buffer, for dialogue
   black=0
   white=255+long(256)*(255+long(256)*255)
   IF iprint eq 0 THEN read,'return to continue',chr1; delay
; Find the index of the current-velocity array corresponding to the
; ferry crossing.
   dummy=min((d1pstcat-tpstcr(irun))^2,icrs); find icrs, the index
	                     ; of the velocity point nearest the crossing.
;   offset=double(fix(d1pstcat(0))); integer day number (floated)
   dummy=label_date(date_format='%h:%i'); set up for axis labels in hh:mm
   plot, d1pstcat+dtj70, $
      background=white,color=black, $
      cvss1cat-total(cvss1cat(icrs-10:icrs+10))/21., $
      xrange=[dw0+dtj70,dw1+dtj70], $
      title=infname,xtitle= $
      'time (PST)', $
      ytitle='Velocity (m/sec)', $
      yrange=[-.2,.2], $
      xstyle=1, $
      xtickformat='label_date', $
      xminor=3
   datestr=label_date(1,1,dw0+dtj70,1,date_format='%m %d, %y'); Generate date
      ;string.
   xyouts,0.15,0.13,datestr,/normal,color=black; write date in LL corner of plot.
; Now mark the crossing time.
   oplot,[-.0001,0,.0001,0,0]+tpstcr(irun)+dtj70, $
      [.14,.13,.14,.13,.16],color=black
; And print the direction.
   dirstr='S -> N'
   IF idir(irun) eq 0 THEN dirstr='N -> S'
   xyouts,tpstcr(irun)+dtj70-.0004,.175,dirstr,charsize=1.5,color=black
;     Show the tidal stage and current velocity at Richmond.
   getrtcw,irmc,rmc,irmw,rmw
   xyouts,0.65,0.20,'current '+string(rmc,format='(f5.2)')+' m/s', $
      /normal,charsize=1.5,color=black
   xyouts,0.65,0.14,'stage   '+string(rmw,format='(f5.2)')+' m', $
      /normal,charsize=1.5,color=black
; The next section marks the predicted times of ferry passage.
   WHILE (imarksf eq 1) and (bday lt d1) DO BEGIN
      IF (bday gt d0) and (bday lt d1) THEN BEGIN
         print,'itype, bday = ', bday,itype
         xyouts,(bday-offset)*24.,0.17, $
         strtrim(string(itype),1),color=black
         oplot,[-.004,0,.004,0,0]+(bday-offset)*24., $
	    [.14,.13,.14,.13,.16],color=black
         print,'(bday-offset)*24. = ',(bday-offset)*24.
      ENDIF
      readf,lun1,'$(f15.5,i8)',bday,itype
      bday=bday-25569.	; conversion from Excel to DOS.
      print,'Just read from file:  bday, itype = ',bday,itype
   ENDWHILE
   IF irun lt nrun-1 THEN BEGIN
      irun=irun+1
      dw0=tpstcr(irun)/86400.-dday0; start time for wake window
      dw1=dw0+dday0+dday1; end time for wake window
   ENDIF ELSE idone=1
   RETURN
END; ****************************************************************

PRO endall; *********************************************************
; close things out.
  COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
  COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
     d0,d1,dw0,dw1
  IF iprint eq 1 THEN BEGIN
    device,/close
    set_plot,'x'
    ; spawn,'lp -o simplex idl.ps' ; this prints on only one side.
    ; spawn,'lp -s idl.ps'; print; -s for long jobs.
    ; spawn,'lp idl.ps'; standard print.
  ENDIF
END; *****************************************************************

PRO getrtcw,irmc,rmc,irmw,rmw; ***************************
   COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
   COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
      d0,d1,dw0,dw1
   COMMON getrtcw,rmcfile,rmwfile,vrmc,trmc,hrmw,trmw
   caldat,dw0+dtj70,im0,id0,iy0; Get calendar date for start of wake.
;   timtx1, dw0,iy0,im0,id0,tday0; starting year, month, day number.
   print,'im0, id0, iy0 = ',im0,id0,iy0
   rmcf='/d4/Darmat/data/naut/tides/'+string(iy0,format='(i4.4)')+ $
      '/'+string(im0,format='(i2.2)')+'/arr'+string(iy0 mod 10, $
      format='(i1.1)')+ $
      string(im0,format='(i2.2)')+'00.rmc'
   rmwf=strmid(rmcf,0,strlen(rmcf)-1)+'w'
   IF jprint ge 1 THEN BEGIN
      print,'rmcf = ',rmcf
      print,'rmwf = ',rmwf
   ENDIF
   irmc=1; 1 means a good value was found; 0 means not.
   irmw=1
   rmc=0.; 0. means no good data.
   rmw=0.

; ***** Start of block to get current value.
   IF rmcf ne rmcfile THEN BEGIN; Need a new data file.
      list=findfile(rmcf,count=nrmcf); Look for the file.
      IF nrmcf eq 0 THEN BEGIN
         IF jprint ge 1 THEN print,'Cannot find rmc file ',rmcf
         irmc = 0; Data not available.
      ENDIF ELSE BEGIN
         rmcfile=rmcf;
         rdsk,vrmc,rmcfile,1; Velocities, projected on 340 deg.
         rdsk,trmc,rmcfile,2; Times in days since start of Jan 1 1970.
      ENDELSE
   ENDIF; End of search for a new file.                   
   IF irmc ne 0 THEN BEGIN; Search for a good value in the data.
      IF jprint ge 2 THEN print,'vrmc(0:4) = ',vrmc(0:4)
      IF jprint ge 2 THEN print,'trmc(0:4) = ',trmc(0:4)
      dtsq=min((trmc-tpstcr(irun))^2,jrmc); jrmc is index of
         ;  closest value to crossing time.
      rmc=vrmc(jrmc)
      IF jprint ge 2 THEN BEGIN
         print,'jrmc,dtsq = ',jrmc,dtsq
         print,'trmc(jrmc-5:jrmc+5) = ',trmc(jrmc-5:jrmc+5)
         print,'vrmc(jrmc-5:jrmc+5) = ',vrmc(jrmc-5:jrmc+5)
      ENDIF
      IF rmc eq 0. THEN BEGIN; bad data; search 5 points on each side.
         diplus=0
         diminus=0
         FOR di = 1,5 DO BEGIN
	    IF diplus eq 0 and trmc(jrmc+di) ne 0. THEN $
	       diplus=di
	    IF diminus eq 0 and trmc(jrmc-di) ne 0. THEN $
	       diminus=di
         ENDFOR
; See if a usable point was found on each side, and if so interpolate.
	 IF diplus ne 0 and diminus ne 0 THEN rmc=(diplus* $
	       vrmc(jrmc-diminus)+diminus*vrmc(jrmc+diplus))/ $
	       (diplus+diminus)
      ENDIF
      IF rmc eq 0. THEN irmc = 0; we didn't find a usable value.
   print,'rmc = ',rmc
   ENDIF; ***** end of search for current value.

; ***** Start of block to get stage value.
   IF rmwf ne rmwfile THEN BEGIN; Need a new data file.
      list=findfile(rmwf,count=nrmwf); Look for the file.
      IF nrmwf eq 0 THEN BEGIN
         IF jprint ge 1 THEN print,'Cannot find rmw file ',rmwf
         irmw = 0; Data not available.
      ENDIF ELSE BEGIN
         rmwfile=rmwf;
         rdsk,hrmw,rmwfile,1; Stage.
         rdsk,trmw,rmwfile,2; Times in days since start of Jan 1 1970.
      ENDELSE
   ENDIF; End of search for a new file.                   
   IF irmw ne 0 THEN BEGIN; Search for a good value in the data.
      IF jprint ge 2 THEN print,'hrmw(0:4) = ',hrmw(0:4)
      IF jprint ge 2 THEN print,'trmw(0:4) = ',trmw(0:4)
      dtsq=min((trmw-tpstcr(irun))^2,jrmw); jrmc is index of
         ;  closest value to crossing time.
      IF jprint ge 2 THEN BEGIN
         print,'jrmw,dtsq = ',jrmw,dtsq
         print,'trmw(jrmw-5:jrmw+5) = ',trmw(jrmw-5:jrmw+5)
         print,'hrmw(jrmw-5:jrmw+5) = ',hrmw(jrmw-5:jrmw+5)
      ENDIF
      rmw=hrmw(jrmw)
      IF rmw eq 0. THEN BEGIN; bad data; search 5 points on each side.
         diplus=0
         diminus=0
         FOR di = 1,5 DO BEGIN
	    IF diplus eq 0 and trmw(jrmw+di) ne 0. THEN $
	       diplus=di
	    IF diminus eq 0 and trmw(jrmw-di) ne 0. THEN $
	       diminus=di
         ENDFOR
; See if a usable point was found on each side, and if so interpolate.
	 IF diplus ne 0 and diminus ne 0 THEN rmw=(diplus* $
	       hrmw(jrmw-diminus)+diminus*hrmw(jrmw+diplus))/ $
	       (diplus+diminus)
      ENDIF
      IF rmw eq 0. THEN irmw = 0; we didn't find a usable value.
   print,'rmw = ',rmw
   ENDIF; ***** end of search for stage value.
END; *****************************************************************

PRO getfs; ***********************************************************
; This routine gets scheduled ferries.
  COMMON control,nfile1,nfile2,ifile,iprint,jprint,idone,imarksf, $ 
     dday0,dday1,itype,bday,vtmfile,irun,nrun,dtj70,infname
  COMMON data,cvss1old,d1pstold,cvss1cat,d1pstcat,tpstcr,idir, $
     d0,d1,dw0,dw1
  fname='/d4/Darmat/data/naut/ferry/fer9711.prn'
  openr,lun1,fname,/get_lun
END


PRO timtx1, day,iyear,imonth,iday,tday
; This procedure takes a double-precision day number and calculates the
; year, month, day, and seconds in the day.  The day number is assumed to
; start at 0.0 at the start of January 1, 1970.  (This is the MS-DOS
; convention for system time.)
  idays=long(day)
  FOR iy=2010,1969,-1 DO BEGIN
    id=dosday(iy,1,1)
    IF id le idays THEN BREAK
  ENDFOR
  iyear=iy
  FOR im=12,1,-1 DO BEGIN
    id=dosday(iyear,im,1)
    IF id le idays THEN BREAK
  ENDFOR
  imonth=im
  iday=idays-id+1
  tday=day-idays
;  print,day,idays,iyear,imonth,iday,tday
END