home *** CD-ROM | disk | FTP | other *** search
/ Antennas / Antennas_CD-ROM_Walnut_Creek_September_1996.iso / thinwire / fortran / tw4.for < prev    next >
Text File  |  1996-06-30  |  35KB  |  1,188 lines
  1. C        MAIN PROGRAM THINWIRE
  2. C         ORIGINAL PROGRAM BY J. H. RICHMOND, 1974
  3. C          MODIFIED 1992, 1994, R. P. HAVILAND
  4.       COMPLEX EP2,EP3,ETA,GAM,Y11,Z11,ZS
  5.       COMPLEX EPPS,EPTS,ETPS,ETTS,EX,EY,EZ
  6.       COMPLEX C(1830),CJ(60),EP(60),ET(60),EPP(60),ETT(60)
  7.       DIMENSION I1(60),I2(60),I3(60),JA(60),JB(60)
  8.       COMPLEX CGD(50),SGD(50),CG(100),VG(100),ZLD(100)
  9.       DIMENSION D(50),IA(50),IB(50),ISC(50),MD(50,4),ND(50)
  10.       DIMENSION X(55),Y(55),Z(55)
  11.       CHARACTER * 31 INFILE
  12.       CHARACTER * 31 OUTFIL
  13.       CHARACTER * 31 NAS
  14.       CHARACTER * 31 DAS
  15.       CHARACTER * 31 FPS
  16.       CHARACTER * 2 CR,LF,SI
  17.       CHARACTER * 3 ESC
  18.       DATA PI,TP/3.14159,6.28318/
  19.       DATA E0,U0/8.854E-12,1.2566E-6/
  20.    2  FORMAT ( 1X,8F15.7)
  21.    3  FORMAT ( 5F15.7)
  22.    4  FORMAT ( 1X,1I5,8F14.6)
  23.    5  FORMAT ( 1H0)
  24.    6  FORMAT ( 1X,6F15.7)
  25.    7  FORMAT ( 8F10.5)
  26.    8  FORMAT ( 1X,1I4,13I5)
  27.    9  FORMAT ( 3X,'MAX = ',I5,3X,'MIN = ',I5,3X,'N = ',I5)
  28. C  10  FORMAT ( 4Z2.0)
  29.       CR=CHAR(13)
  30.       LF=CHAR(10)
  31.       ESC=CHAR(27)
  32.       SI=CHAR(15)
  33.       ICJ=60
  34.       INM=50
  35.       DO 15 J=1,INM
  36.   15  ISC(J)=0
  37.  6020 PRINT *, ' ENTER FILEPATH+FILENAME FOR INPUT, INCLUDE QUOTES  '
  38.       READ (*,*)INFILE
  39. C      SETUP FOR FILE INPUT
  40.       OPEN (UNIT=5, FILE=INFILE,ERR=7020)
  41.       GOTO 6021
  42.  7020 GOTO 6020
  43. C      GET FILEOUT DATA
  44.  6021 PRINT *,' ENTER 1=OUTPUT TO PRINTER, 2=OUTPUT TO FILE  '
  45.       READ (*,*)IFILOU
  46.       IF (IFILOU.EQ. 1) GOTO 6022
  47.       IF (IFILOU.NE. 2) GOTO 6021
  48.  6050 PRINT *,' ENTER FILEPATH+FILENAME FOR OUTPUT, INCLUDE QUOTES  '
  49.       READ (*,*)OUTFIL
  50. C       SETUP FOR FILE OUTPUT
  51.       OPEN (UNIT=6,FILE=OUTFIL,ERR=7050,STATUS='NEW')
  52.       GOTO 6023
  53.  7050 GOTO 6050
  54. C      SETUP FOR PRINTER OUTPUT
  55.  6022 OPEN (UNIT=6,FILE='PRN')
  56.       WRITE (6,*) SI
  57. C      NOW GET DATA
  58.  6023 READ (5,*)NAS,DAS,FPS
  59.       CALL GETDAT(IYR,IMON,IDAY)
  60.       CALL GETTIM(IHR,IMIN)
  61.       WRITE (6,*)'       ',IDAY,IMON,IYR,'  AT   ',IHR,IMIN
  62.       WRITE (6,*) ' THIN-WIRE ANALYSIS PROGRAM '
  63.       WRITE (6,*) '  ANTENNA      ',NAS
  64.       WRITE (6,*) '   OF FILE     ',FPS
  65.       WRITE (6,*) '    DATED      ',DAS
  66.       WRITE (6,*) '    REFERENCE  ',INFILE
  67.       WRITE (6,*) ' '
  68.       READ (5,*)BM,ER2,SIG2,TD2
  69.       WRITE (6,*)'INSULATION DIA.           K           LOSS'
  70.       WRITE (6,2)BM,ER2,SIG2,TD2
  71.       READ (5,*)AM,CMM,ER3,SIG3,TD3
  72.       WRITE (6,*)'WIRE     DIA            RHO              MEDIUM'
  73.       WRITE (6,2)AM,CMM,ER3,SIG3,TD3
  74.       READ (5,*)IBISC,IGAIN,INEAR,ISCAT,IWR,NGEN,NM,NP,NLOAD
  75.       WRITE (6,*)' BIS GAIN NEAR SCAT CUR :GENS SEGS PTS LOADS'
  76.       WRITE (6,8)IBISC,IGAIN,INEAR,ISCAT,IWR,NGEN,NM,NP,NLOAD
  77.       READ (5,*)FMC,PHA,PHAI,PHAN,THA,THAI,THAN
  78.       WRITE (6,*)'FREQUENCY, MHZ'
  79.       WRITE (6,2)FMC
  80.       WRITE (6,*)'THETA INIT., INC., STEPS     PHI INIT., INCR.,  STEPS'
  81.       WRITE (6,*)'FAR FIELD ANGLES'
  82.       WRITE (6,2)PHA,PHAI,PHAN,THA,THAI,THAN
  83.       READ (5,*)PHI,PHII,PHIN,THI,THII,THIN
  84.       WRITE (6,*)'BACK SCATTER ANGLES'
  85.       WRITE (6,2)PHI,PHII,PHIN,THI,THII,THIN
  86.       READ (5,*)PHS,PHSI,PHSN,THS,THSI,THSN
  87.       WRITE (6,*)'BISTATIC SCATTER ANGLES'
  88.       WRITE (6,2)PHS,PHSI,PHSN,THS,THSI,THSN
  89.       WRITE (6,*)'SEGMENT, END1, END2'
  90.       DO 22 J=1,NM
  91.       READ (5,*)IA(J),IB(J)
  92.   22  WRITE (6,8)J,IA(J),IB(J)
  93.       WRITE (6,*)'POINT       LOCATIONS, XYZ  METERS'
  94.       DO 40 I=1,NP
  95.       READ (5,*)X(I),Y(I),Z(I)
  96.   40  WRITE (6,4)I,X(I),Y(I),Z(I)
  97.   41  WRITE (6,*)' NEAR FIELD   X Y Z'
  98.       READ (5,*) XP,YP,ZP
  99.   42  WRITE (6,*) XP,YP,ZP
  100.   43  FHZ=FMC*1.E6
  101.       OMEGA=TP*FHZ
  102.       DO 44 I=1,NP
  103.       IF (BM .GT. AM) ISC(I)=1
  104.   44  CONTINUE 
  105.       IF (SIG2.LT.0.) EP2=ER2*E0*CMPLX(1.,-TD2)
  106.       IF (TD2.LT.0.) EP2=CMPLX(ER2*E0,-SIG2/OMEGA)
  107.       IF (SIG3.LT.0.) EP3=ER3*E0*CMPLX(1.,-TD3)
  108.       IF (TD3.LT.0.) EP3=CMPLX(ER3*E0,-SIG3/OMEGA)
  109.       ETA=CSQRT(U0/EP3)
  110.       GAM=OMEGA*CSQRT(-U0*EP3)
  111.   59  CALL SORT(IA,IB,I1,I2,I3,JA,JB,MD,ND,NM,NP,N,MAX,MIN,ICJ,INM)
  112.       WRITE (6,*)'   SEGMENT COUPLING '
  113.       WRITE (6,9)MAX,MIN,N
  114.       WRITE (6,*)' '
  115.       IF (MIN.LT.1) WRITE (6,*)'***    UNCONNECTED WIRE   ***'
  116.       IF (MAX.GT.4 .OR. MIN.LT.1 .OR. N.GT.ICJ) GO TO 800
  117.       INT=4
  118.       I12=1
  119.       DO 60 J=1,NM
  120.       VG(J)=(.0,.0)
  121.       ZLD(J)=(.0,.0)
  122.       JJ=J+NM
  123.       VG(JJ)=(.0,.0)
  124.   60  ZLD(JJ)=(.0,.0)
  125.       WRITE (6,*)'LOAD           LR, LI                     POINT#'
  126.       NFIN=NLOAD
  127.       IF (NFIN .EQ. 0) NFIN=1
  128.       DO 61 I=1,NFIN
  129.       READ (5,*)LOADR,LOADI,JLOAD
  130.       ZLD(JLOAD)=CMPLX(LOADR,LOADI)
  131.   61  WRITE (6,*)ZLD(I),JLOAD
  132.       WRITE (6,*)'GENERATOR      VR, VI                     POINT#'
  133.       NFIN=NGEN
  134.       IF (NFIN .EQ. 0) NFIN=1
  135.       DO 62 I=1,NFIN
  136.       READ (5,*)GENR,GENI,JGEN
  137.       JPT=JGEN
  138.       VG(JPT)=CMPLX(GENR,GENI)
  139.   62  WRITE (6,*)VG(JPT),JGEN
  140.       PRINT *, 'SET INTEGRATION, 0=NO CHANGE, -1=EXACT'
  141.       PRINT *, ' <4=FAST, LOW, >4=SLOW, HIGH ACCURACY'
  142.       READ (*,*) INTSET
  143.       IF (INTSET .LE. -1) INT=0
  144.       IF (INTSET .GT. 0) INT=INTSET
  145.       WRITE (6,*) 'INT VARIABLE =',INT      
  146.       CALL SGANT(IA,IB,INM,INT,ISC,I1,I2,I3,JA,JB,MD,N,ND,NM,NP
  147.      2,AM,BM,C,CGD,CMM,D,EP2,EP3,ETA,FHZ,GAM,SGD,X,Y,Z,ZLD,ZS)
  148.       IF (N.LE.0) GO TO 800
  149.       IF (NGEN.LE.0) GO TO 400
  150.   63  CALL GANT1(IA,IB,INM,IWR,I1,I2,I3,I12,JA,JB,MD,N,ND,NM,AM
  151.      2,C,CJ,CG,CMM,D,EFF,GAM,GG,CGD,SGD,VG,Y11,Z11,ZLD,ZS)
  152.       WRITE (6,*)'EFFICIENCY             INPOWER         ZIN'
  153.       WRITE (6,3)EFF,GG,Z11
  154.   200 IF (INEAR.LE.0) GO TO 300
  155.       CALL  GNFLD(IA,IB,INM,I1,I2,I3,MD,N,ND,NM,AM,CGD,SGD,ETA,GAM
  156.      2,CJ,D,X,Y,Z,XP,YP,ZP,EX,EY,EZ)
  157.       WRITE (6,*)'NEAR FIELD VALUES'
  158.       WRITE (6,*)'X      Y      Z'
  159.       WRITE (6,3)XP,YP,ZP
  160.       WRITE (6,*)'ELECTRIC FIELD'
  161.       WRITE (6,3)EX,EY,EZ
  162.   300 IF (IGAIN.LE.0) GO TO 400
  163.       INC=0
  164.       WRITE (6,*)'FAR FIELD VALUES'
  165.       WRITE (6,*)'AZMUITH   ELEVATION   H-FIELD   V-FIELD  TOTAL'     
  166.       DO 310 I=PHA,PHAI*PHAN,PHAI
  167.       DO 310 J=THA,THAI*THAN,THAI
  168.       PH=I
  169.       TH=J
  170.       CALL   GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
  171.      2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
  172.      3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
  173.      4,X,Y,Z,ZLD,ZS,ETA,GAM)
  174.       TOT=GPP*GPP+GTT*GTT
  175.       IF (TOT .LE. 0) DTOT=-999
  176.       IF (TOT .GT. 0) DTOT=5*LOG10(TOT)
  177.       IF (TOT .GT. 1E35) DTOT=999
  178.       WRITE (6,3)PH,TH,GPP,GTT,DTOT
  179.   310 CONTINUE
  180.   400 IF (ISCAT.LE.0) GO TO 600
  181.       INC=1
  182.       WRITE (6,*)'BACKSCATTER VALUES,   ECHO AREAS,  CROSS SECTIONS'
  183.       WRITE (6,*)'AZMUITH   ELEVATION   APP    APT    ATP    ATT'
  184.       WRITE (6,*)'ABSORB P,T   EXTINCT P,T   SCATTER P,T'
  185.       DO 410 I=PHI,PHII*PHIN,PHII
  186.       DO 410  J=THI,THII*THIN,THII
  187.       PH=I
  188.       TH=J
  189.       CALL   GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
  190.      2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
  191.      3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
  192.      4,X,Y,Z,ZLD,ZS,ETA,GAM)
  193.       WRITE (6,6)PH,TH,SPPM,SPTM,STPM,STTM
  194.       WRITE (6,6)ACSP,ACST,ECSP,ECST,SCSP,SCST
  195.   410 CONTINUE
  196.   500 IF (IBISC.LE.0) GO TO 600
  197.       INC=2
  198.       WRITE (6,*)'BI-STAT SCATTER  VALUES'
  199.       WRITE (6,*)' WITH INCIDENT WAVE PHI, THETA OF'
  200.       WRITE (6,*) PHI, THI
  201.       WRITE (6,*)'AZMUITH     ELEVATION     H-FIELD     V-FIELD'
  202.       DO 510 I=PHS,PHSI*PHSN,PHSI
  203.       DO 510 J=THS,THSI*THSN,THSI
  204.       PH=I
  205.       TH=J    
  206.       CALL   GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
  207.      2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
  208.      3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
  209.      4,X,Y,Z,ZLD,ZS,ETA,GAM)
  210.       WRITE (6,6)PH,TH,SPPM,SPTM,STPM,STTM
  211.   510 CONTINUE
  212.   600 CONTINUE
  213.   800 CLOSE (5)
  214.       CLOSE (2)
  215.       CLOSE (6)
  216.       WRITE (6,*) ESC,'@'
  217.       PRINT *,'    RUN COMPLETED, RERUN FOR OTHER VALUES'
  218.       STOP 
  219.       END
  220.       SUBROUTINE SORT(IA,IB,I1,I2,I3,JA,JB,MD,ND,NM,NP,N,MAX,MIN
  221.      2,ICJ,INM)
  222.       DIMENSION JSP(20)
  223.       DIMENSION I1(1),I2(1),I3(1),JA(1),JB(1)
  224.       DIMENSION IA(1),IB(1),ND(1),MD(INM,4)
  225.       I=0
  226.       DO 24 K=1,NP
  227.       NJK=0
  228.       DO 20,J=1,NM
  229.       IND=(IA(J)-K)*(IB(J)-K)
  230.       IF (IND.NE.0) GO TO 20
  231.       NJK=NJK+1
  232.       JSP(NJK)=J
  233.    20 CONTINUE
  234.       MOD=NJK-1
  235.       IF (MOD.LE.0) GO TO 24
  236.       DO 22 IMD=1,MOD
  237.       I=I+1
  238.       IF (I.GT.ICJ) GO TO 22
  239.       IPD=IMD+1
  240.       JAI=JSP(IMD)
  241.       JA(I)=JAI
  242.       JBI=JSP(IPD)
  243.       JB(I)=JBI
  244.       I1(I)=IA(JAI)
  245.       IF (IA(JAI).EQ.K) I1(I)=IB(JAI)
  246.       I2(I)=K
  247.       I3(I)=IA(JBI)
  248.       IF (IA(JBI).EQ.K) I3(I)=IB(JBI)
  249.    22 CONTINUE
  250.    24 CONTINUE
  251.       N=I
  252.       DO 30 J=1,NM
  253.       ND(J)=0
  254.       DO 30 K=1,4
  255.    30 MD(J,K)=0
  256.       III=N
  257.       IF (N.GT.ICJ) III=ICJ
  258.       DO 40 I=1,III
  259.       J=JA(I)
  260.       DO 38 L=1,2
  261.       ND(J)=ND(J)+1
  262.       K=1
  263.       M=0
  264.    32 MJK=MD(J,K)
  265.       IF (MJK.NE.0) GO TO 34
  266.       M=1
  267.       MD(J,K)=I
  268.    34 K=K+1
  269.       IF (K.GT.4) GO TO 38
  270.       IF (M.EQ.0) GO TO 32
  271.    38 J=JB(I)
  272.    40 CONTINUE
  273.       MIN=100
  274.       MAX=0
  275.       DO 46 J=1,NM
  276.       NDJ=ND(J)
  277.       IF (NDJ.GT.MAX) MAX=NDJ
  278.    46 IF (NDJ.LT.MIN) MIN=NDJ
  279.       RETURN
  280.       END
  281.       SUBROUTINE SGANT(IA,IB,INM,INT,ISC,I1,I2,I3,JA,JB,MD,N,ND,NM,NP
  282.      2,AM,BM,C,CGD,CMM,D,EP2,EP3,ETA,FHZ,GAM,SGD,X,Y,Z,ZLD,ZS)
  283.       COMPLEX ZG,ZH,ZS,EGD,GD,CGDS,SGDS,SGDT,B01
  284.       COMPLEX P11,P12,P21,P22,Q11,Q12,Q21,Q22,EP2,EP,ETA,GAM,EP3
  285.       COMPLEX EPSILA,CWEA,BETA,ZARG
  286.       COMPLEX P(2,2),Q(2,2),CGD(1),SGD(1),C(1),ZLD(1)
  287.       DIMENSION X(1),Y(1),Z(1),D(1),IA(1),IB(1),MD(INM,4)
  288.       DIMENSION I1(1),I2(1),I3(1),JA(1),JB(1),ND(1),ISC(1)
  289.       DATA E0,TP,U0/8.854E-12,6.28318,1.256E-6/
  290.     2 FORMAT (3X,'AM = ',E10.3,3X,'DMAX = ',E10.3,3X,'DMIN = ',E10.3)
  291.       EP=EP3
  292.       ICC=(N*N+N)/2
  293.       DO 10 I=1,ICC
  294.    10 C(I)=(.0,.0)
  295.       ZS=(.0,.0)
  296.       IF (CMM.LE.0.) GO TO 12
  297.       OMEGA=TP*FHZ
  298.       EPSILA=CMPLX(E0,-CMM*1.E6/OMEGA)
  299.       CWEA=(.0,1.)*OMEGA*EPSILA
  300.       BETA=OMEGA*SQRT(U0)*CSQRT(EPSILA-EP)
  301.       ZARG=BETA*AM
  302.       CALL CBES(ZARG,B01)
  303.       ZS=BETA*B01/CWEA
  304.    12 ZH=ZS/(TP*AM*GAM)
  305.       DMIN=1.E30
  306.       DMAX=.0
  307.       DO 20 J=1,NM
  308.       K=IA(J)
  309.       L=IB(J)
  310.       D(J)=SQRT((X(K)-X(L))**2+(Y(K)-Y(L))**2+(Z(K)-Z(L))**2) 
  311.       IF (D(J).LT.DMIN) DMIN=D(J)
  312.       IF (D(J).GT.DMAX) DMAX=D(J)
  313.       EGD= CEXP (GAM*D(J))
  314.       CGD(J)=(EGD+1./EGD)/2.
  315.    20 SGD(J)=(EGD-1./EGD)/2.
  316.       IF (DMIN.LT.2.*AM) GO TO 25
  317.       IF (CABS(GAM*AM).GT.0.06) GO TO 25
  318.       IF (CABS(GAM*DMAX).GT.3.) GO TO 25                 
  319.       IF (AM.GT.0.) GO TO 30
  320.    25 N=0
  321.       WRITE (6,*)'SEGMENT VALUES'
  322.       WRITE (6,2)AM,DMAX,DMIN
  323.       RETURN
  324.    30 DO 200 K=1,NM
  325.       NDK=ND(K)
  326.       KA=IA(K)
  327.       KB=IB(K)
  328.       DK=D(K)
  329.       CGDS=CGD(K)
  330.       SGDS=SGD(K)
  331.       DO 200 L=1,NM
  332.       NDL=ND(L)
  333.       LA=IA(L)
  334.       LB=IB(L)
  335.       DL=D(L)
  336.       SGDT=SGD(L)
  337.       NIL=0
  338.       DO 200 II=1,NDK
  339.       I=MD(K,II)
  340.       MM=(I-1)*N-(I*I-I)/2
  341.       FI=1.
  342.       IF (KB.EQ.I2(I)) GO TO 36
  343.       IF (KB.EQ.I1(I)) FI=-1.
  344.       IS=1
  345.       GO TO 40
  346.    36 IF (KA.EQ.I3(I)) FI=-1.
  347.       IS=2
  348.    40 DO 200 JJ=1,NDL
  349.       J=MD(L,JJ)
  350.       MMM=MM+J
  351.       IF (I.GT.J) GO TO 200
  352.       FJ=1.
  353.       IF (LB.EQ.I2(J)) GO TO 46
  354.       IF (LB.EQ.I1(J)) FJ=-1.
  355.       JS=1
  356.       GOTO 50
  357.    46 IF (LA.EQ.I3(J)) FJ=-1.
  358.       JS=2
  359.    50 IF (NIL.NE.0) GO TO 168
  360.       NIL=1
  361.       IF (K.EQ.L) GO TO 120
  362.       IND=(LA-KA)*(LB-KA)*(LA-KB)*(LB-KB)
  363.       IF (IND.EQ.0) GO TO 80
  364. C     SEGMENTS K AND L SHARE NO POINTS
  365.       CALL GGS(X(KA),Y(KA),Z(KA),X(KB),Y(KB),Z(KB),X(LA),Y(LA),Z(LA)
  366.      2,X(LB),Y(LB),Z(LB),AM,DK,CGDS,SGDS,DL,SGDT,INT,ETA,GAM
  367.      3,P(1,1),P(1,2),P(2,1),P(2,2))
  368.       GO TO 168
  369. C     SEGMENTS K AND L SHARE ONE POINT (THEY INTERSECT)
  370.    80 KG=0
  371.       JM=KB
  372.       JC=KA
  373.       KF=1
  374.       IND=(KB-LA)*(KB-LB)
  375.       IF (IND.NE.0) GO TO 82
  376.       JC=KB
  377.       KF=-1
  378.       JM=KA
  379.       KG=3
  380.    82 LG=3
  381.       JP=LA
  382.       LF=-1
  383.       IF (LB.EQ.JC) GO TO 83
  384.       JP=LB
  385.       LF=1
  386.       LG=0
  387.    83 SGN=KF*LF
  388.       CPSI=((X(JP)-X(JC))*(X(JM)-X(JC))+(Y(JP)-Y(JC))*(Y(JM)-Y(JC))
  389.      2+(Z(JP)-Z(JC))*(Z(JM)-Z(JC)))/(DK*DL)
  390.       CALL GGMM(.0,DK,.0,DL,AM,CGDS,SGDS,SGDT,CPSI,ETA,GAM
  391.      2,Q(1,1),Q(1,2),Q(2,1),Q(2,2))
  392.       DO 98 KK=1,2
  393.       KP=IABS(KK-KG)
  394.       DO 98 LL=1,2
  395.       LP=IABS(LL-LG)
  396.       P(KP,LP)=SGN*Q(KK,LL)
  397.    98 CONTINUE
  398.       GOTO 168
  399. C     K=L (SELF REACTION OF SEGMENT K)
  400.   120 Q11=(.0,.0)
  401.       Q12=(.0,.0)
  402.       IF (CMM.LE.0.) GO TO 150
  403.       GD=GAM*DK
  404.       ZG=ZH/(SGDS**2)
  405.       Q11=ZG*(SGDS*CGDS-GD)/2.
  406.       Q12=ZG*(GD*CGDS-SGDS)/2.
  407.   150 ISCK=ISC(K)
  408.       P11=(.0,.0)
  409.       P12=(.0,.0)
  410.       IF (ISCK.EQ.0) GO TO 155
  411.       IF (BM.LE.AM) GO TO 155
  412.       CALL DSHELL(AM,BM,DK,CGDS,SGDS,EP2,EP,ETA,GAM,P11,P12)
  413.   155 Q11=P11+Q11
  414.       Q12=P12+Q12
  415.       CALL GGMM(.0,DK,.0,DK,AM,CGDS,SGDS,SGDS,1.
  416.      2,ETA,GAM,P11,P12,P21,P22)
  417.       Q11=P11+Q11
  418.       Q12=P12+Q12
  419.       P(1,1)=Q11
  420.       P(1,2)=Q12
  421.       P(2,1)=Q12
  422.       P(2,2)=Q11
  423.       IF (KA.NE.LA) GO TO 160
  424.       GOTO 168
  425.   160 P(1,1)=-Q12
  426.       P(1,2)=-Q11
  427.       P(2,1)=-Q11
  428.       P(2,2)=-Q12
  429.   168 C(MMM)=C(MMM)+FI*FJ*P(IS,JS)
  430.   200 CONTINUE
  431.       DO 220 I=1,N
  432.       IJ=(I-1)*N-(I*I-I)/2+I
  433.       J1=JA(I)
  434.       IF (I2(I).EQ.IB(J1)) J1=J1+NM
  435.       J2=JB(I)
  436.       IF (I2(I).EQ.IB(J2)) J2=J2+NM
  437.   220 C(IJ)=C(IJ)+ZLD(J1)+ZLD(J2)
  438.       RETURN
  439.       END
  440.       SUBROUTINE DSHELL(AM,BM,DK,CGDS,SGDS,EP2,EP,ETA,GAM,P11,P12)
  441.       COMPLEX CGDS,SGDS,EP2,EP,ETA,GAM,P11,P12,GD,CST
  442.       DATA PI/3.14159/
  443.       GD=GAM*DK
  444.       CST=(EP2-EP)*ETA*ALOG(BM/AM)/(4.*PI*EP2*SGDS*SGDS)
  445.       P11=-CST*(GD+SGDS*CGDS)
  446.       P12=CST*(GD*CGDS+SGDS)
  447.       RETURN
  448.       END
  449.       SUBROUTINE CBES(Z,B01)
  450.       COMPLEX ARG,CC,CS,EX
  451.       COMPLEX B01,Z,TERMJ,TERMN,MZ24,JN(2)
  452.       DATA PI/3.14159/
  453.       IF (CABS(Z).GE.12.0) GO TO 10
  454.       FACTOR=0.0
  455.       TERMN=(0.,0.)
  456.       MZ24=-0.25*Z*Z
  457.       TERMJ=(1.0,0.0)
  458.       DO 1 NP=1,2
  459.       N=NP-1
  460.       JN(NP)=TERMJ
  461.       M=0
  462.     2 M=M+1
  463.       TERMJ=TERMJ*MZ24/FLOAT(M*(N+M))
  464.       JN(NP)=JN(NP)+TERMJ
  465.       IF (NP.NE.1) GO TO 3
  466.       FACTOR=FACTOR+1.0/FLOAT(M)
  467.       TERMN=TERMN+TERMJ*FACTOR
  468.     3 ERROR=CABS(TERMJ)
  469.       IF (ERROR.GT.1.0E-10) GO TO 2
  470.     1 TERMJ=0.5*Z
  471.       B01=JN(1)/JN(2)
  472.       RETURN
  473.    10 Y=AIMAG(Z)
  474.       IF (ABS(Y).GT.20.) GO TO 20
  475.       ARG=(.0,1.)*Z
  476.       EX=CEXP(ARG)
  477.       CC=EX+1./EX
  478.       CS=(.0,-1.)*(EX-1./EX)
  479.       B01=(CS+CC)/(CS-CC)  
  480.       RETURN
  481.    20 B01=(.0,-1.0)
  482.       IF (Y.LT.0.) B01=(0.,1.)
  483.       RETURN
  484.       END
  485.       SUBROUTINE GGS(XA,YA,ZA,XB,YB,ZB,X1,Y1,Z1,X2,Y2,Z2,AM   
  486.      2,DS,CGDS,SGDS,DT,SGDT,INT,ETA,GAM,P11,P12,P21,P22)
  487.       COMPLEX P11,P12,P21,P22,EJA,EJB,EJ1,EJ2,ETA,GAM,C1,C2,CST
  488.       COMPLEX EGD,CGDS,SGDS,SGDT,ER1,ER2,ET1,ET2
  489.       DATA FP/12.56637/
  490.       CA=(X2-X1)/DT
  491.       CB=(Y2-Y1)/DT
  492.       CG=(Z2-Z1)/DT
  493.       CAS=(XB-XA)/DS
  494.       CBS=(YB-YA)/DS
  495.       CGS=(ZB-ZA)/DS
  496.       CC=CA*CAS+CB*CBS+CG*CGS
  497.       IF (ABS(CC).GT..997) GO TO 200
  498.    20 SZ=(X1-XA)*CAS+(Y1-YA)*CBS+(Z1-ZA)*CGS
  499.       IF (INT.LE.0) GO TO 300
  500.       INS=2*(INT/2)
  501.       IF (INS.LT.2) INS=2
  502.       IP=INS+1
  503.       DELT=DT/INS
  504.       T=.0
  505.       DSZ=CC*DELT
  506.       P11=(.0,.0)
  507.       P12=(.0,.0)
  508.       P21=(.0,.0)
  509.       P22=(.0,.0)
  510.       AMS=AM*AM
  511.       SGN=-1.
  512.       DO 100,IN=1,IP
  513.       ZZ1=SZ
  514.       ZZ2=SZ-DS
  515.       XXZ=X1+T*CA-XA-SZ*CAS
  516.       YYZ=Y1+T*CB-YA-SZ*CBS
  517.       ZZZ=Z1+T*CG-ZA-SZ*CGS
  518.       RS=XXZ**2+YYZ**2+ZZZ**2
  519.       R1=SQRT(RS+ZZ1**2)
  520.       EJA=CEXP(-GAM*R1)
  521.       EJ1=EJA/R1
  522.       R2=SQRT(RS+ZZ2**2)
  523.       EJB=CEXP(-GAM*R2)
  524.       EJ2=EJB/R2
  525.       ER1=EJA*SGDS+ZZ1*EJ1*CGDS-ZZ2*EJ2
  526.       ER2=-EJB*SGDS+ZZ2*EJ2*CGDS-ZZ1*EJ1
  527.       FAC=.0
  528.       IF (RS.GT.AMS) FAC=(CA*XXZ+CB*YYZ+CG*ZZZ)/RS
  529.       ET1=CC*(EJ2-EJ1*CGDS)+FAC*ER1
  530.       ET2=CC*(EJ1-EJ2*CGDS)+FAC*ER2
  531.       C=3.+SGN
  532.       IF (IN.EQ.1 .OR. IN.EQ.IP) C=1.
  533.       EGD=CEXP(GAM*(DT-T))
  534.       C1=C*(EGD-1./EGD)/2.
  535.       EGD=CEXP(GAM*T)
  536.       C2=C*(EGD-1./EGD)/2.
  537.       P11=P11+ET1*C1
  538.       P12=P12+ET1*C2
  539.       P21=P21+ET2*C1
  540.       P22=P22+ET2*C2
  541.       T=T+DELT
  542.       SZ=SZ+DSZ
  543.   100 SGN=-SGN
  544.       CST=-ETA*DELT/(3.*FP*SGDS*SGDT)
  545.       P11=CST*P11
  546.       P12=CST*P12
  547.       P21=CST*P21
  548.       P22=CST*P22
  549.       RETURN
  550.   200 SZ1=(X1-XA)*CAS+(Y1-YA)*CBS+(Z1-ZA)*CGS
  551.       RH1=SQRT((X1-XA-SZ1*CAS)**2+(Y1-YA-SZ1*CBS)**2+(Z1-ZA-SZ1*CGS)**2)
  552.       SZ2=SZ1+DT*CC
  553.       RH2=SQRT((X2-XA-SZ2*CAS)**2+(Y2-YA-SZ2*CBS)**2+(Z2-ZA-SZ2*CGS)**2)
  554.       DDD=(RH1+RH2)/2.
  555.       IF (DDD.GT.20.*AM .AND. INT.GT.0) GO TO 20
  556.       IF (DDD.LT.AM) DDD=AM
  557.       CALL GGMM(.0,DS,SZ1,SZ2,DDD,CGDS,SGDS,SGDT,1.
  558.      2,ETA,GAM,P11,P12,P21,P22)
  559.       RETURN
  560.   300 SS=SQRT(1.-CC*CC)
  561.       CAD=(CGS*CB-CBS*CG)/SS
  562.       CBD=(CAS*CG-CGS*CA)/SS
  563.       CGD=(CBS*CA-CAS*CB)/SS
  564.       DK=(X1-XA)*CAD+(Y1-YA)*CBD+(Z1-ZA)*CGD
  565.       DK=ABS(DK)
  566.       IF (DK.LT.AM) DK=AM
  567.       XZ=XA+SZ*CAS
  568.       YZ=YA+SZ*CBS
  569.       ZZ=ZA+SZ*CGS
  570.       XP1=X1-DK*CAD
  571.       YP1=Y1-DK*CBD
  572.       ZP1=Z1-DK*CGD
  573.       CAP=CBS*CGD-CGS*CBD
  574.       CBP=CGS*CAD-CAS*CGD
  575.       CGP=CAS*CBD-CBS*CAD
  576.       P1=CAP*(XP1-XZ)+CBP*(YP1-YZ)+CGP*(ZP1-ZZ)
  577.       T1=P1/SS
  578.       S1=T1*CC-SZ
  579.       CALL GGMM(S1,S1+DS,T1,T1+DT,DK,CGDS,SGDS,SGDT,CC,ETA,GAM
  580.      2,P11,P12,P21,P22)
  581.       RETURN
  582.       END
  583.       SUBROUTINE GGMM(S1,S2,T1,T2,D,CGDS,SGD1,SGD2,CPSI,ETA,GAM
  584.      2,P11,P12,P21,P22)
  585.       DOUBLE PRECISION R1,R2,DPQ,SIS,TS1,TS2,ST1,ST2,CD,BD,CPSS,SK
  586.      2,TL1,TL2,TD1,TD2,SDI,DPSI,DD,ZD
  587.       COMPLEX CGDS,SGDS,SGDT,SGD1,SGD2,ETA,GAM,P11,P12,P21,P22
  588.       COMPLEX CST,EB,EC,EK,EL,EKL,EGZI,ES1,ES2,ET1,ET2,EXPA,EXPB
  589.       COMPLEX E(2,2),F(2,2)
  590.       COMPLEX EGZ(2,2),GM(2),GP(2)
  591.       DATA PI/3.14159/
  592.       DSQ=D*D
  593.       SGDS=SGD1
  594.       IF (S2.LT.S1) SGDS=-SGD1
  595.       SGDT=SGD2
  596.       IF (T2.LT.T1) SGDT=-SGD2
  597.       IF (ABS(CPSI).GT..997) GO TO 110
  598.       ES1=CEXP(GAM*S1)
  599.       ES2=CEXP(GAM*S2)
  600.       ET1=CEXP(GAM*T1)
  601.       ET2=CEXP(GAM*T2)
  602.       DD=D
  603.       DPSI=CPSI
  604.       TD1=T1
  605.       TD2=T2
  606.       CPSS=DPSI*DPSI
  607.       CD=DD/DSQRT(1.D0-CPSS)
  608.       C=CD
  609.       BD=CD*DPSI
  610.       B=BD
  611.       EB=CEXP(GAM*CMPLX(.0,B))
  612.       EC=CEXP(GAM*CMPLX(.0,C))
  613.       DO 10 K=1,2
  614.       DO 10 L=1,2
  615.    10 E(K,L)=(.0,.0)
  616.       TS1=TD1*TD1
  617.       TS2=TD2*TD2
  618.       DPQ=DD*DD
  619.       SI=S1
  620.       DO 100 I=1,2
  621.       FI=(-1)**I
  622.       SDI=SI
  623.       SIS=SDI*SDI
  624.       ST1=2.*SDI*TD1*DPSI
  625.       ST2=2.*SDI*TD2*DPSI
  626.       R1=DSQRT(DPQ+SIS+TS1-ST1)
  627.       R2=DSQRT(DPQ+SIS+TS2-ST2)
  628.       EK=EB
  629.       DO 50 K=1,2
  630.       FK=(-1)**K
  631.       SK=FK*SDI
  632.       EL=EC
  633.       DO 40 L=1,2
  634.       FL=(-1)**L
  635.       EKL=EK*EL
  636.       XX=FK*BD+FL*CD
  637.       TL1=FL*TD1
  638.       TL2=FL*TD2
  639.       RR1=R1+SK+TL1
  640.       RR2=R2+SK+TL2
  641.       CALL EXPJ(GAM*CMPLX(RR1,-XX),GAM*CMPLX(RR2,-XX),EXPA)
  642.       CALL EXPJ(GAM*CMPLX(RR1,XX),GAM*CMPLX(RR2,XX),EXPB)
  643.       E(K,L)=E(K,L)+FI*(EXPA*EKL+EXPB/EKL)
  644.    40 EL=1./EC
  645.    50 EK=1./EB
  646.       ZD=SDI*DPSI
  647.       ZC=ZD
  648.       EGZI=CEXP(GAM*ZC)
  649.       RR1=R1+ZD-TD1
  650.       RR2=R2+ZD-TD2
  651.       CALL EXPJ(GAM*RR1,GAM*RR2,EXPB)
  652.       RR1=R1-ZD+TD1
  653.       RR2=R2-ZD+TD2
  654.       CALL EXPJ(GAM*RR1,GAM*RR2,EXPA)
  655.       F(I,1)=2.*SGDS*EXPA/EGZI
  656.       F(I,2)=2.*SGDS*EXPB*EGZI
  657.   100 SI=S2
  658.       CST=ETA/(16.*PI*SGDS*SGDT)
  659.       P11=CST*(( F(1,1)+E(2,2)*ES2-E(1,2)/ES2)*ET2
  660.      A        +(-F(1,2)-E(2,1)*ES2+E(1,1)/ES2)/ET2)
  661.       P12=CST*((-F(1,1)-E(2,2)*ES2+E(1,2)/ES2)*ET1
  662.      B       + ( F(1,2)+E(2,1)*ES2-E(1,1)/ES2)/ET1)
  663.       P21=CST*((-F(2,1)-E(2,2)*ES1+E(1,2)/ES1)*ET2
  664.      C        +( F(2,2)+E(2,1)*ES1-E(1,1)/ES1)/ET2)
  665.       P22=CST*(( F(2,1)+E(2,2)*ES1-E(1,2)/ES1)*ET1
  666.      D        +(-F(2,2)-E(2,1)*ES1+E(1,1)/ES1)/ET1)
  667.       RETURN
  668.   110 IF (CPSI.LT.0.) GO TO 120
  669.       TA=T1
  670.       TB=T2
  671.       GO TO 130
  672.   120 TA=-T1
  673.       TB=-T2
  674.       SGDT=-SGDT
  675.   130 SI=S1
  676.       DO 150 I=1,2
  677.       TJ=TA
  678.       DO 140 J=1,2
  679.       ZIJ=TJ-SI
  680.       R=SQRT(DSQ+ZIJ*ZIJ)
  681.       W=R+ZIJ
  682.       IF (ZIJ.LT.0.) W=DSQ/(R-ZIJ)
  683.       V=R-ZIJ
  684.       IF (ZIJ.GT.0.) V=DSQ/(R+ZIJ)
  685.       IF (J.EQ.1)V1=V
  686.       IF (J.EQ.1)W1=W
  687.       EGZ(I,J)=CEXP(GAM*ZIJ)
  688.   140 TJ=TB
  689.       CALL EXPJ(GAM*V1,GAM*V,GP(I))
  690.       CALL EXPJ(GAM*W1,GAM*W,GM(I))
  691.   150 SI=S2
  692.       CST=-ETA/(8.*PI*SGDS*SGDT)
  693.       P11=CST*(GM(2)*EGZ(2,2)+GP(2)/EGZ(2,2)
  694.      2-CGDS*(GM(1)*EGZ(1,2)+GP(1)/EGZ(1,2)))
  695.       P12=CST*(-GM(2)*EGZ(2,1)-GP(2)/EGZ(2,1)
  696.      2+CGDS*(GM(1)*EGZ(1,1)+GP(1)/EGZ(1,1)))
  697.       P21=CST*(GM(1)*EGZ(1,2)+GP(1)/EGZ(1,2)
  698.      2-CGDS*(GM(2)*EGZ(2,2)+GP(2)/EGZ(2,2)))
  699.       P22=CST*(-GM(1)*EGZ(1,1)-GP(1)/EGZ(1,1)
  700.      2+CGDS*(GM(2)*EGZ(2,1)+GP(2)/EGZ(2,1)))
  701.       RETURN
  702.       END
  703.       SUBROUTINE EXPJ(V1,V2,W12)
  704.       COMPLEX EC,E15,S,T,UC,VC,V1,V2,W12,Z
  705.       DIMENSION V(21),W(21),D(16),E(16)
  706.       DATA V/   0.22284667E 00,
  707.      20.11889321E 01,0.29927363E 01,0.57751436E 01,0.98374674E 01,
  708.      20.15982874E 02,0.93307812E-01,0.49269174E 00,0.12155954E 01,
  709.      20.22699495E 01,0.36676227E 01,0.54253366E 01,0.75659162E 01,
  710.      20.10120228E 02,0.13130282E 02,0.16654408E 02,0.20776479E 02,
  711.      20.25623894E 02,0.31407519E 02,0.38530683E 02,0.48026086E 02/
  712.       DATA W/   0.45896460E 00,
  713.      20.41700083E 00,0.11337338E 00,0.10399197E-01,0.26101720E-03,
  714.      20.89854791E-06,0.21823487E 00,0.34221017E 00,0.26302758E 00,
  715.      20.12642582E 00,0.40206865E-01,0.85638778E-02,0.12124361E-02,
  716.      20.11167440E-03,0.64599267E-05,0.22263169E-06,0.42274304E-08,
  717.      20.39218973E-10,0.14565152E-12,0.14830270E-15,0.16005949E-19/
  718.       DATA D/   0.22495842E 02,
  719.      2 0.74411568E 02,-0.41431576E 03,-0.78754339E 02, 0.11254744E 02,
  720.      2 0.16021761E 03,-0.23862195E 03,-0.50094687E 03,-0.68487854E 02,
  721.      2 0.12254778E 02,-0.10161976E 02,-0.47219591E 01, 0.79729681E 01,
  722.      2-0.21069574E 02, 0.22046490E 01, 0.89728244E 01/
  723.       DATA E/   0.21103107E 02,
  724.      2-0.37959787E 03,-0.97489220E 02, 0.12900672E 03, 0.17949226E 02,
  725.      2-0.12910931E 03,-0.55705574E 03, 0.13524801E 02, 0.14696721E 03,
  726.      2 0.17949528E 02,-0.32981014E 00, 0.31028836E 02, 0.81657657E 01,
  727.      2 0.22236961E 02, 0.39124892E 02, 0.81636799E 01/
  728.       Z=V1
  729.       DO 100 JIM=1,2
  730.       X=REAL(Z)
  731.       Y=AIMAG(Z)
  732.       E15=(.0,.0)
  733.       AB=CABS(Z)
  734.       IF (AB.EQ.0.) GO TO 90
  735.       IF (X.GE.0. .AND. AB.GT.10.) GO TO 80
  736.       YA=ABS(Y)
  737.       IF (X.LE.0. .AND. YA.GT.10.) GO TO 80
  738.       IF(YA-X.GE.17.5.OR.YA.GE.6.5.OR.X+YA.GE.5.5.OR.X.GE.3.)GO TO 20
  739.       IF (X.LE.-9.) GO TO 40
  740.       IF (YA-X.GE.2.5) GO TO 50
  741.       IF (X+YA.GE.1.5) GO TO 30
  742.    10 N=6.+3.*AB
  743.       E15=1./(N-1.)-Z/N**2
  744.    15 N=N-1
  745.       E15=1./(N-1.)-Z*E15/N
  746.       IF (N.GE.3) GO TO 15
  747.       E15=Z*E15-CMPLX(.577216+ALOG(AB),ATAN2(Y,X))
  748.       GO TO 90
  749.    20 J1=1
  750.       J2=6
  751.       GOTO 31
  752.    30 J1=7
  753.       J2=21
  754.    31 S=(.0,.0)
  755.       YS=Y*Y
  756.       DO 32 I=J1,J2
  757.       XI=V(I)+X
  758.       CF=W(I)/(XI*XI+YS)
  759.    32 S=S+CMPLX(XI*CF,-YA*CF)
  760.       GO TO 54
  761.    40 T3=X*X-Y*Y
  762.       T4=2.*X*YA
  763.       T5=X*T3-YA*T4
  764.       T6=X*T4+YA*T3
  765.       UC=CMPLX(D(11)+D(12)*X+D(13)*T3+T5-E(12)*YA-E(13)*T4,
  766.      2         E(11)+E(12)*X+E(13)*T3+T6+D(12)*YA+D(13)*T4)
  767.       VC=CMPLX(D(14)+D(15)*X+D(16)*T3+T5-E(15)*YA-E(16)*T4,
  768.      2         E(14)+E(15)*X+E(16)*T3+T6+D(15)*YA+D(16)*T4)
  769.       GO TO 52
  770.    50 T3=X*X-Y*Y
  771.       T4=2.*X*YA
  772.       T5=X*T3-YA*T4
  773.       T6=X*T4+YA*T3
  774.       T7=X*T5-YA*T6
  775.       T8=X*T6+YA*T5
  776.       T9=X*T7-YA*T8
  777.       T10=X*T8+YA*T7
  778.       UC=CMPLX(D(1)+D(2)*X+D(3)*T3+D(4)*T5+D(5)*T7+T9-(E(2)*YA+E(3)*T4
  779.      2+E(4)*T6+E(5)*T8),E(1)+E(2)*X+E(3)*T3+E(4)*T5+E(5)*T7+T10+
  780.      3(D(2)*YA+D(3)*T4+D(4)*T6+D(5)*T8))
  781.       VC=CMPLX(D(6)+D(7)*X+D(8)*T3+D(9)*T5+D(10)*T7+T9-(E(7)*YA+E(8)*T4
  782.      2+E(9)*T6+E(10)*T8),E(6)+E(7)*X+E(8)*T3+E(9)*T5+E(10)*T7+T10+
  783.      3(D(7)*YA+D(8)*T4+D(9)*T6+D(10)*T8))
  784.    52 EC=UC/VC
  785.       S=EC/CMPLX(X,YA)
  786.    54 EX=EXP(-X)
  787.       T=EX*CMPLX(COS(YA),-SIN(YA))
  788.       E15=S*T
  789.    56 IF (Y.LT.0.) E15=CONJG(E15)
  790.       GOTO 90
  791.    80 E15=.409319/(Z+.193044)+.421831/(Z+1.02666)+.147126/(Z+2.56788)+
  792.      2.206335E-1/(Z+4.90035)+.107401E-2/(Z+8.18215)+.158654E-4/(Z+
  793.      312.7342)+.317031E-7/(Z+19.3957)
  794.       E15=E15*CEXP(-Z)
  795.    90 IF (JIM.EQ.1) W12=E15
  796.   100 Z=V2
  797.       Z=V2/V1
  798.       TH=ATAN2(AIMAG(Z),REAL(Z))-ATAN2(AIMAG(V2),REAL(V2))
  799.      2+ATAN2(AIMAG(V1),REAL(V1))
  800.       AB=ABS(TH)
  801.       IF(AB.LT.1.) TH=.0
  802.       IF (TH.GT.1.) TH=6.2831853
  803.       IF (TH.LT.-1.) TH=-6.2831853
  804.       W12=W12-E15+CMPLX(.0,TH)
  805.       RETURN
  806.       END
  807.       SUBROUTINE GANT1(IA,IB,INM,IWR,I1,I2,I3,I12,JA,JB,MD,N,ND,NM,AM
  808.      2,C,CJ,CG,CMM,D,EFF,GAM,GG,CGD,SGD,VG,Y11,Z11,ZLD,ZS)
  809.       COMPLEX C(1),CJ(1),CGD(1),SGD(1),VG(1),ZLD(1),Y11,Z11,ZS,GAM,CG(1)
  810.       DIMENSION D(1),IA(1),IB(1),JA(1),JB(1)
  811.       DIMENSION I1(1),I2(1),I3(1),MD(INM,4),ND(1)
  812.     2 FORMAT (1X,I15,8F10.2)
  813.     5 FORMAT (1H0)
  814.       DO 50 I=1,N
  815.       CJ(I)=(.0,.0)
  816.       K=JA(I)
  817.       DO 40 KK=1,2
  818.       KA=IA(K)
  819.       KB=IB(K)
  820.       JJ=K
  821.       FI=1.
  822.       IF (KB.EQ.I2(I)) GO TO 36
  823.       IF (KB.EQ.I1(I)) FI=-1.
  824.       CJ(I)=CJ(I)+FI*VG(JJ)
  825.       GO TO 40
  826.    36 IF (KA.EQ.I3(I)) FI=-1.
  827.       JJ=K+NM
  828.       CJ(I)=CJ(I)+FI*VG(JJ)
  829.    40 K=JB(I)
  830.    50 CONTINUE
  831.       DO 55 I=1,N
  832.    55 CG(I)=CJ(I)
  833.       CALL SQROT(C,CJ,0,I12,N)
  834.       I12=2
  835.       Y11=(.0,.0)
  836.       DO 80 I=1,N
  837.    80 Y11=Y11+CJ(I)*CONJG(CG(I))
  838.       CALL RITE (IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,CJ,CG)
  839.       GG=REAL(Y11)
  840.       Z11=1./Y11
  841.       PIN=GG
  842.       CALL       GDISS(AM,CG,CMM,D,DISS,GAM,NM,SGD,ZLD,ZS)
  843.       PRAD=PIN-DISS
  844.       EFF=100.*PRAD/PIN
  845.       RETURN
  846.       END
  847.       SUBROUTINE SQROT(C,S,IWR,I12,NEQ)
  848.       COMPLEX C(1),S(1),SS
  849.     2 FORMAT (1X,I15,1F10.3,1F15.7,1F10.0,2F15.6)
  850.     3 FORMAT(1H0)
  851.       N=NEQ
  852.       IF (I12.EQ.2) GO TO 20
  853.       C(1)=CSQRT(C(1))
  854.       DO 4 K=2,N
  855.     4 C(K)=C(K)/C(1)
  856.       DO 10 I=2,N
  857.       IMO=I-1
  858.       IPO=I+1
  859.       ID=(I-1)*N-(I*I-I)/2
  860.       II=ID+I
  861.       DO 5 L=1,IMO
  862.       LI=(L-1)*N-(L*L-L)/2+I
  863.     5 C(II)=C(II)-C(LI)*C(LI)
  864.       C(II)=CSQRT(C(II))
  865.       IF (IPO.GT.N) GO TO 10
  866.       DO 8 J=IPO,N
  867.       IJ=ID+J
  868.       DO 6 M=1,IMO
  869.       MD=(M-1)*N-(M*M-M)/2
  870.       MI=MD+I
  871.       MJ=MD+J
  872.     6 C(IJ)=C(IJ)-C(MJ)*C(MI)
  873.     8 C(IJ)=C(IJ)/C(II)
  874.    10 CONTINUE
  875.    20 S(1)=S(1)/C(1)
  876.       DO 30 I=2,N
  877.       IMO=I-1
  878.       DO 25 L=1,IMO
  879.       LI=(L-1)*N-(L*L-L)/2+I
  880.    25 S(I)=S(I)-C(LI)*S(L)
  881.       II=(I-1)*N-(I*I-I)/2+I
  882.    30 S(I)=S(I)/C(II)
  883.       NN=((N+1)*N)/2
  884.       S(N)=S(N)/C(NN)
  885.       NMO=N-1
  886.       DO 40 I=1,NMO
  887.       K=N-I
  888.       KPO=K+1
  889.       KD=(K-1)*N-(K*K-K)/2
  890.       DO 35 L=KPO,N
  891.       KL=KD+L
  892.    35 S(K)=S(K)-C(KL)*S(L)
  893.       KK=KD+K
  894.    40 S(K)=S(K)/C(KK)
  895.       IF (IWR.LE.0) GO TO 100
  896.       CNOR=.0
  897.       DO 50 I=1,N
  898.       SA=CABS(S(I))
  899.    50 IF (SA.GT.CNOR) CNOR=SA
  900.       IF (CNOR.LE.0.) CNOR=1.
  901.       DO 60 I=1,N
  902.       SS=S(I)
  903.       SA=CABS(SS)
  904.       SNOR=SA/CNOR
  905.       PH=0.
  906.       IF (SA.GT.0.) PH=57.29578*ATAN2(AIMAG(SS),REAL(SS))
  907.    60 WRITE(6,2)I,SNOR,SA,PH,SS
  908.       WRITE (6,*)'MATRIX VALUES'
  909.       WRITE(6,3)
  910.   100 RETURN
  911.       END
  912.       SUBROUTINE RITE(IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,CJ,CG)
  913.       COMPLEX CJ(1),CG(1),CJA,CJB
  914.       DIMENSION IA(1),IB(1),I1(1),I2(1),I3(1),MD(INM,4),ND(1)
  915.     2 FORMAT (1X,1I5,2F10.3,2F10.0,4F15.6)
  916.     5 FORMAT (1H0)
  917.       AMAX=.0
  918.       DO 100 K=1,NM
  919.       KA=IA(K)
  920.       KB=IB(K)
  921.       CJA=(.0,.0)
  922.       CJB=(.0,.0)
  923.       NDK=ND(K)
  924.       DO 40 II=1,NDK
  925.       I=MD(K,II)
  926.       FI=1.
  927.       IF (KB.EQ.I2(I)) GO TO 36
  928.       IF (KB.EQ.I1(I)) FI=-1.
  929.       CJA=CJA+FI*CJ(I)
  930.       GO TO 40
  931.    36 IF (KA.EQ.I3(I)) FI=-1.
  932.       CJB=CJB+FI*CJ(I)
  933.    40 CONTINUE
  934.       CG(K)=CJA
  935.       KK=K+NM
  936.       CG(KK)=CJB
  937.       ACJ=CABS(CJA)
  938.       BCJ=CABS(CJB)
  939.       IF (ACJ.GT.AMAX) AMAX=ACJ
  940.       IF (BCJ.GT.AMAX) AMAX=BCJ
  941.   100 CONTINUE
  942.       IF (IWR.GT.0) GO TO 110
  943.       RETURN
  944.   110 IF (AMAX.LE.0.) AMAX=1.
  945.       WRITE (6,*)'BRANCH CURRENTS'
  946.       WRITE (6,*)'SSEG iA  iB  aA  aB  cA      cB'
  947.       DO 200 K=1,NM
  948.       CJA=CG(K)
  949.       KK=K+NM
  950.       CJB=CG(KK)
  951.       ACJ=CABS(CJA)/AMAX
  952.       BCJ=CABS(CJB)/AMAX
  953.       PA=57.29578*ATAN2(AIMAG(CJA),REAL(CJA))
  954.       PB=57.29578*ATAN2(AIMAG(CJB),REAL(CJB))
  955.   200 WRITE (6,2)K,ACJ,BCJ,PA,PB,CJA,CJB
  956.       WRITE (6,5)
  957.       RETURN
  958.       END
  959.       SUBROUTINE GDISS(AM,CG,CMM,D,DISS,GAM,NM,SGD,ZLD,ZS)
  960.       COMPLEX CG(1),SGD(1),ZLD(1),CJA,CJB,GAM,ZS
  961.       DIMENSION D(1)
  962.       DATA PI/3.14159/
  963.       DISS=.0
  964.       IF (CMM.LE.0.) GO TO 120
  965.       ALPH=REAL(GAM)
  966.       BETA=AIMAG(GAM)
  967.       RH=REAL(ZS)/(4.*PI*AM)
  968.       DO 100 K=1,NM
  969.       DK=D(K)
  970.       DEN=CABS(SGD(K))**2
  971.       EAD=EXP(ALPH*DK)
  972.       CAD=(EAD+1./EAD)/2.
  973.       CBD=COS(BETA*DK)
  974.       SAD=DK
  975.       IF (ALPH.NE.0.) SAD=(EAD-1./EAD)/(2.*ALPH)
  976.       SBD=DK
  977.       IF (BETA.NE.0.) SBD=SIN(BETA*DK)/BETA
  978.       FA=RH*(SAD*CAD-SBD*CBD)/DEN
  979.       FB=2.*RH*(CAD*SBD-SAD*CBD)/DEN
  980.       CJA=CG(K)
  981.       L=K+NM
  982.       CJB=CG(L)
  983.   100 DISS=DISS+FA*(CABS(CJA)**2+CABS(CJB)**2)
  984.      2+FB*(REAL(CJA)*REAL(CJB)+AIMAG(CJA)*AIMAG(CJB))
  985.   120 DO 140 J=1,NM
  986.       K=J+NM
  987.   140 DISS=DISS+REAL(ZLD(J))*(CABS(CG(J))**2)
  988.      2+REAL(ZLD(K))*(CABS(CG(K))**2)
  989.       RETURN
  990.       END
  991.       SUBROUTINE GNFLD(IA,IB,INM,I1,I2,I3,MD,N,ND,NM,AM,CGD,SGD,ETA,GAM
  992.      2,CJ,D,X,Y,Z,XP,YP,ZP,EX,EY,EZ)
  993.       COMPLEX EX,EY,EZ,EX1,EY1,EZ1,EX2,EY2,EZ2,ETA,GAM
  994.       COMPLEX CJ(1),CGD(1),SGD(1)
  995.       DIMENSION IA(1),IB(1),I1(1),I2(1),I3(1),D(1),X(1),Y(1),Z(1)
  996.       DIMENSION MD(INM,4),ND(1)
  997.       DATA PI,TP/3.14159,6.28318/
  998.       EX=(.0,.0)
  999.       EY=(.0,.0)
  1000.       EZ=(.0,.0)
  1001.       DO 140 K=1,NM
  1002.       KA=IA(K)
  1003.       KB=IB(K)
  1004.       CALL GNF(X(KA),Y(KA),Z(KA),X(KB),Y(KB),Z(KB),XP,YP,ZP,AM,D(K)
  1005.      2,CGD(K),SGD(K),ETA,GAM,EX1,EY1,EZ1,EX2,EY2,EZ2)
  1006.       NDK=ND(K)
  1007.       DO 140 II=1,NDK
  1008.       I=MD(K,II)
  1009.       FI=1.
  1010.       IF (KB.EQ.I2(I)) GO TO 136
  1011.       IF (KB.EQ.I1(I)) FI=-1.
  1012.       EX=EX+FI*EX1*CJ(I)
  1013.       EY=EY+FI*EY1*CJ(I)
  1014.       EZ=EZ+FI*EZ1*CJ(I)
  1015.       GO TO 140
  1016.   136 IF (KA.EQ.I3(I)) FI=-1.
  1017.       EX=EX+FI*EX2*CJ(I)
  1018.       EY=EY+FI*EY2*CJ(I)
  1019.       EZ=EZ+FI*EZ2*CJ(I)
  1020.   140 CONTINUE
  1021.       RETURN
  1022.       END     
  1023.       SUBROUTINE GNF(XA,YA,ZA,XB,YB,ZB,X,Y,Z,AM,DS,CGDS,SGDS,ETA,GAM
  1024.      2,EX1,EY1,EZ1,EX2,EY2,EZ2)
  1025.       COMPLEX EJA,EJB,EJ1,EJ2,ER1,ER2,ES1,ES2,SGDS,GAM,CST,CGDS,ETA
  1026.       COMPLEX EX1,EY1,EZ1,EX2,EY2,EZ2
  1027.       DATA PI/3.14159/
  1028.       CAS=(XB-XA)/DS
  1029.       CBS=(YB-YA)/DS
  1030.       CGS=(ZB-ZA)/DS
  1031.       SZ=(X-XA)*CAS+(Y-YA)*CBS+(Z-ZA)*CGS
  1032.       ZZ1=SZ
  1033.       ZZ2=SZ-DS
  1034.       XXZ=X-XA-SZ*CAS
  1035.       YYZ=Y-YA-SZ*CBS
  1036.       ZZZ=Z-ZA-SZ*CGS
  1037.       RS=XXZ**2+YYZ**2+ZZZ**2
  1038.       R1=SQRT(RS+ZZ1**2)
  1039.       EJA=CEXP(-GAM*R1)
  1040.       EJ1=EJA/R1
  1041.       R2=SQRT(RS+ZZ2**2)
  1042.       EJB=CEXP(-GAM*R2)
  1043.       EJ2=EJB/R2
  1044.       ES1=EJ2-EJ1*CGDS
  1045.       ES2=EJ1-EJ2*CGDS
  1046.       ER1=(.0,.0)
  1047.       ER2=(.0,.0)
  1048.       AMS=AM*AM
  1049.       IF (RS.LT.AMS) GO TO 80
  1050.       CTH1=ZZ1/R1
  1051.       CTH2=ZZ2/R2
  1052.       ER1=( EJA*SGDS+EJA*CGDS*CTH1-EJB*CTH2)/RS
  1053.       ER2=(-EJB*SGDS+EJB*CGDS*CTH2-EJA*CTH1)/RS
  1054.    80 CST=ETA/(4.*PI*SGDS)
  1055.       EX1=CST*(ES1*CAS+ER1*XXZ)
  1056.       EY1=CST*(ES1*CBS+ER1*YYZ)      
  1057.       EZ1=CST*(ES1*CGS+ER1*ZZZ)
  1058.       EX2=CST*(ES2*CAS+ER2*XXZ)
  1059.       EY2=CST*(ES2*CBS+ER2*YYZ)
  1060.       EZ2=CST*(ES2*CGS+ER2*XXZ)
  1061.       RETURN
  1062.       END
  1063.       SUBROUTINE GFFLD(IA,IB,INC,INM,IWR,I1,I2,I3,I12,MD,N,ND,NM,AM
  1064.      2,ACSP,ACST,C,CGD,CG,CJ,CMM,D,ECSP,ECST,EP,ET,EPP,ETT,EPPS,EPTS
  1065.      3,ETPS,ETTS,GG,GPP,GTT,PH,SGD,SCSP,SCST,SPPM,SPTM,STPM,STTM,TH
  1066.      4,X,Y,Z,ZLD,ZS,ETA,GAM)
  1067.       COMPLEX CJI,ET1,ET2,EP1,EP2,EPPS,ETTS,EPTS,ETPS,ZS,VP,VT
  1068.       COMPLEX C(1),CJ(1),EP(1),ET(1),EPP(1),ETT(1),ZLD(1)
  1069.       COMPLEX ETA,GAM,CGD(1),SGD(1),CG(1)
  1070.       DIMENSION IA(1),IB(1),I1(1),I2(1),I3(1),ND(1),MD(INM,4)
  1071.       DIMENSION D(1),X(1),Y(1),Z(1)
  1072.       DATA PI,TP/3.14159,6.28318/
  1073.       CJI=-4.*PI/(ETA*GAM)
  1074.       GGG=REAL(1./ETA)
  1075.       THR=.0174533*TH
  1076.       CTH=COS(THR)
  1077.       STH=SIN(THR)
  1078.       PHR=.0174533*PH
  1079.       CPH=COS(PHR)
  1080.       SPH=SIN(PHR)
  1081.       DO 130 I=1,N
  1082.       ETT(I)=(.0,.0)
  1083.   130 EPP(I)=(.0,.0)
  1084.       DO 140 K=1,NM
  1085.       KA=IA(K)
  1086.       KB=IB(K)
  1087.       CALL GFF(X(KA),Y(KA),Z(KA),X(KB),Y(KB),Z(KB),D(K)
  1088.      2,CGD(K),SGD(K),CTH,STH,CPH,SPH,GAM,ETA,ET1,ET2,EP1,EP2)
  1089.       NDK=ND(K)
  1090.       DO 140 II=1,NDK
  1091.       I=MD(K,II)
  1092.       FI=1.
  1093.       IF (KB.EQ.I2(I)) GO TO 136
  1094.       IF (KB.EQ.I1(I)) FI=-1.
  1095.       EPP(I)=EPP(I)+FI*EP1
  1096.       ETT(I)=ETT(I)+FI*ET1
  1097.       GOTO 140
  1098.   136 IF (KA.EQ.I3(I)) FI=-1.
  1099.       EPP(I)=EPP(I)+FI*EP2
  1100.       ETT(I)=ETT(I)+FI*ET2
  1101.   140 CONTINUE
  1102.       EPPS=(.0,.0)
  1103.       ETTS=(.0,.0)
  1104.       IF (INC.EQ.0) GO TO 200
  1105.       IF (INC.EQ.2) GO TO 170
  1106.       DO 150 I=1,N
  1107.       ET(I)=ETT(I)*CJI
  1108.   150 EP(I)=EPP(I)*CJI
  1109.       CALL SQROT(C,EP,0,I12,N)
  1110.       I12=2
  1111.       CALL SQROT(C,ET,0,I12,N)
  1112.       CALL RITE(IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,EP,CG)
  1113.       CALL GDISS(AM,CG,CMM,D,PDIS,GAM,NM,SGD,ZLD,ZS)
  1114.       CALL RITE(IA,IB,INM,IWR,I1,I2,I3,MD,ND,NM,ET,CG)
  1115.       CALL GDISS(AM,CG,CMM,D,TDIS,GAM,NM,SGD,ZLD,ZS)
  1116.       ACSP=PDIS/GGG
  1117.       ACST=TDIS/GGG
  1118.       PIN=.0
  1119.       TIN=.0
  1120.       DO 164 I=1,N
  1121.       VP=CJI*EPP(I)
  1122.       VT=CJI*ETT(I)
  1123.       PIN=PIN+REAL(VP*CONJG(EP(I)))
  1124.   164 TIN=TIN+REAL(VT*CONJG(ET(I)))
  1125.       ECSP=PIN/GGG
  1126.       ECST=TIN/GGG
  1127.       SCSP=ECSP-ACSP
  1128.       SCST=ECST-ACST
  1129.   170 EPTS=(.0,.0)
  1130.       ETPS=(.0,.0)
  1131.       DO 180 I=1,N
  1132.       EPPS=EPPS+EP(I)*EPP(I)
  1133.       EPTS=EPTS+EP(I)*ETT(I)
  1134.       ETTS=ETTS+ET(I)*ETT(I)
  1135.   180 ETPS=ETPS+ET(I)*EPP(I)
  1136.       SPPM=2.*TP*(CABS(EPPS)**2)
  1137.       SPTM=2.*TP*(CABS(EPTS)**2)
  1138.       STPM=2.*TP*(CABS(ETPS)**2)
  1139.       STTM=2.*TP*(CABS(ETTS)**2)
  1140.       RETURN
  1141.   200 DO 260 I=1,N
  1142.       ETTS=ETTS+CJ(I)*ETT(I)
  1143.   260 EPPS=EPPS+CJ(I)*EPP(I)
  1144.       APP=CABS(EPPS)
  1145.       ATT=CABS(ETTS)
  1146.       GPP=4.*PI*APP*APP*GGG/GG
  1147.       GTT=4.*PI*ATT*ATT*GGG/GG
  1148.       RETURN
  1149.       END
  1150.       SUBROUTINE GFF(XA,YA,ZA,XB,YB,ZB,D
  1151.      2,CGD,SGD,CTH,STH,CPH,SPH
  1152.      2,GAM,ETA,ET1,ET2,EP1,EP2)
  1153.       COMPLEX ET1,ET2,EP1,EP2,GAM,ETA
  1154.       COMPLEX GD,CGD,SGD,EGD
  1155.       COMPLEX EGFA,EGFB,EGGD,ESA,ESB
  1156.       COMPLEX CST
  1157.       FP=12.56637
  1158.       XAB=XB-XA
  1159.       YAB=YB-YA
  1160.       ZAB=ZB-ZA
  1161.       CA=XAB/D
  1162.       CB=YAB/D
  1163.       CG=ZAB/D
  1164.       G=(CA*CPH+CB*SPH)*STH+CG*CTH
  1165.       GK=1.-G*G
  1166.       ET1=(.0,.0)
  1167.       ET2=(.0,.0)
  1168.       EP1=(.0,.0)
  1169.       EP2=(.0,.0)
  1170.       IF (GK.LT..001) GO TO 200
  1171.       FA=(XA*CPH+YA*SPH)*STH+ZA*CTH
  1172.       FB=(XB*CPH+YB*SPH)*STH+ZB*CTH
  1173.       EGFA=CEXP(GAM*FA)
  1174.       EGFB=CEXP(GAM*FB)
  1175.       EGGD=CEXP(GAM*G*D)
  1176.       CST=ETA/(GK*SGD*FP)
  1177.       ESA=CST*EGFA*(EGGD-G*SGD-CGD)
  1178.       ESB=CST*EGFB*(1./EGGD+G*SGD-CGD)
  1179.       T=(CA*CPH+CB*SPH)*CTH-CG*STH
  1180.       P=-CA*SPH+CB*CPH
  1181.       ET1=T*ESA
  1182.       ET2=T*ESB
  1183.       EP1=P*ESA
  1184.       EP2=P*ESB
  1185.   200 CONTINUE
  1186.       RETURN
  1187.       END
  1188.