📄
字号:
call fft(xt5,yt5,n0,m0)
call fft(xt6,yt6,n0,m0)
open(1,file='e:\work\matlabm\x56.m',action='write',form='formatted',access='sequential')
write(1,14) 'df=',1.0/(n0*dt),';','a=',v_a,';','dt=',dt,';'
!maxx=sqrt(xt5(1)**2+yt5(1)**2)
!do i=1,8192,1
!if(sqrt(xt5(i)**2+yt5(i)**2).gt.maxx) maxx=sqrt(xt5(i)**2+yt5(i)**2)
!end do
!maxx=1.
write(1,16) (xt5(i)/xt6(i),i=1,8192,1)
write(1,*) '];'
close(1)
call fft(xt7,yt7,n0,m0)
call fft(xt8,yt8,n0,m0)
open(1,file='e:\work\matlabm\x78.m',action='write',form='formatted',access='sequential')
write(1,14) 'df=',1.0/(n0*dt),';','a=',v_a,';','dt=',dt,';'
!maxx=sqrt(xt6(1)**2+yt6(1)**2)
!do i=1,8192,1
!if(sqrt(xt6(i)**2+yt6(i)**2).gt.maxx) maxx=sqrt(xt6(i)**2+yt6(i)**2)
!end do
!maxx=1.
write(1,16) (xt7(i)/xt8(i),i=1,8192,1)
write(1,*) '];'
close(1)
call fft(xt9,yt9,n0,m0)
open(1,file='e:\work\matlabm\x96.m',action='write',form='formatted',access='sequential')
write(1,14) 'df=',1.0/(n0*dt),';','a=',v_a,';','dt=',dt,';'
!maxx=sqrt(xt6(1)**2+yt6(1)**2)
!do i=1,8192,1
!if(sqrt(xt6(i)**2+yt6(i)**2).gt.maxx) maxx=sqrt(xt6(i)**2+yt6(i)**2)
!end do
!maxx=1.
write(1,16) (xt9(i)/xt6(i),i=1,8192,1)
write(1,*) '];'
close(1)
!!!!!!!!!!!!!!!
! call fft(xt7,yt7,n0,m0)
!open(1,file='e:\work\matlabm\fftzzh7.m',action='write',form='formatted',access='sequential')
!write(1,14) 'df=',1.0/(n0*dt),';','a=',v_a,';','dt=',dt,';'
!maxx=sqrt(xt7(1)**2+yt7(1)**2)
!do i=1,8192,1
!if(sqrt(xt7(i)**2+yt7(i)**2).gt.maxx) maxx=sqrt(xt7(i)**2+yt7(i)**2)
!end do
!maxx=1.
!write(1,16) (sqrt(xt7(i)**2+yt7(i)**2)/maxx,i=1,8192,1)
!write(1,*) '];'
!close(1)
!call fft(xt8,yt8,n0,m0)
!open(1,file='e:\work\matlabm\fftzzh8.m',action='write',form='formatted',access='sequential')
!write(1,14) 'df=',1.0/(n0*dt),';','a=',v_a,';','dt=',dt,';'
!maxx=sqrt(xt8(1)**2+yt8(1)**2)
!do i=1,8192,1
!if(sqrt(xt8(i)**2+yt8(i)**2).gt.maxx) maxx=sqrt(xt8(i)**2+yt8(i)**2)
!end do
!maxx=1.
!write(1,16) (sqrt(xt8(i)**2+yt8(i)**2)/maxx,i=1,8192,1)
!write(1,*) '];'
!close(1)
! call fft(xt9,yt9,n0,m0)
!open(1,file='e:\work\matlabm\fftzzh9.m',action='write',form='formatted',access='sequential')
!write(1,14) 'df=',1.0/(n0*dt),';','a=',v_a,';','dt=',dt,';'
!maxx=sqrt(xt9(1)**2+yt9(1)**2)
!do i=1,8192,1
!if(sqrt(xt9(i)**2+yt9(i)**2).gt.maxx) maxx=sqrt(xt9(i)**2+yt9(i)**2)
!end do
!maxx=1.
!write(1,16) (sqrt(xt9(i)**2+yt9(i)**2)/maxx,i=1,8192,1)
!write(1,*) '];'
!close(1)
!write(*,*) '####################################################'
17 format(1x,a)
11111 end do
end do
write(*,*) '计算已经完成'
stop
contains
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!###############################################################
subroutine ezfield
implicit none
integer i,j
do i=1,nx-1
do j=1,ny-1
ez(i,j)=ez(i,j)+dt/(jiedian0*jiedian11(i,j)*dx)*(hy(i,j)-hy(i-1,j))&
-dt/(jiedian0*jiedian11(i,j)*dy)*(hx(i,j)-hx(i,j-1))&
+source(i,j)
!ez(i,j)=ez(i,j)+1.0/(2.0*jiedian11(i,j))*(hy(i,j)-hy(i-1,j))&
! -1.0/(2.0*jiedian11(i,j))*(hx(i,j)-hx(i,j-1))&
! +source(i,j)
end do
end do
end subroutine
subroutine hxfield
implicit none
integer i,j
do i=0,nx,1
do j=0,ny-1,1
hx(i,j)=hx(i,j)-(dt/(cidao0*dy))*(ez(i,j+1)-ez(i,j))
!+source(i,j)*sin(30.*pi/180.)
end do
end do
end subroutine
subroutine hyfield
implicit none
integer i,j
do i=0,nx-1,1
do j=0,ny,1
hy(i,j)=hy(i,j)+(dt/(cidao0*dx))*(ez(i+1,j)-ez(i,j))
!-source(i,j)*cos(30.*pi/180.)
!hy(i,j)=hy(i,j)+(1.0/2.0)*(ez(i+1,j)-ez(i,j))
end do
end do
end subroutine
subroutine boundary
implicit none
integer i,j
do j=1,ny-1
ez(0,j)=-ez2(1,j)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(1,j)+ez2(0,j))&
+(2.0*dx)/(light_speed*dt+dx)*(ez1(0,j)+ez1(1,j))&
+(light_speed*dt*light_speed*dt*dx)/(2.0*dy*dy*(light_speed*dt+dx))&
*(ez1(0,j+1)-2.0*ez1(0,j)+ez1(0,j-1)+ez1(1,j+1)-2.0*ez1(1,j)+ez1(1,j-1))
end do
do j=1,ny-1
ez(nx,j)=-ez2(nx-1,j)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(nx-1,j)+ez2(nx,j))&
+(2.0*dx)/(light_speed*dt+dx)*(ez1(nx,j)+ez1(nx-1,j))&
+(light_speed*dt*light_speed*dt*dx)/(2.0*dy*dy*(light_speed*dt+dx))&
*(ez1(nx,j+1)-2.0*ez1(nx,j)+ez1(nx,j-1)+ez1(nx-1,j+1)-2*ez1(nx-1,j)+ez1(nx-1,j-1))
end do
do i=1,nx-1
ez(i,0)=-ez2(i,1)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(i,1)+ez2(i,0))&
+(2.0*dx)/(light_speed*dt+dx)*(ez1(i,0)+ez1(i,1))&
+(light_speed*dt*light_speed*dt*dx)/(2.0*dy*dy*(light_speed*dt+dx))&
*(ez1(i+1,0)-2.0*ez1(i,0)+ez1(i-1,0)+ez1(i+1,1)-2.0*ez1(i,1)+ez1(i-1,1))
end do
do i=1,nx-1
ez(i,ny)=-ez2(i,ny-1)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(i,ny-1)+ez2(i,ny))&
+(2.0*dx)/(light_speed*dt+dx)*(ez1(i,ny)+ez1(i,ny-1))&
+(light_speed*dt*light_speed*dt*dx)/(2.0*dy*dy*(light_speed*dt+dx))&
*(ez1(i+1,ny)-2.0*ez1(i,ny)+ez1(i-1,ny)+ez1(i+1,ny-1)-2*ez1(i,ny-1)+ez1(i-1,ny-1))
end do
!ez(0,0)=ez1(1,0)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(1,0)-ez1(0,0))
!ez(0,ny)=ez1(1,ny)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(1,ny)-ez1(0,ny))
!ez(nx,0)=ez1(nx-1,0)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(nx-1,0)-ez1(nx,0))
!ez(nx,ny)=ez1(nx-1,ny)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ez(nx-1,ny)-ez1(nx,ny))
ez(0,0)=frad*((1-sina1)*(1-cosa1)*ez2(0,0)+(1-sina1)*cosa1*ez2(1,0)+sina1*(1-cosa1)*ez2(0,1)&
+sina1*cosa1*ez2(1,1))
ez(nx,0)=frad*((1-sina1)*(1-cosa1)*ez2(nx,0)+(1-sina1)*cosa1*ez2(nx-1,0)+sina1*(1-cosa1)*ez2(nx,1)&
+sina1*cosa1*ez2(nx-1,1))
ez(nx,ny)=frad*((1-sina1)*(1-cosa1)*ez2(nx,ny)+(1-sina1)*cosa1*ez2(nx-1,ny)+sina1*(1-cosa1)*ez2(nx,ny-1)&
+sina1*cosa1*ez2(nx-1,ny-1))
ez(0,ny)=frad*((1-sina1)*(1-cosa1)*ez2(0,ny)+(1-sina1)*cosa1*ez2(1,ny)+sina1*(1-cosa1)*ez2(0,ny-1)&
+sina1*cosa1*ez2(1,ny-1))
do i=0,nx,1
do j=0,ny,1
ez2(i,j)=ez1(i,j)
ez1(i,j)=ez(i,j)
end do
end do
return
end subroutine
function source(i,j)
implicit none
integer i,j
real*8 source
if(wave=='sinee') then
source=sine(i,j)
else
if(wave=='pulse') then
source=pulse(i,j)
else if(wave/='sine'.or.wave=='pulse') then
write(*,*) ' 你在进行波源的选择时有错误，请核查'
end if
end if
return
end function
function sine(i,j)
implicit none
integer i,j
real*8 sine
if(nsx1.le.i.and.i.le.nsx2.and.nsy1.le.j.and.j.le.nsy2) then
sine=1.0e4*sin(w*t)*exp(real(-(j-j0)**2/jw**2))!+&
!1.0e4*sin(0.395*(2*pi*light_speed/v_a)*t)*exp(real(-(j-j0+100)**2/jw**2))+&
! 1.0e4*sin(0.445*(2*pi*light_speed/v_a)*t)*exp(real(-(j-j0+200)**2/jw**2))
else
sine=0.0
end if
return
end function
function pulse(i,j)
implicit none
integer i,j
real*8 t0,pt
real*8 pulse
t0=21.0*dt
pt=2.0*dt
if(nsx1.le.i.and.i.le.nsx2.and.nsy1.le.j.and.j.le.nsy2) then
pulse=1.0e4*exp(-(t-t0)**2/(pt)**2)*exp(real(-(j-j0)**2/jw**2))
!pulse=1.0*exp(-(nt-20.0)**2/2.0*2)
else
pulse=0.0
end if
return
end function
!Hx chuzhi
! 给磁场加源
function source1(i,j)
implicit none
integer i,j
real*8 source1
if(wave=='sinee') then
source1=cose(i,j)*0.
else
if(wave=='pulse') then
source1=pulse1(i,j)
else if(wave/='sine'.or.wave=='pulse') then
write(*,*) ' 你在进行波源的选择时有错误，请核查'
end if
end if
return
end function
function cose(i,j)
implicit none
integer i,j
real*8 cose
if(nsx1.le.i.and.i.le.nsx2.and.nsy1.le.j.and.j.le.nsy2) then
cose=1.0e3*cos(w*t)
else
cose=0.0
end if
return
end function
function pulse1(i,j)
implicit none
integer i,j
real*8 t0,pt
real*8 pulse1
t0=21.0*dt
pt=2.0*dt
if(nsx1.le.i.and.i.le.nsx2.and.nsy1.le.j.and.j.le.nsy2) then
pulse1=1.0e4*exp(-(t-t0)**2/(pt)**2)
!pulse=1.0*exp(-(nt-20.0)**2/2.0*2)
else
pulse1=0.0
end if
return
end function
subroutine eee
implicit none
integer i,j
integer nxmin,nxmax,nymin,nymax
integer nx1,nx2,ny1,ny2
nx1=1
nx2=10
ny1=1
ny2=1
szhong=0.0
!进入的
do i=0,nx
do j=0,ny
if(i==nsx1-nx1.and.(nsy1-ny1.lt.j.and.j.lt.nsy2+ny2)) then
sru=sru+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*dt
!write(*,*) 'sru1',ez(i,j),hy(i,j),hy(i-1,j), sru
else
end if
if(i==nsx1-nx1.and.(j==nsy1-ny1.or.j==nsy2+ny2)) then
sru=sru+0.25*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*dt
end if
if(i==nsx1+nx2.and.(nsy1-ny1.lt.j.and.j.lt.nsy2+ny2)) then
sru=sru+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
!write(*,*) 'sru2'
else
end if
if(i==nsx1+nx2.and.(j==nsy1-ny1.or.j==nsy2+ny2)) then
sru=sru+0.25*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
end if
if(j==nsy1-ny1.and.(nsx1-nx1.lt.i.and.i.lt.nsx1+nx2)) then
sru=sru+0.5*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*dt
!write(*,*) 'sru3'
else
end if
if(j==nsy1-ny1.and.(nsx1-nx1==i.or.i==nsx1+nx2)) then
sru=sru+0.25*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*dt
!write(*,*) 'sru3'
else
end if
if(j==nsy2+ny2.and.(nsx1-nx1.lt.i.and.i.lt.nsx1+nx2)) then
sru=sru+0.5*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*(-1.)*dt
!write(*,*) 'sru4'
else
end if
if(j==nsy2+ny2.and.(nsx1-nx1==i.or.i==nsx1+nx2)) then
sru=sru+0.25*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*(-1.)*dt
!write(*,*) 'sru4'
else
end if
!外面层出
if(i==1.and.(1.lt.j.and.j.lt.ny-1)) then
schu=schu+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
!write(*,*) 'schu1',ez(i,j),hy(i,j),hy(i-1,j), schu
else
end if
if(i==1.and.(1==j.or.j==ny-1)) then
schu=schu+0.25*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*(-1.)*dt
!write(*,*) 'schu1',ez(i,j),hy(i,j),hy(i-1,j), schu
else
end if
if(i==nx-1.and.(1.lt.j.and.j.lt.ny-1)) then
schu=schu+0.5*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*dt
!write(*,*) 'schu2'
else
end if
if(i==nx-1.and.(1==j.or.j==ny-1)) then
schu=schu+0.25*(-1.)*ez(i,j)*(hy(i,j)+hy(i-1,j))*dy*dt
!write(*,*) 'schu2'
else
end if
if(j==1.and.(1.lt.i.and.i.lt.nx-1)) then
schu=schu+0.5*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*(-1.)*dt
!write(*,*) 'schu3'
else
end if
if(j==1.and.(1==i.or.i==nx-1)) then
schu=schu+0.25*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*(-1.)*dt
!write(*,*) 'schu3'
else
end if
if(j==ny-1.and.(1.lt.i.and.i.lt.nx-1)) then
schu=schu+0.5*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*dt
!write(*,*) 'schu4'
else
end if
if(j==ny-1.and.(1==i.or.i==nx-1)) then
schu=schu+0.25*ez(i,j)*(hx(i,j)+hx(i,j-1))*dx*dt
!write(*,*) 'schu4'
else
end if
!内部存储的
!N1部分
goto 1000
if((1.lt.i.and.i.lt.(nsx1-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(*,*) 'szhong1',ez(i,j),hy(i,j),hx(i,j), szhong
else
end if
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -