📄
字号:
!N2部分
if(((nsx1+1).lt.i.and.i.lt.nx-1).and.(1.lt.j.and.j.lt.ny-1)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))*dx*dy
!write(*,*) 'szhong2'
else
end if
!N3部分
if(((nsx1-1).lt.i.and.i.lt.(nsx1+1)).and.(1.lt.j.and.j.lt.(nsy1-1))) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))*dx*dy
else
end if
!N4部分
if(((nsx1-1).lt.i.and.i.lt.(nsx1+1)).and.((nsy2+1).lt.j.and.j.lt.ny-1)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))*dx*dy
else
end if
!可以综合得：
!szhong=0.0
!if((1.lt.i.and.i.lt.nx-1).and.(1.lt.j.and.j.lt.ny-1).and.(.not.((nsx1-2.lt.i.and.i.lt.nsx1+2).and.(nsy1-2.lt.j.and.j.lt.nsy2+2)))) then
!szhong=szhong+0.5*(jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))
!else
!end if
!N1部分(gaihou)
!goto 2000
1000 nxmin=1
nxmax=nsx1-nx1
nymin=1
nymax=ny-1
if((nxmin.lt.i.and.i.lt.nxmax).and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))*dx*dy
else
end if
!角
if(i==nxmin.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
if(i==nxmin.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
!边
if(i==nxmin.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i,j)**2))*dx*(dy/2.)
end if
if(i==nxmax.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i-1,j)**2))*dx*(dy/2.)
end if
if(j==nymin.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j)**2))*dx*(dy/2.)
end if
if(j==nymax.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j-1)**2))*dx*(dy/2.)
end if
nxmin=nsx1+nx2
nxmax=nx-1
nymin=1
nymax=ny-1
if((nxmin.lt.i.and.i.lt.nxmax).and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))*dx*dy
else
end if
!角
if(i==nxmin.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
if(i==nxmin.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
!边
if(i==nxmin.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i,j)**2))*dx*(dy/2.)
end if
if(i==nxmax.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i-1,j)**2))*dx*(dy/2.)
end if
if(j==nymin.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j)**2))*dx*(dy/2.)
end if
if(j==nymax.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j-1)**2))*dx*(dy/2.)
end if
nxmin=nsx1-nx1
nxmax=nsx1+nx2
nymin=1
nymax=nsy1-ny1
if((nxmin.lt.i.and.i.lt.nxmax).and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))*dx*dy
else
end if
!角
if(i==nxmin.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
if(i==nxmin.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
!边
if(i==nxmin.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i,j)**2))*dx*(dy/2.)
end if
if(i==nxmax.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i-1,j)**2))*dx*(dy/2.)
end if
if(j==nymin.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j)**2))*dx*(dy/2.)
end if
if(j==nymax.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j-1)**2))*dx*(dy/2.)
end if
nxmin=nsx1-nx1
nxmax=nsx1+nx2
nymin=nsy2+ny2
nymax=ny-1
if((nxmin.lt.i.and.i.lt.nxmax).and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+((hy(i,j)+hy(i-1,j))/2.)**2))*dx/dy
else
end if
!角
if(i==nxmin.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymin) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
if(i==nxmin.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i,j)**2))*dx*(dy/4.)
end if
if(i==nxmax.and.j==nymax) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(hx(i,j-1)**2+hy(i-1,j)**2))*dx*(dy/4.)
end if
!边
if(i==nxmin.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i,j)**2))*dx*(dy/2.)
end if
if(i==nxmax.and.(nymin.lt.j.and.j.lt.nymax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hx(i,j)+hx(i,j-1))/2.)**2+hy(i-1,j)**2))*dx*(dy/2.)
end if
if(j==nymin.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j)**2))*dx*(dy/2.)
end if
if(j==nymax.and.(nxmin.lt.i.and.i.lt.nxmax)) then
szhong=szhong+0.5*(jiedian0*jiedian11(i,j)*(ez(i,j)**2)+cidao0*(((hy(i,j)+hy(i-1,j))/2.)**2+hx(i,j-1)**2))*dx*(dy/2.)
end if
2000 end do
end do
if (0.le.nt.and.nt.le.2000) then
write(*,*) '时间步:',nt ,'能量守衡判断：'
write(*,*) 'sru:', sru,'schu:',schu,'szhong:',szhong
if(szhong.ne.0.0) write(*,*) '能量误差为',100-abs(100*(schu+sru)/szhong),'%'
! 188 format(1x,a, es13.6e2,5x,a, es13.6e2,5x,a,es13.6e2,5x,a,es13.6e2)
else
end if
return
end subroutine
subroutine touguolv1
implicit none
integer i,j
if (3000.le.nt.and.nt.le.5999) then
do i=0,nx
do j=0,ny
if(i==nx-30.and.(235.lt.j.and.j.lt.275)) then
schu0=schu0+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
else
end if
if(i==10.and.(135.lt.j.and.j.lt.175)) then
sru0=sru0+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
!sru0=sru0+(-1.)*source(i,j)*source(i,j)*sqrt(jiedian0/cidao0)*dy*(-1.)*dt
else
end if
end do
end do
else
end if
return
end subroutine
subroutine touguolv2
implicit none
integer i,j
if (3000.le.nt.and.nt.le.5999) then
do i=0,nx
do j=0,ny
if(i==nx-3.and.(30.lt.j.and.j.lt.80)) then
!schu01=schu01+0.5*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*(-1.)*dt
schu01=schu01+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
!write(*,*) 'schu3'
else
end if
if(i==3.and.(130.lt.j.and.j.lt.180)) then
sru01=sru01+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
!sru0=sru0+(-1.)*source(i,j)*source(i,j)*sqrt(jiedian0/cidao0)*dy*(-1.)*dt
else
end if
end do
end do
else
end if
return
end subroutine
end
function even_jiedian(i,j,dx,dy,r,v_a,nx,ny,grid_vs_colomnx,grid_vs_colomny,jiedian1,jiedian2,nx_colomn)
implicit none
integer grid_vs_colomnx,grid_vs_colomny
integer i,j
integer n,m
integer nx,ny ,nx_colomn
real*8 dx, dy,r
real*8 v_a
real*8 h1,h2
real*8 s1,s2
real*8 s0,s
real*8 angle_a,angle_b
real*8 jiedian1,jiedian2
real*8 even_jiedian
real*8 l,lmin
real*8 d
real*8, parameter::pi=3.141592654d0
d=(dx+dy)/4
lmin=2*v_a
s0=pi*(d**2)
nx_colomn=50
if(r<=d) then
write(*,*) 'The choice of photonic crystals parameter has errors,please check&
perhaps it is because of spacestep too large or the radius of colomn&
too small.please debug!'
return
else
end if
do n=0,nx,grid_vs_colomnx
do m=0,ny,grid_vs_colomny
l=sqrt((i*dx-n*dx)**2+(j*dy-m*dy)**2)
if(lmin>l) then
lmin=l
else
end if
10 end do
end do
if(lmin>=r+d) then
even_jiedian=jiedian1
else
if(lmin<=r-d) then
even_jiedian=jiedian2
else
if (r-d<lmin.and.lmin<r+d) then
if(lmin.lt.r.and.d.le.sqrt(r**2-lmin**2)) then
h1=-(d**2+lmin**2-r**2)/(2*lmin)
h2=(r**2+lmin**2-d**2)/(2*lmin)
s=sqrt(d**2-h1**2)
angle_a=asin(s/d)
angle_b=asin(s/r)
s1=(d**2)*angle_a-h1*s
s2=(r**2)*angle_b-h2*s
s=pi*d**2-(s1-s2)
even_jiedian=((s0-s)*jiedian1+s*jiedian2)/s0
else
h1=(d**2+lmin**2-r**2)/(2*lmin)
h2=(r**2+lmin**2-d**2)/(2*lmin)
s=sqrt(d**2-h1**2)
angle_a=asin(s/d)
angle_b=asin(s/r)
s1=(d**2)*angle_a-h1*s
s2=(r**2)*angle_b-h2*s
s=s1+s2
even_jiedian=((s0-s)*jiedian1+s*jiedian2)/s0
end if
end if
end if
end if
end function even_jiedian
SUBROUTINE fft(xt,PI,N,K)
REAL*8,DIMENSION(N):: xt,PI,FR,FI
REAL*8:: P,Q,S,VR,VI,PODDR,PODDI
REAL*8,PARAMETER::PAI=3.1415926
integer L,IL
L=0
IL=1
DO IT=0,N-1
M=IT
IS=0
DO I=0,K-1
J=M/2
IS=2*IS+(M-2*J)
M=J
END DO
FR(IT+1)=xt(IS+1)
FI(IT+1)=PI(IS+1)
END DO
xt(1)=1.
PI(1)=0.
xt(2)=COS(2*PAI/N)
PI(2)=-SIN(2*PAI/N)
IF(L.NE.0)PI(2)=-PI(2)
DO I=3,N
P=xt(I-1)*xt(2)
Q=PI(I-1)*PI(2)
S=(xt(I-1)+PI(I-1))*(xt(2)+PI(2))
xt(I)=P-Q
PI(I)=S-P-Q
END DO
DO IT=0,N-2,2
VR=FR(IT+1)
VI=FI(IT+1)
FR(IT+1)=VR+FR(IT+2)
FI(IT+1)=VI+FI(IT+2)
FR(IT+2)=VR-FR(IT+2)
FI(IT+2)=VI-FI(IT+2)
END DO
M=N/2
NV=2
DO L0=K-2,0,-1
M=M/2
NV=2*NV
DO IT=0,(M-1)*NV,NV
DO J=0,(NV/2)-1
P=xt(M*J+1)*FR(IT+J+1+NV/2)
Q=PI(M*J+1)*FI(IT+J+1+NV/2)
S=xt(M*J+1)+PI(M*J+1)
S=S*(FR(IT+J+1+NV/2)+FI(IT+J+1+NV/2))
PODDR=P-Q
PODDI=S-P-Q
FR(IT+J+1+NV/2)=FR(IT+J+1)-PODDR
FI(IT+J+1+NV/2)=FI(IT+J+1)-PODDI
FR(IT+J+1)=FR(IT+J+1)+PODDR
FI(IT+J+1)=FI(IT+J+1)+PODDI
END DO
END DO
END DO
IF(L.NE.0)THEN
DO I=1,N
FR(I)=FR(I)/N
FI(I)=FI(I)/N
END DO
ENDIF
IF(IL.NE.0)THEN
DO I=1,N
xt(I)=SQRT(FR(I)*FR(I)+FI(I)*FI(I))
PI(I)=ATAN(FI(I)/FR(I))*180/PAI
END DO
END IF
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -