function [ output_args ] = Calculation( input_args )
%UNTITLED1 Summary of this function goes here
%  Detailed explanation goes here

clear global M Gx Gy EV kx ky EVS k dbg_reci dbg_ndx1;
global M Gx Gy EV kx ky k EVS k dbg_reci dbg_ndx1;

a=100.0;
b=a*0.8;

na=3.5;
nb=1.0;

% determines number of reciprocal lattice vectors to use as
% Gx=(-Range:Range)*pi/a and Gy=(-Range:Range)*pi/a, giving a total of
% (Range*2+1)^2 vectors. This also determines the number of elements in the
% matrix.
Range=7;

% First we compute all the needed reciprocal lattice vectors. G = (q,p)pi/a
for p= -Range:Range
    index=p+Range+1;
    Gx(index)=p*pi/a;
    Gy(index)=p*pi/a;
end

% number of points to calculate, this determines the smothness of the curve
Resolution=30;


% Main loop, in this loop we setup the matrix 'Resolution' number of times,
% and compute eigenvalues for each loop.
mx=1;
my=1;
for i=1:Resolution
    i %write out status
    kx(i) = 0;
    ky(i) = ((i-1)/(Resolution-1))*(pi/(a));
    k(i)  = sqrt(kx(i)^2 + ky(i)^2);

    % Setup matrix
    for p=-Range:Range % Gx
        for q=-Range:Range % Gy, Outer loops, matrix rows (G loop)

            for m=-Range:Range % Gx'
                for n=-Range:Range % Gy', Inner loops, matrix columns (G' loop)
                    ndx_n=n+Range+1; ndx_m=m+Range+1;
                    ndx_q=q+Range+1; ndx_p=p+Range+1; % indices

                    g1x = Gx(ndx_p); g1y = Gy(ndx_q); % G
                    g2x = Gx(ndx_m); g2y = Gy(ndx_n); % G'

                    ggx = g1x - g2x; % X-component of (G-G')
                    ggy = g1y - g2y; % Y-component of (G-G')

                    % Dotp[(G1+k),(G2+k)] written out
                    g1kg2k = g1x*g2x + g1y*g2y + g1x*kx(i) + g1y*ky(i) + g2x*kx(i) + g2y*ky(i) + k(i)^2;

                    % Determine the coefficient Eta, when ggx or ggy (or
                    % both) is equal zero we have to take the limit of
                    % Sin(x*a)/x as x->0.
                    if(ggx ~= 0 && ggy ~= 0)
                        Eta = 1/(a^2*ggx*ggy)*((1/na^2)*sin(ggx*a)*sin(ggy*a) + (1/nb^2 - 1/na^2)*sin(ggx*b)*sin(ggy*b));
                    elseif(ggx == 0 && ggy ~= 0) % gxx->0
                        Eta = 1/(a*na^2*ggy)*sin(ggy*a) + b/(a^2*ggy)*(1/nb^2 - 1/na^2)*sin(ggy*b);
                    elseif(ggy == 0 && ggx ~= 0) % gyy->0
                        Eta = 1/(a*na^2*ggx)*sin(ggx*a) + b/(a^2*ggx)*(1/nb^2 - 1/na^2)*sin(ggx*b);
                    else % gxx->0 and gyy->0
                        Eta = 1/(na^2) + b^2/a^2*(1/nb^2 - 1/na^2);
                    end

                    M(my,mx) = Eta*g1kg2k;
                    mx=mx+1;
                end
            end
            mx=1;
            my=my+1;
        end
    end

    EV=eig(M);
    EV=sort(EV);
    EVSize=size(EV);
    for p=1:EVSize(1,1)
        % exclude negative eigenvalues as these gives rise to (physically meaningless) complex
        % frequencies.
        if(EV(p)>0)
            EVS(i,p)=sqrt(EV(p))*a;
        else
            EVS(i,p)=0;
        end
    end

    my=1;
end

plot(k, EVS);
axis([0,k(Resolution-1)/2,0,2.5]);
