function interp_eng(engin,engout,mode,dt,flagplot)
% interp_eng Temporal interpolation of CMS-Wave spectra
%
% Usage:
%  interp_eng(engin,engout,mode,dt)
%
% Description:
%   Interpolates the energy spectra in the CMS-Wave *.eng file to a
%   specified output interval which may be specified in intervals or
%   time units (hours). For example, an output interval of 0.5 index
%   units, interpolates the spectra at half intervals. This subroutine
%   is useful for filling small gaps in CMS.
%
% Input:
%   engin = Input eng file
%   engout = Output eng file
%   mode = 1 - Index, 2 - Time
%   dtout = Output steering interval in hours for mode=2
%         or index for mode 1 (can be less than 1)
%
% written by Alex Sanchez, USACE-CHL


%------ For Testing ------------
if nargin==0
    fclose all;
    engin = 'Wave_AugSep05.eng';
    engout = 'Wave_AugSep05_30min.eng';
    mode = 2;
    dt = 0.5;
    interp_eng(engin,engout,mode,dt,1)
    return
end
%------------------------------

if nargin<5
    flagplot = 0;
end

fid = fopen(engin,'r');
fid2 = fopen(engout,'w');

%Formats
fmt21 = '%8.5f';
fmt22 = '%11.5f'; %format for spectra
fmt23 = '%d %8.5f %8.5f %8.5f %8.5f'; %header format

%Read spectra size
nf = fscanf(fid,'%d',1);   %Number of frequencies
ndir = fscanf(fid,'%d',1); %Number of directional bins
nfdir = nf*ndir;           %Size of matrix
dir = (-(ndir-1)/2:(ndir-1)/2)*180/(ndir+1);

%Frequencies
f = fscanf(fid,'%f',nf);

%Read first header and spectra
idate1 = fscanf(fid,'%d',1);
disp(['Reading spectra ',int2str(idate1)])
ws1 = fscanf(fid,'%f',1); %wind speed, m/s
wd1 = fscanf(fid,'%f',1); %wind direction, deg
fp1 = fscanf(fid,'%f',1); %peak frequency, 1/s
tide1 = fscanf(fid,'%f',1); %tide level, m
E1 = fscanf(fid,'%f',nfdir); %Energy spectra

if flagplot==1
    figure
    EE = reshape(E1,ndir,nf);
    him = imagesc(f,dir,EE);
    set(gca,'ydir','normal')
    colorbar
    xlabel('Frequency, 1/s')
    ylabel('Direction, deg')
end

%Timing
if mode==1 %Index
    t1 = idate1;
    t = t1 + dt;   %Next index
    idate = idate1 + 1; %used for output only
else      %Time
    [yyyy,mm,dd,hh] = index1calender(idate1);
    t1 = datenum(yyyy,mm,dd,hh,0,0);
    t = t1 + dt/24; %units are in days, Next time
    adate = datestr(t,'yymmddhh'); %used for output only
    idate = str2double(adate);
end

%Write to output file
fprintf(fid2,'%d %d',nf,ndir);
fprintf(fid2,'%s',char(10));
ncol = 10;
for i=1:ceil(nf/ncol)
    istart = 1 + (i-1)*ncol;
    iend = min(istart+ncol-1,nf);
    fprintf(fid2,fmt21,f(istart:iend));
    fprintf(fid2,'%s',char(10));
end
fprintf(fid2,fmt23,idate1,ws1,wd1,fp1,tide1);
fprintf(fid2,'%s',char(10));
ncol = 7;
for i=1:ceil(nfdir/ncol)
    istart = 1 + (i-1)*ncol;
    iend = min(istart+ncol-1,nfdir);
    fprintf(fid2,fmt22,E1(istart:iend));
    fprintf(fid2,'%s',char(10));
end

while 1
    idate2 = fscanf(fid,'%d',1);
    if isempty(idate2), break, end
    disp(['Reading spectra ',int2str(idate2)])
    ws2 = fscanf(fid,'%f',1);
    wd2 = fscanf(fid,'%f',1);
    fp2 = fscanf(fid,'%f',1);
    tide2 = fscanf(fid,'%f',1);
    if mode==1 %Index
        t2 = idate2;
    else       %Time
        [yyyy,mm,dd,hh] = index1calender(idate2);
        t2 = datenum(yyyy,mm,dd,hh,0,0);
    end
    E2 = fscanf(fid,'%f',nfdir);
    if flagplot==1
        EE = reshape(E2,ndir,nf);
        set(him,'CData',EE)
        drawnow
        title(num2str(idate))
    end
    while t<=t2
        fac = (t-t1)/(t2-t1);
        E = (1.0-fac)*E1 + fac*E2;
        ws = (1.0-fac)*ws1 + fac*ws2;
        wd = (1.0-fac)*wd1 + fac*wd2;
        fp = (1.0-fac)*fp1 + fac*fp2;
        tide = (1.0-fac)*tide1 + fac*tide2;
        fprintf(fid2,fmt23,idate,ws,wd,fp,tide);
        fprintf(fid2,'%s',char(10));
        for i=1:ceil(nfdir/ncol)
            istart = 1 + (i-1)*ncol;
            iend = min(istart+ncol-1,nfdir);
            fprintf(fid2,fmt22,E(istart:iend));
            fprintf(fid2,'%s',char(10));
        end
        if mode==1 %Index
            t = t + dt;
            idate = idate + 1; %used for output only
        else       %Time
            t = t + dt/24;
            adate = datestr(t,'yymmddhh'); %used for output only
            idate = str2double(adate);
        end
    end %while
end %while

fclose(fid);
fclose(fid2);
return
%
% 30 35
%  0.04000 0.05000 0.06000 0.07000 0.08000 0.09000 0.10000 0.11000 0.12000 0.13000
%  0.14000 0.15000 0.16000 0.17000 0.18000 0.19000 0.20000 0.21000 0.22000 0.23000
%  0.24000 0.25000 0.26000 0.27000 0.28000 0.29000 0.30000 0.31000 0.32000 0.33000
%    1  0.000000  0.000000  0.190000  0.000000
%     0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000
%     0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000


function [yyyy,mm,dd,hh] = index1calender(ind)
% Converts the CMS-Wave index to a calender date
%

str = int2str(ind);

if(ind<99999) %mddhh
    yyyy = 2000;
    mm = str2double(str(1));
    dd = str2double(str(2:3));
    hh = str2double(str(4:5));
elseif(ind<999999) %mmddhh
    yyyy = 2000;
    mm = str2double(str(1:2));
    dd = str2double(str(3:4));
    hh = str2double(str(5:6));
elseif(ind<9999999)
    yyyy = 2000+str2double(str(1));
    mm = str2double(str(2:3));
    dd = str2double(str(4:5));
    hh = str2double(str(6:7));
elseif(ind<99999999)
    yyyy = 2000+str2double(str(1:2));
    if(yyyy>2030), yyyy = yyyy - 100; end
    mm = str2double(str(3:4));
    dd = str2double(str(5:6));
    hh = str2double(str(7:8));
else %(ind>99999999)
    yyyy = str2double(str(1:4));
    mm = str2double(str(5:6));
    dd = str2double(str(7:8));
    hh = str2double(str(9:10));
end


return
