PROGRAM SET4
C=======================================================================
C DATA GENERATING PROGRAM FOR HELM4Q.FOR
C PROJECT: DOUBLE-LET PROJECT NAME: DOMAIN
C DOMAIN: INFINITE
C ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C EIJI FUKUMORI DECEMBER 29, 1993
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=4,MXE=5000000,MXN=5200000, MXI=1000, NCIR=5000 )
PARAMETER (MXD=10)
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN),FJ(MXE)
DIMENSION XI(MXI), YI(MXI),ICIR(NCIR)
COMPLEX*16 BASEVEC , UNITVEC
C=======================================================================
C RADIUS=RADIUS OF ROD ELEMENT
C AREA = CROSS SECTIONAL AREA OF THE ROD
C XC1 AND YC1 ARE THE CENTER COORDINATES OF THE UPPER conductor
C XC2 AND YC2 ARE THE CENTER COORDINATES OF THE LOWER conductor
C GAP = SPACING BETWEEN THE UPPER AND THE LOWER conductorS
C FJ(I) = CHARGE PER UNIT AREA
C GROSS CHARGE IN THE UPPER ROD = +1
C GROSS CHARGE IN THE LOWER ROD = -1
C Q = GROSS CHARGE VALUE FED INTO RODS
C CHARGE = Q PER UNIT AREA
C=======================================================================
PI = 4.D0* ATAN( 1.D0)
Q = 1.D0
RADIUS = 1.0D0
AREA = PI * RADIUS*RADIUS
C=======================================================================
CHARGE = Q / AREA
C=======================================================================
WRITE (*,*) 'CHARGE=',Q
WRITE (*,*) 'RADIUS=',RADIUS
WRITE (*,*) 'Type in a value of GAP'
READ (*,*) GAP
WRITE (*,*) 'GAP=',GAP
IF ( GAP .LE. 0.D0 ) STOP 'GAP MUST BE MORE THAN ZERO!!!!'
WRITE (*,*) 'D/a=', (GAP+2.D0*RADIUS)/RADIUS
C=======================================================================
XC1 = 0.
YC1 = RADIUS + GAP/2.D0
XC2 = XC1
YC2 = -YC1
NEX = 14
WRITE (*,*) 'NEX=', NEX
C=======================================================================
C ELEMENT CREATION
NEX = NEX/2*2
NEY = NEX
NEX = NEX/2*2
NEY = NEY/2*2
NDX = NEX + 1
NDY = NEY + 1
NE = 0
DO I = 1 , NEY
DO J = 1 , NEX
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
NODEX(NE,1) = NDX*(I-1) + J
NODEX(NE,2) = NODEX(NE,1) + 1
NODEX(NE,3) = NODEX(NE,2) +NDX
NODEX(NE,4) = NODEX(NE,1)+NDX
FJ(NE) = CHARGE
END DO
END DO
C
DO I = 1 , NEY
DO J = 1 , NEX
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
NODEX(NE,1) = NDX*(I-1) + J + NDX*NDY
NODEX(NE,2) = NODEX(NE,1) + 1
NODEX(NE,3) = NODEX(NE,2) +NDX
NODEX(NE,4) = NODEX(NE,1)+NDX
FJ(NE) = -CHARGE
END DO
END DO
C=======================================================================
DL = RADIUS/(NEX/2)
NNODE = 0
C=========================== UPPER CONDUCTOR ===========================
DO J = 1 , NDY
DO I = 1 , NDX
NNODE = NNODE + 1
XCOORD(NNODE) = (I-(NEX/2+1))*DL
YCOORD(NNODE) = (J-(NEY/2+1))*DL
END DO
END DO
C--------- CIRCULARIZATION ----------
C -- CENTER NODE = JNODE --
JNODE = NDX*(NEX/2)+NEX/2+1
XCOORD(JNODE) = 0.D0
YCOORD(JNODE) = 0.D0
DO I = 1 , NNODE
IF ( I .NE. JNODE ) THEN
VECTOR = DSQRT (XCOORD(I)**2 + YCOORD(I)**2)
AX = XCOORD(I)/VECTOR
AY = YCOORD(I)/VECTOR
R = DMAX1 ( DABS(XCOORD(I)), DABS(YCOORD(I)) )
XCOORD(I) = R * AX
YCOORD(I) = R * AY
IF ( IDINT(DABS(R-RADIUS)*100000.D0) .EQ. 0 ) THEN
END IF
END IF
END DO
C------------ NODES ON SURFACE OF UPPER CONDUCTOR -------------
NC = NEX*2 + NEY*2
C ------- JSTART = STARTING NODE -------
JSTART = ( NEY/2 +2 )*NDX
ICIR(1) = JSTART
DO J = 2, NC
AX = XCOORD(ICIR(J-1))
AY = YCOORD(ICIR(J-1))
DO IEL = 1 , NE/2
DO K = 1 , ND
IF ( NODEX(IEL,K) .EQ. ICIR(J-1) ) THEN
NEXT = K + 1
IF ( K .EQ. ND ) NEXT = 1
NI = NODEX(IEL,NEXT)
R = DSQRT ( XCOORD(NI)**2 + YCOORD(NI)**2 )
IF ( IDINT(DABS(R-RADIUS)*100000.D0) .EQ. 0 ) THEN
BX = XCOORD(NI)
BY = YCOORD(NI)
IF ( AX*BY-BX*AY .GT. 0.D0 ) THEN
ICIR(J) = NI
END IF
END IF
END IF
END DO
END DO
END DO
C------------ EVEN SPACING BETWEEN OUTMOST MODES -------------
BASEVEC = DCMPLX( RADIUS, 0.D0 )
UNITANG = 2.D0*PI / (NEX*2 + NEY*2)
UNITVEC = DCMPLX( DCOS(UNITANG), DSIN(UNITANG) )
DO I = 1 , NC
BASEVEC = BASEVEC * UNITVEC
J = ICIR(I)
XCOORD (J) = DREAL( BASEVEC )
YCOORD (J) = DIMAG( BASEVEC )
END DO
C========================= LOWER CONDUCTOR =========================
DO J = 1 , NDY
DO I = 1 , NDX
NNODE = NNODE + 1
XCOORD(NNODE) = (I-(NEX/2+1))*DL
YCOORD(NNODE) = (J-(NEY/2+1))*DL
END DO
END DO
C--------- CIRCULARIZATION -----
C CENTER NODE = JNODE
JNODE = NDX*(NEX/2)+NEX/2+1 + NDX*NDY
XCOORD(JNODE) = 0.D0
YCOORD(JNODE) = 0.D0
DO I = NDX*NDY+1 , NNODE
IF ( I .NE. JNODE ) THEN
VECTOR = DSQRT (XCOORD(I)**2 + YCOORD(I)**2)
AX = XCOORD(I)/VECTOR
AY = YCOORD(I)/VECTOR
R = DMAX1 ( DABS(XCOORD(I)), DABS(YCOORD(I)) )
XCOORD(I) = R * AX
YCOORD(I) = R * AY
END IF
END DO
C------------ NODES ON SURFACE OF UPPER CONDUCTOR -------------
NC = 2*(NEX*2 + NEY*2)
C ------- JSTART = STARTING NODE -------
JSTART = ( NEY/2 +2 )*NDX + NDX*NDY
ICIR(NC/2+1) = JSTART
DO J = NC/2+2, NC
AX = XCOORD(ICIR(J-1))
AY = YCOORD(ICIR(J-1))
DO IEL = NE/2+1 , NE
DO K = 1 , ND
IF ( NODEX(IEL,K) .EQ. ICIR(J-1) ) THEN
NEXT = K + 1
IF ( K .EQ. ND ) NEXT = 1
NI = NODEX(IEL,NEXT)
R = DSQRT ( XCOORD(NI)**2 + YCOORD(NI)**2 )
IF ( IDINT(DABS(R-RADIUS)*100000.D0) .EQ. 0 ) THEN
BX = XCOORD(NI)
BY = YCOORD(NI)
IF ( AX*BY-BX*AY .GT. 0 ) THEN
ICIR(J) = NI
END IF
END IF
END IF
END DO
END DO
END DO
C------------ EVEN SPACING BETWEEN OUTMOST MODES -------------
BASEVEC = DCMPLX( RADIUS, 0.D0 )
UNITANG = 2.D0*PI / (NEX*2 + NEY*2)
UNITVEC = DCMPLX( DCOS(UNITANG), DSIN(UNITANG) )
DO I = NC/2+1 , NC
BASEVEC = BASEVEC * UNITVEC
J = ICIR(I)
XCOORD (J) = DREAL( BASEVEC )
YCOORD (J) = DIMAG( BASEVEC )
END DO
C-------- RELOACTION OF NODAL COORDINATES -------
DO I = 1 , NNODE/2
YCOORD (I) = YCOORD(I) + YC1
END DO
DO I = NNODE/2+1, NNODE
YCOORD (I) = YCOORD(I) + YC2
END DO
C=======================================================================
C================INTERNAL POINTS WHERE U(X,Y) IS EVALUATED==============
C=======================================================================
EPS = 1.0D-6
DELADD = RADIUS*PI*EPS
NIP = 0
DIA = 2.D0*RADIUS
C------- UPPER OUTER SPACE ------
NDIV = 10
DY = DIA
DO I = 1 , NDIV
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = (NDIV+1-I)*DY + 2.D0*DIA + GAP/2.D0
END DO
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = 2.D0*DIA + GAP/2.D0 + DELADD
C-------- ABOVE THE UPPER CONDUCTOR -------
NDIV = 10
Y1 = DIA + GAP/2.D0 + EPS
Y2 = 2.D0*DIA + GAP/2.D0 - EPS
DY = ( Y2 - Y1 ) / NDIV
DO I = 1 , NDIV-1
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = Y2 - DY*I
END DO
C------- UPPER CONDUTOR ------
NDIV = 11
Y1 = GAP/2.D0 + DELADD
Y2 = DIA + GAP/2.D0 - DELADD
DY = ( Y2 - Y1 ) / NDIV
DO I = 2 , NDIV
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = Y2 - DY*(I-1)
END DO
C-------- BETWEEN THE ORIGIN AND THE UPPER CONDUTOR
Y1 = EPS
Y2 = GAP/2 - EPS
NDIV = 10
DY = ( Y2 - Y1 ) / NDIV
DO I = 1 , NDIV-1
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = Y2 - DY*I
END DO
C-------- ORIGIN OF THE COORDINATE SYSTEM
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = 0.00D0
C-------- BETWEEN THE ORIGIN AND THE LOWER CONDUTOR
Y1 = -EPS
Y2 = -(GAP/2 - EPS)
NDIV = 10
DY = ( Y2 - Y1 ) / NDIV
DO I = 1 , NDIV-1
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = Y1 + DY*I
END DO
C------- LOWER CONDUTOR ------
NDIV = 11
Y1 = -(GAP/2.D0 + DELADD)
Y2 = -(DIA + GAP/2.D0 - DELADD)
DY = ( Y2 - Y1 ) / NDIV
DO I = 2 , NDIV
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = Y1 + DY*(I-1)
END DO
C-------- BELOW THE LOWER CONDUCTOR -------
NDIV = 10
Y1 = -(DIA + GAP/2.D0 + EPS)
Y2 = -(2.D0*DIA + GAP/2.D0 - EPS)
DY = ( Y2 - Y1 ) / NDIV
DO I = 1 , NDIV
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = Y1 + DY*I
END DO
C------- LOWER OUTER SPACE ------
NDIV = 10
DY = DIA
Y1 = -(2.D0*DIA + GAP/2.D0 +DELADD)
DO I = 1 , NDIV+1
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = Y1 - DY*(I-1)
END DO
C=======================================================================
C MAKING DATA FILES
OPEN ( 1, FILE='DOMAIN.DAT', STATUS = 'UNKNOWN' )
WRITE(1,*) NE
C------ ELEMENTS CONECTIVITY INFORMATION
DO I = 1 , NE
WRITE (1,*) I,(NODEX(I,J),J=1,ND), FJ(I)
END DO
C------ NODAL COORDINATES
WRITE(1,*) NNODE
DO I = 1 , NNODE
WRITE(1,*) I, XCOORD(I), YCOORD(I)
END DO
C------ INTERNAL POINT
WRITE(1,*) NIP
IF ( NIP .GE. 1 ) THEN
DO I = 1 , NIP
WRITE(1,*) XI(I), YI(I)
END DO
END IF
WRITE (1,*) NC
DO I = 1 , NC
WRITE (1,*) ICIR(I)
END DO
CLOSE (1)
WRITE(*,*) 'NE=',NE
WRITE(*,*) 'NNODE=',NNODE
OPEN ( 2, FILE='ELEMENTCG.DAT', STATUS = 'UNKNOWN' )
DO i = 1, NE
DO j = 1, ND
WRITE (2,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (2,*) XCOORD(NODEX(I,1)), YCOORD(NODEX(I,1))
WRITE (2,*)
END DO
CLOSE (2)
C=======================================================================
OPEN ( 4, FILE='SURFACE.DAT', STATUS = 'UNKNOWN' )
DO I = 1 , NC/2
J = I + 1
IF ( I .EQ. NC/2 ) J = 1
WRITE (4,*) XCOORD (ICIR(I)), YCOORD(ICIR(I))
WRITE (4,*) XCOORD (ICIR(J)), YCOORD(ICIR(J))
WRITE (4,*)
END DO
DO I = NC/2+1 , NC
J = I + 1
IF ( I .EQ. NC ) J = NC/2+1
WRITE (4,*) XCOORD (ICIR(I)), YCOORD(ICIR(I))
WRITE (4,*) XCOORD (ICIR(J)), YCOORD(ICIR(J))
WRITE (4,*)
END DO
CLOSE (4)
C=======================================================================
C----- COORDINATES FOR COMPUTATION OF MAGNETIC FIELD DENSITY AROUND WIRE
OPEN ( 3, FILE='DPDN.DAT', STATUS = 'UNKNOWN' )
WRITE (3,*) NC, XC1, YC1
DN = 0.001D0
DO J = 1 , NC/2
I = J - 1
IF ( I .EQ. 0 ) I = NC/2
X1 = XCOORD(ICIR(I)) - XC1
Y1 = YCOORD(ICIR(I)) - YC1
X2 = XCOORD(ICIR(J)) - XC1
Y2 = YCOORD(ICIR(J)) - YC1
XP = (X1+X2)/2.D0
YP = (Y1+Y2)/2.D0
C***************************************************
CCCCCCCC XP = XP * 2.D0
CCCCCCCC YP = YP * 2.D0
C***************************************************
DX = X2 - X1
DY = Y2 - Y1
DS = DSQRT(DX**2 + DY**2)
FNX = DY/DS
FNY = -DX/DS
WRITE (3,*) XP+XC1, YP+YC1
WRITE (3,*) XP+XC1+FNX*DN, YP+YC1+FNY*DN
END DO
CLOSE (3)
C=======================================================================
C
C ELEMENT CREATION FOR VECTOR PLOT GRAPHICS
C
C=======================================================================
DL = RADIUS / 2.6D0
XLENGTH = 8*RADIUS
YLENGTH = 6*RADIUS + GAP
NEX = XLENGTH / DL
NEY = YLENGTH / DL
NEX = NEX/2*2
NEY = NEY/2*2
DX = XLENGTH / NEX
DY = YLENGTH / NEY
NDX = NEX + 1
NDY = NEY + 1
NE = 0
DO I = 1 , NEY
DO J = 1 , NEX
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
NODEX(NE,1) = NDX*(I-1) + J
NODEX(NE,2) = NODEX(NE,1) + 1
NODEX(NE,3) = NODEX(NE,2) +NDX
NODEX(NE,4) = NODEX(NE,1)+NDX
END DO
END DO
C--------------- NODAL INFORMATION -------------------
NNODE = 0
DO J = 1 , NDY
DO I = 1 , NDX
NNODE = NNODE + 1
XCOORD(NNODE) = (I-(NEX/2+1))*DX
YCOORD(NNODE) = (J-(NEY/2+1))*DY
END DO
END DO
C--------------MAKING INPUT FILE------------------------
OPEN ( 5, FILE='VECTORCG.DAT', STATUS = 'UNKNOWN' )
WRITE(5,*) NE
DO I = 1 , NE
WRITE (5,*) I,(NODEX(I,J),J=1,ND)
END DO
WRITE(5,*) NNODE
DO I = 1 , NNODE
WRITE(5,*) I, XCOORD(I), YCOORD(I)
END DO
CLOSE (5)
C=======================================================================
STOP "NORMAL TERMINATION"
END