heksapolin_analyysi.m

% Analysoidaan heksapolimagneettikenttää. Lasketaan keskimääräinen
% kentän intensiteetti rlim-säteisen pinta-alan sisällä.
%
% Tehdään sama analyysi käyttäen kahta for-silmukkaa ja toisella
% tavalla hyödyntäen matriiseja. Näillä kahdella eri tavalla on huomattava
% nopeusero Octavessa ja jonkinlainen nopeusero Matlabissa. Mitataan
% nopeus tic/toc -kellolla.

function heksapolin_analyysi()
    fprintf( 'Analyysi käyttäen silmukoita\n' )
    tic
    analyysi_silmukka()
    toc

    fprintf( '\nAnalyysi ilman silmukoita, käyttäen matriiseja\n' )
    tic
    analyysi_matriisi()
    toc
end


function analyysi_silmukka()
    rlim = 39;
    Bmax = 1.07;
    A = -2*Bmax/rlim^2;

    N = 200;
    
    Bsum = 0;
    BN = 0;
    
    % Lasketaan käyttäen kahta sisäkkäistä for-silmukkaa. Indeksi i
    % viittaa x-koordinaatteihin ja j y-koordinaatteihin
    for i=1:N
        
        % Lasketaan x-koordinaatti indeksistä i
        x = -rlim + 2*rlim*(i-1)/(N-1);
        
        for j=1:N
            
            % Lasketaan y-koordinaatti indeksistä j
            y = -rlim + 2*rlim*(j-1)/(N-1);
        
            % Lasketaan r-koordinaatti ja testataan ollaanko ulkona
            % Hypätään seuraavaan pisteeseen jos ollaan
            r = sqrt(x^2 + y^2);
            if r > rlim
                continue
            end
            
            % Lasketaan magneettikentän intensiteetti
            Bx = A*x*y;
            By = 0.5*A*(x^2-y^2);
            Bint = sqrt(Bx^2+By^2);
            
            % Kontribuutio summaan
            Bsum = Bsum + Bint;
            BN = BN + 1;
            
        end
    end
    
    % Lasketaan keskimääräinen kentän intensiteetti
    Bave = Bsum/BN
end



function analyysi_matriisi()
    rlim = 39;
    Bmax = 1.07;
    A = -2*Bmax/rlim^2;

    % Tuotetaan x- ja y-koordinaatit joissa lasketaan
    N = 200;
    [x, y] = meshgrid(linspace(-rlim,rlim,N));
            
    % Lasketaan magneettikentän intensiteetti näissä pisteissä
    Bx = A*x.*y;
    By = 0.5*A*(x.^2-y.^2);
    Bint = sqrt(Bx.^2+By.^2);

    % Lasketaan totuusarvotaulukko (tai maski) missä kaikki pisteet,
    % joissa säde on pienempi kuin rlim on merkitty arvolla tosi ja
    % muut pisteet ovat epätosi
    rmaski = sqrt(x.^2 + y.^2) < rlim;

    % Poimitaan totuusarvomaskia käyttäen magneettikentän
    % intensiteettitaulukosta pisteet jotka ovat pinta-alan sisällä
    Bsisalla = Bint(rmaski);

    % Lasketaan keskimääräinen kentän intensiteetti
    Bave = sum(Bsisalla)/length(Bsisalla)
end