/VF/parallel_of_all_3/ParallelPole.for
FORTRAN Legacy | 498 lines | 303 code | 125 blank | 70 comment | 16 complexity | 834be950f8d7745becbd09fc721d6caa MD5 | raw file
Possible License(s): GPL-3.0
- !!=========================================================================
- module init
-
- DOUBLE PRECISION Freq,OMega,Mu1,Mu2,Mu0,Mur1,Mur2,Pi
- DOUBLE PRECISION Epsilonr1,Epsilonr2,Epsilon0,Zn1,Zn2,d
- DOUBLE PRECISION tanDelta1, tanDelta2, Error,lambda
- COMPLEX*16 Epsilon1, Epsilon2,Kp, Kz1,Kz2,j, Diff, K1,K2,DDa, DDb
- INTEGER Iteration
-
- end module init
-
- !!=========================================================================
- subroutine initValue
-
- use init
- implicit none
-
- j =(0.0d0,1.0d0)
- !Pi = 3.141592653589793d0
- Pi = acos(0.d0)*2
- Mu0 = 4*Pi*1d-7
- Mur1 = 1.00d0
- Mur2 = 1.00d0
- Mu1 = Mur1*Mu0
- Mu2 = Mur2*Mu0
- Epsilon0 = 8.854d-12
- Epsilonr1 = 1.000d0
- Epsilonr2 = 4.000d0
- tanDelta1 = 0.05d0
- tanDelta2 = 0.05d0
- ! Epsilon1 = Epsilon0*Epsilonr1*((1.0d0,0.0d0)-tanDelta1*j)
- ! Epsilon2 = Epsilon0*Epsilonr2*((1.0d0,0.0d0)-tanDelta2*j)
- Epsilon1 = Epsilon0*(Epsilonr1*(1.0d0,0.0d0)-1.d-6*j)
- Epsilon2 = Epsilon0*(Epsilonr2*(1.0d0,0.0d0)-1.d-4*j)
-
- Freq = 12.0e+9
- Lambda = 3.e8/Freq
- Omega = 2.0d0*Pi*Freq
- Diff = (Omega**2)*(Mu2*Epsilon2-Mu1*Epsilon1)
-
- K1=(Omega**2)*Mu1*Epsilon1
- K2=(Omega**2)*Mu2*Epsilon2
-
- ! Zn1=0.0007620d0
- ! Zn2=0.0000762d0
- zn1=Lambda/10
- zn2=10*Zn1
- d =ZN1+ZN2
-
-
- Error = 1e-4
- Iteration = 100 !1e+4
-
-
- end subroutine initValue
-
- !!=========================================================================
-
- module Parallel_POLE
-
- !use init_module
- use init
- implicit none
-
- contains
-
- !!!***********************************************************
- FUNCTION TheMerge(A,B)
- DOUBLE PRECISION A(:),B(:)
- DOUBLE PRECISION TheMerge(1:size(A)+size(B))
-
- TheMerge(1:size(A))=A
- TheMerge(size(A)+1:size(A)+size(B))=B
- end FUNCTION TheMerge
-
- !!!***********************************************************
- FUNCTION TheSort (A)
- DOUBLE PRECISION A(:)
- DOUBLE PRECISION TheSort(1:size(A))
- INTEGER sizeA, nn,mm
- DOUBLE PRECISION tmp
-
- sizeA = size (A)
- TheSort(1:sizeA)=A
-
- do nn=1, sizeA
- do mm=nn+1, sizeA
- if (TheSort(nn) > TheSort(mm)) then
- tmp=TheSort(nn)
- TheSort(nn)=TheSort(mm)
- TheSort(mm)=tmp
- end if
- end do
- end do
- end FUNCTION TheSort
-
- !!!***********************************************************
-
- subroutine TEtheDb(theMu1, theMu2,theEpsilon1,theEpsilon2,theZn1,
- & theZn2, theKz2, theFreq,thedb,thedfdr)
-
-
- DOUBLE PRECISION theMu1, theMu2,thePi, theFreq, theOmega,
- & theZn1, theZn2
- complex*16,intent(inout) :: thedb,thedfdr
-
- COMPLEX*16 theEpsilon1, theEpsilon2,theKz2,theKZ1,theDiff
- & ,tan_kz1,tan_kz2
- thePi = PI
- !theOmega = 2.0d0*thePi*theFreq
- theOmega = Omega
- theDiff = diff
- theKZ1 = sqrt(theKz2**2-theDiff)
- tan_kz1 = tan(thekz1*theZn2)
- tan_kz2 = tan(thekz2*theZn1)
-
- !TEtheDb = tan_kz1*theKz2+tan_kz2*theKz1
-
- thedb = tan_kz1*theKz2+tan_kz2*theKz1
- thedfdr = tan_kz2*thekz2/thekz1+thekz1*(1+tan_kz2*tan_kz2)*theZn1
- &+tan_kz1+(1+tan_kz1*tan_kz1)*thezn2*thekz2*thekz2/thekz1
-
-
-
- end subroutine TEtheDb
-
- !!!***********************************************************
-
- COMPLEX*16 FUNCTION TEtheSearchRoot(tMu1, tMu2, tEpsilon1
- & , tEpsilon2, tZn1,tZn2, tFreq, start_val, end_val, iter, err)
-
- DOUBLE PRECISION tMu1, tMu2, tZn1, tZn2 , InitGuess
- & tPi, tFreq, tOmega, err, mag, magD,
- & start_val, end_val,tKz2r, tKz2i,InitGuess
- & ,startVal,endVal
-
- COMPLEX*16 tEpsilon1, tEpsilon2, tkz1,tKz2, tDiff,
- & tDb, tDfdr
- INTEGER iter, n
-
- !!! tzn1=d2,tzn2=d1 !*** Important ***
- !theDiff = theOmega**2*(theMu2*theEpsilon2-theMu1*theEpsilon1)
- !!!tdb=FU, tdfdr=DFU
-
- startVal=start_val
- endval =end_val
-
- tdiff=Diff
- n=0
- InitGuess= (startVal+endval)/2.0d0
-
- tKz2r = InitGuess
-
- tKz2i =0.0d0
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- tkz1 = sqrt(tKz2**2-tDiff)
-
- CALL TEtheDb (tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDb,tDfdr)
-
- !mag= sqrt(REAL(tDb)**2+AIMAG(tDb)**2)
- mag= abs(tDb)
- magD= REAL(tDfdr)**2-AIMAG(tDfdr)**2
-
- do while ( (n<=iter) .AND. (mag>err*1.e-4*abs(tan(tkz1*tZn2)*tkz2)) )
-
- if(REAL(tKz2)<startVal) then
-
- endval=InitGuess
- InitGuess= (startVal+endval)/2.0d0
-
- tKz2r = InitGuess
- tKz2i =0.0d0
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- tkz1 = sqrt(tKz2**2-tDiff)
-
- CALL TEtheDb (tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDb,tDfdr)
-
- mag= abs(tDb)
- magD= REAL(tDfdr)**2-AIMAG(tDfdr)**2
- n=n+1
-
- else if(REAL(tKz2)>endval) then
-
- startVal=InitGuess
- InitGuess= (startVal+endval)/2.0d0
- tKz2r = InitGuess
- tKz2i =0.0d0
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- tkz1 = sqrt(tKz2**2-tDiff)
-
- CALL TEtheDb (tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDb,tDfdr)
-
- mag= abs(tDb)
- magD= REAL(tDfdr)**2-AIMAG(tDfdr)**2
-
- n=n+1
-
- else
-
- tKz2r=tKz2r-(REAL(tDb)*REAL(tDfdr)-j*AIMAG(tDb)*AIMAG(tDfdr))/magD
- tKz2i=tKz2i-(j*REAL(tDb)*AIMAG(tDfdr)+AIMAG(tDb)*REAL(tDfdr))/magD
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- !
- !evaluate again
- CALL TEtheDb (tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDb,tDfdr)
-
- !mag= sqrt(REAL(tDb)**2+AIMAG(tDb)**2)
- mag= abs(tDb)
- magD= REAL(tDfdr)**2-AIMAG(tDfdr)**2
- n=n+1
-
- end if ! if-else-if
- end do ! while
-
- TEtheSearchRoot=tKz2
-
- end function TEtheSearchRoot
-
- !!!***********************************************************
- Subroutine TMtheDa(theMu1, theMu2, theEpsilon1, theEpsilon2
- & ,theZn1,theZn2, theKz2, theFreq,theda,thedfdr)
-
- DOUBLE PRECISION theMu1, theMu2,
- & thePi, theFreq, theOmega,theZn1, theZn2
- complex*16,intent(inout)::theda,thedfdr
- COMPLEX*16 theEpsilon1, theEpsilon2, theKz2, theDiff,
- & thekz1,tan_kz1,tan_kz2,Er
- thePi = PI
- !theOmega = 2.0d0*thePi*theFreq
- theOmega = Omega
- !theDiff = theOmega**2*(theMu2*theEpsilon2-theMu1*theEpsilon1)
- theDiff = Diff
- theKZ1 = sqrt(theKz2**2-theDiff)
- tan_kz1 = tan(thekz1*theZn2)
- tan_kz2 = tan(thekz2*theZn1)
- Er = theEpsilon2/theEpsilon1
-
- theda = Er*thekz1*tan_kz1+thekz2*tan_kz2
- thedfdr= Er*tan_kz1*thekz2/thekz1+Er*(1+tan_kz1**2)*thekz2*thezn2
- & +tan_kz2+(1+tan_kz2**2)*thekz2*thezn1
-
-
- end Subroutine TMtheDa
-
-
- !!!***********************************************************
-
- COMPLEX*16 FUNCTION TMtheSearchRoot(tMu1, tMu2, tEpsilon1,
- & tEpsilon2, tZn1,tZn2, tFreq, start_Val, end_Val, iter, err)
-
- DOUBLE PRECISION tMu1, tMu2, tZn1, tZn2 , InitGuess
- & tPi, tFreq, tOmega, err, mag, magD,
- & start_Val, end_Val,
- & startVal, endVal,tKz2r, tKz2i,InitGuess
-
- COMPLEX*16:: tEpsilon1, tEpsilon2, tKz2,tkz1,tDiff,
- & tDa, tMdfdr
- INTEGER iter, n
-
- tdiff=Diff
- startVal=start_val
- endval =end_val
- n=0
- InitGuess= (startVal+endVal)/2.0d0
-
- tKz2r = InitGuess
-
- tKz2i =0.0d0
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- tkz1 = sqrt(tKz2**2-tDiff)
-
- CALL TMtheDa (tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDa,tMdfdr)
-
- !mag= sqrt(REAL(tDb)**2+AIMAG(tDb)**2)
- mag= abs(tDa)
- magD= REAL(tMdfdr)**2-AIMAG(tMdfdr)**2
-
- do while ( (n<=iter) .AND. (mag>err*1.e-4*abs(tan(tkz1*tZn2)*tkz2)) )
-
- if(REAL(tKz2)<startVal) then
-
- endVal=InitGuess
- InitGuess= (startVal+endVal)/2.0d0
- tKz2r = InitGuess
- tKz2i =0.0d0
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- tkz1 = sqrt(tKz2**2-tDiff)
-
- CALL TMtheDa(tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDa,tMdfdr)
-
- mag= abs(tDa)
- magD= REAL(tMdfdr)**2-AIMAG(tMdfdr)**2
- n=n+1
-
- else if(REAL(tKz2)>endVal) then
-
- startVal=InitGuess
- InitGuess= (startVal+endVal)/2.0d0
- tKz2r = InitGuess
- tKz2i =0.0d0
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- tkz1 = sqrt(tKz2**2-tDiff)
-
- CALL TMtheDa (tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDa,tMdfdr)
-
- mag= abs(tDa)
- magD= REAL(tMdfdr)**2-AIMAG(tMdfdr)**2
-
- n=n+1
-
- else
-
- tKz2r=tKz2r-(REAL(tDa)*REAL(tMdfdr)-j*AIMAG(tDa)
- & *AIMAG(tMdfdr))/magD
- tKz2i=tKz2i-(j*REAL(tDa)*AIMAG(tMdfdr)+AIMAG(tDa)
- & *REAL(tMdfdr))/magD
- tKz2 =tKz2r*(1.0d0,0.0d0)+tKz2i*(0.0d0,1.0d0)
- !
- !evaluate again
- CALL TMtheDa (tMu1,tMu2,tEpsilon1,tEpsilon2,
- & tZn1,tZn2,tKz2,tFreq,tDa,tMdfdr)
-
- !mag= sqrt(REAL(tDa)**2+AIMAG(tDa)**2)
- mag= abs(tDa)
- magD= REAL(tMdfdr)**2-AIMAG(tMdfdr)**2
- n=n+1
-
- end if
- end do
-
- TMtheSearchRoot=tKz2
-
- end function TMtheSearchRoot
- !!!***********************************************************
-
- FUNCTION foundKpforTE(N)
-
- ! For calculating the poles of TE mode for Horizontal Dipole in stripline
-
- !use init
- ***********************************
- INTEGER N,i,ii
- DOUBLE PRECISION tmag, testStart, testEnd
-
- COMPLEX*16 Db, Dfdr,foundKpforTE(1:N), Kp1,Kp2
- DOUBLE PRECISION, ALLOCATABLE :: Asys1(:), Asys2(:),
- & Asys(:), SortedAys(:)
-
-
- ***********************************
- call initValue
- *********************************
- allocate (Asys1(1:(N+5)))
- allocate (Asys2(1:(N+5)))
- allocate (Asys(1:2*(N+5)))
- allocate (SortedAys(1:2*(N+5)))
- ****************************************
- Do i=1, N+5
- Asys1(i)=(i-0.5d0)*Pi/Zn1
- Asys2(i)=sqrt(((i-0.5d0)*Pi/Zn2)**2+Diff)
- end Do
- Asys=TheMerge(Asys1,Asys2)
- SortedAys=TheSort(Asys)
- c print *, SortedAys
-
- open (2, FILE="./results/theResultTE.txt",STATUS='replace')
-
- ii=1
- Do i=1,2*(N+5)-1
- testStart= SortedAys(i)
- testEnd = SortedAys(i+1)
-
- Kz2= TEtheSearchRoot(Mu1, Mu2, Epsilon1, Epsilon2, Zn1,
- & Zn2, Freq, testStart, testEnd, Iteration, Error)
-
- Call TEtheDb(Mu1,Mu2,Epsilon1,Epsilon2,
- & Zn1,Zn2,Kz2,Freq,db,Dfdr)
-
- tmag= sqrt(REAL(Db)**2+AIMAG(Db)**2)
- Kz1= sqrt(Kz2**2-Diff)
- Kp1=sqrt(Omega**2*Mu1*Epsilon1-Kz1**2)
- Kp2=sqrt(Omega**2*Mu2*Epsilon2-Kz2**2)
-
- write (2,"(A20)" ) ""
- write (2,"('the start point=',F25.15)") (testStart)
- write (2,"('the end point=',F25.15)") (testEnd)
- write (2,"('Kz2= ', F24.15, F24.15,'j')") (REAL(Kz2),AIMAG(Kz2))
- write (2,"('Kz1= ', F24.15, F24.15,'j')") (REAL(Kz1),AIMAG(Kz1))
- write (2,"('Db =', F25.15, F25.15,'j')") (REAL(Db),AIMAG(Db))
- write (2,"('the error=',F25.15)") (tmag)
-
- write (2,"('Kp1(derived from K_z1)=', F22.15, F22.15,'j')")
- & (REAL(Kp1),AIMAG(Kp1))
- write (2,"('Kp2(derived from K_z2)=', F22.15, F22.15,'j')")
- & (REAL(Kp2),AIMAG(Kp2))
-
- if (tmag<Error) then
- write (2, "(A10)") "ROOT FOUND"
- if (ii<=N) then
- foundKpforTE(ii)=Kp1
- ii=ii+1
- end if
- end if
-
- if (tmag>=Error) then
- write (2, "(A19)") "ROOT NOT FOUND !!!"
- end if
-
- end Do
- end FUNCTION foundKpforTE
-
- !!!***********************************************************
-
- FUNCTION foundKpforTM(N)
- ! For calculating the poles of TM mode for Horizontal Dipole in stripline
-
- !use init
- ***********************************
- INTEGER N,i,ii
- DOUBLE PRECISION tmag, testStart, testEnd
-
- COMPLEX*16 Da, Dfdr,foundKpforTM(1:N), Kp1,Kp2
- DOUBLE PRECISION, ALLOCATABLE :: Asys1(:), Asys2(:),
- & Asys(:), SortedAys(:)
-
-
- ***********************************
- call initValue
- *********************************
- allocate (Asys1(1:(N+5)))
- allocate (Asys2(1:(N+5)))
- allocate (Asys(1:2*(N+5)))
- allocate (SortedAys(1:2*(N+5)))
- ****************************************
- Do i=1, N+5
- Asys1(i)=(i-0.5)*Pi/Zn1
- Asys2(i)=real(sqrt(((i-0.5)*Pi/Zn2)**2+Diff))
- end Do
- Asys=TheMerge(Asys1,Asys2)
- SortedAys=TheSort(Asys)
- c print *, SortedAys
-
- open (2, FILE="./results/theResultTM.txt",STATUS='replace')
-
- ii=1
- Do i=1,2*(N+5)-1
- testStart= SortedAys(i)
- testEnd = SortedAys(i+1)
-
- Kz2= TMtheSearchRoot(Mu1, Mu2, Epsilon1, Epsilon2, Zn1,
- & Zn2, Freq, testStart, testEnd, Iteration, Error)
-
- Call TMtheDa(Mu1,Mu2,Epsilon1,Epsilon2,
- & Zn1,Zn2,Kz2,Freq,Da,Dfdr)
-
- tmag= sqrt(REAL(Da)**2+AIMAG(Da)**2)
- Kz1= sqrt(Kz2**2-Diff)
- Kp1=sqrt(Omega**2*Mu1*Epsilon1-Kz1**2)
- Kp2=sqrt(Omega**2*Mu2*Epsilon2-Kz2**2)
-
- write (2,"(A20)" ) ""
- write (2,"('the start point=',F25.15)") (testStart)
- write (2,"('the end point=',F25.15)") (testEnd)
- write (2,"('Kz2= ', F24.15, F24.15,'j')") (REAL(Kz2),AIMAG(Kz2))
- write (2,"('Kz1= ', F24.15, F24.15,'j')") (REAL(Kz1),AIMAG(Kz1))
- write (2,"('Da =', F25.15, F25.15,'j')") (REAL(Da),AIMAG(Da))
- write (2,"('the error=',F25.15)") (tmag)
-
- write (2,"('Kp1(derived from Kz1)=', F22.15, F22.15,'j')")
- & (REAL(Kp1),AIMAG(Kp1))
- write (2,"('Kp2(derived from Kz2)=', F22.15, F22.15,'j')")
- & (REAL(Kp2),AIMAG(Kp2))
-
- if (tmag<Error) then
- write (2, "(A10)") "ROOT FOUND"
- if (ii<=N) then
- foundKpforTM(ii)=Kp1
- ii=ii+1
- end if
- end if
-
- if (tmag>=Error) then
- write (2, "(A19)") "ROOT NOT FOUND !!!"
- end if
-
- end Do
- end FUNCTION foundKpforTM
-
-
-
- end module Parallel_POLE