function [VNMO2D VNMO2DS]=vrmsmodel
% 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);

% NMO or Stacking velocity
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];

% 
VNMO=zeros(nfunc,nt);
for icdp=1:nfunc
    count=1;
    twt =eval(['tnmo_',num2str(icdp)]);
    vnmo=eval(['vnmo_',num2str(icdp)]);
    [dummy ndata]=size(twt);
    %vertical interpolation
    for itwt=1:nt
        if count < ndata
            if  itwt*dt > twt(count+1)
                count=count+1;
            end
            VNMO(icdp,itwt)=vnmo(count);
        else
            VNMO(icdp,itwt)=vnmo(count);
        end
    end
end

% Horizontal Interporation
VNMO2D=zeros(nx,nt);
for itwt=1:nt
    count=1;
    for ix=1:nx
        icdp=ix+cdps;
        if icdp < cdp(1)
            VNMO2D(ix,itwt)=VNMO(1,itwt);
        elseif icdp > cdp(nfunc)
            VNMO2D(ix,itwt)=VNMO(nfunc,itwt);
        else
            if  icdp > cdp(count+1)
                count=count+1;
            end
            VNMO2D(ix,itwt)=(icdp-cdp(count))*VNMO(count+1,itwt)+(cdp(count+1)-icdp)*VNMO(count,itwt);
            VNMO2D(ix,itwt)=VNMO2D(ix,itwt)/(cdp(count+1)-cdp(count));
        end
    end
end



% Smoothing
VNMO2DS=VNMO2D;
ap=100;
for it=1:nt
    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;
        VNMO2DS(ix,it)=mean(VNMO2D(istart:iend,it));
    end
end
