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 ( 4F15.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,*)'  FLAGS                    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
      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'     
      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)
      WRITE (6,3)PH,TH,GPP,GTT
  310 CONTINUE
  400 IF (ISCAT.LE.0) GO TO 600
      INC=1
      WRITE (6,*)'BACKSCATTER VALUES, ECHO AREAS A--'
      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
  700 PRINT *,' ENTER NUMBER OF ROUTINE TO EXECUTE',CR,LF
      PRINT *,'    1= CHANGE GENERATOR VOLTAGES',CR,LF
      PRINT *,'    2= CHANGE FREQUENCY',CR,LF
      PRINT *,'    3= ADD WIRE INSULATION',CR,LF
      PRINT *,'    4= ADD, CHANGE NEAR FIELD X,Y,Z',CR,LF
      PRINT *,'    5= CHANGE LOADS',CR,LF
      PRINT *,'    6= EXIT THINWIRE',CR,LF,CR,LF
      READ (*,*)IEXIT1
      GOTO (710,720,730,740,750,800), IEXIT1
  710 PRINT *,' INPUT NEW REAL, IMAGINARY VOLTAGE FOR GENERATOR-  '
      DO 711 I=1, NGEN
      PRINT *,'SEGMENT ',I,'  '
      READ (*,* )VG(I)
  711 CONTINUE
      GOTO 63
  720 PRINT *,' INPUT NEW FREQUENCY, MHZ  '
      READ (*,*)FMC
      WRITE (6,*)' NEW FREQUENCY=',FMC
      GOTO 43
  730 PRINT *, ' INPUT INSULATION RADIUS, METERS, 0= NO INSULATION  '
      READ (*,*)BM
      WRITE (6,*)' NEW INSULATION DIAMETER, METERS=',BM
      GOTO 59
  740 PRINT *,' INPUT NEAR FIELD POINT, X,Y,Z, METERS  '
      READ (*,*)XP,YP,ZP
      GOTO 41
  750 PRINT *,' INPUT REAL, IMAGINARY VALUES FOR LOAD-  '
      DO 751 I=1, NGEN
      PRINT *,'SEGMENT ',I
      READ (*,*)ZLD(I)
  751 CONTINUE
  800 CLOSE (5)
      CLOSE (2)
      CLOSE (6)
      WRITE (6,*) ESC,'@'
      STOP 'RUN COMPLETED'
      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

