(* Mathematica program for Problem 4.2. For details see the sugestions 
   in the homework.

  date : 22/09/2013
  author: Jorge C. Romão
  email : jorge.romao@tecnico.ulisboa.pt
*)

(* Result for M1 and M2 *)

M1={{1/2 (1+I k/K) E^(K a) E^(-I k a),
     1/2 (1-I k/K) E^(K a) E^(I k a)},
    {1/2 (1-I k/K) E^(-K a) E^(-I k a),
     1/2 (1+I k/K) E^(-K a) E^(I k a)}};

M2={{1/2 (1-I K/k) E^(K a) E^(-I k a),
     1/2 (1+I K/k) E^(-K a) E^(-I k a)},
    {1/2 (1+I K/k) E^(K a) E^(I k a),
     1/2 (1-I K/k) E^(-K a) E^(I k a)}};

(* Obtain Q *)

Qaux= Simplify[M2 . M1] ;

(* Simpler form for Q *)

Q=Function[{a},Q11=I (Sinh[2 K a] (k^2-K^2)- 2 I k K Cosh[2 K a])
	/(2 k K) E^(-2 I k a); 
Q12=-I Sinh[2 K a] (k^2 + K^2)/(2 k K);Q21=-Q12; 
Q22=-I (Sinh[2 K a] (k^2-K^2)+ 2 I k K Cosh[2 K a])/(2 k K) E^(2 I k a); 
{{Q11,Q12},{Q21,Q22}}];

(* Check that the two forms are equivalent *)

TestQ=Simplify[Qaux - Q[a]];

Print["Test Q=M2*M1 (If matrix is null OK).   Q - M2 . M1=",TestQ]

(* Get the Reflection and Transmission Coefficients for
 waves incident from the left *)

T=Function[M,M[[1,1]]-M[[1,2]]*M[[2,1]]/M[[2,2]]];
R=Function[M,-M[[2,1]]/M[[2,2]]];

ReT=ComplexExpand[Re[Simplify[T[Q[a]]]]];
ImT=ComplexExpand[Im[Simplify[T[Q[a]]]]];
T2=Simplify[ReT^2+ImT^2];

ReR=ComplexExpand[Re[Simplify[R[Q[a]]]]];
ImR=ComplexExpand[Im[Simplify[R[Q[a]]]]];
R2=Simplify[ReR^2+ImR^2];

Print["|T|^2 = ",T2]
Print["|R|^2 = ", R2]
Print["|T|^2 + |R|^2 = ",Simplify[R2+T2]]