

      COMPLEX A(401),B(401),C(401),E(401),SIGMA(401) 
      COMPLEX ZFILL
      REAL VEL,FSTART,FSTOP,FREQ(401),X(401),T,PI,TH
      REAL TEST,AMP(401),PHASE(401)
      CHARACTER*128 FILVAR
      INTEGER NPOINTS
                     
      PARAMETER(PI=3.141593)
      PARAMETER(TH=1.0E8)
      PARAMETER(ZFILL=(0.0,0.0)) 
      filvar='sph_14.x'
c      filvar='sph_12'
      VEL = 2.998E+8 
C     FREQUENCIES IN GHz, SPHERE DIAMETER IN INCHES      
       FSTART = 9.3
       FSTOP =  9.7 
c      FSTART = 1.1
c      FSTOP =  1.4 
      DIAMETER = 14.
      RADIUS = DIAMETER/2.0
      RADIUS = RADIUS*25.4/1000.0
      NPOINTS = 21
      FDELTA = (FSTOP-FSTART)/(NPOINTS-1)
      DO 11, I=1,NPOINTS
        FREQ(I)=FSTART+(I-1)*FDELTA
        X(I)=2*PI*FREQ(I)*1E9*RADIUS/VEL
 11   CONTINUE
      DO 14, I=1,NPOINTS
      B(I)=CMPLX(1.0,0.0)
      C(I)=CMPLX(1.0,0.0)
      E(I)=CMPLX(0.0,0.0)
      T=0.0
      DO 12, J=1,1000
        T=T+1
        A(I)=B(I)
        B(I)=C(I)
        C(I)=A(I)+B(I)*CMPLX(0.0,(1.0-2.0*T)/X(I))
        E(I)=E(I)+CMPLX(0.0,T+0.5)/(C(I)*(B(I)+C(I)*CMPLX(0.0,-T/X(I))))
        TEST=ABS(REAL(C(I)))
        IF(TEST.GT.TH) GOTO 13
 12   CONTINUE
 13   CONTINUE
      AMP(I)=(CABS(E(I))**2)*4.0*PI*(RADIUS**2)/X(I)**2
      AMP(I)=SQRT(AMP(I))
      PHASE(I)=ATAN2(AIMAG(E(I)),REAL(E(I)))
      SIGMA(I)=AMP(I)*CMPLX(COS(PHASE(I)),SIN(PHASE(I)))/sqrt(4.*pi)
      PHASE(I)=180*PHASE(I)/PI
      WRITE(*,*) FREQ(I),'   ',20*LOG10(AMP(I)),'     ',PHASE(I)
C     WRITE(*,*) FREQ(I),'   ',SIGMA(I)
 14   CONTINUE
      OPEN(UNIT=1,FILE=FILVAR)
      DO 16, J=1,NPOINTS
        WRITE(1,*) FREQ(J),SIGMA(J)
 16   CONTINUE
      CLOSE(UNIT=1)            
      STOP 
 310  FORMAT(I3,4X,F11.7)
 311  FORMAT(I3,4X,E12.5,4X,E12.5)
      END

