N-layer model. [F. Abeles, Ann. Phys. (Paris) 5, (1950) 596. W. N. Hansen, J. Opt. Soc. Amer. 58 (1968) 380.]
The following equations have not been checked by myself, but it works!!

図 12: N-layer model for SPR measurement.

The tangential fields at the first boundary $z = z_1 = 0$ are related to those at the final boundary $z = z_{N-1}$ by

$$\left[ \begin{array}{c} U_1 \\ V_1 \end{array}\right] = M_... ...left[ \begin{array}{c} U_{N-1} \\ V_{N-1} \end{array}\right]$$ (215)

For p-wave at boundary $k$,

$$U_k = H_y^T + H_y^R$$ (216)

$$V_k = E_x^T + E_x^R$$ (217)

and

$$M_k = \left[ \begin{array}{cc} \cos\beta_k & -i\sin\beta_k / q_k \\ -i q_k \sin\beta_k & \cos \beta_k \end{array}\right]$$ (218)

$${\rm Here} \ \ q_k = (\mu_k /\tilde{\epsilon_k})^{1/2} \cos\theta_k$$ (219)

$$\mu_k \cong 1$$

$$q_k \cong (1 /\tilde{\epsilon_k})^{1/2} \cos\theta_k = \frac{ (\tilde{\epsilon_k} - n_1^2\sin^2\theta_1 )^{1/2} } {\tilde{\epsilon}_k}$$ (220)

$$\beta_k = \frac{2\pi}{\lambda_0} \tilde{n}_k \cos\theta_k (z_k - z_{k-1}) =... ...k-1})\frac{2\pi}{\lambda_0} (\tilde{\epsilon }_k - n_1^2 \sin^2\theta_1 )^{1/2}$$ (221)

The reflection and transmission coefficient for p-wave is

$$r^p = \frac{ (M_{11} + M_{12} q_N )q_1 - ( M_{21} + M_{22} q_N )} {(M_{11} + M_{12} q_N )q_1 + ( M_{21} + M_{22} q_N )}$$ (222)

$$M_{ij} = \left( \prod_{k=2}^{N-1} M_k \right)_{ij}, \ \ \ \ \ i,j = 1, 2$$ (223)

$$R_p = \vert r_p \vert^2$$ (224)

$$r_p = R_p^{1/2} e^{i\phi_p^r}$$ (225)

$$\phi_p^r = {\rm arg}( r^p )$$ (226)

$$t_{H}^p = \frac{2q_1} {(M_{11} + M_{12} q_N )q_1 + ( M_{21} + M_{22} q_N )}$$ (227)

$$t_{E}^p = \frac{\mu_N n_1}{\mu_1 \tilde{n}_N} t_{H}^p$$ (228)

$$T_p = \frac{\mu_N {\rm Re}(\tilde{n}_N \cos\theta_N /\tilde{n}_N^2 )} {\mu_1 n_1 \cos\theta_1 /n_1^2 } \vert t_{H}^p \vert^2$$ (229)

$$\phi_p^t = {\rm arg}( t_{E}^p )$$ (230)

図 13: SPR resonance curves for the SF10(n=1.723)$\vert$Au( $n + ik = 0.1726+ i3.4218$, 50nm)$\vert$Air(n=1.0) and SF10$\vert$Au(50nm)$\vert$SAM(n=1.61245, 1nm)$\vert$Air systems.

The SPR resonance calculation FORTRAN90 programs for N-layer system are

c234567-- spr_angle_Nlayer.f ---
c     1 | 2 | 3 ...  N-2|N-1|N :N layer system
c     complex calculation
      implicit real*8 (a-h,o-z)

      parameter (nlay=10)
      complex*16 e(nlay) , em(nlay,2,2) ,emtot(2,2)
      complex*16 emtot1(2,2)
      dimension en(nlay), ek(nlay), d(nlay)
      complex*16 beta,q,rp,q1,qn,ref,tp,tra
      complex*16 fukso

      c=2.99792458d8
      hbar=6.5822d-16
      pi=acos(-1.0d0)
c     ------- N layer ----
      nlayer=4
c     ------ air ----
      en(4)=1.d0
      ek(4)=0.0d0
      er=en(4)**2-ek(4)**2
      ei=2.0d0*en(4)*ek(4)
      e(4)=dcmplx(er,ei)
c     ----- SF10 633 nm
      en(1)=1.723d0
      ek(1)=0.0d0
      er=en(1)**2-ek(1)**2
      ei=2.0d0*en(1)*ek(1)
      e(1)=dcmplx(er,ei)
c     ----- gold 633 nm
      en(2)=0.1726d0
      ek(2)=3.4218d0
      er=en(2)**2-ek(2)**2
      ei=2.0d0*en(2)*ek(2)
      e(2)=dcmplx(er,ei)
c     ---- gold thickness (m)
      d(2)=50.0d-9
c     ------ SAM -----------
      en(3)=1.61245
      ek(3)=0.0d0
      er=en(3)**2-ek(3)**2
      ei=2.0d0*en(3)*ek(3)
      e(3)=dcmplx(er,ei)
c     ----- SAM thickness ---
      d(3)=1.0d-9
c
      ramd=633.0d-9
      omega=2.0d0*pi/ramd*c
      fukso=dcmplx(0.0d0,1.0d0)
      write (6,*) ramd,omega,hbar*omega
c     --------- angle scan ----------
      ang0=35.0d0
      ang1=45.0d0
      do i=1, 1001
       theta=(ang0+dble(i-1)/1000.0d0*(ang1-ang0))/180.0d0*pi
       q1=sqrt(e(1)-en(1)**2*sin(theta)**2)/e(1)
       qn=sqrt(e(nlayer)-en(1)**2*sin(theta)**2)/e(nlayer)
       do j=2, nlayer-1
        beta=d(j)*2.0d0*pi/ramd*sqrt(e(j)-en(1)**2*sin(theta)**2)
        q=sqrt(e(j)-en(1)**2*sin(theta)**2)/e(j)
        em(j,1,1)=cos(beta)
        em(j,1,2)=-fukso*sin(beta)/q
        em(j,2,1)=-fukso*sin(beta)*q
        em(j,2,2)=cos(beta)
       enddo
        emtot(1,1)=dcmplx(1.0d0,0.0d0)
        emtot(2,2)=dcmplx(1.0d0,0.0d0)
        emtot(1,2)=dcmplx(0.0d0,0.0d0)
        emtot(2,1)=dcmplx(0.0d0,0.0d0)
       do j=2, nlayer-1
        emtot1(1,1)=em(j,1,1)
        emtot1(1,2)=em(j,1,2)
        emtot1(2,1)=em(j,2,1)
        emtot1(2,2)=em(j,2,2)
        emtot=matmul(emtot,emtot1)
       enddo

       rp=( (emtot(1,1)+emtot(1,2)*qn)*q1 -
     &      (emtot(2,1)+emtot(2,2)*qn) )
     &    / ( (emtot(1,1)+emtot(1,2)*qn)*q1 +
     &      (emtot(2,1)+emtot(2,2)*qn) )
       tp=2.0d0*q1/( (emtot(1,1)+emtot(1,2)*qn)*q1 +
     &      (emtot(2,1)+emtot(2,2)*qn) )

       ref=rp*conjg(rp)
       tra=tp*conjg(tp)/cos(theta)*en(1)*dble(qn)
       write (6,*) theta/pi*180.0d0,dble(ref),dble(tra)
      enddo
      end