antenna.for

fortran code of dipole antenna

cccccccccccccccccccccccccccccccccccccccc     
      implicit real*8 (a-h,o-z)
	integer xs,ys,zs
	parameter(xs=60)
	parameter(ys=60)
	parameter(zs=60)
ccccccccccccccccccccccccccccccccccccccccc
      common /ia/hxi(0:xs,0:ys-1,0:zs-1)
	common /ib/hyi(0:xs-1,0:ys,0:zs-1)
	common /ic/hzi(0:xs-1,0:ys-1,0:zs)
	common /id/exi(0:xs-1,0:ys,0:zs)
	common /ie/eyi(0:xs,0:ys-1,0:zs)
	common /if/ezi(0:xs,0:ys,0:zs-1)
ccccccccccccccccccccccccccccccccccccccccc
	common /ha1/exiya1(0:xs-1,0:1,0:zs)
      common /hb1/exiyb1(0:xs-1,ys-1:ys,0:zs)
      common /hc1/exiza1(0:xs-1,0:ys,0:1)
	common /hd1/exizb1(0:xs-1,0:ys,zs-1:zs)
	common /ha2/exiya2(0:xs-1,0:1,0:zs)
	common /hb2/exiyb2(0:xs-1,0:ys-1,0:zs)
	common /hc2/exiza2(0:xs-1,0:ys,0:1)
	common /hd2/exizb2(0:xs-1,0:ys,zs-1:zs)
cccccccccccccccccccccccccccccccccccccccccccccc
      common /he1/eyiza1(0:xs,0:ys-1,0:1)
      common /hf1/eyizb1(0:xs,0:ys-1,zs-1:zs)
      common /hg1/eyixa1(0:1,0:ys-1,0:zs)
	common /hh1/eyixb1(xs-1:xs,0:ys-1,0:zs)
	common /he2/eyiza2(0:xs,0:ys-1,0:1)
	common /hf2/eyizb2(0:xs,0:ys-1,zs-1:zs)
	common /hg2/eyixa2(0:1,0:ys-1,0:zs)
	common /hh2/eyixb2(xs-1:xs,0:ys-1,0:zs)
ccccccccccccccccccccccccccccccccccccccccccccccc
      common /hi1/ezixa1(0:1,0:ys,0:zs-1)
      common /hj1/ezixb1(xs-1:xs,0:ys,0:zs-1)
      common /hk1/eziya1(0:xs,0:1,0:zs-1)
	common /hl1/eziyb1(0:xs,ys-1:ys,0:zs-1)
	common /hi2/ezixa2(0:1,0:ys,0:zs-1)
	common /hj2/ezixb2(xs-1:xs,0:ys,0:zs-1)
	common /hk2/eziya2(0:xs,0:1,0:zs-1)
	common /hl2/eziyb2(0:xs,ys-1:ys,0:zs-1)
ccccccccccccccccccccccccccccccccccccccccccccccc
	common /a/rh1,rh2,dx,dy,dz,ts,nt
	common /b/vz,pi
	common /c/ux1,ux2,ux3,ux4,uy1,uy2,uy3,uy4,uz1,uz2,uz3,uz4
cccccccccccccccccccccccccccccccccccccccccccccccccc

cccccccc Calcution ccccccccccccccccccccccc
	call shoki
	open(1,file='v-i.dat')

	do n=0,nt
	 
	call inc(n)

	hxn=hxi(30,30-1,30)-hxi(30,30,30)
	hyn=hyi(30,30,30)-hyi(30-1,30,30)
	cur=(hxn*dx+hyn*dy)
	write(1,*)n,vz,cur

	end do

	close(1)

	call furie

	stop
	end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine shoki

	implicit real*8 (a-h,o-z)

	common /a/rh1,rh2,dx,dy,dz,ts,nt
	common /b/vz,pi
	common /c/ux1,ux2,ux3,ux4,uy1,uy2,uy3,uy4,uz1,uz2,uz3,uz4

	dx=0.02
	dy=0.02
	dz=0.02

	nt=2000

	PI=4.*atan(1.0)
	c=2.998e8
	e0=8.854e-12
	u0=PI*4.e-7

	dxt=c/dx
	dyt=c/dy
	dzt=c/dz
	ts=1./sqrt(dxt**2+dyt**2+dzt**2)
	rh1=ts/u0
	rh2=ts/e0

	vvx=c*ts-dx
	wwx=c*ts+dx
	ux1=vvx/wwx
	ux2=2.*dx/wwx
	ux3=dx*c*c*ts*ts/2./dy/dy/wwx
	ux4=dx*c*c*ts*ts/2./dz/dz/wwx

	vvy=c*ts-dy
	wwy=c*ts+dy
	uy1=vvy/wwy
	uy2=2.*dy/wwy
	uy3=dy*c*c*ts*ts/2./dz/dz/wwy
	uy4=dy*c*c*ts*ts/2./dx/dx/wwy
	
	vvz=c*ts-dz
	wwz=c*ts+dz
	uz1=vvz/wwz
	uz2=2.*dz/wwz
	uz3=dz*c*c*ts*ts/2./dx/dx/wwz
	uz4=dz*c*c*ts*ts/2./dy/dy/wwz

      return
	end

ccccccccccc Forier Trans cccccccccccccccccccccc
      subroutine furie

	implicit real*8 (a-h,o-z)
	real*8 vz(0:2000),cur(0:2000)
	complex an1(0:500),an2(0:500),ans,ci

	common /a/rh1,rh2,dx,dy,dz,ts,nt
	common /b/vz,pi

	do i=0,500
		an1(i)=0
		an2(i)=0
	end do
	

	open(2,file='v-i.dat')
		do i=0,nt
			read (2,*),vz(i),cur(i)
		end do
	close(2)

	do n=0,nt
		if(n.lt.65) then
	     do i=0,500
			fre =1000e+6/500 * i
			omega = 2.*PI*fre
			ci=cmplx(0.,-1.)
			an2(i)=an2(i)
     &			   +real(vz(n))*cexp(ci*real(omega*real(n)*ts))
		 end do
		end if
		
	     do i=0,500
		 fre =1000e+6/500 * i
		 omega = 2.*PI*fre
		 ci=cmplx(0.,-1.)
		 an1(i)=an1(i)
     &	        +real(cur(n))*cexp(ci*real(omega*real(n)*ts))
     		 end do
		 
      end do
	
	open(3,file='z.dat')
		do i=0,500
			ans=an2(i)/an1(i)
			s=real(i)/500
			write(3,200)s,ans
		end do
	close(3)
200	format(1h,e12.5,2e20.5)
	
	return
	end
	
ccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine inc(n)
	
	implicit real*8 (a-h,o-z)
	integer xs,ys,zs
	parameter(xs=60)
	parameter(ys=60)
	parameter(zs=60)
	
	common /ia/hxi(0:xs,0:ys-1,0:zs-1)
	common /ib/hyi(0:xs-1,0:ys,0:zs-1)
	common /ic/hzi(0:xs-1,0:ys-1,0:zs)
	common /id/exi(0:xs-1,0:ys,0:zs)
	common /ie/eyi(0:xs,0:ys-1,0:zs)
	common /if/ezi(0:xs,0:ys,0:zs-1)
ccccccccccccccccccccccccccccccccccccccccc
	common /ha1/exiya1(0:xs-1,0:1,0:zs)
      common /hb1/exiyb1(0:xs-1,ys-1:ys,0:zs)
      common /hc1/exiza1(0:xs-1,0:ys,0:1)
	common /hd1/exizb1(0:xs-1,0:ys,zs-1:zs)
	common /ha2/exiya2(0:xs-1,0:1,0:zs)
	common /hb2/exiyb2(0:xs-1,0:ys-1,0:zs)
	common /hc2/exiza2(0:xs-1,0:ys,0:1)
	common /hd2/exizb2(0:xs-1,0:ys,zs-1:zs)
cccccccccccccccccccccccccccccccccccccccccccccc
      common /he1/eyiza1(0:xs,0:ys-1,0:1)
      common /hf1/eyizb1(0:xs,0:ys-1,zs-1:zs)
      common /hg1/eyixa1(0:1,0:ys-1,0:zs)
	common /hh1/eyixb1(xs-1:xs,0:ys-1,0:zs)
	common /he2/eyiza2(0:xs,0:ys-1,0:1)
	common /hf2/eyizb2(0:xs,0:ys-1,zs-1:zs)
	common /hg2/eyixa2(0:1,0:ys-1,0:zs)
	common /hh2/eyixb2(xs-1:xs,0:ys-1,0:zs)
ccccccccccccccccccccccccccccccccccccccccccccccc
      common /hi1/ezixa1(0:1,0:ys,0:zs-1)
      common /hj1/ezixb1(xs-1:xs,0:ys,0:zs-1)
      common /hk1/eziya1(0:xs,0:1,0:zs-1)
	common /hl1/eziyb1(0:xs,ys-1:ys,0:zs-1)
	common /hi2/ezixa2(0:1,0:ys,0:zs-1)
	common /hj2/ezixb2(xs-1:xs,0:ys,0:zs-1)
	common /hk2/eziya2(0:xs,0:1,0:zs-1)
	common /hl2/eziyb2(0:xs,ys-1:ys,0:zs-1)
ccccccccccccccccccccccccccccccccccccccccccccccc
	common /a/rh1,rh2,dx,dy,dz,ts,nt
	common /b/vz
	common /c/ux1,ux2,ux3,ux4,uy1,uy2,uy3,uy4,uz1,uz2,uz3,uz4
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	if (n.eq.0) then

	do i=0,xs-1
		do j=0,ys
			do k=0,zs
				exi(i,j,k)=0
			end do
		end do
	end do

	do i=0,xs
		do j=0,ys-1
			do k=0,zs
				eyi(i,j,k)=0
			end do
		end do
	end do

	do i=0,xs
		do j=0,ys
			do k=0,zs-1
				ezi(i,j,k)=0
			end do
		end do
	end do
	
	do i=0,xs
		do j=0,ys-1
			do k=0,zs-1
				hxi(i,j,k)=0
			end do
		end do
	end do

	do i=0,xs-1
		do j=0,ys
			do k=0,zs-1
				hyi(i,j,k)=0
			end do
		end do
	end do

	do i=0,xs-1
		do j=0,ys-1
			do k=0,zs
				hzi(i,j,k)=0
			end do
		end do
	end do

	end if

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	do k=-12,-1
	ezi(30,30,k+30)=0
	end do
	do k=1,12
	ezi(30,30,k+30)=0
	end do
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	alpha=1./64.
	if(n.lt.65) then
		vz=exp(-alpha*(real(32-n))**2)
	else
		vz=0
	end if
		ezi(30,30,30)=-vz/dz
cccccccccccc H_FEILD ccccccccccccccccccccccccccccc
	do i=0,xs
		do j=0,ys-1
			do k=0,zs-1
				aax=((eyi(i,j,k+1)-eyi(i,j,k))/dz)
     &				-((ezi(i,j+1,k)-ezi(i,j,k))/dy)
				hxi(i,j,k)=hxi(i,j,k)+rh1*aax
			end do
		end do
	end do
	
	do i=0,xs-1
		do j=0,ys
			do k=0,zs-1
				aay=((ezi(i+1,j,k)-ezi(i,j,k))/dx)
     &				-((exi(i,j,k+1)-exi(i,j,k))/dz)
				hyi(i,j,k)=hyi(i,j,k)+rh1*aay
			end do
		end do
	end do

	do i=0,xs-1
		do j=0,ys-1
			do k=0,zs
				aaz=((exi(i,j+1,k)-exi(i,j,k))/dy)
     &				-((eyi(i+1,j,k)-eyi(i,j,k))/dx)
				hzi(i,j,k)=hzi(i,j,k)+rh1*aaz
			end do
		end do
	end do
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	do i=0,xs-1
		do k=0,zs
			exiya2(i,0,k)=exiya1(i,0,k)
			exiya2(i,1,k)=exiya1(i,1,k)
			exiyb2(i,ys,k)=exiyb1(i,ys,k)
			exiyb2(i,ys-1,k)=exiyb1(i,ys-1,k)
			exiya1(i,0,k)=exi(i,0,k)
			exiya1(i,1,k)=exi(i,1,k)