	SUBROUTINE ZANTEN(F)
C
C	CALCULATES SELF AND MUTUAL IMPEDANCES OF ANTENNAS
C
	COMPLEX YF,ZA,TBAR,TEMP,ZT,CUR(50),Z11,Z12,Z21,Z22
	COMMON /MATR/ YF(50,50),ZA(50,50),TBAR(50,50),TEMP(50,50),CUR
	COMMON /DESIGN/ NELEM,X(50),S(50),H(50),VEL,XL,CLT,XSPAC,XDIA,ZT
	COMMON /DIPOLES/ H1,H2,SESP4,FREQMC,Z11,Z12,Z21,Z22
	DATA TWOPI /6.2831853/
C
	W = TWOPI*F
C	PRINT*,'FREQ',F
	FREQMC = 1.E-6*F
C	DO N = 1,NELEM
C	  DO M = 1,NELEM
C	  ZA(N,M) = 0.
C	  ENDDO
C	ENDDO
	DO I = 1,NELEM-1
	  DO J = I+1,NELEM 
	    SESP4 = X(J) - X(I)
	    H1 = 2.*H(I)
	    H2 = 2.*H(J)
	    CALL ESP4SR
	    ZA(I,I) = Z11
	    ZA(J,J) = Z22
C		ALTERNATE CONNECTION TO FEEDER
	    ISGN = (J-I).AND.1
	    IF(ISGN.GT.0) THEN
	      ZA(I,J) = -Z12
	      ZA(J,I) = -Z21
	    ELSE
	      ZA(I,J) = Z12
	      ZA(J,I) = Z21
	    ENDIF
C	    PRINT*,' ',ISGN,ZA(I,J)
C	    PRINT*,I,J,Z11,Z12
C	    PRINT*,I,J,Z21,Z22
	  ENDDO
	ENDDO
C	FAKE VERY LARGE IMPEDANCE
C	DO N = 1,NELEM
C	  ZA(N,N) = 1.E3
C	ENDDO
	RETURN
	END
	FUNCTION AGAIN(F,THETA,PHI)
C
C	CALCULATES GAIN OF SYSTEM IN DIRECTION THETA,PHI AT FREQ. F
C	THETA AND PHI ARE IN DEGREES
C	THETA IS MEASURED FROM THE WIRE AXIS DIRECTION, PHI FROM DIRECTION
C	OF ANTENNA RADIATION AXIS, A WIERD CONVENTION, BUT APPARENTLY USUAL IN
C	ANTENNA CALCULATIONS
C
	COMPLEX YF,ZA,TBAR,TEMP,ZT,AGAIN,E,B,CUR(50)
	COMMON /MATR/ YF(50,50),ZA(50,50),TBAR(50,50),TEMP(50,50),CUR
	COMMON /DESIGN/ NELEM,X(50),S(50),H(50),VEL,XL,CLT,XSPAC,XDIA,ZT
	DATA TWOPI /6.2831853/
C
	W = TWOPI*F
	XK = W/CLT
C
C	DISTANCES ARE MEASURED FROM THE CENTER OF THE SHORTEST ELEMENT
C	THE POSITIVE DIRECTION OF CURRENT IS OUT OF THE POSITIVE SIDE OF
C	THE FEEDER, HENCE MUTUAL IMPEDANCE IS REVERSED.
C
C	CALCULATE ANGLE ANTENNA RADIATION AXIS TO VIEW DIRECTION
C	RADIATION AXIS IS  X DIRECTION
C	X(N) WAS POSITIVE, BUT IS COUNTED NEGATIVE IN THIS SYSTEM
C	INCREASING /X/ MEANS PHASE DELAY, NEGATIVE PHASE SHIFT
C
	SINTH = SIND(THETA)
	COSTH = COSD(THETA)
	COSAXD = SINTH*COSD(PHI)
	SINAXD = SQRT(1. - COSAXD**2)
C
	E = CMPLX(0.,0.)
	B = CMPLX(0.,0.)
	SG = 1.
	DO N = 1,NELEM
	  HKN = XK*H(N)
	  DPHASE = -X(N)*XK*COSAXD
	  IF(SINTH.NE.0..AND.SIN(HKN).NE.0.)
     1	  E = E + SG*CUR(N)*((COS(HKN*COSTH) - COS(HKN))/SIN(HKN))*
     1	  CEXP(CMPLX(0.,DPHASE))/SINTH
	  SG = -SG
	ENDDO
	AGAIN = E*CONJG(E)
	RETURN
	END

