%% Övning 1
% Table of Contents
% Uppgift 1: Find solutions of  using Newton's method.
% Uppgift 2: Solve for  in
% Uppgift 3: Heat conduction. A long bar is momentarily heated at . How long is above 20C?
% Uppgift 4: Compute the integral
% Uppgift 5: Draw the Fresnel integrals
%
%% Uppgift 1: Find solutions of f\left(x\right)=\tan \left(x\right)-x=0 using Newton's method.
% The nature of the function means that initial x_0 must be chosen very carefully, since the roots are near the asymptotes.
% Below, for instance, x_0 =2\ldotp 5 finds something around 10.9 instead of 4.3!
% This can be mitigated by using a "friendlier" function with the same solution set that is continuous.
% We can rearrange f(x) to obtain such a function g(x).
% \tan \left(x\right)-x=0\to \frac{\textrm{sinx}}{\textrm{cosx}}-x=0\to \textrm{sinx}-\textrm{xcosx}=0
% Thus g\left(x\right)=\textrm{sinx}-\textrm{xcosx}=0

clear all
close all
f=@(x) tan(x)-x;
fprim = @(x) sec(x).^2-1; % or tan(x).^2
fprim = @(x) tan(x).^2;
x=linspace(0,12,500);
plot(x, f(x)), hold on, grid on, axis([0 12 -10 10])
plot(x, fprim(x), 'r--')
title('f(x)=tan(x)-x')
legend('f(x)', 'f''(x)')

kmax = 10; 
tol  = 0.5e-4; % tol = 'tolerance'
x0   = 2.5;    % Bad guess!
 
for k=1:kmax
    d=f(x0)/fprim(x0);
    x0=x0-d;
    if abs(d)<tol; break; end
end

%
% Man ser tydligt ovan att vissa nollställen hamnar nära kurvans asymptoter.
% För att lättare tillämpa Newtons metod så kan man hitta en funktion med samma lösningsmängd: tanx - x = 0 => sinx/cosx - x = 0 => sinx - xcosx = 0

close all

% The analytical functions.
f     = @(x) sin(x) - x.*cos(x);
fprim = @(x) x.*sin(x);

plot(x, f(x), x, fprim(x), '--')
hold on
title('f(x)=sin(x)-x cos(x) Roots shown')
xlabel('x')
grid on

x0=4;
for k=1:kmax
    d=f(x0)/fprim(x0);
    x0=x0-d;
    if abs(d)<tol; break; end
end
k,x0, plot(x0, f(x0), 'ro')

x0=8;
for k=1:kmax
    d=f(x0)/fprim(x0);
    x0=x0-d;
    if abs(d)<tol; break; end
end
k,x0, plot(x0, f(x0), 'ro')

x0=10;
for k=1:kmax
    d=f(x0)/fprim(x0);
    x0=x0-d;
    if abs(d)<tol; break; end
end
k,x0, plot(x0, f(x0), 'ro')
legend('f(x)', 'f''(x)')

%
%% Uppgift 2: Solve for r in L=L_0 -20\log_{10} \left(r\right)-\beta r=20\textrm{dB}
% Sound volume [dB] r meters from a source is given by L(r), where L_0 is the volume 1m from the source and \beta is an attenuation factor.

close all

L0   = 80;      % dB 1m from source.
beta = 1.15e-3; % attenuation factor
L    = @(r) L0 - 20.*log10(r) - beta*r;
r    = linspace(1,1000, 10000); % Start from 1m due to L0 being defined at 1m.

plot(r, L(r)), grid on, hold on
xlabel('r'); ylabel('L[dB]')

% Since fzero() finds roots, and we are trying to find a solution at y=20, we must offset the function by that amount.

res = fzero(@(r)L(r)-20, 1000) % In order to insert the 20dB offset, the function variables must be explicitly called.
plot(res, L(res), 'ro')

%
%% Uppgift 3: Heat conduction. A long bar is momentarily heated at x=0. How long is x=1above 20C?

close all
k = 5.17;
C = 28.3;

% T\left(x,t\right)=\frac{C}{\sqrt{t}}\exp \left(-\frac{x^2 }{\textrm{kt}}\right),t>0

T=@(x,t) C./sqrt(t).*exp(-(x.^2)./(k.*t));
t=linspace(0.1,2, 1000);

% If we plot the resulting temperature at x=1, we get the following curve.

plot(t, T(1,t)), hold on, grid on
xlabel('t'); ylabel('T'); title('T at x=1')

% We use the same trick as before to find the interval T>20.

x1 = fzero(@(t)T(1,t)-20, 0.1);
x2 = fzero(@(t)T(1,t)-20, 1.8);

duration = x2 - x1;
plot([x1 x2], [20 20])

% Just for fun, we can animate the temperature of the bar over time (t>0).

xvec = (-1:0.1:1);
Tvec = T(xvec,1);
figure
for t = t
    Tvec = T(xvec,t);
    plot(xvec, Tvec)
    axis([-1 1 0 100])
    xlabel('x'); ylabel('T')
    pause(0.01)
    title(['t= ' num2str(t)])
end

%
%% Uppgift 4: Compute the integral \int_0^1 \frac{\textrm{sinx}}{\sqrt{1-x^2 }}\textrm{dx}
% The integrand contains singularities. Find a function that has the same solution and compute the integral of that instead.
% Through variable substitution, we find an equivalent expression \sin \left(\sin \left(x\right)\right)

% Integralen kan uttryckas sin(sin(t)) genom variabelsubstitution
% (x=sin(t))

f=@(t)sin(sin(t));
integral(f, 0, pi/2)

%
%% Uppgift 5: Draw the Fresnel integrals 0\le x\le 5
% To draw the integral, really we need the primitive function F.
% The first fundamental theorem of calculus states that a primitive function may be obtained as the integral of f with variable bounds of integration.
% F\left(x\right)=\int_a^x f\left(t\right)\mathrm{dt}
% The primitive function may change by a constant depending on the choice of a, but in all cases the second fundamental theorem of calculus states that  F^{\prime } \left(x\right)=f\left(x\right) for all of the infinitely many primitive functions.
% Thus, the Fresnel Integrals are specific primitive functions.
%
% C\left(x\right)=\int_0^x \cos \left(\pi \frac{t^2 }{2}\right)\mathrm{dt}, S\left(x\right)=same for sine.

close all
a=-5;
b=5;    
n=400;

x=linspace(a,b,n);
C=zeros(size(x));
S=zeros(size(x));

fresnelIntegrandC = @(t) cos(pi*t.^2./2);
fresnelIntegrandS = @(t) sin(pi*t.^2./2);

C=zeros(size(x));
S=zeros(size(x));

% More efficient as it computes only each interval
for k=1:n
    C(k) =  integral(fresnelIntegrandC, 0, x(k));
    S(k) =  integral(fresnelIntegrandS, 0, x(k)); 
end

subplot(4,2,1:2)
plot(x, C, x, S)
hold on
title('Fresnel''s Integrals C(x)=\int_0^x cos(\pi t^2/2)')
legend('C(\chi)', 'S(\chi)')
xlabel('x')
grid on

subplot(4,2,3:4)
plot(x, cos(pi/2.*x.^2))
hold on
plot(x, sin(pi/2.*x.^2))
legend('C integrand', 'S integrand')
grid on
xlabel('x')
title('Integrands')

subplot(4,2,5:8)
plot(C,S), title('Cornu''s Spiral'), axis square
xlabel('C(x)'); ylabel('S(x)');
grid on