% Testprogram for å løse karuselloppgaven (oppgave 7) i Coriolis-
% oppgavesettet i uke 15, når ADAPTIV 4. ordens Runge-Kutta 
% benyttes, og vi bruker 3D vektornotasjoner (nesten) over alt.
% (Det ville vært mest naturlig HER å bruke 2D-beregninger.)
% Programmet er skrevet av Arnt Inge Vistnes. Denne versjonen
% er fra 1. mai 2005.
% Selve rutinene for adaptiv RK (rk4a.m) og basis RK (rk4.m) er
% hentet fra Garcia's bok: Numerical Methods for Physics, 2.ed. 2000.
% Disse må finnes i samme "current directory" som dette hovedprog-
% rammet når dette kjøres. 
% I tillegg må vi samme sted ha et program for beregning av akselerasjon
% (og hastighet) for det systemet som behandles. Vi bruker her 
% rutinen corioliskraft3d.m. Det må være overensstemmelse mellom 
% hovedprogram og dette siste tilleggsprogrammet mhp om vi arbeider
% i 2D eller 3D, mens rk4a og rk4 tåler vilkårlige dimensjoner 
% uten omskriving.

clear all;
N = 100;                % Ca antall punkter vi antar å bruke i beregningene
                        % Brukes her bare for å dimensjonere allokering
                        % av hukommelse til arrayer for plotting o.l.
err = 1.0e-12;          % Feilmargin brukt ved adaptiv Runge-Kutta.
ts = 0.0;               % Velger å starte klokka (t=0) ved forsøkets start
tn = 1.0;               % Omtrentlig tidsrom vi ønsker å studere bevegelsen i.
                        % Brukes HER bare for å gi et første tips til
                        % delta_t for beregningene.
                        % Enheter må vi selv holde nøye rede med!!!
r = zeros(N,3);        % Allokerer plass i hukommelsen for
v = zeros(N,3);        % posisjon og hastighet,

t = zeros(N,1);        % Allokerer plass for faktisk brukt tid i bane

% Initialbetingelser (posisjon og hastighet i SI-enheter)
r1(1) =  0.0;
r1(2) = -3.0;
r1(3) =  0.0;
v1(1) =  0.0;
v1(2) =  6.0;
v1(3) =  0.0;

% Hele initialbetingelsene samles i én vektor og paramtre i en annen
status = [r1 v1]; 
omegaparam = 2*pi*12/60;  % Vinkelhastighet

%tid=linspace(ts,tn,N);   % Oppretter hele tidsarrayen i en jafs
                          % Blir overskrevet ved adaptiv metode
delta_t=(tn-ts)/(N-1);    % Lengde på tidsstegene Bare startverdi. 
                          % Endres underveis ved adaptiv metode.
n=0;
avst=norm(r);             % Avstand origo til punkt
while (avst<3.0001)       % Løkke for å beregne de resterende punktene              
    n = n+1;
    r(n,1) = status(1); r(n,2) = status(2); r(n,3) = status(3);
    v(n,1) = status(4); v(n,2) = status(5); v(n,3) = status(6);
    t(n)=ts;
    [status ts delta_t]=rk4a(status,ts,delta_t,err,'corioliskraft3d', ...
        omegaparam);
    avst = norm(r(n,:));
                        % 4. ordens adaptiv Runge-Kutta
                        % hvor likningene for kraften (akselerasjonen)
                        % ligger i funksjonen coriolsikraft3d.m. 
end;
N = n;                  % Virkelig antall beregnede punkter
t(n)                    % Utskrift av antall sek. før ball når randen

% Lager en sirkel for å visualisere resultatet bedre
for n=1:100
    phi=2*pi/100;
    rsirkel(n,1)=3.0*cos(phi*n);
    rsirkel(n,2)=3.0*sin(phi*n);
end;

% Diverse plot. Valgt 2D 
plot(r(:,1),r(:,2),'-r');
hold;
plot(rsirkel(:,1),rsirkel(:,2),'-b');
daspect([1 1 1]);       % Ønsker å ha lik skalering på alle aksene
xlabel('X');
ylabel('Y');
