📄 3dpml.f90

📁 计算同轴谐振腔的一个小软件
💻 F90
字号:
!A spot at (nx/2,ny/2,0) with gauss impulse as a stimulus
!zero coordinence at left lower
!trip_d 
!gedeny pml as abc
!f=900M,dif_t=2.5E-11,delx=0.015,dely=0.015,delz=0.015
PROGRAM gedney_pml
 IMPLICIT NONE
 REAL, PARAMETER:: PI=3.14159
 REAL, PARAMETER:: miu=PI*4.E-7  																			 
 REAL, PARAMETER:: ipu=1/(4.*PI*9.E9)
 REAL, PARAMETER:: dif_t=11.11e-12,delx=0.005,dely=0.005,delz=0.01
 REAL, PARAMETER:: freq=1E9
 REAL, PARAMETER:: T=16*dif_t
 REAL, PARAMETER:: V0=10.
 INTEGER, PARAMETER:: vol=3e8
 INTEGER, PARAMETER:: Nt=700
 INTEGER, PARAMETER:: nx=6,ny=6,nz=6
 
 
!define variables
 DOUBLE PRECISION V
!WHEN X LT -10+10
 DOUBLE PRECISION dx_e(-9:nx+9)
 DOUBLE PRECISION dx_h(-9:nx+9)
 DOUBLE PRECISION dx_bi_e(-10:nx+9)
 DOUBLE PRECISION dx_bi_h(-10:nx+9)
 DOUBLE PRECISION dy_e(-9:ny+9)
 DOUBLE PRECISION dy_h(-9:ny+9)
 DOUBLE PRECISION dy_bi_e(-10:ny+9)
 DOUBLE PRECISION dy_bi_h(-10:ny+9)
 DOUBLE PRECISION dz_e(-9:nz+9)
 DOUBLE PRECISION dz_h(-9:nz+9)
 DOUBLE PRECISION dz_bi_e(-10:nz+9)
 DOUBLE PRECISION dz_bi_h(-10:nz+9)
 !(orient,x,y,z)
 DOUBLE PRECISION Dx(-10:nx+9,-10:ny+10,-10:nz+10)
 DOUBLE PRECISION Dy(-10:nx+10,-10:ny+9,-10:nz+10)
 DOUBLE PRECISION Dz(-10:nx+10,-10:ny+10,-10:nz+9)
 DOUBLE PRECISION Bx(-10:nx+10,-10:ny+9,-10:nz+9)
 DOUBLE PRECISION By(-10:nx+9,-10:ny+10,-10:nz+9)
 DOUBLE PRECISION Bz(-10:nx+9,-10:ny+9,-10:nz+10)

 DOUBLE PRECISION Ex(-10:nx+9,-10:ny+10,-10:nz+10)
 DOUBLE PRECISION Ey(-10:nx+10,-10:ny+9,-10:nz+10)
 DOUBLE PRECISION Ez(-10:nx+10,-10:ny+10,-10:nz+9)
 DOUBLE PRECISION Hx(-10:nx+10,-10:ny+9,-10:nz+9)
 DOUBLE PRECISION Hy(-10:nx+9,-10:ny+10,-10:nz+9)
 DOUBLE PRECISION Hz(-10:nx+9,-10:ny+9,-10:nz+10)


 DOUBLE PRECISION TEMP
 INTEGER x,y,z
 INTEGER N		            !N OF REIT
 INTEGER I                  !I AS INTERMIT
 REAL *8 :: dmex=1/(30*pi*delx),dmey=1/(30*pi*dely),dmez=1/(30*pi*delz)	 !max diele
 REAL *8  dmhx,dmhy,dmhz
 dmhx=dmex*miu/ipu;dmhy=dmey*miu/ipu;dmhz=dmez*miu/ipu
 V=100
 Ex=0.; Ey=0.; Ez=0.; Hx=0.; Hy=0.; Hz=0.
 Dx=0.; Dy=0.; Dz=0.; Bx=0.; By=0.; Bz=0.
 OPEN(UNIT=100,FILE='Exn.DAT')
 OPEN(UNIT=201,FILE='Ex.DAT')
 OPEN(UNIT=202,FILE='Ey.DAT')
 OPEN(UNIT=203,FILE='Ez.DAT')
 OPEN(UNIT=204,FILE='Hx.DAT')
 OPEN(UNIT=205,FILE='Hy.DAT')
 OPEN(UNIT=206,FILE='Hz.DAT')
 DO x=-9,nx+9
    if(x.LE.0)then
    dx_e(x)=dmex*((x)/10.)**4 
    dx_h(x)=dmhx*((x)/10.)**4 
	else if(x.LE.nx-1) then
	dx_e(x)=0.
	dx_h(x)=0.
	else
	dx_e(x)=dmex*((x-nx)/10.)**4 
    dx_h(x)=dmhx*((x-nx)/10.)**4 
	endif				   
 ENDDO
 DO x=-10,nx+9
    if(x.LE.-1)then
	dx_bi_e(x)=dmex*((x+0.5)/10.)**4 
	dx_bi_h(x)=dmhx*((x+0.5)/10.)**4 
    elseif(x.LE.nx-1)then
	dx_bi_e(x)=0.
	dx_bi_h(x)=0.
	else
    dx_bi_e(x)=dmex*((x-nx+0.5)/10.)**4 
	dx_bi_h(x)=dmhx*((x-nx+0.5)/10.)**4 
	endif
 ENDDO
 DO y=-9,ny+9
    if(y.LE.0.)then
    dy_e(y)=dmey*((y)/10.)**4 
    dy_h(y)=dmhy*((y)/10.)**4 
	elseif(y.LE.ny-1)then
	dy_e(y)=0.
	dy_h(y)=0.
	else
	dy_e(y)=dmey*((y-ny)/10.)**4 
	dy_h(y)=dmhy*((y-ny)/10.)**4 
	endif
 ENDDO
 DO y=-10,ny+9
    if(y.LE.-1.)then
    dy_bi_e(y)=dmey*((y+0.5)/10.)**4 
	dy_bi_h(y)=dmhy*((y+0.5)/10.)**4 
	elseif(y.LE.ny-1)then
	dy_bi_e(y)=0.
	dy_bi_h(y)=0.
	else
	dy_bi_e(y)=dmey*((y-ny+0.5)/10.)**4 
	dy_bi_h(y)=dmhy*((y-ny+0.5)/10.)**4 
	endif
 ENDDO
 DO z=-10+1,nz+9
    if(z.LE.0.)then
    dz_e(z)=dmez*((z)/10.)**4 
    dz_h(z)=dmhz*((z)/10.)**4 
	elseif(z.LE.nz-1)then
	dz_e(z)=0.
	dz_e(z)=0.
 	else
	dz_e(z)=dmez*((z-nz)/10.)**4 
    dz_h(z)=dmhz*((z-nz)/10.)**4 
    endif
 ENDDO
 DO z=-10,nz+9
    if(z.LE.-1)then
    dz_bi_e(z)=dmez*((z+0.5)/10.)**4 
	dz_bi_h(z)=dmhz*((z+0.5)/10.)**4 
	elseif(z.LE.nz-1)then
	dz_bi_e(z)=0.
	dz_bi_e(z)=0.
	else
	dz_bi_e(z)=dmez*((z-nz+0.5)/10.)**4 
	dz_bi_h(z)=dmhz*((z-nz+0.5)/10.)**4 
    endif
 ENDDO
 DO N=1,Nt
 write(*,*)N
!galvanize source
!  V=V0*SIN(2*PI*freq*N*dif_t)
  v=v0*EXP(-((N-1)*dif_t-3*T)**2/T**2)
  Ey(nx/2,ny/2,nz/2)=V/delz
!refurbish Hx at space domain
   DO x=-10+1,nx+9
     DO y=-10,ny+9
	   DO z=-10,nz+9
TEMP=Bx(x,y,z)
Bx(x,y,z)=Bx(x,y,z)*(2*miu-dz_bi_h(z)*dif_t)/(2*miu+dz_bi_h(z)*dif_t)+&
      2*dif_t/(2*miu+dif_t*dz_bi_h(z))*((Ez(x,y,z)-Ez(x,y+1,z))/dely+(Ey(x,y,z+1)-Ey(x,y,z))/delz)
Hx(x,y,z)=Hx(x,y,z)*(2*miu-dy_bi_h(y)*dif_t)/(2*miu+dy_bi_h(y)*dif_t)+&
	  Bx(x,y,z)*(2*miu+dx_h(x)*dif_t)/(2*miu+dy_bi_h(y)*dif_t)+&
	  TEMP*(dx_h(x)*dif_t-2*miu)/(2*miu+dy_bi_h(y)*dif_t)
       ENDDO
     ENDDO
   ENDDO
!refurbish Hy at space domain
   DO x=-10,nx+9
     DO y=-10+1,ny+9
	   DO z=-10,nz+9
TEMP=By(x,y,z)		    
By(x,y,z)=By(x,y,z)*(2*miu-dz_bi_h(z)*dif_t)/(2*miu+dz_bi_h(z)*dif_t)+&
      2*dif_t/(2*miu+dif_t*dz_bi_h(z))*((Ex(x,y,z)-Ex(x,y,z+1))/delz+(Ez(x+1,y,z)-Ez(x,y,z))/delx)
Hy(x,y,z)=Hy(x,y,z)*(2*miu-dx_bi_h(x)*dif_t)/(2*miu+dx_bi_h(x)*dif_t)+&
	  By(x,y,z)*(2*miu+dy_h(y)*dif_t)/(2*miu+dx_bi_h(x)*dif_t)+&
      TEMP*(dy_h(y)*dif_t-2*miu)/(2*miu+dx_bi_h(x)*dif_t)
 	   ENDDO
	 ENDDO
   ENDDO
!refurbish Hz at space domain
   DO x=-10,nx+9
     DO y=-10,ny+9
	   DO z=-10+1,nz+9
TEMP=Bz(x,y,z)		    
Bz(x,y,z)=Bz(x,y,z)*(2*miu-dy_bi_h(y)*dif_t)/(2*miu+dy_bi_h(y)*dif_t)+&
      2*dif_t/(2*miu+dif_t*dy_bi_h(y))*((Ey(x,y,z)-Ey(x+1,y,z))/delx+(Ex(x,y+1,z)-Ex(x,y,z))/dely)
Hz(x,y,z)=Hz(x,y,z)*(2*miu-dx_bi_h(x)*dif_t)/(2*miu+dx_bi_h(x)*dif_t)+&
      Bz(x,y,z)*(2*miu+dz_h(z)*dif_t)/(2*miu+dx_bi_h(x)*dif_t)+&
	  TEMP*(dz_h(z)*dif_t-2*miu)/(2*miu+dx_bi_h(x)*dif_t)
	    ENDDO
     ENDDO
   ENDDO
!refurbish Ex at space domain
   DO x=-10,nx+9
     DO y=-10+1,ny+9
	   DO z=-10+1,nz+9
	  
TEMP=Dx(x,y,z)		    
Dx(x,y,z)=Dx(x,y,z)*(2*ipu-dz_e(z)*dif_t)/(2*ipu+dz_e(z)*dif_t)+&
      2*dif_t/(2*ipu+dif_t*dz_e(z))*((Hz(x,y,z)-Hz(x,y-1,z))/dely+(Hy(x,y,z-1)-Hy(x,y,z))/delz)
Ex(x,y,z)=Ex(x,y,z)*(2*ipu-dy_e(y)*dif_t)/(2*ipu+dy_e(y)*dif_t)+&
      Dx(x,y,z)*(2*ipu+dx_bi_e(x)*dif_t)/(2*ipu+dy_e(y)*dif_t)+&
      TEMP*(dx_bi_e(x)*dif_t-2*ipu)/(2*ipu+dy_e(y)*dif_t)
	  
		ENDDO
	  ENDDO
   ENDDO
!refurbish Ey at space domain
   DO x=-10+1,nx+9
     DO y=-10,ny+9
	   DO z=-10+1,nz+9
TEMP=Dy(x,y,z)		    
Dy(x,y,z)=Dy(x,y,z)*(2*ipu-dz_e(z)*dif_t)/(2*ipu+dz_e(z)*dif_t)+&
      2*dif_t/(2*ipu+dif_t*dz_e(z))*((Hx(x,y,z)-Hx(x,y,z-1))/delz+(Hz(x-1,y,z)-Hz(x,y,z))/delx)
Ey(x,y,z)=Ey(x,y,z)*(2*ipu-dx_e(x)*dif_t)/(2*ipu+dx_e(x)*dif_t)+&
      Dy(x,y,z)*(2*ipu+dy_bi_e(y)*dif_t)/(2*ipu+dx_e(x)*dif_t)+&
      TEMP*(dy_bi_e(y)*dif_t-2*ipu)/(2*ipu+dx_e(x)*dif_t)
	   ENDDO
	 ENDDO
   ENDDO
!refurbish Ez at space domain
   DO x=-10+1,nx+9
     DO y=-10+1,ny+9
	   DO z=-10,nz+9
TEMP=Dz(x,y,z)		    
Dz(x,y,z)=Dz(x,y,z)*(2*ipu-dy_e(y)*dif_t)/(2*ipu+dy_e(y)*dif_t)+&
      2*dif_t/(2*ipu+dif_t*dy_e(y))*((Hy(x,y,z)-Hy(x-1,y,z))/delx+(Hx(x,y-1,z)-Hx(x,y,z))/dely)
Ez(x,y,z)=Ez(x,y,z)*(2*ipu-dx_e(x)*dif_t)/(2*ipu+dx_e(x)*dif_t)+&
      Dz(x,y,z)*(2*ipu+dz_bi_e(z)*dif_t)/(2*ipu+dx_e(x)*dif_t)+&
	  TEMP*(dz_bi_e(z)*dif_t-2*ipu)/(2*ipu+dx_e(x)*dif_t)
		ENDDO
     ENDDO
  ENDDO
  
  WRITE(100,100) Ey(nx/2+2,ny/2,nz/2)
  IF(N.EQ.20) THEN	
 !**********************************need being modified to get a satisfying result
  DO I=0,nx
  WRITE(201,100) Ex(I,ny/2,0:nz)
  WRITE(202,100) Ey(I,ny/2,0:nz)
  WRITE(203,100) Ez(I,ny/2,0:nz)
  WRITE(204,100) Hx(I,ny/2,0:nz)
  WRITE(205,100) Hy(I,ny/2,0:nz)
  WRITE(206,100) Hz(I,ny/2,0:nz)
!**********************************need being modified,too
100 FORMAT(1X,81E16.5e4)

  ENDDO
  ENDIF 
 ENDDO
CLOSE(100)
CLOSE(200)
END
💿 文件大小 3 K
👤 上传用户 scorpioll
📂 所属分类数学计算
🏷️ 相关标签

#计算 #同轴 #谐振腔 #软件
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -