用于计算矩形波导本征模式的小程序.txt

用于计算矩形波导本征模式的小程序,直接运行就可以得到

!  WaveGuideMode.f90 
!
!  FUNCTIONS:
!	WaveGuideMode      - Entry point of console application.
!
!	Example of displaying 'Hello World' at execution time.
!

!****************************************************************************
!
!  PROGRAM: WaveGuideMode
!
!  PURPOSE:  Entry point for 'Hello World' sample console application.
!
!****************************************************************************

	program WaveGuideMode

	implicit none

	
	integer :: i, j, iid=1, maxnum=10, cal_type
	integer, parameter :: m=20, n=20
	integer, dimension(1:(m+1)*(n+1)) :: m_ind, n_ind

	real :: miu=1, epsilon=1, a=0.07214, b=0.03404, f_0=6e9, k_0, omega
	real, dimension(1:(m+1)*(n+1)) :: alpha, beta, v_p, v_g, lambda_g
	real, dimension(1:(m+1)*(n+1)) :: k_c, lambda_c, f_c
	real, parameter :: PI=3.1415926, c=3e8

	print *, "Input waveguide size, a and b"
	read *, a, b
	print *, "Input miu_r and epsilon_r"
	read *, miu, epsilon
	print *, "Choose calculate type, 1 for lambda_c and 2 for beta, 3 for all"
	read *, cal_type
	if (cal_type==2 .or. cal_type==3) then
		print *, "Input working frequency (Hz)"
		read *, f_0
	end if		
	print *, "Input maxnumber of modes to calculate:"
	read "(I2)", maxnum

	do i=0,m
		do j=0,n
			k_c(iid)=sqrt((i*PI/a)**2+(j*PI/b)**2)
			lambda_c(iid)=2/sqrt((i/a)**2+(j/b)**2)
			f_c(iid)=c/sqrt(miu*epsilon)/lambda_c(iid)
			m_ind(iid)=i
			n_ind(iid)=j
			!print *, m_ind(iid), n_ind(iid), k_c(iid), lambda_c(iid), f_c(iid)
			iid=iid+1
		end do
	end do

	do j=(m+1)*(n+1)-1, 1, -1
		do i=1, j
			if (lambda_c(i)<lambda_c(i+1)) then
				call Swap(lambda_c(i), lambda_c(i+1))
				call Swap(k_c(i), k_c(i+1))
				call Swap(f_c(i), f_c(i+1))
				call Swap(m_ind(i), m_ind(i+1))
				call Swap(n_ind(i), n_ind(i+1))
			end if
		end do
	end do
	
	omega=2*PI*f_0
	k_0=2*PI*f_0/c*sqrt(miu*epsilon)
	do i=1, (m+1)*(n+1)
		if (k_c(i) < k_0) then
			alpha(i)=0
			beta(i)=sqrt(k_0**2-k_c(i)**2)
			v_p(i)=omega/beta(i)
			v_g(i)=(c/sqrt(miu*epsilon))**2/v_p(i)
			lambda_g(i)=v_p(i)/f_0
		else
			alpha(i)=sqrt(k_c(i)**2-k_0**2)
			beta(i)=0
			v_p(i)=0
			v_g(i)=0
			lambda_g(i)=0
		end if
	end do

	if (cal_type==1 .or. cal_type==3) then
		print *
		print *, "*********************"
		print *, "***Eigenvalue Mode***"
		print *, "*********************"
		print *
		print *, "Parameters:"
		print "(3X, A, EN10.2, 3X, A, EN10.2)", "a =", a, "b =", b
		print "(3X, A, F5.2, 4X, A, F5.2)", "miu_r =", miu, "epsilon_r =", epsilon
		print *
		print *, "Results:"
		print *, " m  n        k_c          lambda_c          f_c"
		print *, " -  -    ------------   ------------    -----------"
		print "(2I3, 3X, EN12.3, 3X, EN12.3, 3X, EN12.3)", &
			  (m_ind(iid), n_ind(iid), k_c(iid), lambda_c(iid), f_c(iid), iid=1, maxnum)
	end if
	if (cal_type==2 .or. cal_type==3) then
		print *
		print *, "**************************"
		print *, "***Given Frequency Mode***"
		print *, "**************************"
		print *
		print *, "Parameters:"
		print "(3X, A, EN10.2, 3X, A, EN10.2)", "a =", a, "b =", b
		print "(3X, A, F5.2, 4X, A, F5.2)", "miu_r =", miu, "epsilon_r =", epsilon
		print "(3X, A, EN12.2, 3X, A, EN12.2, 3X, A, EN12.2, 3X, A, EN12.2)", &
		      "f_0", f_0, "k_0", k_0, "v_0", c/sqrt(miu*epsilon), "lambda_0", 2*PI/k_0
		print *
		print *, "Results:"
		print *, " m  n        alpha         beta          v_p           v_g         lambda_g"
		print *, " -  -     -----------  -----------   -----------   -----------    ----------" 
  		print "(2I3, 1X, 5EN14.3)", &
		      (m_ind(iid), n_ind(iid), alpha(iid), beta(iid), &
			   v_p(iid), v_g(iid), lambda_g(iid), iid=1, maxnum)
	end if
	end program WaveGuideMode

	subroutine Swap(a, b)
		implicit none
		real, intent(inout) :: a, b
		real temp
		temp=a
		a=b
		b=temp
	end subroutine Swap