!+ Subroutine to calculate single scattering after Mitchell and Macke.
!
! Method:
!	The cross-sections of extinction and scattering and the 
!	asymmetry of  a varietu of shapes are calculated using 
!	the enhancements to ADT intriduced by Mitchell with the
!       asymmetry taken from Macke's ray-tracing programs.
!
! Special Note:
!	This is not for general release and so is removed on
!	configuration if private is not set.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		26-06-96		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE ENHANCED_ADT(WAVELENGTH, RADIUS, ASPECT
     &   , N_RELATIVE
     &   , XS_SCATTER, XS_EXT, ASYMMETRY
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/math_cnst.ccf'
!
!     DUMMY ARGUMENTS
      REAL*8	!, INTENT(IN)
     &     WAVELENGTH
!             WAVELENGTH
     &   , RADIUS
!             RADIUS OF CIRCUMSCRIPTION
     &   , ASPECT
!             ASPECT RATIO OF CRYSTAL
      COMPLEX*16	!, INTENT(IN)
     &     N_RELATIVE
!             RELATIVE REFRACTIVE INDEX
      REAL*8	!, INTENT(OUT)
     &     XS_SCATTER
!             SCATTERING CROSS_SECTION
     &   , XS_EXT
!             EXTINCTIVE CROSS_SECTION
     &   , ASYMMETRY
!             ASYMMTRY OF SCATTERING
!
!     LOCAL VARIABLES
      REAL*8
     &     AREA
!             SURFACE AREA OF CRYSTAL
     &   , RADIUS_EQV
!             EQUIVALENT RADIUS OF CRYSTAL
     &   , MI_SIZE
!             IM(REF. INDEX)* SIZE PARAMETER
     &   , OMEGA
!             SINGLE SCATTERING ALBEDO
     &   , INV_REFRACT
!             DIFFERENCE OF INVERSE
!                REFRACTIVE INDICES
!
      REAL*8
     &     DA
!             EFFECTIVE DIMENSION FOR THE ASYMMETRY
     &   , ASYMMETRY_LOW
!             ASYMMETRY FACTOR AT LOWER END OF WAVELENGTH INTERVAL
     &   , ASYMMETRY_HIGH
!             ASYMMETRY FACTOR AT HIGHER END OF WAVELENGTH INTERVAL
!
      REAL*8
     &     SQRT_3_BY_2
!             CONSTANT
!
!
!
      SQRT_3_BY_2=SQRT(3.0D+00)/2.0D+00
!
      IF (I_SHAPE.NE.IP_SHAPE_HEXCYL) THEN
         WRITE(IU_ERR, '(/A)') 
     &      '*** ERROR: THIS ROUTINE HAS NOT YET BEEN EXTENDED '
     &      //'TO DEAL WITH NON-HEXAGONAL COLUMNS.'
         IERR=I_ERR_FATAL
         RETURN
      ENDIF
!
      IF (I_SHAPE.EQ.IP_SHAPE_HEXCYL) THEN
!
!        SET DA.
         AREA=6.0D+00*RADIUS*RADIUS*(SQRT_3_BY_2+2.0D+00*ASPECT)
         RADIUS_EQV=SQRT(AREA/(4.0D+00*PI))
         MI_SIZE=2.0D+00*PI*RADIUS*DIMAG(N_RELATIVE)/WAVELENGTH
         DA=RADIUS_EQV
!
!        CALCULATE ASYMMETRIES USING LINEAR INTERPOLATION BETWEEN 
!        THE FITS TO MACKE'S RAY-TRACING RESULTS.
         IF (WAVELENGTH.LT.5.5D-07) THEN
            ASYMMETRY=9.712D-01*DA**4.473D-02
         ELSE IF ( (WAVELENGTH.GE.5.5D-07)
     &      .AND.(WAVELENGTH.LT.1.0D-06) ) THEN
            ASYMMETRY_LOW=9.712D-01*DA**4.473D-02
            ASYMMETRY_HIGH=9.780D-01*DA**4.473D-02
            ASYMMETRY=ASYMMETRY_LOW+(ASYMMETRY_HIGH-ASYMMETRY_LOW)
     &         *(WAVELENGTH-5.5D-07)/4.5D-07
         ELSE IF ( (WAVELENGTH.GE.1.0D-06)
     &      .AND.(WAVELENGTH.LT.1.6D-06) ) THEN
            ASYMMETRY_LOW=9.780D-01*DA**4.473D-02
            ASYMMETRY_HIGH=1.007D+00*DA**4.473D-02
            ASYMMETRY=ASYMMETRY_LOW+(ASYMMETRY_HIGH-ASYMMETRY_LOW)
     &         *(WAVELENGTH-1.0D-06)/6.0D-07
         ELSE IF ( (WAVELENGTH.GE.1.6D-06)
     &      .AND.(WAVELENGTH.LT.2.2D-06) ) THEN
            ASYMMETRY_LOW=1.007D+00*DA**4.473D-02
            ASYMMETRY_HIGH=1.046D+00*DA**4.473D-02
            ASYMMETRY=ASYMMETRY_LOW+(ASYMMETRY_HIGH-ASYMMETRY_LOW)
     &         *(WAVELENGTH-1.0D-06)/6.0D-07
         ELSE IF ( (WAVELENGTH.GE.2.2D-06)
     &      .AND.(WAVELENGTH.LT.3.0D-06) ) THEN
            ASYMMETRY_LOW=1.046D+00*DA**4.473D-02
            ASYMMETRY_HIGH=9.595D-01+1.076D+00*DA
            ASYMMETRY=ASYMMETRY_LOW+(ASYMMETRY_HIGH-ASYMMETRY_LOW)
     &         *(WAVELENGTH-1.0D-06)/6.0D-07
         ELSE IF (WAVELENGTH.GE.3.0D-06) THEN
            ASYMMETRY=9.595D-01+1.076D+00*DA
         ENDIF
!
      ENDIF
!
!
!
      RETURN
      END
