%% Ziege-Weide Aufgabe
% <https://de.wikipedia.org/wiki/Ziegenproblem_(Geometrie)> (meine Loesung schaut aehnlich aus)
% 
% Das symbolischen Paket von Matlab wird benutzt so weit es geht und zum 
% Schluss muss numerisch geloest werden.

%%  Variablen
% symbolisch (gesuchte Groessen)
clc; clear all; %clf;

% syms l real positive;
l = sym('l','positive');           % Leinenlaenge
% syms x real positive;        
x = sym('x','positive');          % x-wert des Schnittpunktes Wiesenkreis mit Leinenkreis
% gegebene Groessen, hier als Zahlenwert
R = 10;                      % Wiesenradius
assumeAlso(l<2*R);

%% Die beiden Kurven (symbolisch)
yR = R - sqrt(R^2-x^2);
yL = sqrt(l^2-x^2);

%% bestimme den Schnittpunkt
X = solve(yL - yR, x);

% Ausgabe
disp('Schnittpunkt der beiden Kurven')
pretty(yR); pretty(yL)
disp('ist')
pretty(X)

%% unbestimmtes Integrieren 
I = int(yL - yR, x);

%% Einsetzen der Integrationsgrenzen
%  ==> abgegraste Flaeche, symbolische Funktion A(l)  
A = 2*int(yL - yR, x, 0,X); 

% Ausgabe
disp('abgegrasten Flaeche, Funktion A(l)')
pretty(A)
%%  Aufstellen der Gleichung, sodass die halbe Flaeche abgegrast wurde
% Hier <http://de.mathworks.com/help/symbolic/matlabfunction.html konvertieren>
%       wir eine symbolische Funktion in eine numerische Funktion.

Fs= A - pi*R^2/2;          % symbolische Funktion
F = matlabFunction(Fs);    % numerische Funktion
                           
% Grafik
% xx=linspace(0,2*R,1001);
% plot(xx,F(xx),[xx(1) xx(end)],[0 0 ],'b')
% title('A(L)-pi*R^2/2'); xlabel('L')

fplot(Fs,[0,2*R])
hold on 
plot([0,2*R],[0 0 ],'b')
hold off
stringFs=strjoin(arrayfun(@char, Fs, 'uniform', 0));
title('Area = A(L)-pi*R^2/2'); xlabel('L')

%% numerisches Bestimmen der Nullstelle (= Loesung)
%  Hier wird nicht exakt gelöst !!
options = optimset('Display','iter');
L = fzero(F, R, options)                     % only with optimization package
%options = optimoptions(@fsolve,'Display','iter') % Wie geht das genau?
% options = optimoptions(optimoptions(fsolve),'Display','iter')
% L = fsolve(F,R)


%% Grafik 
n  = 100;
% Leinenkurve
xxL = linspace(-L, L, n);
yyL = subs(subs(yL,l,L),x,xxL);
% Wiesenkurve
xxR = linspace(-R, R, n);
yyR = subs(yR,xxR);

plot(xxL,yyL,'r');
hold on;
plot(xxR,yyR,'b');
plot(xxR,2*R-yyR);
plot(xxL,R);
axis equal
title(['Weidenradius ',num2str(R),' m und Leinenlaenge ',num2str(L),' m']);
legend('Leine','Weide');

% publish('ziege_wiese.m');





