MATLAB 3D plot on dispersion curve

I am working on wave propagation in multilayer structures. I wrote a code for 2-D plotting dispersion curves at specific frequency (j=1:1 condition) which is presented at the end of the question. My question is that how I can plot those curves at different frequencies (for example, j=1:3) on the same plotting. I think that it should be 3-D plotting and third axes should be frequency. Can anyone give me an idea or source for this condition ?

Thank you,

clc; 
  clear all;

    % System Definition

    N=3; %Number of layers including half spaces N:1 top half space N:bottom half space
    G=zeros((N-1)*4,(N-1)*4);

    %Thickness
    th(1)=1000;%mm thickness of the top half space arbitrary meters
    th(2)=1E-3;%mm thickness of the layer meters
    th(N)=1000;%mm thickness of the layer meters
    %th(4)=1E-3;%mm thickness of the layer meters
    %th(N)=0;%mm thickness of the bottom half space arbitrary meters

    poi=0.3;
    coeff=[0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3];
    mass_coeff=[0.2 0.3 .4 .5 .6 .7 .8 .9 1 2];

    for ii=1:10
    for jj=1:11

    %Material Props
    %Layer 1

    rho(1)=1000; %(if vacuum 0)
    alfa(1)=1480; %(if vacuum 1000)%m/s    longitidunal velocity
    beta(1)=1; %(if vacuum 1000)%m/s  shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition

    %Layer 2

    rho(2)=rho(1)*mass_coeff(ii); %(if vacuum 0)
    alfa(2)=alfa(1)*coeff(jj); %(if vacuum 1000)%m/s    longitidunal velocity
    alfa22(jj)=alfa(2);
    beta(2)=alfa(2)/(sqrt((1-2*poi)/(2*(1-poi))))^-1; %(if vacuum 1000)%m/s  shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition

    %Layer 2

    rho(3)=rho(1); %(if vacuum 0)
    alfa(3)=alfa(1); %(if vacuum 1000)%m/s    longitidunal velocity
    beta(3)=beta(1); %(if vacuum 1000)%m/s  shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition

    % %Layer N
    % rho(4)=rho(2); %(if vacuum 0)
    % alfa(4)=alfa(2); %(if vacuum 1000)%m/s    longitidunal velocity
    % beta(4)=beta(2); %(if vacuum 1000)%m/s  shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition
    % 
    % %Layer N
    % rho(5)=rho(1); %(if vacuum 0)
    % alfa(5)=alfa(1); %(if vacuum 1000)%m/s    longitidunal velocity
    % beta(5)=beta(1); %(if vacuum 1000)%m/s  shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition

    %Velocity Range

    Vstop=2000; %Top Vel
    Vstart=10; %Bottom Vel
    f_step=5E5; %frequency step

    for j=1:1

    %Function Vstart
    f(j)=j*f_step;
    w=2*pi*f(j);

    for i=1:(Vstop-Vstart);

    cp=Vstart+i-1;%km/s
    k=w/cp;

    for r=1:N
    Ca=sqrt(((w^2)/(alfa(r)^2))-k^2);
    Cb=sqrt(((w^2)/(beta(r)^2))-k^2);
    B=(w^2)-(2*(beta(r)^2)*(k^2));
    ga=exp(1i*th(r)*sqrt(((w^2)/(alfa(r)^2))-k^2));
    gb=exp(1i*th(r)*sqrt(((w^2)/(beta(r)^2))-k^2));

    %Wave vectors

    L_plus=[k; Ca; B*1i*rho(r); 2*1i*rho(r)*k*(beta(r)^2)*Ca];
    S_plus=[Cb; -k; -2*1i*rho(r)*k*(beta(r)^2)*Cb; B*1i*rho(r)];
    L_minus=[k; -Ca; B*1i*rho(r); -2*1i*rho(r)*k*(beta(r)^2)*Ca];
    S_minus=[-Cb; -k; 2*1i*rho(r)*k*(beta(r)^2)*Cb; B*1i*rho(r)];

    a=-mod(r,2)*2+1;
    if (r==1)
        G(1:4,1:2)=[-L_minus -S_minus];
    elseif (r==N)
        G((N-1)*4-3:(N-1)*4,(N-1)*4-1:(N-1)*4)=[a*L_plus a*S_plus];
    else
       G((r-1)*4-3:(r-1)*4,(r-1)*4-1:(r-1)*4+2)=[a*L_plus a*S_plus a*L_minus a*S_minus];
       G((r-1)*4+1:(r-1)*4+4,(r-1)*4-1:(r-1)*4+2)=[(a*ga)*L_plus (a*gb)*S_plus (a/ga)*L_minus (a/gb)*S_minus];
    end
    end

    Deter_G_mag(i)=abs(det(G));
    end

    v_range=Vstart:(Vstop-1);
    M(j,:)=Deter_G_mag;
    n=1;

    for m=2:(Vstop-Vstart-1)

        if (m==2)
        der_det(m,j)=(Deter_G_mag(m+1)-Deter_G_mag(1));
        else if (m==(Vstop-Vstart-1))    
        der_det(m,j)=(Deter_G_mag(Vstop-Vstart-1)-Deter_G_mag(m-1));
            else
        der_det(m,j)=(Deter_G_mag(m+1)-Deter_G_mag(m-1));
            end
        end

        if (sign(der_det(m-1,j))==-1 && sign(der_det(m,j))==1 && abs(v_range(m)-alfa(2))>1 && abs(v_range(m)-alfa(1))>1 && abs(v_range(m)-alfa(3))>1 && abs(v_range(m)-beta(1))>1 && abs(v_range(m)-beta(2))>1 && abs(v_range(m)-beta(3))>1);
            V_p(n,j)=v_range(m); 
             n=n+1;
        end

    end
    end
    Vpp(jj)=V_p(1,1);
    end

    ratio=Vpp./alfa22;

    hold on 
    plot(coeff,ratio,'*');

    end

[Ghanem Alatteili]

Hi,

I am working something similar. I am trying to plot a dispersion relation by computiong dipole intaeraction matrix. Then I will insert the translation vectors to get the dipersion relation. I would like to know how do you define your k as a number or a function. I might help.

Thank you,

Answer 1

[MANYINGA MARC]

Hi everyone i can bring you some solution about your problem. My code is base of Stiffness matrix method of Rokhlin and Wang. Here i try to calculate reflexion and transmission coefficient for a multilayer a a function of frequency/incudent angle. My code didn't give me a result what i expected but it can help you for your Problem. And you can also help me to improve it

if true
  % code
end