%=========================================
% box_plot_theory_modal_vibes.m
% written by prof. s. errede 03/05/2009
%=========================================
close all;
clear all;

double          xs(10000);
double          ys(10000);
double          zs(10000);

double displ_amplx(10000);
double displ_amply(10000);
double displ_amplz(10000);

double press_amplx(10000);
double press_amply(10000);
double press_amplz(10000);

double        xloc(32,32);
double        yloc(32,32);
double        zloc(32,32);

double  displ_ampl(32,32);

double  press_ampl(32,32);


%====================================================
% Specify 3D box dimensions: 
% n.b. must be within +_15.5 units x-y-z scan
% n.b. must have ordering Lx > Ly > Lz
%====================================================
Lx       =  30.0;
Ly       =  25.0;
Lz       =  20.0;

x_lo     = -15.5;
x_hi     =  15.5;
y_lo     = -15.5;
y_hi     =  15.5;
z_lo     = -15.5;
z_hi     =  15.5;

pi = 4.0*atan(1.0);

% ===============================================================
% n.b. need to change the following items for each mode:
%
%   1: Specify the mode indices l, m, n (= integers: 0,1,2,3...)
%   2: Specify the A_l,m,n mode in all figure titles (below)
%
% ===============================================================
l = 3; % x mode index
m = 2; % y mode index
n = 0; % z mode index

l_odd = ODD(l); % odd(x) rounds x to nearest odd integer
l_parity = 0;
if (l_odd == l); % then l is odd
    l_parity = 1; % set l-parity flag
end;

m_odd = ODD(m); % odd(x) rounds x to nearest odd integer
m_parity = 0;
if (m_odd == m); % then m is odd
    m_parity = 1; % set m-parity flag
end;

n_odd = ODD(l); % odd(x) rounds x to nearest odd integer
n_parity = 0;
if (n_odd == n); % then n is odd
    n_parity = 1; % set n-parity flag
end;

% wavenumbers:
k_l = l*pi/Lx;
k_m = m*pi/Ly;
k_n = n*pi/Lz;

% ==================================================================
% 1-D scan in x of D & P amplitudes at fixed y = 0 and t = 0:
% ==================================================================
for i = 1:3100;
    xs(i) = 0.01*(i-1) + x_lo;
    ys(i) = 0.01*(i-1) + y_lo;
    zs(i) = 0.01*(i-1) + z_lo;
    
    displ_amplx(i) = 0.0;
    press_amplx(i) = 0.0;
    if (abs(xs(i)) < Lx/2.0);
        if ((l > 0) && (m == 0) && (n == 0)); % 1-D l00 axial modes
            if (l_parity == 0);
                displ_amplx(i) = 1.0*sin(k_l*xs(i));
                press_amplx(i) = 1.0*cos(k_l*xs(i));
            end;
            if (l_parity == 1);
                displ_amplx(i) = 1.0*cos(k_l*xs(i));
                press_amplx(i) = 1.0*sin(k_l*xs(i));
            end;
        end;
        if ((l > 0) && (m > 0)  && (n == 0)); % 2-D lm0 tangential modes
            if (l_parity == 0);
                displ_amplx(i) = 1.0*sin(k_l*xs(i));
                press_amplx(i) = 1.0*cos(k_l*xs(i));
            end;
            if (l_parity == 1);
                displ_amplx(i) = 1.0*cos(k_l*xs(i));
                press_amplx(i) = 1.0*sin(k_l*xs(i));
            end;
        end;
        if ((l > 0) && (m == 0)  && (n > 0)); % 2-D l0n tangential modes
            if (l_parity == 0);
                displ_amplx(i) = 1.0*sin(k_l*xs(i));
                press_amplx(i) = 1.0*cos(k_l*xs(i));
            end;
            if (l_parity == 1);
                displ_amplx(i) = 1.0*cos(k_l*xs(i));
                press_amplx(i) = 1.0*sin(k_l*xs(i));
            end;
       end;
    end;

    displ_amply(i) = 0.0;
    press_amply(i) = 0.0;
    if (abs(ys(i)) < Ly/2.0);
        if ((l == 0) && (m > 0) && (n == 0)); % 1-D 0l0 axial modes
            if (m_parity == 0);
                displ_amply(i) = 1.0*sin(k_m*ys(i));
                press_amply(i) = 1.0*cos(k_m*ys(i));
            end;
            if (m_parity == 1);
                displ_amply(i) = 1.0*cos(k_m*ys(i));
                press_amply(i) = 1.0*sin(k_m*ys(i));
            end;
        end;
        if ((l > 0) && (m > 0)  && (n == 0)); % 2-D lm0 tangential modes
            if (m_parity == 0);
                displ_amply(i) = 1.0*sin(k_m*ys(i));
                press_amply(i) = 1.0*cos(k_m*ys(i));
            end;
            if (m_parity == 1);
                displ_amply(i) = 1.0*cos(k_m*ys(i));
                press_amply(i) = 1.0*sin(k_m*ys(i));
            end;
        end;
        if ((l == 0) && (m > 0)  && (n > 0)); % 2-D 0mn tangential modes
            if (m_parity == 0);
                displ_amply(i) = 1.0*sin(k_m*ys(i));
                press_amply(i) = 1.0*cos(k_m*ys(i));
            end;
            if (m_parity == 1);
                displ_amply(i) = 1.0*cos(k_m*ys(i));
                press_amply(i) = 1.0*sin(k_m*ys(i));
            end;
       end;
    end;

    displ_amplz(i) = 0.0;
    press_amplz(i) = 0.0;
    if (abs(zs(i)) < Lz/2.0);
        if ((l == 0) && (m == 0) && (n > 0)); % 1-D 00l axial modes
            if (n_parity == 0);
                displ_amplz(i) = 1.0*sin(k_n*zs(i));
                press_amplz(i) = 1.0*cos(k_n*zs(i));
            end;
            if (n_parity == 1);
                displ_amplz(i) = 1.0*cos(k_n*zs(i));
                press_amplz(i) = 1.0*sin(k_n*zs(i));
            end;
        end;
        if ((l == 0) && (m > 0)  && (n > 0)); % 2-D 0mn tangential modes
            if (n_parity == 0);
                displ_amplz(i) = 1.0*sin(k_n*zs(i));
                press_amplz(i) = 1.0*cos(k_n*zs(i));
            end;
            if (n_parity == 1);
                displ_amplz(i) = 1.0*cos(k_n*zs(i));
                press_amplz(i) = 1.0*sin(k_n*zs(i));
            end;
        end;
        if ((l > 0) && (m == 0)  && (n > 0)); % 2-D l0n tangential modes
            if (n_parity == 0);
                displ_amplz(i) = 1.0*sin(k_n*zs(i));
                press_amplz(i) = 1.0*cos(k_n*zs(i));
            end;
            if (n_parity == 1);
                displ_amplz(i) = 1.0*cos(k_n*zs(i));
                press_amplz(i) = 1.0*sin(k_n*zs(i));
            end;
       end;
    end;
end;

% ====================================================================
% 2-D scan in x-y of D & P amplitudes at fixed z = 0 and t = 0:
% ====================================================================
for j=1:32;
    y = 1.0*(j-1) + y_lo;
    for i = 1:32;
        x    = 1.0*(i-1) + x_lo;
        
        xloc(i,j) = x;
        yloc(i,j) = y;

        displ_ampl(i,j) = 0.0;
        press_ampl(i,j) = 0.0;
        if ((abs(x) < Lx/2.0) && (abs(y) < Ly/2.0));
            if ((l > 0) && (m == 0) && (n == 0)); % l00 axial modes
                if (l_parity == 0);
                    displ_ampl(i,j) = 1.0*sin(k_l*x);
                    press_ampl(i,j) = 1.0*cos(k_l*x);
                end;
                if (l_parity == 1);
                    displ_ampl(i,j) = 1.0*cos(k_l*x);
                    press_ampl(i,j) = 1.0*sin(k_l*x);
                end;
            end;
            if ((l == 0) && (m > 0) && (n == 0)); % 0m0 axial modes
                if (m_parity == 0);
                    displ_ampl(i,j) = 1.0*sin(k_m*y);
                    press_ampl(i,j) = 1.0*cos(k_m*y);
                end;
                if (m_parity == 1);
                    displ_ampl(i,j) = 1.0*cos(k_m*y);
                    press_ampl(i,j) = 1.0*sin(k_m*y);
                end;
            end;
            if ((l > 0) && (m > 0)  && (n == 0)); % lm0 tangential modes
                if ((l_parity == 0) && (m_parity == 0));
                    displ_ampl(i,j) = 1.0*sin(k_l*x)*sin(k_m*y);
                    press_ampl(i,j) = 1.0*cos(k_l*x)*cos(k_m*y);
                end;
                if ((l_parity == 1) && (m_parity == 0));
                    displ_ampl(i,j) = 1.0*cos(k_l*x)*sin(k_m*y);
                    press_ampl(i,j) = 1.0*sin(k_l*x)*cos(k_m*y);
                end;
                if ((l_parity == 0) && (m_parity == 1));
                    displ_ampl(i,j) = 1.0*sin(k_l*x)*cos(k_m*y);
                    press_ampl(i,j) = 1.0*cos(k_l*x)*sin(k_m*y);
                end;
                if ((l_parity == 1) && (m_parity == 1));
                    displ_ampl(i,j) = 1.0*cos(k_l*x)*cos(k_m*y);
                    press_ampl(i,j) = 1.0*sin(k_l*x)*sin(k_m*y);
                end;
            end;
        end;
    end;
end;

xg(1) = -20.0;
xg(2) =  20.0;

yg(1) =   0.0;
yg(2) =   0.0;

figure(01);
plot(xs,displ_amplx,'b',xg,yg,'k');
axis tight
grid on
xlabel('X')
ylabel('Amplitude (Arb. Units)')
title('320 Mode X-Displacement Amplitude @(x,0,0)');

figure(02);
plot(ys,displ_amply,'b',xg,yg,'k');
axis tight
grid on
xlabel('Y')
ylabel('Amplitude (Arb. Units)')
title('320 Mode Y-Displacement Amplitude @(x,0,0)');

figure(03);
plot(zs,displ_amplz,'b',xg,yg,'k');
axis tight
grid on
xlabel('Z')
ylabel('Amplitude (Arb. Units)')
title('320 Mode Z-Displacement Amplitude @(x,0,0)');

figure(04);
plot(xs,press_amplx,'b',xg,yg,'k');
axis tight
grid on
xlabel('X')
ylabel('Amplitude (Arb. Units)')
title('320 Mode Pressure Amplitude @(x,0,0)');

figure(05);
plot(ys,press_amply,'b',xg,yg,'k');
axis tight
grid on
xlabel('Y')
ylabel('Amplitude (Arb. Units)')
title('320 Mode Pressure Amplitude @(x,0,0)');

figure(06);
plot(zs,press_amplz,'b',xg,yg,'k');
axis tight
grid on
xlabel('Z')
ylabel('Amplitude (Arb. Units)')
title('320 Mode Pressure Amplitude @(x,0,0)');

figure(11);
surf(xloc,yloc,displ_ampl);
shading interp;
xlabel('x');
ylabel('y');
zlabel('Amplitude (arb. units)');
title ('320 Mode Displacement Amplitude @(x,y,0)');

figure(12);
surf(xloc,yloc,press_ampl);
shading interp;
xlabel('x');
ylabel('y');
zlabel('Amplitude (arb. units)');
title ('320 Mode Pressure Amplitude @(x,y,0)');

%==============================================================