% 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)