% Define topography
scl=0.1;
xpft=[-4 -3 -2 -1 0 +1 +2 +3 +4];
zpft=scl*[2 -1 -1 -1 0 -1 -1 -0.4 2];

PF=polyfit(xpft,zpft,4);

xmn=min(xpft); xmx=max(xpft);
xi=0.1; xex=4;
x=xmn:xi:xmx;
xfull=[xmn-xex x xmx+xex];
zfull=polyval(PF,xfull);
z=zfull(2:end-1);
zrn=max(z)-min(z);
rb=0.02*zrn; rb2=rb*2; rb10=rb*10;
zbot=min(z)-0.05*zrn; ztop=max(z)+rb10;

% Slope and acceleration
g=9.8;
s=gradient(zfull,xfull);
ang=atan(s);
sang=sin(ang); cang=cos(ang);
a=-g*sang.*cang;
afm=cang.^2;

rbm=1.5/scl;
rbv=rb./cos(atan(rbm*s));

% Divisions
sfir=s(1:end-1); ssec=s(2:end);
ieq2=find((ssec-sfir)<0 & (ssec.*sfir)<0);
xeq2=mean(x(ieq2:ieq2+1));

x1=[x(1:ieq2) x(ieq2) x(1)];
z1=[z(1:ieq2) zbot zbot];
x2=[x(ieq2:end) x(end) x(ieq2)];
z2=[z(ieq2:end) zbot zbot];


% Plot topography
gcfp=[100 100 960 600];
fgn=figure(1); clf, set(gcf,'position',gcfp)
ax1=axes('position',[0.1 0.02 0.7 0.96],'color',[1 1 1]); box on, hold on
set(ax1,'xtick',[],'ytick',[])

ax=[x(1) x(end) zbot ztop];
axis(ax)

coleq(1,:)=0.6*ones(1,3);
coleq(2,:)=0.3*ones(1,3);

fill(x1,z1,coleq(1,:),'linestyle','none')
fill(x2,z2,coleq(2,:),'linestyle','none')
%plot(x,0.2*a)