subroutine phasmat(Z,S11,S12,S21,S22)
      real Z(4,4)
      COMPLEX*16 S11,S12,S21,S22
Cf2py intent(out) Z
      Z(1,1)=0.5D0*(S11*DCONJG(S11)+S12*DCONJG(S12)
     &          +S21*DCONJG(S21)+S22*DCONJG(S22))
      Z(1,2)=0.5D0*(S11*DCONJG(S11)-S12*DCONJG(S12)
     &          +S21*DCONJG(S21)-S22*DCONJG(S22))
      Z(1,3)=-S11*DCONJG(S12)-S22*DCONJG(S21)
      Z(1,4)=(0D0,1D0)*(S11*DCONJG(S12)-S22*DCONJG(S21))
      Z(2,1)=0.5D0*(S11*DCONJG(S11)+S12*DCONJG(S12)
     &          -S21*DCONJG(S21)-S22*DCONJG(S22))
      Z(2,2)=0.5D0*(S11*DCONJG(S11)-S12*DCONJG(S12)
     &          -S21*DCONJG(S21)+S22*DCONJG(S22))
      Z(2,3)=-S11*DCONJG(S12)+S22*DCONJG(S21)
      Z(2,4)=(0D0,1D0)*(S11*DCONJG(S12)+S22*DCONJG(S21))
      Z(3,1)=-S11*DCONJG(S21)-S22*DCONJG(S12)
      Z(3,2)=-S11*DCONJG(S21)+S22*DCONJG(S12)
      Z(3,3)=S11*DCONJG(S22)+S12*DCONJG(S21)
      Z(3,4)=(0D0,-1D0)*(S11*DCONJG(S22)+S21*DCONJG(S12))
      Z(4,1)=(0D0,1D0)*(S21*DCONJG(S11)+S22*DCONJG(S12))
      Z(4,2)=(0D0,1D0)*(S21*DCONJG(S11)-S22*DCONJG(S12))
      Z(4,3)=(0D0,-1D0)*(S22*DCONJG(S11)-S12*DCONJG(S21))
      Z(4,4)=S22*DCONJG(S11)-S12*DCONJG(S21)
      return
      end subroutine phasmat