clc
clear all
close all
syms x 
rho=1.2;%air density
f=10000%frequency Hz
w=2*pi*f;%frequency
c=343;%speed of sound
k=w/c;
z1=1.2*343;%impedance on left
z4=z1;%impedance on right
xo=.02;
po=5;
l=.073;
%equations
general_form1=[exp(-i*k*x),exp(i*k*x),0,0];
 general_form2=[0,0,exp(-i*k*x),exp(i*k*x)];
 u1=1/(rho*i*w).*(-1*diff(general_form1,x));
u2=1/(rho*i*w).*(-1*diff(general_form2,x));
x=0;
l1=[eval(general_form1-z1.*(u1))];
x=xo;
l2=[eval(diff(general_form1,x)-diff(general_form2,x))];
l3=[eval(general_form1)];
x=l;
l4=[eval(general_form2-z4.*(u2))];

L=[l1;l2;l3;l4]
R=[0 0 po 0]'
P=inv(L)*R
num=100;
x=linspace(0,l,num);
for ii=1:num
    if x(ii)<=xo
        Pressure(ii)=P(1)*exp(-i*k*x(ii))+P(2)*exp(i*k*x(ii));
    else
        Pressure(ii)=P(3)*exp(-i*k*x(ii))+P(4)*exp(i*k*x(ii));
    end
end
PL=P(1).*exp(-i*k.*x)+P(2).*exp(i*k.*x);
PR=P(3).*exp(-i*k.*x)+P(4).*exp(i*k.*x);
plot(x,(Pressure),x,imag(Pressure)),legend real imag
% figure
% plot(x,PL,x,PR)
% legend 1 2
kk=1;
jj=20;
  L=[exp(-i*k*x(kk)) exp(i*k*x(kk));
        exp(-i*k*x(jj)) exp(i*k*x(jj))];
    R=[Pressure(kk);Pressure(jj)];
    L^-1*R
    
%  Q=[x' PRES'];
%  save abspress10000.txt Q -ascii