function n = EPIC_vector_average(in_name,out_name,instart,outstep,dirvar,spdvar)

% EPIC_vector_average:  time average variables in an EPIC 
%    netCDF file, including a current vector expressed 
%    as (direction, speed).
%
%    in_name is the name of an epic netcdf file.
%    out_name is the name of an epic netcdf file which will be
%        written, and will contain all the variables and 
%        attributes of in_name.
%    instart is the start time of the in_name data to be used, a 
%        Matlab datevec, [year,month,day,hour,minute,second].
%    outstep is the time step of out_name, in seconds.
%    n is the length of new variables created.
%    dirvar is the name of a current direction variable.
%    spdvar is the name of a current speed variable.
%
%         Fran Hotchkiss,  09 Mar 1999.

%Investigate in_name time base.
nc = netcdf(in_name);
ncatt = att(nc);
time = nc{'time'}(:);
time2 = nc{'time2'}(:);
tsecs = time*86400 + time2/1000;
tdifs = diff(tsecs);
[theMode, theCount]=modecount(tdifs)%;
incount = length(time);

%Create a netcdf variable with a regular time base which 
%  is close to the inc interval but a simple fraction of 
%  the out_name interval.
avcount = round(outstep/mean(theMode))%;
regstep = outstep/avcount/86400%;
regfirst = datenum(instart(1),instart(2),instart(3),instart(4),instart(5),instart(6))%;
secsfirst = regfirst*86400;
secsmax = 86400*ep_datenum([time(incount) time2(incount)]);
regcount = avcount * fix((secsmax-secsfirst)/outstep);
reglast = regfirst + regcount*regstep - regstep;
regtime = [regfirst:regstep:reglast].';
regeptime = EP_Time(regtime);

regc = netcdf('scratch.nc','clobber');
regc('time') = 0;
regc{'time'} = nclong('time');
regc{'time'}.units = ncchar('True Julian Day');
regc{'time'}.type = ncchar('EVEN');
regc{'time'}.epic_code = nclong(624);
regc{'time2'} = nclong('time') ;
regc{'time2'}.units = ncchar('msec since 0:00 GMT');
regc{'time2'}.type = ncchar('EVEN');
regc{'time2'}.epic_code = nclong(624);
regc{'count'} = nclong('time') ;
endef(regc)
nreg = length(regtime);
regc{'time'}(1:nreg) = regeptime (:,1);
regc{'time2'}(1:nreg) = regeptime (:,2);
regc{'count'}(1:nreg) = [1:nreg];
n = nreg/avcount;

%Copy and update inc attributes for out_name.
outc = netcdf(out_name,'noclobber');
ncatt = att(nc);
for i = [1:length(ncatt)]
   copy (ncatt{i}, outc);
end
outc.CREATION_DATE = ncchar(datestr(now,0));
history = ['Averaged, input file is ' in_name ':  ' nc.history(:)];
outc.history = ncchar(history);
outc.DELTA_T = ncchar(num2str(outstep));

% Create out_name variables and their attributes.
% Dimension definitions and attributes.
ncdim = dim(nc);
for i = [1:length(ncdim)]
   copy (ncdim{i}, outc);
end
% Variable definitions and attributes.
ncvar = var(nc);
for i = [1:length(ncvar)]
   ivarname = name(ncvar{i})
   copy (ncvar{i}, outc, 0, 1);
end
endef(outc)
%Cycle through in_name variables.
% Coordinate variables. 
nccoord = coord(nc);
for i = [1:length(nccoord)]
   ivarname = name(nccoord{i})
   if ~(strncmp('time',ivarname,4)) %exclude time and time2
      copy (nccoord{i}, outc, 1, 0);
   end
end

% Time variables.
time = regc{'time'}(:);
time2 = regc{'time2'}(:);
invar = ep_datenum([time time2]);
varray = reshape(invar,avcount,n);
outvar = mean((varray)).';
eptime = EP_Time(outvar);
outc{'time'}(1:n) = eptime (:,1);
outc{'time2'}(1:n) = eptime (:,2);
outc.start_time = datestr(outvar(1));
outc.stop_time = datestr(outvar(n));

% Scalar average record variables.
outcrecvar = recvar(nc);
for i = [1:length(outcrecvar)]
   ivarname = name(outcrecvar{i})
   if ~(strncmp('time',ivarname,4)) %exclude time and time2
      % Use timterp to make in_name variables have regular time base.
      invar = timterp(nc{ivarname},regc{'count'});
      %Reshape and average.
      varray = reshape(invar,avcount,n);
      outvar = mean((varray)).';
      %Update min/max attributes.
      outc{ivarname}.minimum = ncfloat(min(outvar));
      outc{ivarname}.maximum = ncfloat(max(outvar));
      %Put data in out_name.
      outc{ivarname}(1:n) = outvar(1:n);
   end
end

% Vector average speed and direction      
spd = timterp(nc{spdvar},regc{'count'});
dir = timterp(nc{dirvar},regc{'count'});
[u,v] = polar2uv(dir,spd);
%Reshape and average.
varray = reshape(v,avcount,n);
outv = mean((varray)).';
varray = reshape(u,avcount,n);
outu = mean((varray)).';
[dir,spd] = uv2polar(outu,outv);
%Update min/max attributes.
outc{dirvar}.minimum = ncfloat(min(dir));
outc{dirvar}.maximum = ncfloat(max(dir));
outc{spdvar}.minimum = ncfloat(min(spd));
outc{spdvar}.maximum = ncfloat(max(spd));
%Put data in out_name.
outc{dirvar}(1:n) = dir(1:n);
outc{spdvar}(1:n) = spd(1:n);

ncclose