% Lösning till ENM 3.14 -- rita upp kurvorna
%
% Skriv om ekvationen i polära koordinater, låt vinkeln theta gå
% från 0 till 2*pi och lös ut avståndet r för varje theta med
% Newtons metod (för en skalär ekvation).
%
% Joar Bagge 2016-11-24

clear

nth = 500; % antal theta-punkter vi vill ha längs kurvan
Theta = linspace(0, 2*pi, nth); % vektor med alla theta

% Funktionen som implicit definierar r som funktion av theta
% och dess derivata med avseende på r.
f = @(r,th,a) r^2 - 1 - a*sin(r^2*sin(2*th)/2);
dfdr = @(r,th,a) 2*r - a*cos(r^2*sin(2*th)/2)*r*sin(2*th);

clf
for a=[0.8, 1.6, 2.4, 3.2]
  R = []; % vektor med alla r
  r = 1; % startgissning (exakt för theta=0)
  for j=1:nth
    th = Theta(j);
    % Newtons metod
    dr = Inf;
    while abs(dr) >= 1e-10
      dr = -f(r,th,a) / dfdr(r,th,a);
      r = r + dr;
    end
    R = [R r];
  end
  % Polära koordinater (transformera tillbaka)
  x = R.*cos(Theta);
  y = R.*sin(Theta);
  plot(x, y, 'LineWidth', 1.5), hold on
end
axis equal
