%% Ziege-Weide Aufgabe
% <http://www.bigbandi.de/dokus/ziege/ Problemstellung und Loseung> (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;

%% Lade das symbolische Paket fuer Octave
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if (isOctave)
   setenv PYTHON python
   pkg load symbolic
end

% 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);
X=X(2)       % Hack for octave to get the positive solution

% Ausgabe
disp('Schnittpunkt der beiden Kurven')
%pretty(yR); pretty(yL)
yR
yL
disp('ist')
%pretty(X)
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)
A
%%  Aufstellen der Gleichung, dass 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])     % Matlab
fplot(F,[0,2*R])        % Octave
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 = double(subs(subs(yL,l,L),x,xxL));
%yyL = double(subs(subs(yL,l,vpa(L)),x,vpa(xxL))); % avoids warnings but takes forever
% Wiesenkurve
xxR = linspace(-R, R, n);
yyR = double(subs(yR,xxR));
%yyR = double(subs(yR,vpa(xxR))); % avoids warnings but takes forever

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');





