function Bogenbruecke

% http://geometrie.hs-nb.de/MATLAB-Programme/Bogenbruecke.m   Bearbeitungsstand: 4.2.2004
% (C) Martin Nitschke, nitschke@hs-nb.de

close all
clc
Dat= [ 0.0  7.5 15.0 22.5 30.0 37.5 45.0 52.5 60.0 67.5 75.0 80.0 90.0
      69.3 68.7 66.8 63.6 59.3 53.6 46.8 38.8 29.3 18.8  7.6  0.0  NaN
      73.3 73.0 71.7 69.5 66.4 62.2 57.0 50.7 43.6 35.2 26.0  NaN  7.1]; % Originaldaten Muengstener Bruecke

x0=[-Dat(1,end:-1:2) Dat(1,:)]+Dat(1,end);
y0=Dat(2:3,[end:-1:2 1:end]);

I=4:length(x0)-3;
xI=x0(I);
xi=union(xI,linspace(xI(1),xI(end),100));
P=[polyfit(xI,y0(1,I),2)
   polyfit(xI,y0(2,I),2)];
yI=[polyval(P(1,:),xI)
    polyval(P(2,:),xI)];

dP=[polyder(P(1,:))
    polyder(P(2,:))];

[xLU,yLU]=Gerade(x0([2 3 4]),[y0(1,[2 3]) yI(1,1)],x0(4),yI(1,1),polyval(dP(1,:),x0(4)));
[xRU,yRU]=Gerade(x0([end-3 end-2 end-1]),[yI(1,end) y0(1,[end-2 end-1])],x0(end-3),yI(1,end),polyval(dP(1,:),x0(end-3)));
xLO=x0([1 3 4]);
gLO=polyfit(xLO,[y0(2,[1 3]) yI(2,1)],1);
yLO=polyval(gLO,xLO);
xRO=x0([end-3 end-2 end]);
gRO=polyfit(xRO,[yI(2,end) y0(2,[end-2 end])],1);
yRO=polyval(gRO,xRO);
alpha=atan(gLO(1))*180/pi;
beta=atan(polyval(dP(2,:),xLO(end)))*180/pi;
delta=beta-alpha;

set(plot([xI;xI],yI,'.k',[xi;xi]',[polyval(P(1,:),xi);polyval(P(2,:),xi)]','k',...
    xLU(2:3),yLU(2:3),'.k',xLU,yLU,'k',xLO(2:3),yLO(2:3),'.k',xLO,yLO,'k',...
    xRU(1:2),yRU(1:2),'.k',xRU,yRU,'k',xRO(1:2),yRO(1:2),'.k',xRO,yRO,'k',...
    [xLU(end-1) xI xRU(2);xLO(end-1) xI xRO(2)],[yLU(end-1) yI(1,:) yRU(2);yLO(end-1) yI(2,:) yRO(2)],'k',...
    [xLU(end-1) xI(1:ceil(end/2)-1);xI(1:ceil(end/2))],[yLO(end-1) yI(2,1:ceil(end/2)-1);yI(1,1:ceil(end/2))],'k',...
    [xI(ceil(end/2):end);xI(ceil(end/2)+1:end) xRU(2)],[yI(1,ceil(end/2):end);yI(2,ceil(end/2)+1:end) yRO(2)],'k',...
    x0([1 end]),max(yI(2,:))*[1 1],'k'),'MarkerSize',15);
axis('equal')
axis([x0([1 end]) 0 1.1*max(y0(2,:))]);
axis('off')
print('-deps','Bogenbruecke')

function [xG,yG]=Gerade(x,y,x0,y0,m);
xG=x;
yG=y0+m*(xG-x0);
