      COMPLEX REFREL,S1(400),S2(400)
      REAL PIE
      real RAD
      real freq
      real QSCA,QEXT,QBACK
      
      REFREL=CMPLX(8.15107, 1.76359)
      RAD=1.e+0 !1.e-6
      freq=10.000000

      PIE = 4.0*ATAN(1.0)
      wavel = 0.3/freq
      X = 2.0*PIE*RAD/WAVEL

      CALL BHMIE(X,REFREL,181,S1,S2,QEXT,QSCA,QBACK)
      print *,'qsca=',QSCA,', qext=',QEXT,', qback=',QBACK,
     &        ',pirsq=',PIE*RAD*RAD

      STOP
      END

      SUBROUTINE OLD_MAIN()
C*******************************C
C                               C
C            MIE.FTN            C
C   (Bohren and Huffman,1983)   C
C*******************************C
C                               C
C  MIE.FTN calculates the size  C
C  parameter (kD) and relative  C
C  refractive index (n) for a   C
C  given sphere refractive      C
C  index, medium refractive     C
C  index, radius and free space C
C  wavelength.  The amplitude   C
C  scattering matrix elements   C
C  and efficiencies are also    C
C  computed.                    C
C*******************************C
C
      COMPLEX REFREL,S1(400),S2(400)
      REAL PIE, phase
      character*70 name
      WRITE(6,11)
C
      PIE = 4.0*ATAN(1.0)
C
C**** (REAL) REFRACTIVE INDEX OF SURROUNDING MEDIUM = REFMED
C
      PRINT *, 'ENTER refractive index of medium (real)'
      READ *, REFMED
C
C**** REFRACTIVE INDEX OF SPHERE = REFRE + iREFIM
C
      PRINT *, 'ENTER refractive index of SPHERE (real)'
      READ *, REFRE
      PRINT *, 'ENTER refractive index of SPHERE (Imag)'
      READ *, REFIM
      REFREL = CMPLX(REFRE,REFIM)/REFMED
      WRITE(6,12) REFMED,REFRE,REFIM
C
C**** INPUT SPHERE RADIUS (RAD) AND WAVELENGTH (WAVEL)
C
      PRINT *, 'Enter radius of  SPHERE (m), and frequency in GHz'
      READ *, RAD, freq

      wavel = 0.3/freq
      X = 2.0*PIE*RAD*REFMED/WAVEL
      xx = 2.0*pie/ wavel
      WRITE(6,13) RAD,WAVEL
      WRITE(6,14) X
C
C**** COMPUTE MATRIX ELEMENTS AT OBSERVATION ANGLES
C
c      PRINT *, 'NUMBER OF OBSERVATION ANGLES'
c      PRINT *, 'BETWEEN THETA = 0 AND 90 DEGREES'
c      READ *, NANG

      nang = 181
      DANG = (PIE/2.0)/FLOAT(NANG-1)
      CALL BHMIE(X,REFREL,NANG,S1,S2,QEXT,QSCA,QBACK)
      WRITE(6,65) QSCA,QEXT,QBACK
c      WRITE(6,15)
c      WRITE(6,16)
c      WRITE(6,17)
      print*,'Enter the name of the O/P file'
      read(*,500) name
  500 format(a70)
      open(11,file=name)
C
C**** S33 AND S34 MATRIX ELEMENTS ARE NORMALIZED BY S11
C**** S11 IS NORMALIZED TO 1.0 IN THE FORWARD DIRECTION
C**** POL = DEGREE OF POLARIZATION (UNPOLARIZED INCIDENT
C**** FIELD
C
      S11NOR = 0.5*(CABS(S2(1))**2 + CABS(S1(1))**2)
      NAN = 2*NANG-1
      DO 355 J=1,NAN
        s2(j) = -conjg(s2(j) * cmplx(0.,1.0))
        s1(j) =  conjg(s1(j) * cmplx(0.,1.0))
        write(11,*) (j-1)*dang*180.0/pie,
     $        10.*alog10(4.0*pie*(cabs(s1(j)/xx))**2),
     $        10.*alog10(4.0*pie*(cabs(s2(j)/xx))**2)
     $       , phase(s2(j)/s1(j))   
 355  CONTINUE
C
C **** PRINT FORMATS
C
  65  FORMAT(//,1X,'Qsca = ',E13.6,3X,'Qext =',E13.6,3X,
     &      'Qback = ',E13.6)
  75  FORMAT(1X,F6.2,2X,E13.6,2X,E13.6,2X,E13.6,2X,E13.6)
  11  FORMAT(/,'MIE SCATTERING PROGRAM -- BOHREN AND HUFFMAN'//)
  12  FORMAT(2X,'n(Medium) = ',F6.4,3X,'Real n(Sphere) = ',E10.4,3X,
     &      'Imag n(Sphere)= ',E10.4)
  13  FORMAT(2X,'SPHERE RADIUS = ',F7.3,3X,'WAVELENGTH = ',F7.4)
  14  FORMAT(2X,'SIZE PARAMETER (kD) =',F8.3)
  15  FORMAT(/,21X,'Stokes Matrix Elements')
  16  FORMAT(/,20X,'I  S11  S12   0   0    I',
     &       /,20X,'I  S12  S11   0   0    I',
     &       /,20X,'I   0    0   S33  S34  I',
     &       /,20X,'I   0    0  -S34  S33  I')
  17  FORMAT(//,2X,'ANGLE',7X,'S11',8X,'POL(-S12/S11)'
     &        ,6X,'S33',12X,'S34'//)
      STOP
      END
C
C**** SUBROUTINE BHMIE CALCULATES THE AMPLITUDE SCATTERING
C**** MATRIX ELEMENTS AND EFFICIENCIES FOR EXTINCTION, THE
C**** TOTAL SCATTERING AND BACKSCATTERING FOR A GIVEN SIZE
C**** PARAMETER (kD) AND RELATIVE INDEX OF REFRACTION
C
      SUBROUTINE BHMIE(X,REFREL,NANG,S1,S2,QEXT,QSCA,QBACK)
C
      DIMENSION AMU(400),THETA(400),PI(400),TAU(400),PI0(400)
      DIMENSION PI1(400)
      COMPLEX D(12000),Y,REFREL,XI,XI0,XI1,AN,BN,S1(400),S2(400)
      REAL*8 PSI0,PSI1,PSI,DN,DX,PIE,ONE
C
      DX = X
      Y = X*REFREL
      ONE = 1.0000000000
      PIE = 4.0*DATAN(ONE)
C
C**** SERIES SUMMATION TERMINATED AFTER NSTOP TERMS
C
      XSTOP = X + 4.0*X**0.3333 + 2.0
      NSTOP = XSTOP
      YMOD = CABS(Y)
      NMX = AMAX1(XSTOP,YMOD) + 15
      DANG = (PIE/2.0)/FLOAT(NANG - 1)
C
           DO 555 J=1,NANG
              THETA(J) = (FLOAT(J) - 1)*DANG
              AMU(J) = COS(THETA(J))
  555      CONTINUE
C
C**** LOGARITHMIC DERIVATIVE D(J) CALCULATED BY DOWNWARD
C**** RECURRENCE BEGINNING WITH THE INITIAL VALUE 0 + i0
C**** AT J = NMX
C
      D(NMX) = CMPLX(0.0,0.0)
      NN = NMX - 1
      DO 120 N=1,NN
         RN = NMX - N + 1
         D(NMX - N) = (RN/Y) - (1./(D(NMX-N+1)+RN/Y))
  120 CONTINUE
C
      DO 140 J=1,NANG
         PI0(J) = 0.0
         PI1(J) = 1.0
  140 CONTINUE
C
      NN = 2*NANG -1
C
      DO 160 J=1,NN
         S1(J) = CMPLX(0.0,0.0)
         S2(J) = CMPLX(0.0,0.0)
  160 CONTINUE
C
C**** RICCATI-BESSEL FUNCTION WITH REAL ARGUMENT X
C**** AND CALCULATED BY UPWARD RECURRENCE
C
      PSI0 = DCOS(DX)
      PSI1 = DSIN(DX)
      CHI0 = -SIN(X)
      CHI1 = COS(X)
      APSI0 = PSI0
      APSI1 = PSI1
      XI0 = CMPLX(APSI0,-CHI0)
      XI1 = CMPLX(APSI1,-CHI1)
      QSCA = 0.0
      N = 1
 200  DN = N
      RN = N
      FN =(2.0*RN+1.0)/(RN*(RN+1.0))
      PSI = (2.0*DN-1.)*PSI1/DX-PSI0
      APSI = PSI
      CHI = (2.0*RN-1.0)*CHI1/X - CHI0
      XI = CMPLX(APSI,-CHI)
      AN = (D(N)/REFREL+RN/X)*APSI - APSI1
      AN = AN/((D(N)/REFREL+RN/X)*XI - XI1)
      BN = (REFREL*D(N)+RN/X)*APSI - APSI1
      BN = BN/((REFREL*D(N)+RN/X)*XI - XI1)
      QSCA = QSCA+(2.0*RN+1.)*(CABS(AN)*CABS(AN)+CABS(BN)*CABS(BN))
C
      DO 210 J=1,NANG
         JJ= 2*NANG-J
         PI(J) = PI1(J)
         TAU(J) = RN*AMU(J)*PI(J) - (RN+1.)*PI0(J)
         P = (-1.)**(N-1)
         S1(J) = S1(J)+FN*(AN*PI(J)+BN*TAU(J))
         T = (-1.)**N
         S2(J) = S2(J)+FN*(AN*TAU(J)+BN*PI(J))
         IF(J .EQ. JJ) GO TO 210
         S1(JJ) = S1(JJ) + FN*(AN*PI(J)*P + BN*TAU(J)*T)
         S2(JJ) = S2(JJ) + FN*(AN*TAU(J)*T + BN*PI(J)*P)
  210 CONTINUE
C
      PSI0 = PSI1
      PSI1 = PSI
      APSI1 = PSI1
      CHI0 = CHI1
      CHI1 = CHI
      XI1 = CMPLX(APSI1,-CHI1)
      N = N+1
      RN = N
C
      DO 999 J=1,NANG
         PI1(J) = ((2.*RN-1.)/(RN-1.))*AMU(J)*PI(J)
         PI1(J) = PI1(J)-RN*PI0(J)/(RN-1.)
         PI0(J) = PI(J)
  999 CONTINUE
C
      IF (N-1-NSTOP) 200,300,300
  300 QSCA = (2./(X*X))*QSCA
      QEXT = (4./(X*X))*REAL(S1(1))
      QBACK = (4./(X*X))*CABS(S1(2*NANG-1))*CABS(S1(2*NANG-1))
      RETURN
      END
c---------------------------------------------------------------------
c---------------------------------------------------------------------
c
        FUNCTION PHASE(Z)
c
c---------------------------------------------------------------------
        COMPLEX Z 
        PI=4.*ATAN(1.)
        X=REAL(Z)
        Y=AIMAG(Z)
        PHASE=(180./PI)*ATAN2(Y,X)
c        if(phase.lt.-150.0) phase=phase+360
        RETURN
        END
