pcf_1.m

来自「one dimensional PCF band simulation」· M 代码 · 共 59 行

a=1.55+2.8
h1=1.55
step=.02
abegin=0
bend=1
h2=2.8
n1=1.55
n2=2.8
e1=n1^2
e2=n2^2
pi=3.14159
i=sqrt(-1)
for omiga=0:step:bend
for kx=-2 :.02:0
   
kz1=sqrt(  (omiga*2*pi/a)^2*e1-(kx*pi/a)^2     )
   kz2=sqrt(  (omiga*2*pi/a)^2*e2-(kx*pi/a)^2     )
   delta1=e2/e1*kz1/kz2
   p1=kz1*h1
   p2=kz2*h2
T11=(1+delta1)^2/(4*delta1)*exp(i*(p1+p2))-(1-delta1)^2/(4*delta1)*exp(-i*(p2-p1))
T22=(1+delta1)^2/(4*delta1)*exp(-i*(p1+p2))-(1-delta1)^2/(4*delta1)*exp(i*(p2-p1))
pqr=.5*(T11+T22)-sqrt( (.5*(T11+T22) )^2-1 )
K1=i/a*log(.5*(T11+T22)-sqrt( (.5*(T11+T22) )^2-1 ) )
K2=i/a*log(.5*(T11+T22)+sqrt( (.5*(T11+T22) )^2-1 ) )
   if (   (K1>=-pi/a )&(K1<=pi/a)|(K2>=-pi/a )&(K2<=pi/a)  )
    plot(kx,omiga,'r')
    hold on
   end
end
end

for omiga=0:step:bend
for kx=0:step:2*bend
   
kz1=sqrt(  (omiga*2*pi/a)^2*e1-(kx*pi/a)^2     )
   kz2=sqrt(  (omiga*2*pi/a)^2*e2-(kx*pi/a)^2     )
   delta2=kz1/kz2
   p1=kz1*h1
   p2=kz2*h2
T11=(1+delta2)^2/(4*delta2)*exp(i*(p1+p2))-(1-delta2)^2/(4*delta2)*exp(-i*(p2-p1))
T22=(1+delta2)^2/(4*delta2)*exp(-i*(p1+p2))-(1-delta2)^2/(4*delta2)*exp(i*(p2-p1))
pqr=.5*(T11+T22)-sqrt( (.5*(T11+T22) )^2-1 )
K1=i/a*log(.5*(T11+T22)-sqrt( (.5*(T11+T22) )^2-1 ) )
K2=i/a*log(.5*(T11+T22)+sqrt( (.5*(T11+T22) )^2-1 ) )

   if (   (K1>=-pi/a )&(K1<=pi/a)|(K2>=-pi/a )&(K2<=pi/a)  )

    plot(kx,omiga,'r')
    hold on
   end
end
end