function [vave vaves]=vavemod(dz,nz)
% NMO velocity interpolation

nt=2750;dt=0.004;twtmax=(nt-1)*dt;
nx=401;dx=16.67;
cdps=900;

cdp=[933,958,983,1008,1033,1058,1083,1108,1133,1158,1183,1208,1233,1258,1283];

[dummy nfunc]=size(cdp);

tnmo_1=[0.0,5.52263,6.44307,6.97119,7.65775,8.50274,10.9849];
vnmo_1=[1500.,1522.81,1460.57,1622.4,1834.01,2331.93,2892.08];
tnmo_2=[0.0,5.53018,6.49588,7.55967,9.19685,10.9774];
vnmo_2=[1500.,1497.92,1497.92,1634.84,2966.77,3016.56];
tnmo_3=[0.0,5.52263,6.35254,7.40878,9.27229,10.9849];
vnmo_3=[1500.,1497.92,1460.57,1659.74,2307.03,2991.67];
tnmo_4=[0.0,5.53772,7.48422,8.93279,10.9925];
vnmo_4=[1500.,1497.92,1634.84,2282.14,3564.27];
tnmo_5=[0.0,5.55281,6.99383,7.5144,8.48011,10.9925];
vnmo_5=[1500.,1522.81,1659.74,1734.43,2319.48,2991.67];
tnmo_6=[0.0,5.53018,6.69959,7.46159,9.45336,10.9774];
vnmo_6=[1500.,1485.47,1572.6,1871.35,2282.14,2929.43];
tnmo_7=[0.0,5.53018,6.79012,7.58985,8.88752,10.9774];
vnmo_7=[1500.,1497.92,1622.4,1771.77,2406.61,2979.22];
tnmo_8=[0.0,5.53772,6.29973,7.3786,8.7668,10.9698];
vnmo_8=[1500.,1497.92,1497.92,1759.32,2842.29,3589.17];
tnmo_9=[0.0,5.53772,7.30315,9.16667,10.9925];
vnmo_9=[1500.,1510.36,1672.19,2630.68,3053.91];
tnmo_10=[0.0,5.52263,7.39369,9.60425,10.9849];
vnmo_10=[1500.,1497.92,1721.98,2443.96,3004.11];
tnmo_11=[0.0,5.53018,6.93347,7.55213,10.9849];
vnmo_11=[1500.,1535.26,1634.84,1734.43,2867.19];
tnmo_12=[0.0,5.52263,7.43896,8.2915,10.9849];
vnmo_12=[1500.,1510.36,1647.29,2282.14,2979.22];
tnmo_13=[0.0,5.53018,7.49931,8.13306,10.9698];
vnmo_13=[1500.,1510.36,1659.74,1921.15,2555.99];
tnmo_14=[0.0,5.54527,7.35597,10.9849];
vnmo_14=[1500.,1510.36,1697.08,3004.11];
tnmo_15=[0.0,5.53772,7.29561,8.19342,9.08368,10.9925];
vnmo_15=[1500.,1522.81,1746.88,2456.41,2879.64,4248.91];
%

%
VINT=zeros(nfunc,nz);
for icdp=1:nfunc
    count=1;
    twt =eval(['tnmo_',num2str(icdp)]);
    vnmo=eval(['vnmo_',num2str(icdp)]);
    [dummy ndata]=size(twt);
    %
    vint=zeros(1,ndata);
    tkns=vint;
    vint(1)=vnmo(1);
    for count=1:ndata
        if count < ndata     
             % Dix equation
            vint(count)=sqrt((vnmo(count+1)^2*twt(count+1)-vnmo(count)^2*twt(count))/(twt(count+1)-twt(count)));
            tkns(count)=(twt(count+1)-twt(count))*vint(count)/2;
        else
            vint(ndata)=vint(ndata-1);
            tkns(ndata)=(nt*dt-twt(ndata))*vint(ndata-1)/2;
        end
    end

    %Vertical interpolation in depth domain
    depth=tkns(1);
    count=1;
    for iz=1:nz
        if count < ndata
            if  dz*iz > depth
                count=count+1;
                depth=depth+tkns(count);
            end
%             vint(count)=sqrt((vnmo(count+1)^2*twt(count+1)-vnmo(count)^2*twt(count))/(twt(count+1)-twt(count)));
            VINT(icdp,iz)=vint(count);
        else
            VINT(icdp,iz)=vint(ndata-1);
        end
    end
end

% Interpolation in x direction
VINT2D=zeros(nx,nz);
for iz=1:nz
    count=1;
    for ix=1:nx
        icdp=ix+cdps;
        if icdp < cdp(1)
            VINT2D(ix,iz)=VINT(1,iz);
        elseif icdp > cdp(nfunc)
            VINT2D(ix,iz)=VINT(nfunc,iz);
        else
            if  icdp > cdp(count+1)
                count=count+1;
            end
            %             TWT(icdp,itwt)=round(twt(count)/dt)+1;
            VINT2D(ix,iz)=(icdp-cdp(count))*VINT(count+1,iz)+(cdp(count+1)-icdp)*VINT(count,iz);
            VINT2D(ix,iz)=VINT2D(ix,iz)/(cdp(count+1)-cdp(count));
        end
    end
end

%Plot Average Velocity Model
figure;
imagesc([1:nx]*dx,dz*[1:nz],VINT2D',[1500 6000]);
title('Interval Velocity from Dix Inversion','FontSize',18);
xlabel('X (m)','FontSize',16);
ylabel('Depth (m)','FontSize',16);
colorbar;

% vave=VINT2D;
% for iz=1:nz
%     for ix=1:nx
%        vave(ix,iz)=mean(VINT2D(ix,1:iz));
%     end
% end
vave=VINT2D;
for ix=1:nx
    t=0;
    for iz=1:nz 
       delt=dz/VINT2D(ix,iz);
       t=t+delt;
       vave(ix,iz)=iz*dz/t;
    end
end

% Smoothing
VINT2DS=VINT2D;
ap=100;
for iz=1:nz
    for ix=1:nx
        istart=1 + (ix-1)-ap;iend=1+ (ix-1)+ap;
        if istart<1;istart=1;end;
        if iend>nx;iend=nx;end;
        VINT2DS(ix,iz)=mean(VINT2D(istart:iend,iz));
    end
end

vaves=VINT2DS;
for ix=1:nx
    t=0;
    for iz=1:nz
       delt=dz/VINT2DS(ix,iz);
       t=t+delt;
       vaves(ix,iz)=iz*dz/t;
    end
end





