home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Antennas
/
Antennas_CD-ROM_Walnut_Creek_September_1996.iso
/
thinwire
/
fortran
/
tw3.for
< prev
next >
Wrap
Text File
|
1996-06-30
|
35KB
|
1,188 lines
C MAIN PROGRAM THINWIRE
C ORIGINAL PROGRAM BY J. H. RICHMOND, 1974
C MODIFIED 1992, 1994, R. P. HAVILAND
COMPLEX EP2,EP3,ETA,GAM,Y11,Z11,ZS
COMPLEX EPPS,EPTS,ETPS,ETTS,EX,EY,EZ
COMPLEX C(1830),CJ(60),EP(60),ET(60),EPP(60),ETT(60)
DIMENSION I1(60),I2(60),I3(60),JA(60),JB(60)
COMPLEX CGD(50),SGD(50),CG(100),VG(100),ZLD(100)
DIMENSION D(50),IA(50),IB(50),ISC(50),MD(50,4),ND(50)
DIMENSION X(55),Y(55),Z(55)
CHARACTER * 31 INFILE
CHARACTER * 31 OUTFIL
CHARACTER * 31 NAS
CHARACTER * 31 DAS
CHARACTER * 31 FPS
CHARACTER * 2 CR,LF,SI
CHARACTER * 3 ESC
DATA PI,TP/3.14159,6.28318/
DATA E0,U0/8.854E-12,1.2566E-6/
2 FORMAT ( 1X,8F15.7)
3 FORMAT ( 5F15.7)
4 FORMAT ( 1X,1I5,8F14.6)
5 FORMAT ( 1H0)
6 FORMAT ( 1X,6F15.7)
7 FORMAT ( 8F10.5)
8 FORMAT ( 1X,1I4,13I5)
9 FORMAT ( 3X,'MAX = ',I5,3X,'MIN = ',I5,3X,'N = ',I5)
C 10 FORMAT ( 4Z2.0)
CR=CHAR(13)
LF=CHAR(10)
ESC=CHAR(27)
SI=CHAR(15)
ICJ=60
INM=50
DO 15 J=1,INM
15 ISC(J)=0
6020 PRINT *, ' ENTER FILEPATH+FILENAME FOR INPUT, INCLUDE QUOTES '
READ (*,*)INFILE
C SETUP FOR FILE INPUT
OPEN (UNIT=5, FILE=INFILE,ERR=7020)
GOTO 6021
7020 GOTO 6020
C GET FILEOUT DATA
6021 PRINT *,' ENTER 1=OUTPUT TO PRINTER, 2=OUTPUT TO FILE '
READ (*,*)IFILOU
IF (IFILOU.EQ. 1) GOTO 6022
IF (IFILOU.NE. 2) GOTO 6021
6050 PRINT *,' ENTER FILEPATH+FILENAME FOR OUTPUT, INCLUDE QUOTES '
READ (*,*)OUTFIL
C SETUP FOR FILE OUTPUT
OPEN (UNIT=6,FILE=OUTFIL,ERR=7050,STATUS='NEW')
GOTO 6023
7050 GOTO 6050
C SETUP FOR PRINTER OUTPUT
6022 OPEN (UNIT=6,FILE='PRN')
WRITE (6,*) SI
C NOW GET DATA
6023 READ (5,*)NAS,DAS,FPS
CALL GETDAT(IYR,IMON,IDAY)
CALL GETTIM(IHR,IMIN)
WRITE (6,*)' ',IDAY,IMON,IYR,' AT ',IHR,IMIN
WRITE (6,*) ' THIN-WIRE ANALYSIS PROGRAM '
WRITE (6,*) ' ANTENNA ',NAS
WRITE (6,*) ' OF FILE ',FPS
WRITE (6,*) ' DATED ',DAS
WRITE (6,*) ' REFERENCE ',INFILE
WRITE (6,*) ' '
READ (5,*)BM,ER2,SIG2,TD2
WRITE (6,*)'INSULATION DIA. K LOSS'
WRITE (6,2)BM,ER2,SIG2,TD2
READ (5,*)AM,CMM,ER3,SIG3,TD3
WRITE (6,*)'WIRE DIA RHO MEDIUM'
WRITE (6,2)AM,CMM,ER3,SIG3,TD3
READ (5,*)IBISC,IGAIN,INEAR,ISCAT,IWR,NGEN,NM,NP,NLOAD
WRITE (6,*)' BIS GAIN NEAR SCAT CUR :GENS SEGS PTS LOADS'
WRITE (6,8)IBISC,IGAIN,INEAR,ISCAT,IWR,NGEN,NM,NP,NLOAD
READ (5,*)FMC,PHA,PHAI,PHAN,THA,THAI,THAN
WRITE (6,*)'FREQUENCY, MHZ'
WRITE (6,2)FMC
WRITE (6,*)'THETA INIT., INC., STEPS PHI INIT., INCR., STEPS'
WRITE (6,*)'FAR FIELD ANGLES'
WRITE (6,2)PHA,PHAI,PHAN,THA,THAI,THAN
READ (5,*)PHI,PHII,PHIN,THI,THII,THIN
WRITE (6,*)'BACK SCATTER ANGLES'
WRITE (6,2)PHI,PHII,PHIN,THI,THII,THIN
READ (5,*)PHS,PHSI,PHSN,THS,THSI,THSN
WRITE (6,*)'BISTATIC SCATTER ANGLES'
WRITE (6,2)PHS,PHSI,PHSN,THS,THSI,THSN
WRITE (6,*)'SEGMENT, END1, END2'
DO 22 J=1,NM
READ (5,*)IA(J),IB(J)
22 WRITE (6,8)J,IA(J),IB(J)
WRITE (6,*)'POINT LOCATIONS, XYZ METERS'
DO 40 I=1,NP
READ (5,*)X(I),Y(I),Z(I)
40 WRITE (6,4)I,X(I),Y(I),Z(I)
41 WRITE (6,*)' NEAR FIELD X Y Z'
READ (5,*) XP,YP,ZP
42 WRITE (6,*) XP,YP,ZP
43 FHZ=FMC*1.E6
OMEGA=TP*FHZ
DO 44 I=1,NP
IF (BM .GT. AM) ISC(I)=1
44 CONTINUE
IF (SIG2.LT.0.) EP2=ER2*E0*CMPLX(1.,-TD2)
IF (TD2.LT.0.) EP2=CMPLX(ER2*E0,-SIG2/OMEGA)
IF (SIG3.LT.0.) EP3=ER3*E0*CMPLX(1.,-TD3)
IF (TD3.LT.0.) EP3=CMPLX(ER3*E0,-SIG3/OMEGA)
ETA=CSQRT(U0/EP3)
GAM=OMEGA*CSQRT(-U0*EP3)
59 CALL SORT(IA,IB,I1,I2,I3,JA,JB,MD,ND,NM,NP,N,MAX,MIN,ICJ,INM)
WRITE (6,*)' SEGMENT COUPLING '
WRITE (6,9)MAX,MIN,N
WRITE (6,*)' '
IF (MIN.LT.1) WRITE (6,*)'*** UNCONNECTED WIRE ***'
IF (MAX.GT.4 .OR. MIN.LT.1 .OR. N.GT.ICJ) GO TO 800
INT=4
I12=1
DO 60 J=1,NM
VG(J)=(.0,.0)
ZLD(J)=(.0,.0)
JJ=J+NM
VG(JJ)=(.0,.0)
60 ZLD(JJ)=(.0,.0)
WRITE (6,*)'LOAD LR, LI POINT#'
NFIN=NLOAD
IF (NFIN .EQ. 0) NFIN=1
DO 61 I=1,NFIN
READ (5,*)LOADR,LOADI,JLOAD
ZLD(JLOAD)=CMPLX(LOADR,LOADI)
61 WRITE (6,*)ZLD(I),JLOAD
WRITE (6,*)'GENERATOR VR, VI POINT#'
NFIN=NGEN
IF (NFIN .EQ. 0) NFIN=1
DO 62 I=1,NFIN
READ (5,*)GENR,GENI,JGEN
JPT=JGEN
VG(JPT)=CMPLX(GENR,GENI)
62 WRITE (6,*)VG(JPT),JGEN
PRINT *, 'SET INTEGRATION, 0=NO CHANGE, -1=EXACT'
PRINT *, ' <4=FAST, LOW, >4=SLOW, HIGH ACCURACY'
READ (*,*) INTSET
IF (INTSET .LE. -1) INT=0
IF (INTSET .GT. 0) INT=INTSET
WRITE (6,*) 'INT VARIABLE =',INT
CALL SGANT(IA,IB,INM,INT,ISC,I1,I2,I3,JA,JB,MD,N,ND,NM,NP
2,AM,BM,C,CGD,CMM,D,EP2,EP3,ETA,FHZ,GAM,SGD,X,Y,Z,ZLD,ZS)
IF (N.LE.0) GO TO 800
IF (NGEN.LE.0) GO TO 400
63 CALL GANT1(IA,IB,INM,IWR,I1,I2,I3,I12,JA,JB,MD,N,ND,NM,AM
2,C,CJ,CG,CMM,D,EFF,GAM,GG,CGD,SGD,VG,Y11,Z11,ZLD,ZS)
WRITE (6,*)'EFFICIENCY INPOWER ZIN'
WRITE (6,3)EFF,GG,Z11
200 IF (INEAR.LE.0) GO TO 300
CALL GNFLD(IA,IB,INM,I1,I2,I3,MD,N,ND,NM,AM,CGD,SGD,ETA,GAM
2,CJ,D,X,Y,Z,XP,YP,ZP,EX,EY,EZ)
WRITE (6,*)'NEAR FIELD VALUES'
WRITE (6,*)'X Y Z'
WRITE (6,3)XP,YP,ZP
WRITE (6,*)'ELECTRIC FIELD'
WRITE (6,3)EX,EY,EZ
300 IF (IGAIN.LE.0) GO TO 400
INC=0
WRITE (6,*)'FAR FIELD VALUES'
WRITE (6,*)'AZMUITH ELEVATION H-FIELD V-FIELD TOTAL'
DO 310 I=PHA,PHAI*PHAN,PHAI
DO 310 J=THA,THAI*THAN,THAI
PH=I
TH=J
CALL GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
4,X,Y,Z,ZLD,ZS,ETA,GAM)
TOT=GPP*GPP+GTT*GTT
IF (TOT .LE. 0) DTOT=-999
IF (TOT .GT. 0) DTOT=5*LOG10(TOT)
IF (TOT .GT. 1E35) DTOT=999
WRITE (6,3)PH,TH,GPP,GTT,DTOT
310 CONTINUE
400 IF (ISCAT.LE.0) GO TO 600
INC=1
WRITE (6,*)'BACKSCATTER VALUES, ECHO AREAS, CROSS SECTIONS'
WRITE (6,*)'AZMUITH ELEVATION APP APT ATP ATT'
WRITE (6,*)'ABSORB P,T EXTINCT P,T SCATTER P,T'
DO 410 I=PHI,PHII*PHIN,PHII
DO 410 J=THI,THII*THIN,THII
PH=I
TH=J
CALL GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
4,X,Y,Z,ZLD,ZS,ETA,GAM)
WRITE (6,6)PH,TH,SPPM,SPTM,STPM,STTM
WRITE (6,6)ACSP,ACST,ECSP,ECST,SCSP,SCST
410 CONTINUE
500 IF (IBISC.LE.0) GO TO 600
INC=2
WRITE (6,*)'BI-STAT SCATTER VALUES'
WRITE (6,*)' WITH INCIDENT WAVE PHI, THETA OF'
WRITE (6,*) PHI, THI
WRITE (6,*)'AZMUITH ELEVATION H-FIELD V-FIELD'
DO 510 I=PHS,PHSI,PHSN
DO 510 J=THS,THSI,THSN
PH=I
TH=J
CALL GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
4,X,Y,Z,ZLD,ZS,ETA,GAM)
WRITE (6,6)PH,TH,SPPM,SPTM,STPM,STTM
510 CONTINUE
600 CONTINUE
800 CLOSE (5)
CLOSE (2)
CLOSE (6)
WRITE (6,*) ESC,'@'
PRINT *,' RUN COMPLETED, RERUN FOR OTHER VALUES'
STOP
END
SUBROUTINE SORT(IA,IB,I1,I2,I3,JA,JB,MD,ND,NM,NP,N,MAX,MIN
2,ICJ,INM)
DIMENSION JSP(20)
DIMENSION I1(1),I2(1),I3(1),JA(1),JB(1)
DIMENSION IA(1),IB(1),ND(1),MD(INM,4)
I=0
DO 24 K=1,NP
NJK=0
DO 20,J=1,NM
IND=(IA(J)-K)*(IB(J)-K)
IF (IND.NE.0) GO TO 20
NJK=NJK+1
JSP(NJK)=J
20 CONTINUE
MOD=NJK-1
IF (MOD.LE.0) GO TO 24
DO 22 IMD=1,MOD
I=I+1
IF (I.GT.ICJ) GO TO 22
IPD=IMD+1
JAI=JSP(IMD)
JA(I)=JAI
JBI=JSP(IPD)
JB(I)=JBI
I1(I)=IA(JAI)
IF (IA(JAI).EQ.K) I1(I)=IB(JAI)
I2(I)=K
I3(I)=IA(JBI)
IF (IA(JBI).EQ.K) I3(I)=IB(JBI)
22 CONTINUE
24 CONTINUE
N=I
DO 30 J=1,NM
ND(J)=0
DO 30 K=1,4
30 MD(J,K)=0
III=N
IF (N.GT.ICJ) III=ICJ
DO 40 I=1,III
J=JA(I)
DO 38 L=1,2
ND(J)=ND(J)+1
K=1
M=0
32 MJK=MD(J,K)
IF (MJK.NE.0) GO TO 34
M=1
MD(J,K)=I
34 K=K+1
IF (K.GT.4) GO TO 38
IF (M.EQ.0) GO TO 32
38 J=JB(I)
40 CONTINUE
MIN=100
MAX=0
DO 46 J=1,NM
NDJ=ND(J)
IF (NDJ.GT.MAX) MAX=NDJ
46 IF (NDJ.LT.MIN) MIN=NDJ
RETURN
END
SUBROUTINE SGANT(IA,IB,INM,INT,ISC,I1,I2,I3,JA,JB,MD,N,ND,NM,NP
2,AM,BM,C,CGD,CMM,D,EP2,EP3,ETA,FHZ,GAM,SGD,X,Y,Z,ZLD,ZS)
COMPLEX ZG,ZH,ZS,EGD,GD,CGDS,SGDS,SGDT,B01
COMPLEX P11,P12,P21,P22,Q11,Q12,Q21,Q22,EP2,EP,ETA,GAM,EP3
COMPLEX EPSILA,CWEA,BETA,ZARG
COMPLEX P(2,2),Q(2,2),CGD(1),SGD(1),C(1),ZLD(1)
DIMENSION X(1),Y(1),Z(1),D(1),IA(1),IB(1),MD(INM,4)
DIMENSION I1(1),I2(1),I3(1),JA(1),JB(1),ND(1),ISC(1)
DATA E0,TP,U0/8.854E-12,6.28318,1.256E-6/
2 FORMAT (3X,'AM = ',E10.3,3X,'DMAX = ',E10.3,3X,'DMIN = ',E10.3)
EP=EP3
ICC=(N*N+N)/2
DO 10 I=1,ICC
10 C(I)=(.0,.0)
ZS=(.0,.0)
IF (CMM.LE.0.) GO TO 12
OMEGA=TP*FHZ
EPSILA=CMPLX(E0,-CMM*1.E6/OMEGA)
CWEA=(.0,1.)*OMEGA*EPSILA
BETA=OMEGA*SQRT(U0)*CSQRT(EPSILA-EP)
ZARG=BETA*AM
CALL CBES(ZARG,B01)
ZS=BETA*B01/CWEA
12 ZH=ZS/(TP*AM*GAM)
DMIN=1.E30
DMAX=.0
DO 20 J=1,NM
K=IA(J)
L=IB(J)
D(J)=SQRT((X(K)-X(L))**2+(Y(K)-Y(L))**2+(Z(K)-Z(L))**2)
IF (D(J).LT.DMIN) DMIN=D(J)
IF (D(J).GT.DMAX) DMAX=D(J)
EGD= CEXP (GAM*D(J))
CGD(J)=(EGD+1./EGD)/2.
20 SGD(J)=(EGD-1./EGD)/2.
IF (DMIN.LT.2.*AM) GO TO 25
IF (CABS(GAM*AM).GT.0.06) GO TO 25
IF (CABS(GAM*DMAX).GT.3.) GO TO 25
IF (AM.GT.0.) GO TO 30
25 N=0
WRITE (6,*)'SEGMENT VALUES'
WRITE (6,2)AM,DMAX,DMIN
RETURN
30 DO 200 K=1,NM
NDK=ND(K)
KA=IA(K)
KB=IB(K)
DK=D(K)
CGDS=CGD(K)
SGDS=SGD(K)
DO 200 L=1,NM
NDL=ND(L)
LA=IA(L)
LB=IB(L)
DL=D(L)
SGDT=SGD(L)
NIL=0
DO 200 II=1,NDK
I=MD(K,II)
MM=(I-1)*N-(I*I-I)/2
FI=1.
IF (KB.EQ.I2(I)) GO TO 36
IF (KB.EQ.I1(I)) FI=-1.
IS=1
GO TO 40
36 IF (KA.EQ.I3(I)) FI=-1.
IS=2
40 DO 200 JJ=1,NDL
J=MD(L,JJ)
MMM=MM+J
IF (I.GT.J) GO TO 200
FJ=1.
IF (LB.EQ.I2(J)) GO TO 46
IF (LB.EQ.I1(J)) FJ=-1.
JS=1
GOTO 50
46 IF (LA.EQ.I3(J)) FJ=-1.
JS=2
50 IF (NIL.NE.0) GO TO 168
NIL=1
IF (K.EQ.L) GO TO 120
IND=(LA-KA)*(LB-KA)*(LA-KB)*(LB-KB)
IF (IND.EQ.0) GO TO 80
C SEGMENTS K AND L SHARE NO POINTS
CALL GGS(X(KA),Y(KA),Z(KA),X(KB),Y(KB),Z(KB),X(LA),Y(LA),Z(LA)
2,X(LB),Y(LB),Z(LB),AM,DK,CGDS,SGDS,DL,SGDT,INT,ETA,GAM
3,P(1,1),P(1,2),P(2,1),P(2,2))
GO TO 168
C SEGMENTS K AND L SHARE ONE POINT (THEY INTERSECT)
80 KG=0
JM=KB
JC=KA
KF=1
IND=(KB-LA)*(KB-LB)
IF (IND.NE.0) GO TO 82
JC=KB
KF=-1
JM=KA
KG=3
82 LG=3
JP=LA
LF=-1
IF (LB.EQ.JC) GO TO 83
JP=LB
LF=1
LG=0
83 SGN=KF*LF
CPSI=((X(JP)-X(JC))*(X(JM)-X(JC))+(Y(JP)-Y(JC))*(Y(JM)-Y(JC))
2+(Z(JP)-Z(JC))*(Z(JM)-Z(JC)))/(DK*DL)
CALL GGMM(.0,DK,.0,DL,AM,CGDS,SGDS,SGDT,CPSI,ETA,GAM
2,Q(1,1),Q(1,2),Q(2,1),Q(2,2))
DO 98 KK=1,2
KP=IABS(KK-KG)
DO 98 LL=1,2
LP=IABS(LL-LG)
P(KP,LP)=SGN*Q(KK,LL)
98 CONTINUE
GOTO 168
C K=L (SELF REACTION OF SEGMENT K)
120 Q11=(.0,.0)
Q12=(.0,.0)
IF (CMM.LE.0.) GO TO 150
GD=GAM*DK
ZG=ZH/(SGDS**2)
Q11=ZG*(SGDS*CGDS-GD)/2.
Q12=ZG*(GD*CGDS-SGDS)/2.
150 ISCK=ISC(K)
P11=(.0,.0)
P12=(.0,.0)
IF (ISCK.EQ.0) GO TO 155
IF (BM.LE.AM) GO TO 155
CALL DSHELL(AM,BM,DK,CGDS,SGDS,EP2,EP,ETA,GAM,P11,P12)
155 Q11=P11+Q11
Q12=P12+Q12
CALL GGMM(.0,DK,.0,DK,AM,CGDS,SGDS,SGDS,1.
2,ETA,GAM,P11,P12,P21,P22)
Q11=P11+Q11
Q12=P12+Q12
P(1,1)=Q11
P(1,2)=Q12
P(2,1)=Q12
P(2,2)=Q11
IF (KA.NE.LA) GO TO 160
GOTO 168
160 P(1,1)=-Q12
P(1,2)=-Q11
P(2,1)=-Q11
P(2,2)=-Q12
168 C(MMM)=C(MMM)+FI*FJ*P(IS,JS)
200 CONTINUE
DO 220 I=1,N
IJ=(I-1)*N-(I*I-I)/2+I
J1=JA(I)
IF (I2(I).EQ.IB(J1)) J1=J1+NM
J2=JB(I)
IF (I2(I).EQ.IB(J2)) J2=J2+NM
220 C(IJ)=C(IJ)+ZLD(J1)+ZLD(J2)
RETURN
END
SUBROUTINE DSHELL(AM,BM,DK,CGDS,SGDS,EP2,EP,ETA,GAM,P11,P12)
COMPLEX CGDS,SGDS,EP2,EP,ETA,GAM,P11,P12,GD,CST
DATA PI/3.14159/
GD=GAM*DK
CST=(EP2-EP)*ETA*ALOG(BM/AM)/(4.*PI*EP2*SGDS*SGDS)
P11=-CST*(GD+SGDS*CGDS)
P12=CST*(GD*CGDS+SGDS)
RETURN
END
SUBROUTINE CBES(Z,B01)
COMPLEX ARG,CC,CS,EX
COMPLEX B01,Z,TERMJ,TERMN,MZ24,JN(2)
DATA PI/3.14159/
IF (CABS(Z).GE.12.0) GO TO 10
FACTOR=0.0
TERMN=(0.,0.)
MZ24=-0.25*Z*Z
TERMJ=(1.0,0.0)
DO 1 NP=1,2
N=NP-1
JN(NP)=TERMJ
M=0
2 M=M+1
TERMJ=TERMJ*MZ24/FLOAT(M*(N+M))
JN(NP)=JN(NP)+TERMJ
IF (NP.NE.1) GO TO 3
FACTOR=FACTOR+1.0/FLOAT(M)
TERMN=TERMN+TERMJ*FACTOR
3 ERROR=CABS(TERMJ)
IF (ERROR.GT.1.0E-10) GO TO 2
1 TERMJ=0.5*Z
B01=JN(1)/JN(2)
RETURN
10 Y=AIMAG(Z)
IF (ABS(Y).GT.20.) GO TO 20
ARG=(.0,1.)*Z
EX=CEXP(ARG)
CC=EX+1./EX
CS=(.0,-1.)*(EX-1./EX)
B01=(CS+CC)/(CS-CC)
RETURN
20 B01=(.0,-1.0)
IF (Y.LT.0.) B01=(0.,1.)
RETURN
END
SUBROUTINE GGS(XA,YA,ZA,XB,YB,ZB,X1,Y1,Z1,X2,Y2,Z2,AM
2,DS,CGDS,SGDS,DT,SGDT,INT,ETA,GAM,P11,P12,P21,P22)
COMPLEX P11,P12,P21,P22,EJA,EJB,EJ1,EJ2,ETA,GAM,C1,C2,CST
COMPLEX EGD,CGDS,SGDS,SGDT,ER1,ER2,ET1,ET2
DATA FP/12.56637/
CA=(X2-X1)/DT
CB=(Y2-Y1)/DT
CG=(Z2-Z1)/DT
CAS=(XB-XA)/DS
CBS=(YB-YA)/DS
CGS=(ZB-ZA)/DS
CC=CA*CAS+CB*CBS+CG*CGS
IF (ABS(CC).GT..997) GO TO 200
20 SZ=(X1-XA)*CAS+(Y1-YA)*CBS+(Z1-ZA)*CGS
IF (INT.LE.0) GO TO 300
INS=2*(INT/2)
IF (INS.LT.2) INS=2
IP=INS+1
DELT=DT/INS
T=.0
DSZ=CC*DELT
P11=(.0,.0)
P12=(.0,.0)
P21=(.0,.0)
P22=(.0,.0)
AMS=AM*AM
SGN=-1.
DO 100,IN=1,IP
ZZ1=SZ
ZZ2=SZ-DS
XXZ=X1+T*CA-XA-SZ*CAS
YYZ=Y1+T*CB-YA-SZ*CBS
ZZZ=Z1+T*CG-ZA-SZ*CGS
RS=XXZ**2+YYZ**2+ZZZ**2
R1=SQRT(RS+ZZ1**2)
EJA=CEXP(-GAM*R1)
EJ1=EJA/R1
R2=SQRT(RS+ZZ2**2)
EJB=CEXP(-GAM*R2)
EJ2=EJB/R2
ER1=EJA*SGDS+ZZ1*EJ1*CGDS-ZZ2*EJ2
ER2=-EJB*SGDS+ZZ2*EJ2*CGDS-ZZ1*EJ1
FAC=.0
IF (RS.GT.AMS) FAC=(CA*XXZ+CB*YYZ+CG*ZZZ)/RS
ET1=CC*(EJ2-EJ1*CGDS)+FAC*ER1
ET2=CC*(EJ1-EJ2*CGDS)+FAC*ER2
C=3.+SGN
IF (IN.EQ.1 .OR. IN.EQ.IP) C=1.
EGD=CEXP(GAM*(DT-T))
C1=C*(EGD-1./EGD)/2.
EGD=CEXP(GAM*T)
C2=C*(EGD-1./EGD)/2.
P11=P11+ET1*C1
P12=P12+ET1*C2
P21=P21+ET2*C1
P22=P22+ET2*C2
T=T+DELT
SZ=SZ+DSZ
100 SGN=-SGN
CST=-ETA*DELT/(3.*FP*SGDS*SGDT)
P11=CST*P11
P12=CST*P12
P21=CST*P21
P22=CST*P22
RETURN
200 SZ1=(X1-XA)*CAS+(Y1-YA)*CBS+(Z1-ZA)*CGS
RH1=SQRT((X1-XA-SZ1*CAS)**2+(Y1-YA-SZ1*CBS)**2+(Z1-ZA-SZ1*CGS)**2)
SZ2=SZ1+DT*CC
RH2=SQRT((X2-XA-SZ2*CAS)**2+(Y2-YA-SZ2*CBS)**2+(Z2-ZA-SZ2*CGS)**2)
DDD=(RH1+RH2)/2.
IF (DDD.GT.20.*AM .AND. INT.GT.0) GO TO 20
IF (DDD.LT.AM) DDD=AM
CALL GGMM(.0,DS,SZ1,SZ2,DDD,CGDS,SGDS,SGDT,1.
2,ETA,GAM,P11,P12,P21,P22)
RETURN
300 SS=SQRT(1.-CC*CC)
CAD=(CGS*CB-CBS*CG)/SS
CBD=(CAS*CG-CGS*CA)/SS
CGD=(CBS*CA-CAS*CB)/SS
DK=(X1-XA)*CAD+(Y1-YA)*CBD+(Z1-ZA)*CGD
DK=ABS(DK)
IF (DK.LT.AM) DK=AM
XZ=XA+SZ*CAS
YZ=YA+SZ*CBS
ZZ=ZA+SZ*CGS
XP1=X1-DK*CAD
YP1=Y1-DK*CBD
ZP1=Z1-DK*CGD
CAP=CBS*CGD-CGS*CBD
CBP=CGS*CAD-CAS*CGD
CGP=CAS*CBD-CBS*CAD
P1=CAP*(XP1-XZ)+CBP*(YP1-YZ)+CGP*(ZP1-ZZ)
T1=P1/SS
S1=T1*CC-SZ
CALL GGMM(S1,S1+DS,T1,T1+DT,DK,CGDS,SGDS,SGDT,CC,ETA,GAM
2,P11,P12,P21,P22)
RETURN
END
SUBROUTINE GGMM(S1,S2,T1,T2,D,CGDS,SGD1,SGD2,CPSI,ETA,GAM
2,P11,P12,P21,P22)
DOUBLE PRECISION R1,R2,DPQ,SIS,TS1,TS2,ST1,ST2,CD,BD,CPSS,SK
2,TL1,TL2,TD1,TD2,SDI,DPSI,DD,ZD
COMPLEX CGDS,SGDS,SGDT,SGD1,SGD2,ETA,GAM,P11,P12,P21,P22
COMPLEX CST,EB,EC,EK,EL,EKL,EGZI,ES1,ES2,ET1,ET2,EXPA,EXPB
COMPLEX E(2,2),F(2,2)
COMPLEX EGZ(2,2),GM(2),GP(2)
DATA PI/3.14159/
DSQ=D*D
SGDS=SGD1
IF (S2.LT.S1) SGDS=-SGD1
SGDT=SGD2
IF (T2.LT.T1) SGDT=-SGD2
IF (ABS(CPSI).GT..997) GO TO 110
ES1=CEXP(GAM*S1)
ES2=CEXP(GAM*S2)
ET1=CEXP(GAM*T1)
ET2=CEXP(GAM*T2)
DD=D
DPSI=CPSI
TD1=T1
TD2=T2
CPSS=DPSI*DPSI
CD=DD/DSQRT(1.D0-CPSS)
C=CD
BD=CD*DPSI
B=BD
EB=CEXP(GAM*CMPLX(.0,B))
EC=CEXP(GAM*CMPLX(.0,C))
DO 10 K=1,2
DO 10 L=1,2
10 E(K,L)=(.0,.0)
TS1=TD1*TD1
TS2=TD2*TD2
DPQ=DD*DD
SI=S1
DO 100 I=1,2
FI=(-1)**I
SDI=SI
SIS=SDI*SDI
ST1=2.*SDI*TD1*DPSI
ST2=2.*SDI*TD2*DPSI
R1=DSQRT(DPQ+SIS+TS1-ST1)
R2=DSQRT(DPQ+SIS+TS2-ST2)
EK=EB
DO 50 K=1,2
FK=(-1)**K
SK=FK*SDI
EL=EC
DO 40 L=1,2
FL=(-1)**L
EKL=EK*EL
XX=FK*BD+FL*CD
TL1=FL*TD1
TL2=FL*TD2
RR1=R1+SK+TL1
RR2=R2+SK+TL2
CALL EXPJ(GAM*CMPLX(RR1,-XX),GAM*CMPLX(RR2,-XX),EXPA)
CALL EXPJ(GAM*CMPLX(RR1,XX),GAM*CMPLX(RR2,XX),EXPB)
E(K,L)=E(K,L)+FI*(EXPA*EKL+EXPB/EKL)
40 EL=1./EC
50 EK=1./EB
ZD=SDI*DPSI
ZC=ZD
EGZI=CEXP(GAM*ZC)
RR1=R1+ZD-TD1
RR2=R2+ZD-TD2
CALL EXPJ(GAM*RR1,GAM*RR2,EXPB)
RR1=R1-ZD+TD1
RR2=R2-ZD+TD2
CALL EXPJ(GAM*RR1,GAM*RR2,EXPA)
F(I,1)=2.*SGDS*EXPA/EGZI
F(I,2)=2.*SGDS*EXPB*EGZI
100 SI=S2
CST=ETA/(16.*PI*SGDS*SGDT)
P11=CST*(( F(1,1)+E(2,2)*ES2-E(1,2)/ES2)*ET2
A +(-F(1,2)-E(2,1)*ES2+E(1,1)/ES2)/ET2)
P12=CST*((-F(1,1)-E(2,2)*ES2+E(1,2)/ES2)*ET1
B + ( F(1,2)+E(2,1)*ES2-E(1,1)/ES2)/ET1)
P21=CST*((-F(2,1)-E(2,2)*ES1+E(1,2)/ES1)*ET2
C +( F(2,2)+E(2,1)*ES1-E(1,1)/ES1)/ET2)
P22=CST*(( F(2,1)+E(2,2)*ES1-E(1,2)/ES1)*ET1
D +(-F(2,2)-E(2,1)*ES1+E(1,1)/ES1)/ET1)
RETURN
110 IF (CPSI.LT.0.) GO TO 120
TA=T1
TB=T2
GO TO 130
120 TA=-T1
TB=-T2
SGDT=-SGDT
130 SI=S1
DO 150 I=1,2
TJ=TA
DO 140 J=1,2
ZIJ=TJ-SI
R=SQRT(DSQ+ZIJ*ZIJ)
W=R+ZIJ
IF (ZIJ.LT.0.) W=DSQ/(R-ZIJ)
V=R-ZIJ
IF (ZIJ.GT.0.) V=DSQ/(R+ZIJ)
IF (J.EQ.1)V1=V
IF (J.EQ.1)W1=W
EGZ(I,J)=CEXP(GAM*ZIJ)
140 TJ=TB
CALL EXPJ(GAM*V1,GAM*V,GP(I))
CALL EXPJ(GAM*W1,GAM*W,GM(I))
150 SI=S2
CST=-ETA/(8.*PI*SGDS*SGDT)
P11=CST*(GM(2)*EGZ(2,2)+GP(2)/EGZ(2,2)
2-CGDS*(GM(1)*EGZ(1,2)+GP(1)/EGZ(1,2)))
P12=CST*(-GM(2)*EGZ(2,1)-GP(2)/EGZ(2,1)
2+CGDS*(GM(1)*EGZ(1,1)+GP(1)/EGZ(1,1)))
P21=CST*(GM(1)*EGZ(1,2)+GP(1)/EGZ(1,2)
2-CGDS*(GM(2)*EGZ(2,2)+GP(2)/EGZ(2,2)))
P22=CST*(-GM(1)*EGZ(1,1)-GP(1)/EGZ(1,1)
2+CGDS*(GM(2)*EGZ(2,1)+GP(2)/EGZ(2,1)))
RETURN
END
SUBROUTINE EXPJ(V1,V2,W12)
COMPLEX EC,E15,S,T,UC,VC,V1,V2,W12,Z
DIMENSION V(21),W(21),D(16),E(16)
DATA V/ 0.22284667E 00,
20.11889321E 01,0.29927363E 01,0.57751436E 01,0.98374674E 01,
20.15982874E 02,0.93307812E-01,0.49269174E 00,0.12155954E 01,
20.22699495E 01,0.36676227E 01,0.54253366E 01,0.75659162E 01,
20.10120228E 02,0.13130282E 02,0.16654408E 02,0.20776479E 02,
20.25623894E 02,0.31407519E 02,0.38530683E 02,0.48026086E 02/
DATA W/ 0.45896460E 00,
20.41700083E 00,0.11337338E 00,0.10399197E-01,0.26101720E-03,
20.89854791E-06,0.21823487E 00,0.34221017E 00,0.26302758E 00,
20.12642582E 00,0.40206865E-01,0.85638778E-02,0.12124361E-02,
20.11167440E-03,0.64599267E-05,0.22263169E-06,0.42274304E-08,
20.39218973E-10,0.14565152E-12,0.14830270E-15,0.16005949E-19/
DATA D/ 0.22495842E 02,
2 0.74411568E 02,-0.41431576E 03,-0.78754339E 02, 0.11254744E 02,
2 0.16021761E 03,-0.23862195E 03,-0.50094687E 03,-0.68487854E 02,
2 0.12254778E 02,-0.10161976E 02,-0.47219591E 01, 0.79729681E 01,
2-0.21069574E 02, 0.22046490E 01, 0.89728244E 01/
DATA E/ 0.21103107E 02,
2-0.37959787E 03,-0.97489220E 02, 0.12900672E 03, 0.17949226E 02,
2-0.12910931E 03,-0.55705574E 03, 0.13524801E 02, 0.14696721E 03,
2 0.17949528E 02,-0.32981014E 00, 0.31028836E 02, 0.81657657E 01,
2 0.22236961E 02, 0.39124892E 02, 0.81636799E 01/
Z=V1
DO 100 JIM=1,2
X=REAL(Z)
Y=AIMAG(Z)
E15=(.0,.0)
AB=CABS(Z)
IF (AB.EQ.0.) GO TO 90
IF (X.GE.0. .AND. AB.GT.10.) GO TO 80
YA=ABS(Y)
IF (X.LE.0. .AND. YA.GT.10.) GO TO 80
IF(YA-X.GE.17.5.OR.YA.GE.6.5.OR.X+YA.GE.5.5.OR.X.GE.3.)GO TO 20
IF (X.LE.-9.) GO TO 40
IF (YA-X.GE.2.5) GO TO 50
IF (X+YA.GE.1.5) GO TO 30
10 N=6.+3.*AB
E15=1./(N-1.)-Z/N**2
15 N=N-1
E15=1./(N-1.)-Z*E15/N
IF (N.GE.3) GO TO 15
E15=Z*E15-CMPLX(.577216+ALOG(AB),ATAN2(Y,X))
GO TO 90
20 J1=1
J2=6
GOTO 31
30 J1=7
J2=21
31 S=(.0,.0)
YS=Y*Y
DO 32 I=J1,J2
XI=V(I)+X
CF=W(I)/(XI*XI+YS)
32 S=S+CMPLX(XI*CF,-YA*CF)
GO TO 54
40 T3=X*X-Y*Y
T4=2.*X*YA
T5=X*T3-YA*T4
T6=X*T4+YA*T3
UC=CMPLX(D(11)+D(12)*X+D(13)*T3+T5-E(12)*YA-E(13)*T4,
2 E(11)+E(12)*X+E(13)*T3+T6+D(12)*YA+D(13)*T4)
VC=CMPLX(D(14)+D(15)*X+D(16)*T3+T5-E(15)*YA-E(16)*T4,
2 E(14)+E(15)*X+E(16)*T3+T6+D(15)*YA+D(16)*T4)
GO TO 52
50 T3=X*X-Y*Y
T4=2.*X*YA
T5=X*T3-YA*T4
T6=X*T4+YA*T3
T7=X*T5-YA*T6
T8=X*T6+YA*T5
T9=X*T7-YA*T8
T10=X*T8+YA*T7
UC=CMPLX(D(1)+D(2)*X+D(3)*T3+D(4)*T5+D(5)*T7+T9-(E(2)*YA+E(3)*T4
2+E(4)*T6+E(5)*T8),E(1)+E(2)*X+E(3)*T3+E(4)*T5+E(5)*T7+T10+
3(D(2)*YA+D(3)*T4+D(4)*T6+D(5)*T8))
VC=CMPLX(D(6)+D(7)*X+D(8)*T3+D(9)*T5+D(10)*T7+T9-(E(7)*YA+E(8)*T4
2+E(9)*T6+E(10)*T8),E(6)+E(7)*X+E(8)*T3+E(9)*T5+E(10)*T7+T10+
3(D(7)*YA+D(8)*T4+D(9)*T6+D(10)*T8))
52 EC=UC/VC
S=EC/CMPLX(X,YA)
54 EX=EXP(-X)
T=EX*CMPLX(COS(YA),-SIN(YA))
E15=S*T
56 IF (Y.LT.0.) E15=CONJG(E15)
GOTO 90
80 E15=.409319/(Z+.193044)+.421831/(Z+1.02666)+.147126/(Z+2.56788)+
2.206335E-1/(Z+4.90035)+.107401E-2/(Z+8.18215)+.158654E-4/(Z+
312.7342)+.317031E-7/(Z+19.3957)
E15=E15*CEXP(-Z)
90 IF (JIM.EQ.1) W12=E15
100 Z=V2
Z=V2/V1
TH=ATAN2(AIMAG(Z),REAL(Z))-ATAN2(AIMAG(V2),REAL(V2))
2+ATAN2(AIMAG(V1),REAL(V1))
AB=ABS(TH)
IF(AB.LT.1.) TH=.0
IF (TH.GT.1.) TH=6.2831853
IF (TH.LT.-1.) TH=-6.2831853
W12=W12-E15+CMPLX(.0,TH)
RETURN
END
SUBROUTINE GANT1(IA,IB,INM,IWR,I1,I2,I3,I12,JA,JB,MD,N,ND,NM,AM
2,C,CJ,CG,CMM,D,EFF,GAM,GG,CGD,SGD,VG,Y11,Z11,ZLD,ZS)
COMPLEX C(1),CJ(1),CGD(1),SGD(1),VG(1),ZLD(1),Y11,Z11,ZS,GAM,CG(1)
DIMENSION D(1),IA(1),IB(1),JA(1),JB(1)
DIMENSION I1(1),I2(1),I3(1),MD(INM,4),ND(1)
2 FORMAT (1X,I15,8F10.2)
5 FORMAT (1H0)
DO 50 I=1,N
CJ(I)=(.0,.0)
K=JA(I)
DO 40 KK=1,2
KA=IA(K)
KB=IB(K)
JJ=K
FI=1.
IF (KB.EQ.I2(I)) GO TO 36
IF (KB.EQ.I1(I)) FI=-1.
CJ(I)=CJ(I)+FI*VG(JJ)
GO TO 40
36 IF (KA.EQ.I3(I)) FI=-1.
JJ=K+NM
CJ(I)=CJ(I)+FI*VG(JJ)
40 K=JB(I)
50 CONTINUE
DO 55 I=1,N
55 CG(I)=CJ(I)
CALL SQROT(C,CJ,0,I12,N)
I12=2
Y11=(.0,.0)
DO 80 I=1,N
80 Y11=Y11+CJ(I)*CONJG(CG(I))
CALL RITE (IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,CJ,CG)
GG=REAL(Y11)
Z11=1./Y11
PIN=GG
CALL GDISS(AM,CG,CMM,D,DISS,GAM,NM,SGD,ZLD,ZS)
PRAD=PIN-DISS
EFF=100.*PRAD/PIN
RETURN
END
SUBROUTINE SQROT(C,S,IWR,I12,NEQ)
COMPLEX C(1),S(1),SS
2 FORMAT (1X,I15,1F10.3,1F15.7,1F10.0,2F15.6)
3 FORMAT(1H0)
N=NEQ
IF (I12.EQ.2) GO TO 20
C(1)=CSQRT(C(1))
DO 4 K=2,N
4 C(K)=C(K)/C(1)
DO 10 I=2,N
IMO=I-1
IPO=I+1
ID=(I-1)*N-(I*I-I)/2
II=ID+I
DO 5 L=1,IMO
LI=(L-1)*N-(L*L-L)/2+I
5 C(II)=C(II)-C(LI)*C(LI)
C(II)=CSQRT(C(II))
IF (IPO.GT.N) GO TO 10
DO 8 J=IPO,N
IJ=ID+J
DO 6 M=1,IMO
MD=(M-1)*N-(M*M-M)/2
MI=MD+I
MJ=MD+J
6 C(IJ)=C(IJ)-C(MJ)*C(MI)
8 C(IJ)=C(IJ)/C(II)
10 CONTINUE
20 S(1)=S(1)/C(1)
DO 30 I=2,N
IMO=I-1
DO 25 L=1,IMO
LI=(L-1)*N-(L*L-L)/2+I
25 S(I)=S(I)-C(LI)*S(L)
II=(I-1)*N-(I*I-I)/2+I
30 S(I)=S(I)/C(II)
NN=((N+1)*N)/2
S(N)=S(N)/C(NN)
NMO=N-1
DO 40 I=1,NMO
K=N-I
KPO=K+1
KD=(K-1)*N-(K*K-K)/2
DO 35 L=KPO,N
KL=KD+L
35 S(K)=S(K)-C(KL)*S(L)
KK=KD+K
40 S(K)=S(K)/C(KK)
IF (IWR.LE.0) GO TO 100
CNOR=.0
DO 50 I=1,N
SA=CABS(S(I))
50 IF (SA.GT.CNOR) CNOR=SA
IF (CNOR.LE.0.) CNOR=1.
DO 60 I=1,N
SS=S(I)
SA=CABS(SS)
SNOR=SA/CNOR
PH=0.
IF (SA.GT.0.) PH=57.29578*ATAN2(AIMAG(SS),REAL(SS))
60 WRITE(6,2)I,SNOR,SA,PH,SS
WRITE (6,*)'MATRIX VALUES'
WRITE(6,3)
100 RETURN
END
SUBROUTINE RITE(IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,CJ,CG)
COMPLEX CJ(1),CG(1),CJA,CJB
DIMENSION IA(1),IB(1),I1(1),I2(1),I3(1),MD(INM,4),ND(1)
2 FORMAT (1X,1I5,2F10.3,2F10.0,4F15.6)
5 FORMAT (1H0)
AMAX=.0
DO 100 K=1,NM
KA=IA(K)
KB=IB(K)
CJA=(.0,.0)
CJB=(.0,.0)
NDK=ND(K)
DO 40 II=1,NDK
I=MD(K,II)
FI=1.
IF (KB.EQ.I2(I)) GO TO 36
IF (KB.EQ.I1(I)) FI=-1.
CJA=CJA+FI*CJ(I)
GO TO 40
36 IF (KA.EQ.I3(I)) FI=-1.
CJB=CJB+FI*CJ(I)
40 CONTINUE
CG(K)=CJA
KK=K+NM
CG(KK)=CJB
ACJ=CABS(CJA)
BCJ=CABS(CJB)
IF (ACJ.GT.AMAX) AMAX=ACJ
IF (BCJ.GT.AMAX) AMAX=BCJ
100 CONTINUE
IF (IWR.GT.0) GO TO 110
RETURN
110 IF (AMAX.LE.0.) AMAX=1.
WRITE (6,*)'BRANCH CURRENTS'
WRITE (6,*)'SSEG iA iB aA aB cA cB'
DO 200 K=1,NM
CJA=CG(K)
KK=K+NM
CJB=CG(KK)
ACJ=CABS(CJA)/AMAX
BCJ=CABS(CJB)/AMAX
PA=57.29578*ATAN2(AIMAG(CJA),REAL(CJA))
PB=57.29578*ATAN2(AIMAG(CJB),REAL(CJB))
200 WRITE (6,2)K,ACJ,BCJ,PA,PB,CJA,CJB
WRITE (6,5)
RETURN
END
SUBROUTINE GDISS(AM,CG,CMM,D,DISS,GAM,NM,SGD,ZLD,ZS)
COMPLEX CG(1),SGD(1),ZLD(1),CJA,CJB,GAM,ZS
DIMENSION D(1)
DATA PI/3.14159/
DISS=.0
IF (CMM.LE.0.) GO TO 120
ALPH=REAL(GAM)
BETA=AIMAG(GAM)
RH=REAL(ZS)/(4.*PI*AM)
DO 100 K=1,NM
DK=D(K)
DEN=CABS(SGD(K))**2
EAD=EXP(ALPH*DK)
CAD=(EAD+1./EAD)/2.
CBD=COS(BETA*DK)
SAD=DK
IF (ALPH.NE.0.) SAD=(EAD-1./EAD)/(2.*ALPH)
SBD=DK
IF (BETA.NE.0.) SBD=SIN(BETA*DK)/BETA
FA=RH*(SAD*CAD-SBD*CBD)/DEN
FB=2.*RH*(CAD*SBD-SAD*CBD)/DEN
CJA=CG(K)
L=K+NM
CJB=CG(L)
100 DISS=DISS+FA*(CABS(CJA)**2+CABS(CJB)**2)
2+FB*(REAL(CJA)*REAL(CJB)+AIMAG(CJA)*AIMAG(CJB))
120 DO 140 J=1,NM
K=J+NM
140 DISS=DISS+REAL(ZLD(J))*(CABS(CG(J))**2)
2+REAL(ZLD(K))*(CABS(CG(K))**2)
RETURN
END
SUBROUTINE GNFLD(IA,IB,INM,I1,I2,I3,MD,N,ND,NM,AM,CGD,SGD,ETA,GAM
2,CJ,D,X,Y,Z,XP,YP,ZP,EX,EY,EZ)
COMPLEX EX,EY,EZ,EX1,EY1,EZ1,EX2,EY2,EZ2,ETA,GAM
COMPLEX CJ(1),CGD(1),SGD(1)
DIMENSION IA(1),IB(1),I1(1),I2(1),I3(1),D(1),X(1),Y(1),Z(1)
DIMENSION MD(INM,4),ND(1)
DATA PI,TP/3.14159,6.28318/
EX=(.0,.0)
EY=(.0,.0)
EZ=(.0,.0)
DO 140 K=1,NM
KA=IA(K)
KB=IB(K)
CALL GNF(X(KA),Y(KA),Z(KA),X(KB),Y(KB),Z(KB),XP,YP,ZP,AM,D(K)
2,CGD(K),SGD(K),ETA,GAM,EX1,EY1,EZ1,EX2,EY2,EZ2)
NDK=ND(K)
DO 140 II=1,NDK
I=MD(K,II)
FI=1.
IF (KB.EQ.I2(I)) GO TO 136
IF (KB.EQ.I1(I)) FI=-1.
EX=EX+FI*EX1*CJ(I)
EY=EY+FI*EY1*CJ(I)
EZ=EZ+FI*EZ1*CJ(I)
GO TO 140
136 IF (KA.EQ.I3(I)) FI=-1.
EX=EX+FI*EX2*CJ(I)
EY=EY+FI*EY2*CJ(I)
EZ=EZ+FI*EZ2*CJ(I)
140 CONTINUE
RETURN
END
SUBROUTINE GNF(XA,YA,ZA,XB,YB,ZB,X,Y,Z,AM,DS,CGDS,SGDS,ETA,GAM
2,EX1,EY1,EZ1,EX2,EY2,EZ2)
COMPLEX EJA,EJB,EJ1,EJ2,ER1,ER2,ES1,ES2,SGDS,GAM,CST,CGDS,ETA
COMPLEX EX1,EY1,EZ1,EX2,EY2,EZ2
DATA PI/3.14159/
CAS=(XB-XA)/DS
CBS=(YB-YA)/DS
CGS=(ZB-ZA)/DS
SZ=(X-XA)*CAS+(Y-YA)*CBS+(Z-ZA)*CGS
ZZ1=SZ
ZZ2=SZ-DS
XXZ=X-XA-SZ*CAS
YYZ=Y-YA-SZ*CBS
ZZZ=Z-ZA-SZ*CGS
RS=XXZ**2+YYZ**2+ZZZ**2
R1=SQRT(RS+ZZ1**2)
EJA=CEXP(-GAM*R1)
EJ1=EJA/R1
R2=SQRT(RS+ZZ2**2)
EJB=CEXP(-GAM*R2)
EJ2=EJB/R2
ES1=EJ2-EJ1*CGDS
ES2=EJ1-EJ2*CGDS
ER1=(.0,.0)
ER2=(.0,.0)
AMS=AM*AM
IF (RS.LT.AMS) GO TO 80
CTH1=ZZ1/R1
CTH2=ZZ2/R2
ER1=( EJA*SGDS+EJA*CGDS*CTH1-EJB*CTH2)/RS
ER2=(-EJB*SGDS+EJB*CGDS*CTH2-EJA*CTH1)/RS
80 CST=ETA/(4.*PI*SGDS)
EX1=CST*(ES1*CAS+ER1*XXZ)
EY1=CST*(ES1*CBS+ER1*YYZ)
EZ1=CST*(ES1*CGS+ER1*ZZZ)
EX2=CST*(ES2*CAS+ER2*XXZ)
EY2=CST*(ES2*CBS+ER2*YYZ)
EZ2=CST*(ES2*CGS+ER2*XXZ)
RETURN
END
SUBROUTINE GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
4,X,Y,Z,ZLD,ZS,ETA,GAM)
COMPLEX CJI,ET1,ET2,EP1,EP2,EPPS,ETTS,EPTS,ETPS,ZS,VP,VT
COMPLEX C(1),CJ(1),EP(1),ET(1),EPP(1),ETT(1),ZLD(1)
COMPLEX ETA,GAM,CGD(1),SGD(1),CG(1)
DIMENSION IA(1),IB(1),I1(1),I2(1),I3(1),ND(1),MD(INM,4)
DIMENSION D(1),X(1),Y(1),Z(1)
DATA PI,TP/3.14159,6.28318/
CJI=-4.*PI/(ETA*GAM)
GGG=REAL(1./ETA)
THR=.0174533*TH
CTH=COS(THR)
STH=SIN(THR)
PHR=.0174533*PH
CPH=COS(PHR)
SPH=SIN(PHR)
DO 130 I=1,N
ETT(I)=(.0,.0)
130 EPP(I)=(.0,.0)
DO 140 K=1,NM
KA=IA(K)
KB=IB(K)
CALL GFF(X(KA),Y(KA),Z(KA),X(KB),Y(KB),Z(KB),D(K)
2,CGD(K),SGD(K),CTH,STH,CPH,SPH,GAM,ETA,ET1,ET2,EP1,EP2)
NDK=ND(K)
DO 140 II=1,NDK
I=MD(K,II)
FI=1.
IF (KB.EQ.I2(I)) GO TO 136
IF (KB.EQ.I1(I)) FI=-1.
EPP(I)=EPP(I)+FI*EP1
ETT(I)=ETT(I)+FI*ET1
GOTO 140
136 IF (KA.EQ.I3(I)) FI=-1.
EPP(I)=EPP(I)+FI*EP2
ETT(I)=ETT(I)+FI*ET2
140 CONTINUE
EPPS=(.0,.0)
ETTS=(.0,.0)
IF (INC.EQ.0) GO TO 200
IF (INC.EQ.2) GO TO 170
DO 150 I=1,N
ET(I)=ETT(I)*CJI
150 EP(I)=EPP(I)*CJI
CALL SQROT(C,EP,0,I12,N)
I12=2
CALL SQROT(C,ET,0,I12,N)
CALL RITE(IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,EP,CG)
CALL GDISS(AM,CG,CMM,D,PDIS,GAM,NM,SGD,ZLD,ZS)
CALL RITE(IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,ET,CG)
CALL GDISS(AM,CG,CMM,D,TDIS,GAM,NM,SGD,ZLD,ZS)
ACSP=PDIS/GGG
ACST=TDIS/GGG
PIN=.0
TIN=.0
DO 164 I=1,N
VP=CJI*EPP(I)
VT=CJI*ETT(I)
PIN=PIN+REAL(VP*CONJG(EP(I)))
164 TIN=TIN+REAL(VT*CONJG(ET(I)))
ECSP=PIN/GGG
ECST=TIN/GGG
SCSP=ECSP-ACSP
SCST=ECST-ACST
170 EPTS=(.0,.0)
ETPS=(.0,.0)
DO 180 I=1,N
EPPS=EPPS+EP(I)*EPP(I)
EPTS=EPTS+EP(I)*ETT(I)
ETTS=ETTS+ET(I)*ETT(I)
180 ETPS=ETPS+ET(I)*EPP(I)
SPPM=2.*TP*(CABS(EPPS)**2)
SPTM=2.*TP*(CABS(EPTS)**2)
STPM=2.*TP*(CABS(ETPS)**2)
STTM=2.*TP*(CABS(ETTS)**2)
RETURN
200 DO 260 I=1,N
ETTS=ETTS+CJ(I)*ETT(I)
260 EPPS=EPPS+CJ(I)*EPP(I)
APP=CABS(EPPS)
ATT=CABS(ETTS)
GPP=4.*PI*APP*APP*GGG/GG
GTT=4.*PI*ATT*ATT*GGG/GG
RETURN
END
SUBROUTINE GFF(XA,YA,ZA,XB,YB,ZB,D
2,CGD,SGD,CTH,STH,CPH,SPH
2,GAM,ETA,ET1,ET2,EP1,EP2)
COMPLEX ET1,ET2,EP1,EP2,GAM,ETA
COMPLEX GD,CGD,SGD,EGD
COMPLEX EGFA,EGFB,EGGD,ESA,ESB
COMPLEX CST
FP=12.56637
XAB=XB-XA
YAB=YB-YA
ZAB=ZB-ZA
CA=XAB/D
CB=YAB/D
CG=ZAB/D
G=(CA*CPH+CB*SPH)*STH+CG*CTH
GK=1.-G*G
ET1=(.0,.0)
ET2=(.0,.0)
EP1=(.0,.0)
EP2=(.0,.0)
IF (GK.LT..001) GO TO 200
FA=(XA*CPH+YA*SPH)*STH+ZA*CTH
FB=(XB*CPH+YB*SPH)*STH+ZB*CTH
EGFA=CEXP(GAM*FA)
EGFB=CEXP(GAM*FB)
EGGD=CEXP(GAM*G*D)
CST=ETA/(GK*SGD*FP)
ESA=CST*EGFA*(EGGD-G*SGD-CGD)
ESB=CST*EGFB*(1./EGGD+G*SGD-CGD)
T=(CA*CPH+CB*SPH)*CTH-CG*STH
P=-CA*SPH+CB*CPH
ET1=T*ESA
ET2=T*ESB
EP1=P*ESA
EP2=P*ESB
200 CONTINUE
RETURN
END