Disclaimer: This software is provided "as is", without any form of support. Use at your own risk.
You can use any part of this code in your own work provided that this source is acknowledged.
PROGRAM TLM3D
*************************************************************
* *
* Transmission-line modelling of electromagnetic fields in *
* 3-dimensions using the symmetrical condensed node. *
* *
* Written by: *
* J L Herring *
* *
* E-mail: tlm@jlherring.uk *
* WWW: http://www.jlherring.uk/tlm/ *
* *
* History: *
* Date started : 10 Dec 1995 *
* Updated : 10 Dec 1995 (HP-UX) 1.00 *
* Current version : 30 Jun 1996 (HP-UX) 1.10 *
* - more efficient scattering *
* - two-sided internal boundaries *
* *
*************************************************************
IMPLICIT NONE
*---- constant declarations
CHARACTER*(*) VERSN
INTEGER INUNT,OUTUNT
INTEGER XMAX0,YMAX0,ZMAX0
INTEGER NMBCS0,NMNEX0,NMNOU0
INTEGER PYNX,PZNX,PXNY,PZNY,PYNZ,PXNZ
INTEGER PYPZ,PZPY,PZPX,PXPZ,PXPY,PYPX
REAL Z0
PARAMETER ( VERSN = '1.10' )
PARAMETER ( INUNT=7,OUTUNT=8 )
PARAMETER ( XMAX0=16,YMAX0=16,ZMAX0=16 )
PARAMETER ( NMBCS0=20,NMNEX0=10,NMNOU0=10 )
PARAMETER ( PYNX= 1,PZNX= 2,PXNY= 3,PZNY= 4,PYNZ= 5,PXNZ= 6 )
PARAMETER ( PYPZ= 7,PZPY= 8,PZPX= 9,PXPZ=10,PXPY=11,PYPX=12 )
PARAMETER ( Z0=376.7 )
*---- variable declarations
CHARACTER*80 LINE,COMM
CHARACTER*40 INNAME,OUTNAM
CHARACTER*2 LORIEN
CHARACTER TYP,DIRN,LSIGN
LOGICAL EX
INTEGER ERR,N,X,Y,Z,I
INTEGER X1,Y1,Z1,D,T
INTEGER XMAX,YMAX,ZMAX,NUMDT
INTEGER NUMBCS,NUMNEX,NUMNOU
REAL DL,VE,VH,OUT
REAL VYNX,VZNX,VXNY,VZNY,VYNZ,VXNZ
REAL VYPZ,VZPY,VZPX,VXPZ,VXPY,VYPX
REAL RHO,V0,VTEMP,VDIFF
*---- array declarations
CHARACTER BDIRN(NMBCS0),BSIGN(NMBCS0)
INTEGER BOX1(NMBCS0),BOY1(NMBCS0),BOZ1(NMBCS0)
INTEGER BOX2(NMBCS0),BOY2(NMBCS0),BOZ2(NMBCS0)
REAL BRHO(NMBCS0)
CHARACTER*3 NEID(NMNEX0)
INTEGER NEX1(NMNEX0),NEY1(NMNEX0),NEZ1(NMNEX0)
INTEGER NEX2(NMNEX0),NEY2(NMNEX0),NEZ2(NMNEX0)
REAL NEAMP(NMNEX0)
CHARACTER*2 NOID(NMNOU0)
INTEGER NOX1(NMNOU0),NOY1(NMNOU0),NOZ1(NMNOU0)
COMMON/COMV/V
REAL V(12,XMAX0,YMAX0,ZMAX0)
*---- display program details
WRITE(*,*)
WRITE(*,*) 'TLM3D (v'//VERSN//') - '//
; 'Transmission-Line Modelling in 3-Dimensions'
WRITE(*,*)
WRITE(*,*) 'Limits :'
WRITE(*,99010) 'mesh ',XMAX0,YMAX0,ZMAX0
WRITE(*,99020) 'boundaries',NMBCS0
WRITE(*,99020) 'excite ',NMNEX0
WRITE(*,99020) 'output ',NMNOU0
WRITE(*,*)
99010 FORMAT(1X,' ',A,' = ',I2,' x ',I2,' x ',I2)
99020 FORMAT(1X,' ',A,' = ',I2)
*---- open input file
WRITE(*,*) 'Name of input file : '
READ(*,'(A)') INNAME
INQUIRE(FILE=INNAME,EXIST=EX)
IF (.NOT.EX) STOP 'File does not exist'
OPEN(INUNT,FILE=INNAME,STATUS='OLD',ACCESS='SEQUENTIAL',
; FORM='FORMATTED',IOSTAT=ERR)
IF (ERR.NE.0) STOP 'File cannot be opened'
REWIND(INUNT)
*---- set default data
OUTNAM = '?'
XMAX = XMAX0
YMAX = YMAX0
ZMAX = ZMAX0
DL = 1.0
NUMDT = 0
NUMBCS = 0
NUMNEX = 0
NUMNOU = 0
*---- read input file
11010 READ(INUNT,'(A)',IOSTAT=ERR) LINE
IF (ERR.NE.0) STOP 'Error reading file'
IF (LINE(1:1).EQ.'!') GOTO 11010
IF (LINE(1:1).NE.':') STOP 'Start of data expected'
COMM = LINE(2:80)
*------ mesh size
IF (COMM.EQ.'MESH') THEN
READ(INUNT,*) XMAX,YMAX,ZMAX
IF (XMAX.GT.XMAX0.OR.YMAX.GT.YMAX0.OR.ZMAX.GT.ZMAX0)
; STOP 'Too big'
*------ node spacing
ELSE IF (COMM.EQ.'DL') THEN
READ(INUNT,*) DL
*------ number of timesteps
ELSE IF (COMM.EQ.'TIMESTEPS') THEN
READ(INUNT,*) NUMDT
*------ output file
ELSE IF (COMM.EQ.'OUTPUT_FILE') THEN
READ(INUNT,*) OUTNAM
*------ boundaries
ELSE IF (COMM.EQ.'BOUNDARIES') THEN
READ(INUNT,*) NUMBCS
IF (NUMBCS.GT.NMBCS0) STOP 'Too many boundaries'
DO 12010, I = 1, NUMBCS
READ(INUNT,*) BRHO(I),LORIEN,BOX1(I),BOY1(I),BOZ1(I),
; BOX2(I),BOY2(I),BOZ2(I)
BDIRN(I) = LORIEN(1:1)
BSIGN(I) = LORIEN(2:2)
IF ((BSIGN(I).NE.'N'.AND.BSIGN(I).NE.'P'.AND.
; BSIGN(I).NE.'0').OR.
; BOX1(I).LT.1.OR.BOX2(I).GT.XMAX.OR.BOX2(I).LT.BOX1(I).OR.
; BOY1(I).LT.1.OR.BOY2(I).GT.YMAX.OR.BOY2(I).LT.BOY1(I).OR.
; BOZ1(I).LT.1.OR.BOZ2(I).GT.ZMAX.OR.BOZ2(I).LT.BOZ1(I).OR.
; BDIRN(I).EQ.'X'.AND.(BOX1(I).NE.BOX2(I).OR.
; (BSIGN(I).EQ.'0'.AND.BOX2(I).EQ.XMAX)).OR.
; BDIRN(I).EQ.'Y'.AND.(BOY1(I).NE.BOY2(I).OR.
; (BSIGN(I).EQ.'0'.AND.BOY2(I).EQ.YMAX)).OR.
; BDIRN(I).EQ.'Z'.AND.(BOZ1(I).NE.BOZ2(I).OR.
; (BSIGN(I).EQ.'0'.AND.BOZ2(I).EQ.ZMAX)))
; STOP 'Bad boundary coordinates'
12010 CONTINUE
*------ excitation
ELSE IF (COMM.EQ.'EXCITE') THEN
READ(INUNT,*) NUMNEX
IF (NUMNEX.GT.NMNEX0) STOP 'Too many excite'
DO 13010, I = 1, NUMNEX
READ(INUNT,*) NEID(I),NEAMP(I),NEX1(I),NEY1(I),NEZ1(I),
; NEX2(I),NEY2(I),NEZ2(I)
IF (NEX1(I).LT.1.OR.NEX2(I).GT.XMAX.OR.NEX2(I).LT.NEX1(I).OR.
; NEY1(I).LT.1.OR.NEY2(I).GT.YMAX.OR.NEY2(I).LT.NEY1(I).OR.
; NEZ1(I).LT.1.OR.NEZ2(I).GT.ZMAX.OR.NEZ2(I).LT.NEZ1(I))
; STOP 'Bad excitation coordinates'
13010 CONTINUE
*------ output
ELSE IF (COMM.EQ.'OUTPUT') THEN
READ(INUNT,*) NUMNOU
IF (NUMNOU.GT.NMNOU0) STOP 'Too many output'
DO 14010, I = 1, NUMNOU
READ(INUNT,*) NOID(I),NOX1(I),NOY1(I),NOZ1(I)
IF (NOX1(I).LT.1.OR.NOX1(I).GT.XMAX.OR.
; NOY1(I).LT.1.OR.NOY1(I).GT.YMAX.OR.
; NOZ1(I).LT.1.OR.NOZ1(I).GT.ZMAX)
; STOP 'Bad output coordinates'
14010 CONTINUE
*---- end of read
ELSE IF (COMM.NE.'END') THEN
STOP 'Data not recognised'
END IF
IF (COMM.NE.'END') GOTO 11010
CLOSE(INUNT)
WRITE(*,*) 'Input file closed'
*---- open output file
IF (OUTNAM.EQ.'?') THEN
WRITE(*,*) 'Name of output file : '
READ(*,'(A)') OUTNAM
END IF
OPEN(OUTUNT,FILE=OUTNAM,STATUS='NEW',ACCESS='SEQUENTIAL',
; FORM='FORMATTED',IOSTAT=ERR)
IF (ERR.NE.0) STOP 'Output file cannot be opened'
*---- clear array
DO 21010, Z = 1, ZMAX
DO 21020, Y = 1, YMAX
DO 21030, X = 1, XMAX
DO 21040, D = 1, 12
V(D,X,Y,Z) = 0.0
21040 CONTINUE
21030 CONTINUE
21020 CONTINUE
21010 CONTINUE
*---- excite
DO 22010, N = 1, NUMNEX
TYP = NEID(N)(1:1)
DIRN = NEID(N)(2:2)
IF (TYP.EQ.'E') THEN
VE = - DL * NEAMP(N) / 2.0
ELSE IF (TYP.EQ.'H') THEN
VH = DL * Z0 * NEAMP(N) / 2.0
ELSE
STOP 'Excitation field error'
END IF
DO 22020, Z = NEZ1(N), NEZ2(N)
DO 22030, Y = NEY1(N), NEY2(N)
DO 22040, X = NEX1(N), NEX2(N)
IF (TYP.EQ.'E') THEN
IF (DIRN.EQ.'X') THEN
V(PYNX,X,Y,Z) = V(PYNX,X,Y,Z) + VE
V(PZNX,X,Y,Z) = V(PZNX,X,Y,Z) + VE
V(PYPX,X,Y,Z) = V(PYPX,X,Y,Z) + VE
V(PZPX,X,Y,Z) = V(PZPX,X,Y,Z) + VE
ELSE IF (DIRN.EQ.'Y') THEN
V(PZNY,X,Y,Z) = V(PZNY,X,Y,Z) + VE
V(PXNY,X,Y,Z) = V(PXNY,X,Y,Z) + VE
V(PZPY,X,Y,Z) = V(PZPY,X,Y,Z) + VE
V(PXPY,X,Y,Z) = V(PXPY,X,Y,Z) + VE
ELSE IF (DIRN.EQ.'Z') THEN
V(PXNZ,X,Y,Z) = V(PXNZ,X,Y,Z) + VE
V(PYNZ,X,Y,Z) = V(PYNZ,X,Y,Z) + VE
V(PXPZ,X,Y,Z) = V(PXPZ,X,Y,Z) + VE
V(PYPZ,X,Y,Z) = V(PYPZ,X,Y,Z) + VE
ELSE
STOP 'E-field excitation error'
END IF
ELSE IF (TYP.EQ.'H') THEN
IF (DIRN.EQ.'X') THEN
V(PZNY,X,Y,Z) = V(PZNY,X,Y,Z) + VH
V(PYNZ,X,Y,Z) = V(PYNZ,X,Y,Z) - VH
V(PZPY,X,Y,Z) = V(PZPY,X,Y,Z) - VH
V(PYPZ,X,Y,Z) = V(PYPZ,X,Y,Z) + VH
ELSE IF (DIRN.EQ.'Y') THEN
V(PXNZ,X,Y,Z) = V(PXNZ,X,Y,Z) + VH
V(PZNX,X,Y,Z) = V(PZNX,X,Y,Z) - VH
V(PXPZ,X,Y,Z) = V(PXPZ,X,Y,Z) - VH
V(PZPX,X,Y,Z) = V(PZPX,X,Y,Z) + VH
ELSE IF (DIRN.EQ.'Z') THEN
V(PYNX,X,Y,Z) = V(PYNX,X,Y,Z) + VH
V(PXNY,X,Y,Z) = V(PXNY,X,Y,Z) - VH
V(PYPX,X,Y,Z) = V(PYPX,X,Y,Z) - VH
V(PXPY,X,Y,Z) = V(PXPY,X,Y,Z) + VH
ELSE
STOP 'H-field excitation error'
END IF
END IF
22040 CONTINUE
22030 CONTINUE
22020 CONTINUE
22010 CONTINUE
*---- start calculation
WRITE(*,*) 'Start of calculation'
DO 30010, T = 1, NUMDT
WRITE(*,99100) 100.0*REAL(T-1)/REAL(NUMDT)
99100 FORMAT(1X,F6.3,' %')
*---- output
DO 41010, N = 1, NUMNOU
TYP = NOID(N)(1:1)
DIRN = NOID(N)(2:2)
Z = NOZ1(N)
Y = NOY1(N)
X = NOX1(N)
IF (TYP.EQ.'E') THEN
IF (DIRN.EQ.'X') THEN
OUT = V(PYNX,X,Y,Z) + V(PZNX,X,Y,Z) +
; V(PYPX,X,Y,Z) + V(PZPX,X,Y,Z)
ELSE IF (DIRN.EQ.'Y') THEN
OUT = V(PZNY,X,Y,Z) + V(PXNY,X,Y,Z) +
; V(PZPY,X,Y,Z) + V(PXPY,X,Y,Z)
ELSE IF (DIRN.EQ.'Z') THEN
OUT = V(PXNZ,X,Y,Z) + V(PYNZ,X,Y,Z) +
; V(PXPZ,X,Y,Z) + V(PYPZ,X,Y,Z)
ELSE
STOP 'E-field output error'
END IF
OUT = - OUT / DL / 2.0
ELSE IF (TYP.EQ.'H') THEN
IF (DIRN.EQ.'X') THEN
OUT = V(PZNY,X,Y,Z) - V(PYNZ,X,Y,Z) -
; V(PZPY,X,Y,Z) + V(PYPZ,X,Y,Z)
ELSE IF (DIRN.EQ.'Y') THEN
OUT = V(PXNZ,X,Y,Z) - V(PZNX,X,Y,Z) -
; V(PXPZ,X,Y,Z) + V(PZPX,X,Y,Z)
ELSE IF (DIRN.EQ.'Z') THEN
OUT = V(PYNX,X,Y,Z) - V(PXNY,X,Y,Z) -
; V(PYPX,X,Y,Z) + V(PXPY,X,Y,Z)
ELSE
STOP 'H-field output error'
END IF
OUT = OUT / Z0 / DL / 2.0
ELSE
STOP 'Output field error'
END IF
WRITE(OUTUNT,99200) OUT
99200 FORMAT(1X,E13.6)
41010 CONTINUE
*---- scatter
DO 42010, Z = 1, ZMAX
DO 42020, Y = 1, YMAX
DO 42030, X = 1, XMAX
VYNX = V(PYNX,X,Y,Z)
VZNX = V(PZNX,X,Y,Z)
VXNY = V(PXNY,X,Y,Z)
VZNY = V(PZNY,X,Y,Z)
VYNZ = V(PYNZ,X,Y,Z)
VXNZ = V(PXNZ,X,Y,Z)
VYPZ = V(PYPZ,X,Y,Z)
VZPY = V(PZPY,X,Y,Z)
VZPX = V(PZPX,X,Y,Z)
VXPZ = V(PXPZ,X,Y,Z)
VXPY = V(PXPY,X,Y,Z)
VYPX = V(PYPX,X,Y,Z)
VDIFF = VXPY - VXNY
VTEMP = 0.5 * ( VZNX + VZPX + VDIFF )
V(PYPX,X,Y,Z) = VTEMP
V(PYNX,X,Y,Z) = VTEMP - VDIFF
VDIFF = VXPZ - VXNZ
VTEMP = 0.5 * ( VYNX + VYPX + VDIFF )
V(PZPX,X,Y,Z) = VTEMP
V(PZNX,X,Y,Z) = VTEMP - VDIFF
VDIFF = VYPZ - VYNZ
VTEMP = 0.5 * ( VXNY + VXPY + VDIFF )
V(PZPY,X,Y,Z) = VTEMP
V(PZNY,X,Y,Z) = VTEMP - VDIFF
VDIFF = VYPX - VYNX
VTEMP = 0.5 * ( VZNY + VZPY + VDIFF )
V(PXPY,X,Y,Z) = VTEMP
V(PXNY,X,Y,Z) = VTEMP - VDIFF
VDIFF = VZPX - VZNX
VTEMP = 0.5 * ( VYNZ + VYPZ + VDIFF )
V(PXPZ,X,Y,Z) = VTEMP
V(PXNZ,X,Y,Z) = VTEMP - VDIFF
VDIFF = VZPY - VZNY
VTEMP = 0.5 * ( VXNZ + VXPZ + VDIFF )
V(PYPZ,X,Y,Z) = VTEMP
V(PYNZ,X,Y,Z) = VTEMP - VDIFF
42030 CONTINUE
42020 CONTINUE
42010 CONTINUE
*---- boundaries
DO 43010, N = 1, NUMBCS
DIRN = BDIRN(N)
LSIGN = BSIGN(N)
RHO = BRHO(N)
IF (DIRN.EQ.'X') THEN
X = BOX1(N)
IF (LSIGN.EQ.'N'.AND.X.EQ.1) THEN
DO 43110, Z = BOZ1(N), BOZ2(N)
DO 43120, Y = BOY1(N), BOY2(N)
V(PXNY,X,Y,Z) = RHO * V(PXNY,X,Y,Z)
V(PXNZ,X,Y,Z) = RHO * V(PXNZ,X,Y,Z)
43120 CONTINUE
43110 CONTINUE
ELSE IF (LSIGN.EQ.'P'.AND.X.EQ.XMAX) THEN
DO 43130, Z = BOZ1(N), BOZ2(N)
DO 43140, Y = BOY1(N), BOY2(N)
V(PXPZ,X,Y,Z) = RHO * V(PXPZ,X,Y,Z)
V(PXPY,X,Y,Z) = RHO * V(PXPY,X,Y,Z)
43140 CONTINUE
43130 CONTINUE
ELSE
IF (LSIGN.EQ.'N') X = X - 1
X1 = X + 1
DO 43150, Z = BOZ1(N), BOZ2(N)
DO 43160, Y = BOY1(N), BOY2(N)
VTEMP = V(PXPY,X,Y,Z)
V(PXPY,X,Y,Z) = RHO * V(PXNY,X1,Y,Z)
V(PXNY,X1,Y,Z) = RHO * VTEMP
VTEMP = V(PXPZ,X,Y,Z)
V(PXPZ,X,Y,Z) = RHO * V(PXNZ,X1,Y,Z)
V(PXNZ,X1,Y,Z) = RHO * VTEMP
43160 CONTINUE
43150 CONTINUE
END IF
ELSE IF (DIRN.EQ.'Y') THEN
Y = BOY1(N)
IF (LSIGN.EQ.'N'.AND.Y.EQ.1) THEN
DO 43210, Z = BOZ1(N), BOZ2(N)
DO 43220, X = BOX1(N), BOX2(N)
V(PYNZ,X,Y,Z) = RHO * V(PYNZ,X,Y,Z)
V(PYNX,X,Y,Z) = RHO * V(PYNX,X,Y,Z)
43220 CONTINUE
43210 CONTINUE
ELSE IF (LSIGN.EQ.'P'.AND.Y.EQ.YMAX) THEN
DO 43230, Z = BOZ1(N), BOZ2(N)
DO 43240, X = BOX1(N), BOX2(N)
V(PYPZ,X,Y,Z) = RHO * V(PYPZ,X,Y,Z)
V(PYPX,X,Y,Z) = RHO * V(PYPX,X,Y,Z)
43240 CONTINUE
43230 CONTINUE
ELSE
IF (LSIGN.EQ.'N') Y = Y - 1
Y1 = Y + 1
DO 43250, Z = BOZ1(N), BOZ2(N)
DO 43260, X = BOX1(N), BOX2(N)
VTEMP = V(PYPZ,X,Y,Z)
V(PYPZ,X,Y,Z) = RHO * V(PYNZ,X,Y1,Z)
V(PYNZ,X,Y1,Z) = RHO * VTEMP
VTEMP = V(PYPX,X,Y,Z)
V(PYPX,X,Y,Z) = RHO * V(PYNX,X,Y1,Z)
V(PYNX,X,Y1,Z) = RHO * VTEMP
43260 CONTINUE
43250 CONTINUE
END IF
ELSE IF (DIRN.EQ.'Z') THEN
Z = BOZ1(N)
IF (LSIGN.EQ.'N'.AND.Z.EQ.1) THEN
DO 43310, Y = BOY1(N), BOY2(N)
DO 43320, X = BOX1(N), BOX2(N)
V(PZNX,X,Y,Z) = RHO * V(PZNX,X,Y,Z)
V(PZNY,X,Y,Z) = RHO * V(PZNY,X,Y,Z)
43320 CONTINUE
43310 CONTINUE
ELSE IF (LSIGN.EQ.'P'.AND.Z.EQ.ZMAX) THEN
DO 43330, Y = BOY1(N), BOY2(N)
DO 43340, X = BOX1(N), BOX2(N)
V(PZPX,X,Y,Z) = RHO * V(PZPX,X,Y,Z)
V(PZPY,X,Y,Z) = RHO * V(PZPY,X,Y,Z)
43340 CONTINUE
43330 CONTINUE
ELSE
IF (LSIGN.EQ.'N') Z = Z - 1
Z1 = Z + 1
DO 43350, Y = BOY1(N), BOY2(N)
DO 43360, X = BOX1(N), BOX2(N)
VTEMP = V(PZPX,X,Y,Z)
V(PZPX,X,Y,Z) = RHO * V(PZNX,X,Y,Z1)
V(PZNX,X,Y,Z1) = RHO * VTEMP
VTEMP = V(PZPY,X,Y,Z)
V(PZPY,X,Y,Z) = RHO * V(PZNY,X,Y,Z1)
V(PZNY,X,Y,Z1) = RHO * VTEMP
43360 CONTINUE
43350 CONTINUE
END IF
ELSE
STOP 'Boundary direction error'
END IF
43010 CONTINUE
*---- connect
DO 44010, Z = 1, ZMAX
Z1 = Z + 1
DO 44020, Y = 1, YMAX
Y1 = Y + 1
DO 44030, X = 1, XMAX
X1 = X + 1
IF (X.LT.XMAX) THEN
V0 = V(PXPZ,X,Y,Z)
V(PXPZ,X,Y,Z) = V(PXNZ,X1,Y,Z)
V(PXNZ,X1,Y,Z) = V0
V0 = V(PXPY,X,Y,Z)
V(PXPY,X,Y,Z) = V(PXNY,X1,Y,Z)
V(PXNY,X1,Y,Z) = V0
END IF
IF (Y.LT.YMAX) THEN
V0 = V(PYPZ,X,Y,Z)
V(PYPZ,X,Y,Z) = V(PYNZ,X,Y1,Z)
V(PYNZ,X,Y1,Z) = V0
V0 = V(PYPX,X,Y,Z)
V(PYPX,X,Y,Z) = V(PYNX,X,Y1,Z)
V(PYNX,X,Y1,Z) = V0
END IF
IF (Z.LT.ZMAX) THEN
V0 = V(PZPY,X,Y,Z)
V(PZPY,X,Y,Z) = V(PZNY,X,Y,Z1)
V(PZNY,X,Y,Z1) = V0
V0 = V(PZPX,X,Y,Z)
V(PZPX,X,Y,Z) = V(PZNX,X,Y,Z1)
V(PZNX,X,Y,Z1) = V0
END IF
44030 CONTINUE
44020 CONTINUE
44010 CONTINUE
*---- end calculation
30010 CONTINUE
*---- close output file
CLOSE(OUTUNT)
WRITE(*,*) 'Output file closed'
WRITE(*,*)
END