nh=2.5;
nl=1.5;
nd=2.5;
nstar= 1;
nzero = 1;
nn = [nh,nl,nh,nl,nd,nl,nh,nl,nh,nstar];
% length scale a=1.
wh = 0.25/nh;
wl = 0.25/nl;
% case 1:
wd = 0.25/nd; 
% case 2:
% wd = 0.5/nd;
ww = [wh,wl,wh,wl,wd,wl,wh,wl,wh];
nnlength = 9;

max = 2;
nmax = 400;
for kk=1:nmax
  k0 = 2*pi*max * (kk/nmax);
  R = 0;
  T= 1;
  for j=nnlength:-1:1
    % transition from z_j positive to z_j negative
    P = nn(j+1)/nn(j)*(1-R)/(1+R);
    Rold = R;
    R = (1-P)/(1+P);
    T = T*(1+R)/(1+Rold);
    % transition from z_j negative to z_{j-1} positive
    Q = exp( i * k0 * nn(j) * ww(j));
    T = T*Q;
    R = R*Q^2;
  end
  % from z_0 positive to z_0 negative
  P = nn(1)/nzero*(1-R)/(1+R);
  Rold = R;
  R = (1-P)/(1+P);
  T = T*(1+R)/(1+Rold);
  savek(kk)= max*kk/nmax;
  saveT(kk) = T;
  saveR(kk) = R;
end
subplot(2,1,1)
plot(savek,abs(saveT).^2)
xlabel('k_0a/(2\pi)')
ylabel('|T|^2')
subplot(2,1,2)
plot(savek,abs(saveR).^2)
xlabel('k_0a/(2\pi)')
ylabel('|R|^2')