% This script reads a series of transects (a run), carries out
% calculations, and writes the results to a run structure.
% RWB, June 28, 1998.
   initrun; % Sets up constants and control variables.

   ens.dtrans=run.par.dtrans;
   ens.thetap=run.par.thetap;
   jtran1=run.par.jtran1;
   jtran2=run.par.jtran2;
   nens=run.par.nensmax; % Number of ensembles to provide for.
   nd=run.par.ndmax; % Number of depth bins to provide for.
   ntran=jtran2-jtran1+1; % Number of transects.
   %%%%%% Initialize run arrays.
   run.vp=zeros(nd,nens,ntran);
   run.WNmax=zeros(nens,ntran);
   run.vflag=zeros(nens,ntran);
   run.dist=zeros(nens,ntran);
   run.pstsec=zeros(nens,ntran);
   run.pstsecav=zeros(ntran,1);
   run.xens=zeros(3,nens,ntran);
   run.BTV=zeros(4,nens,ntran);
   run.verav=zeros(nens,ntran);
   run.verrms=zeros(nens,ntran);
   run.depth=zeros(nd,ntran);
   run.dbot=zeros(nens,ntran);
   run.WNmaxmax=zeros(ntran,1);
   run.nens=zeros(ntran,1);
   run.madegood=zeros(ntran,1);
   itran=0;
   for jtran=jtran1:jtran2
      ens.itran=jtran;
      itran=itran+1;
      %%%%%% Open disk file for this transect.
      jstr=int2str(jtran);
      if length(jstr)==1, jstr=['00',jstr];, end;
      if length(jstr)==2, jstr=['0',jstr];, end;
      filename=[run.par.infile(1:length(run.par.infile)-8),jstr, ...
	 run.par.infile(length(run.par.infile)-4:length(run.par.infile))];
      ens.fidin=fopen(filename,'r');
      %%%%%% Initialize transect arrays.
      tran.vp=zeros(nd,nens);
      tran.WNmax=zeros(nens,1);
      tran.vflag=zeros(nens,1);
      tran.dbot=zeros(nens,1);
      tran.pstsec=zeros(nens,1);
      tran.xens=zeros(3,nens);
      tran.BTV=zeros(4,nens);
      tran.verav=zeros(nens,1);
      tran.verrms=zeros(nens,1);

      iens=0;
      ens.EOF=0;
      while ens.EOF == 0
	 iens=iens+1; % First ensemble is 1, MatLab style.
         ens.iens=iens; % caladens needs to recognize the first
	    % ensemble.
         ens=rdadens(ens);
         ens = caladens(ens);
         tran.vp(1:ens.WNmax,iens)=ens.vw(5,1:ens.WNmax); % transfer one
	    % ensemble to 2-d array.
	 tran.WNmax(iens)=ens.WNmax;
	 tran.vflag(iens)=ens.vflag;
	 tran.dbot(iens)=ens.dbot;
	 tran.pstsec(iens)=ens.tsec;
	 tran.xens(:,iens)=ens.xens(:);
	 tran.BTV(:,iens)=ens.BTV(:);
	 tran.verav(iens)=ens.verav;
	 tran.verrms(iens)=ens.verrms;
	 %%%%%% seek the end of file (work around a bug in Matlab!)
	 junk1=fread(ens.fidin,1,'uint8');
	 % position=ftell(ens.fidin)
	 ens.EOF = feof(ens.fidin);
	 fseek(ens.fidin,-1,0); % go back one byte
	 %%%%%%
      end
      fclose(ens.fidin);
      tran.nens=iens;
% Here we calculate the displacement of each ensemble end-point from
% the west end of the transect.
      xwest=tran.xens(:,1); % Take first ensemble as West end.
      xeast=tran.xens(:,tran.nens);
      if tran.xens(1,tran.nens)<tran.xens(1,1)% But now test.
	 xwest=tran.xens(:,tran.nens); % Last ensemble is West end.
	 xeast=tran.xens(:,1);
      end
      tran.madegood=sqrt((xeast-xwest)'*(xeast-xwest));
      for ien=1:tran.nens
         tran.dist(ien)=(xeast(1:3)-xwest(1:3))' ...
	    *(tran.xens(1:3,ien)-xwest(1:3))./tran.madegood;
      end
      tran.WNmaxmax=max(tran.WNmax);
      if tran.nens>0, tran.pstsecav=mean(tran.pstsec(1:tran.nens));, end
      tran.depth=ens.depth(1:tran.WNmaxmax);
      vplot=tran.vp(1:tran.WNmaxmax,1:tran.nens);

      %%%%%% Now store transect results in run structure.
      run.vp(1:tran.WNmaxmax,1:tran.nens,itran)=tran.vp( ...
	 1:tran.WNmaxmax,1:tran.nens); % velocity profile
      run.WNmax(1:tran.nens,itran)=tran.WNmax(1:tran.nens);
      run.dist(1:tran.nens,itran)=tran.dist(1:tran.nens);
      run.pstsec(1:tran.nens,itran)=tran.pstsec(1:tran.nens);
      run.pstsecav(itran)=tran.pstsecav;
      run.xens(:,1:tran.nens,itran)=tran.xens(:,1:tran.nens);
      run.BTV(:,1:tran.nens,itran)=tran.BTV(:,1:tran.nens);
      run.verav(1:tran.nens,itran)=tran.verav(1:tran.nens);
      run.verrms(1:tran.nens,itran)=tran.verrms(1:tran.nens);
      run.WNmaxmax(itran)=tran.WNmaxmax;
      run.vflag(1:tran.nens,itran)=tran.vflag(1:tran.nens);
      run.depth(1:tran.WNmaxmax,itran)=tran.depth(1:tran.WNmaxmax);
      run.dbot(1:tran.nens,itran)=tran.dbot(1:tran.nens);
      run.nens(itran)=tran.nens;
      run.madegood(itran)=tran.madegood;
   end
   run.ntran=ntran;
   run.WNmxmxmx=max(run.WNmaxmax);
   run.nensmax=max(run.nens);

   % Interpolate for logger values of stage and UVM velocity.
   run.hcs=interp1(run.cs20.pstseccs,run.cs20.hcs,run.pstsecav);
   run.vcs=interp1(run.cs20.pstseccs,run.cs20.vcs,run.pstsecav);

   %%%%%% Trim arrays to minimum size necessary.
   run.vp=run.vp(1:run.WNmxmxmx,1:run.nensmax,1:ntran);
   run.xens=run.xens(:,1:run.nensmax,1:ntran);
   run.BTV=run.BTV(:,1:run.nensmax,1:ntran);
   run.WNmax=run.WNmax(1:run.nensmax,1:ntran);
   run.vflag=run.vflag(1:run.nensmax,1:ntran);
   run.dist=run.dist(1:run.nensmax,1:ntran);
   run.pstsec=run.pstsec(1:run.nensmax,1:ntran);
   run.depth=run.depth(1:run.WNmxmxmx,1:ntran);
   run.verav=run.verav(1:run.nensmax,1:ntran);
   run.verrms=run.verrms(1:run.nensmax,1:ntran);
   run.dbot=run.dbot(1:run.nensmax,1:ntran);
   run.madegdmx=max(run.madegood);
   %%%%%%  Add last ensemble and last transect to run structure.
   run.ens=ens;
   run.tran=tran;
   % Calculate ensemble and transect averages of vp.
   vaver2;