PRO calcvtm, iy,im,id
; This program reads a data file of positions and times for analysis.
  IF n_params() eq 0 THEN BEGIN
    print,'calcvtm scans a file of VTM messages.'
    print,'caldvtm,iy,im,id'
    print,'   iy = year number'
    print,'   im = month number'
    print,'   id = day number'
    print,"   Example:  calcvtm,2000,11,1"
    RETURN
  ENDIF
   theta=!pi/6.;	Angle of line of sight from Romberg sea wall.
   n_e1=-sin(theta); x component of unit vector perp to line of sight
   n_e2=cos(theta); y component of unit vector
   re=double(6380000.); 	Radius of the earth, in meters
   onemin=2.*!pi*re/(double(360.)*60.); Arc length for one minute of lattitude.
   x_0mn=double(33.3);	rough coords for Romberg, in minutes from (123 W, 37 N)
   y_0mn=double(53.2)
; Calculate the coordinates of Romberg in meters.
   x_0=x_0mn*sin(!pi/2-(37.+y_0mn/60.)*!pi/180.)*onemin
   y_0=y_0mn*onemin
   print,'x and y for Romberg: ',x_0,y_0
   string1=''
; Get date of runs to process.
;   iy=2000; year number
;   read,'iy = '+string(iy),string1
;   IF strlen(string1) ne 0 THEN iy=fix(string1)
;   im=11; month number
;   read,'im = '+string(im),string1
;   IF strlen(string1) ne 0 THEN im=fix(string1)
;   id=1; day number
;   read,'id = '+string(id),string1
;   IF strlen(string1) ne 0 THEN id=fix(string1)

; Open file with arrays of data in WDSK format.
;  example:  fname='/d4/Darmat/data/naut/vtm/2000/11/mar01101.dat'
   fname='/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)')+'.dat'
  print,fname
  listf=findfile(fname,count=ninfile)
  IF ninfile eq 0 THEN BEGIN
    print,'cannot find file ',fname
    RETURN
  ENDIF

; Read in arrays
   rdsk, xmn, fname, 4; xmn(irun,ipt)
   rdsk, ymn, fname, 5; ymn(irun,ipt)
   rdsk, tsecpst, fname, 6; tsecpst(irun,ipt)
   rdsk, npt, fname, 7; npt(irun)
   rdsk, idir, fname, 8; idir(irun)
; The array jdir will be +1 for S->N runs and -1 for N->S runs.
   jdir=2*idir-1
   nrun=n_elements(xmn(*,0))
   nptmax=max(npt)
   x=dblarr(nrun,nptmax)
   y=dblarr(nrun,nptmax)
   tpstcr=dblarr(nrun)
   xcross=dblarr(nrun)
   ycross=dblarr(nrun)
   vcross=dblarr(nrun)
   icross=intarr(nrun)
   d_perp=dblarr(nrun,nptmax)
;   d_perpmn=fltarr(nrun,nptmax)
   ds=fltarr(nrun,nptmax)
   dt=fltarr(nrun,nptmax)
   v=fltarr(nrun,nptmax)
   icr=long(0)
; Calculations in matrix form.
   x=xmn*sin(!pi/2.-(37.+ymn/60.)*!pi/180.)*onemin; arc distance in meters
   y=ymn*onemin; arc distance in meters
   d_perp=(x-x_0)*n_e1 + (y-y_0)*n_e2
   FOR irun=0,nrun-1 DO BEGIN 
      ds(irun,0)=0.
      dt(irun,0)=1.
      icr=0; Covers anomalous cases.
      FOR ipt=1, npt(irun)-1 DO BEGIN
;         flrad=(37.+ymn(irun,ipt)/60.)*!pi/ 180.; latitude of point, in radians
;         x(irun,ipt)=xmn(irun,ipt)*sin((!pi/2)-flrad)*onemin; arc length in meters
;         y(irun,ipt)=ymn(irun,ipt)*onemin; arc length in meters
; Delete This         d_perpmn(irun,ipt)=(xmn(irun,ipt)-x_omn)*n_e1+(ymn(irun,ipt)-y_omn)*n_e2
;         d_perp(irun,ipt)=(x(irun,ipt)-x_o)*n_e1+(y(irun,ipt)-y_o)*n_e2
         ds(irun,ipt)=sqrt((x(irun,ipt)-x(irun,ipt-1))^2+(y(irun,ipt)-$
            y(irun,ipt-1))^2)
         dt(irun,ipt)=tsecpst(irun,ipt)-tsecpst(irun,ipt-1)
; For runs from south to north, jdir = 1; for north to south, jdir = -1.
; icr is the number of the point just before crossing the line of
; sight.
         IF d_perp(irun,ipt)*jdir(irun) le 0 THEN icr=ipt 
         v(irun,ipt)=ds(irun,ipt)/dt(irun,ipt)        
      ENDFOR; end of loop over ipt
      icross(irun)=icr; Save the number of the point before crossing.
; Find the point where the boat crossed the line of sight from Romberg
; Center.  (Doesn't make sense if there is only one point.)
      IF npt(irun)ge 2 THEN BEGIN
; INTERPOLATION.  Use the points before and after crossing, except in
; exceptional cases.
         IF icr eq npt(irun)-1 THEN BEGIN
	    print,'problem:  irun, icr, npt(irun) = ',irun,icr,npt(irun)
         ENDIF

         i1=min([icr,npt(irun)-2]); first point for interpolation
	 i2=i1+1; second point
         tpstcr(irun)=tsecpst(irun,i1)+(tsecpst(irun,i2)-tsecpst(irun,i1))*$
            (-d_perp(irun,i1))/(d_perp(irun,i2)-d_perp(irun,i1))
         xcross(irun)=x(irun,i1)+(x(irun,i2)-x(irun,i1))*$
            (-d_perp(irun,i1))/(d_perp(irun,i2)-d_perp(irun,i1))
         ycross(irun)=y(irun,i1)+(y(irun,i2)-y(irun,i1))*$
            (-d_perp(irun,i1))/(d_perp(irun,i2)-d_perp(irun,i1))
         vcross(irun)=v(irun,i1)+(v(irun,i2)-v(irun,i1))*$
            (-d_perp(irun,i1))/(d_perp(irun,i2)-d_perp(irun,i1))
      ENDIF 
   ENDFOR
   wdsk,x,fname,11,/insert
   wdsk,y,fname,12,/insert
   wdsk,v,fname,13,/insert
   wdsk,icross,fname,14,/insert
   wdsk,xcross,fname,15,/insert
   wdsk,ycross,fname,16,/insert
   wdsk,vcross,fname,17,/insert
   wdsk,tpstcr,fname,18,/insert
   SAVE
   print,nrun,'runs
   print,'irun idir icross  tpstcr  xcross  ycross  vcross'
   FOR irun=0,nrun-1 DO BEGIN
      print,irun,idir(irun),icross(irun),tpstcr(irun),xcross(irun), $
	 ycross(irun),vcross(irun),format='(i4,i5,i7,4f8.1)'
   ENDFOR
END