PRO scanvtm,iy,im,id
; This program scans a file of VTM messages.  They are assumed to all
; come from a single ship of interest, the Marin of the Golden
; Gate Ferry line.  The data are broken up into ferry runs, from the SF
; Ferry Building to the Larkspur Ferry terminal, or vice versa.  A
; minimum of calculation is done in this program.
; Note on units.  The positions xm and ym and xmn and ymn are measured
; in minutes of arc from the position 37 deg north, 123 deg west, with
; x positive to the East, and y positive to the North.
  IF n_params() eq 0 THEN BEGIN
    print,'scanvtm scans a file of VTM messages.'
    print,'scanvtm,iy,im,id'
    print,'   iy = year number'
    print,'   im = month number'
    print,'   id = day number'
    print,"   Example:  scanvtm,2000,11,1"
    RETURN
  ENDIF
   string1=''; input line buffer
   nline=5; number of lines to print out, for checking purposes.
   nerr1=0; counter for input lines not for the Marin
   nerr2=0; counter for lines out of order
   nerr3=0; counter for bad first line
   nerr4=0; counter for bad second line
   nerr5=0; counter for zero coordinates

   npos=5000; maximum number of positions expected from the Marin in one day.
; xm and ym are the positions in minutes of all the loggings of the Marin in 
; the day being processed here.  They are measured from the point 123
; degrees west longitude, 37 degrees north lattitude.  Note that x is 
; measured EAST from this point.
   xm=dblarr(npos);	
   ym=dblarr(npos)
   tsp=lonarr(npos)
   nr=20; maximum number of runs expected from the Marin in one day.
   nrpt=400; maximum number of positions expected from a single run.
; The arrays below hold the same numbers as xm, ym and tsp, but
; separated into ferry runs, and only points where the ferry was under
; way.
   xmn=dblarr(nr,nrpt)
   ymn=dblarr(nr,nrpt)
   tsecpst=dblarr(nr,nrpt)
   npt=intarr(nr)
   idirect=intarr(nr)
   ipos1=intarr(nr); position index of first point in the run
   ipos2=intarr(nr); last point
; The positions of the terminals are taken from the actual loggings.
   xmnN=29.487	; Larkspur Landing
   ymnN=56.652	; Larkspur Landing
   rmnN=1.	; Radius of error circle, Larkspur Landing
   xmnS=36.475	; San Francisco Ferry Terminal
   ymnS=47.779	; San Francisco Ferry Terminal
   rmnS=1.	; Radius of error circle, San Francisco Ferry Terminal

; Read in the parameters of the file to analyze.  New values can be entered
; from the console; carriage return only means keep the default value.
;   iy=2000; calendar year
;   read,'iy = '+string(iy),string1
;   IF strlen(string1) ne 0 THEN iy=fix(string1)
;   im=11; calendar month
;   read,'im = '+string(im),string1
;   IF strlen(string1) ne 0 THEN im=fix(string1)
;   id=1; calendar day
;   read,'id = '+string(id),string1
;   IF strlen(string1) ne 0 THEN id=fix(string1)
; Piece together the input and output file names.
   infile='/d4/Darmat/data/naut/vtm/'+string(iy,format='(i4.4)')+ $
      '/'+string(im,format='(i2.2)')+'/mar'+string(iy mod 10,format='(i1.1)')+ $
      string(im,format='(i2.2)')+string(id,format='(i2.2)')+'.log'
   outfile=strmid(infile,0,strlen(infile)-3)+'dat'
   print,'infile  = ',infile
   print,'outfile = ',outfile
;   infile='/d4/Darmat/data/naut/vtm/2000/11/mar01101.log'
  listf=findfile(infile,count=ninfile)
  IF ninfile eq 0 THEN BEGIN
    print,'cannot find file ',infile
    RETURN
  ENDIF
   openr,lunin,infile,/get_lun
   print,'lunin,  eof(lunin ): ',lunin, eof(lunin)
   iline=long(0)
   ipos=long(0)
   icross=long(0)
   irun=-1; index for runs identified
   imove=0; flag which is set if the Ferry is on a run.
   idir=0; direction flag (1 = northward, 0 = southward)

   WHILE eof(lunin) eq 0 DO BEGIN
      readf,lunin,string1
      iline=iline+1
      IF iline lt nline then print,string1
      oneline,string1,nstr,strvar; Parse the comma-delimited strings.
      IF strvar(0) ne '$PRDCV' or strvar(1) ne 'Marin' THEN BEGIN
	 pskip,string1,'Wrong type of record'
	 nerr1=nerr1+1
      ENDIF ELSE IF strvar(3) ne 'SFRDV' THEN BEGIN
         pskip,string1, 'Record out of order'
	 nerr2=nerr2+1
      ENDIF ELSE BEGIN
;	 print,'start of a three-line sequence'
	 IF strvar(4) ne 3 or strvar(5) ne 'Marin' or strvar(8) ne $
	    '0003-7381' THEN BEGIN
	    pskip,string1,'Bad first line'
	    nerr3=nerr3+1
         ENDIF ELSE BEGIN
	    IF eof(lunin) THEN BREAK
            readf,lunin,string1
            iline=iline+1
            IF iline lt 6 THEN print,string1
            oneline,string1,nstr,strvar; Parse the comma-delimited strings.
	    IF strvar(3) ne 'GPGGA' THEN BEGIN
	       pskip,string1,'Bad second line'
	       nerr4=nerr4+1
            ENDIF ELSE IF long(strvar(5)) eq 0 or long(strvar(7)) eq 0 THEN BEGIN
; This test is to identify the source of some spurious zero coordinates.
	       print,'error 5; string = ',string1
	       print,'strvar(5), strvar(7) = ',strvar(5),strvar(7)
	       print,'strlen(5), strlen(7) = ',strlen(strvar(5)), $
		  strlen(strvar(7))
	       nerr5=nerr5+1
; We need to resynchronize, assuming the group of three is still good.
	       IF eof(lunin) THEN BREAK
               readf,lunin,string1; read third line         
               iline=iline+1   
	    ENDIF ELSE BEGIN
; We have identified a good logging of position.
; There are three numbers of interest in the messages:  time (GMT),
; lattitude, and longitude.
               time=long(strvar(4)); hhmmss
	       flat=double(strvar(5)); lattitude in ddmm.mmm
	       flong=double(strvar(7)); longitude in dddmm.mmm
	       IF flat eq 0. or flong eq 0. THEN stop

; Convert the time (GMT, in hhmmss format, into time of the day in
; seconds, in PST.
               i1=time
               isec=i1 mod 100
               i2=i1-isec
               imin=(i2 mod 10000)/100
               i3=i2-imin*100
               ihr=i3/10000
               tsecgmt=isec+60*imin+3600*ihr
               tsp(ipos)=(tsecgmt + 16*long(3600)) mod (24*long(3600))

;JUNK               fmin=flat mod 100.
;JUNK               flrad=((flat-fmin)/100.*60.+fmin)*!pi/180./60.
               xm(ipos)=60.-flong+12200.

               ym(ipos)=flat-3700.  
;***** The position is read in.  Now identify start of run, end of run,
;***** and during run conditions.
	       x=xm(ipos)	; short name for convenience
	       y=ym(ipos)
;	       print,imove,idir,irun,x,y
; imove is 1 if the ferry was between terminals at the last logging.
               IF imove eq 0 THEN BEGIN
; See if the ferry is at one of the terminals or another.  If so, record
; the direction of a run from that terminal.  If not, start a run.
	          IF (x-xmnN)^2+(y-ymnN)^2 lt rmnN^2 THEN idir=0 $; run S
                  ELSE IF (x-xmnS)^2+(y-ymnS)^2 lt rmnS^2 THEN idir=1 $; run N
	          ELSE BEGIN
; The ferry has just left a terminal.  Start a new run and start saving
; data.
		     imove=1; Now the ferry is moving.
		     irun=irun+1
		     ipt=0
		     xmn(irun,ipt)=x
		     ymn(irun,ipt)=y
		     tsecpst(irun,ipt)=tsp(ipos)
		     npt(irun)=ipt+1; accumulating
		     idirect(irun)=idir; accumulating
		     ipos1(irun)=ipos; index in one-dim arrays of starting pt.
		     ipt=ipt+1
                  ENDELSE
	       ENDIF ELSE BEGIN
; A run is in progress.  Check for end of run.  If not, add data and
; increment the run number.
		  imove=0; Will be reset to 1 below if necessary.
	          IF (x-xmnN)^2+(y-ymnN)^2 lt rmnN^2 THEN idir=0 $
                  ELSE IF (x-xmnS)^2+(y-ymnS)^2 lt rmnS^2 THEN idir=1 $
	          ELSE BEGIN
		     imove=1; Shouldn't have been changed.
		     xmn(irun,ipt)=x
		     ymn(irun,ipt)=y
		     tsecpst(irun,ipt)=tsp(ipos)
		     npt(irun)=ipt+1; accumulating
		     idirect(irun)=idir; accumulating
		     ipos2(irun)=ipos; will be the position index of last point.
		     ipt=ipt+1
		  ENDELSE
               ENDELSE

	       IF eof(lunin) THEN BREAK
               readf,lunin,string1; read third line         
               iline=iline+1   
	       ipos=ipos+1
	    ENDELSE
	 ENDELSE
      ENDELSE
   ENDWHILE
   close,lunin
   free_lun,lunin

;***** Final stage; resize the arrays and write them out.
   npos=ipos
; truncate one-dimensional arrays
   xm=xm(0:npos-1)
   ym=ym(0:npos-1)
   tsp=tsp(0:npos-1)
   nrun=irun; number of runs with points
   npt=npt(0:nrun-1)
   idirect=idirect(0:nrun-1)
; Truncate two-dimensional arrays.
   nptmax=max(npt)
   xmn=xmn(0:nrun-1,0:nptmax-1)
   ymn=ymn(0:nrun-1,0:nptmax-1)
   tsecpst=tsecpst(0:nrun-1,0:nptmax-1)
   wdsk,xm,outfile,1,/new
   wdsk,ym,outfile,2,/insert
   wdsk,tsp,outfile,3,/insert
   wdsk,xmn,outfile,4,/insert
   wdsk,ymn,outfile,5,/insert
   wdsk,tsecpst,outfile,6,/insert
   wdsk,npt,outfile,7,/insert
   wdsk,idirect,outfile,8,/insert
   wdsk,ipos1,outfile,9,/insert; position index of first point in the run.
   wdsk,ipos2,outfile,10,/insert; last point
   
   print,iline,' lines read.
   print,npos,' positions
   print,nrun,' runs found.
   print,'errors: ',nerr1,nerr2,nerr3,nerr4,nerr5
   print,'irn idr    tbeg   xmbeg   ymbeg    tend   xmend   ymend'+ $
     ' delta-t  length'
   FOR irun=0,nrun-1 DO BEGIN
      nn=npt(irun)-1
      print,irun,idirect(irun), tsecpst(irun,0)/3600.,xmn(irun,0),ymn(irun,0), $
	 tsecpst(irun,nn)/3600.,xmn(irun,nn),ymn(irun,nn), $
	 tsecpst(irun,nn)/3600.-tsecpst(irun,0)/3600., $
	 sqrt((xmn(irun,nn)-xmn(irun,0))^2+(ymn(irun,nn)-ymn(irun,0))^2), $
	 FORMAT='(i3,i4,8f8.3)'
   ENDFOR

   SAVE
END