      PROGRAM CALLBHCOAT
      IMPLICIT NONE

! Local variables:
      INTEGER NCOMPMX
      PARAMETER(NCOMPMX=40)

      INTEGER IDIELC,IDIELM,IHOW,ININDEX,INPIND,JC,JW,NWAVE
      REAL DW,ENIMC,ENIMM,ENRE1C,ENRE1M,EPS11C,EPS11M,EPS2C,EPS2M,
     &   PI,RADCOR,RADMAN,REFMED,TGR,WAVE,WAVE1,WAVE2
      COMPLEX CXI,EPSCOR,EPSMAN
      CHARACTER*60 COMP(1:NCOMPMX)

! Variables passed to subroutine:

      REAL GSCA,QABS,QBACK,QEXT,QSCA,XCOR,XMAN
      COMPLEX REFCOR,REFMAN

      EXTERNAL INDEX
!***********************************************************************
! Program to compute optical properties of coated spheres using
! subroutine BHCOAT from Bohren & Huffman
!
! History:
! 90.10.11 (BTD) Modified structure of CALLBHCOAT
! 92.11.24 (BTD) Explicit declaration of all variables
! 00.05.16 (BTD) Modified to compile under g77
! 12.06.15 (BTD) v2 
!                * Modified to use new version bhcoat_v2.f that
!                  calculates <cos(theta)> for scattered power
! end history
!
!***********************************************************************
!
! Caution:
!     BHCOAT should not be used for large, highly absorbing coated
!     spheres.
!     x*refim1, x*refim2, and y*refim2 should be less than about 30.
!
! REFMED = (real) refractive index of surrounding medium
!***********************************************************************
      PI=4.*ATAN(1.)
      REFMED=1.0
      CXI=(0.,1.)

      COMP(1)='DL84 Graphite, E parallel to c'
      COMP(2)='DL84 Graphite, E perp. to c'
      COMP(3)='DL84 Astronomical silicate'
      COMP(4)='Diamond, Type Ib, N/C=0.003'
      COMP(5)='Meteoritic diamond, Allende CJ'
      COMP(6)='alpha-SiC'
      COMP(7)='Astronomical silicate, smoothed UV'
      COMP(8)='new astronomical silicate'
      COMP(9)='DL84 graphite, E parallel to c, but sigma=62.5 mho/cm'
      COMP(10)='crystalline H2O ice'
      COMP(11)='Draine 2003 astrosilicate'
      COMP(12)='Draine 2003 graphite E parallel to c, sigma=62.5 mho/cm'
      COMP(13)='Draine 2003 graphite E perp. c'
      COMP(14)='Draine 2003 astrosilicate with porosity(a)'
      COMP(15)='Draine 2003 graphite E parallel c, porosity(a)'
      COMP(16)='Draine 2003 graphite E perp. c, porosity(a)'
      COMP(17)='Djuric & Li 1999 graphite, E parallel c, optical data'
      COMP(18)='Djuric & Li 1999 graphite, E parallel c, EELS data'
      COMP(19)='Djuric & Li 1999 graphite, E perp c'
      COMP(20)='800C pyrolyzed cellulose, Jaeger etal 1998'
      COMP(21)='600C pyrolyzed cellulose, Jaeger etal 1998'
      COMP(22)='amorph. carbon AC1, Rouleau & Martin 1991'
      COMP(23)='amorph. carbon BE1, Rouleau & Martin 1991'
      COMP(24)='amorph. carbon FC21PS, Rouleau & Martin 1991'
      COMP(25)='amorph. carbon HAPS, Rouleau & Martin 1991'
      COMP(26)='amorph. carbon BE, Zubko etal 1996'
      COMP(27)='amorph. carbon ACAR, Zubko etal 1996'
      COMP(28)='amorph. carbon ACH2, Zubko etal 1996'
      COMP(29)='amorph. carbon amcACH2, modified from Zubko etal 1996'
      COMP(30)='D03 graphite, E||c, 2-comp model for free electrons'
      COMP(31)='Fe metal'
      COMP(32)='Fe3O4 magnetite'
      COMP(33)='gamma-Fe2O3 maghemite'
      COMP(34)='Hensley & Draine 2014 astrosilicate'
      COMP(35)='Hensley & Draine 2014 astrosil for 10% of Fe in incl.' 
      COMP(36)='Hensley & Draine 2014 astrosil for 25% of Fe in incl.' 
      COMP(37)='Hensley & Draine 2014 astrosil for 50% of Fe in incl.' 
      COMP(38)='Hensley & Draine 2014 astrosil for 100% of Fe in incl.' 
      COMP(39)='amorphous H2O, T=25K, from Mastrapa et al 2009'

      OPEN(UNIT=7,FILE='callbhcoat.out')

! Specify whether to supply refractive index by hand or to obtain it
! using subroutine INDEX

      WRITE(0,FMT='(A,A)')'Wish to input refractive indices by ',
     &   'hand (0) or obtain them from subroutine INDEX (1)?'
      READ(*,*)ININDEX
      IF(ININDEX.GE.1)THEN
         WRITE(0,FMT='(A)')'known materials:'
         DO JC=1,39
            WRITE(0,FMT='(I2,1X,A)')JC,COMP(JC)
         ENDDO
         WRITE(0,FMT='(A)')'Enter IDIEL for core'
         READ(*,*)IDIELC
         WRITE(0,FMT='(A,A)')' core  = ',COMP(IDIELC)
         WRITE(7,FMT='(A,A)')' core  = ',COMP(IDIELC)
         WRITE(0,FMT='(A)')'Enter IDIEL for mantle'
         READ(*,*)IDIELM
         WRITE(0,FMT='(A,A)')'mantle = ',COMP(IDIELM)
         WRITE(7,FMT='(A,A)')'mantle = ',COMP(IDIELM)
         WRITE(0,FMT='(A)')'Enter TGR = temperature (K)'
         READ(*,*)TGR
      ELSE
         WRITE(0,FMT='(A,A)')'Wish to specify ref. indices (0)',
     &     ' or diel.consts (1)?'
         READ(*,*)INPIND
      ENDIF

 0050 WRITE(*,*)'Enter radius of core (physical units (0 to stop)'
      READ(*,*)RADCOR
      WRITE(0,FMT='(1PE10.3,A)')RADCOR,' = radius of core (um)'
      WRITE(7,FMT='(1PE10.3,A)')RADCOR,' = radius of core (um)'
      IF(RADCOR.LE.0.)STOP
      WRITE(*,*)'Enter radius of mantle (physical units)'
      READ(*,*)RADMAN
      WRITE(0,FMT='(1PE10.3,A)')RADMAN,' = outer radius of mantle (um)'
      WRITE(7,FMT='(1PE10.3,A)')RADMAN,' = outer radius of mantle (um)'

      IF(ININDEX.LE.0)THEN
 0100    WRITE(*,FMT='(A,A)')'Enter wavelength (physical units)',
     &      ' (0 to change size)'
         READ(*,*)WAVE
         IF(WAVE.LE.0.)GOTO 0050
         IF(INPIND.EQ.0)THEN
            WRITE(*,*)'Enter complex refractive index of core'
            READ(*,*)REFCOR
            WRITE(*,*)'Enter complex refractive index of mantle'
            READ(*,*)REFMAN
         ENDIF
         IF(INPIND.EQ.1)THEN
            WRITE(*,*)'Enter complex dielec. constant of core'
            READ(*,*)EPSCOR
            WRITE(*,*)'Enter complex dielec. constant of mantle'
            READ(*,*)EPSMAN
            REFCOR=SQRT(EPSCOR)
            REFMAN=SQRT(EPSMAN)
         ENDIF
         XCOR=2.*PI*RADCOR*REFMED/WAVE
         XMAN=2.*PI*RADMAN*REFMED/WAVE
         CALL BHCOAT(XCOR,XMAN,REFCOR,REFMAN,QEXT,QSCA,QBACK,GSCA)
         QABS=QEXT-QSCA
         WRITE(0,6010)RADCOR,RADMAN,WAVE,XMAN,REFCOR,REFMAN,QEXT,QSCA,
     &                QABS,QBACK,GSCA
         GOTO 0100
      ELSE ! use subr. index, loop over wavelength
         WRITE(0,FMT='(A)')'enter first, last wavelength, how many'
         READ(*,*)WAVE1,WAVE2,NWAVE
         WRITE(0,FMT='(A,A)')'how? 1,2,3 for linear in ',
     &      'wave, 1/wave, log(wave)'
         READ(*,*)IHOW
         WRITE(0,FMT='(1X,A,3X,A,6X,A,5X,A,1X,A)')'wave(um)','Q_ext',
     &      'Q_abs','Q_sca','<cos(theta)>'
         WRITE(7,FMT='(1X,A,3X,A,6X,A,5X,A,1X,A)')'wave(um)','Q_ext',
     &      'Q_abs','Q_sca','<cos(theta)>'
         DO JW=1,NWAVE
            IF(IHOW.EQ.1)THEN
               DW=(WAVE2-WAVE1)/MAX(REAL(NWAVE-1),1.E-6)
               WAVE=WAVE1+DW*(JW-1)
            ELSEIF(IHOW.EQ.2)THEN
               DW=(1./WAVE2-1./WAVE1)/MAX(REAL(NWAVE-1),1.E-6)
               WAVE=1./(1./WAVE2+DW*(JW-1))
            ELSEIF(IHOW.EQ.3)THEN
               DW=LOG(WAVE2/WAVE1)/MAX(REAL(NWAVE-1),1.E-6)
               WAVE=WAVE1*EXP(DW*(JW-1))
            ENDIF
!*** diagnostic
!            write(0,*)'callbhcoat_v2 ckpt 1: wave=',wave
!***
            CALL INDEX(RADCOR,TGR,WAVE,IDIELC,EPS11C,EPS2C,ENRE1C,ENIMC)
            CALL INDEX(RADMAN,TGR,WAVE,IDIELM,EPS11M,EPS2M,ENRE1M,ENIMM)
            REFCOR=1.+ENRE1C+CXI*ENIMC
            REFMAN=1.+ENRE1M+CXI*ENIMM
            XCOR=2.*PI*RADCOR*REFMED/WAVE
            XMAN=2.*PI*RADMAN*REFMED/WAVE
            CALL BHCOAT(XCOR,XMAN,REFCOR,REFMAN,QEXT,QSCA,QBACK,GSCA)
            QABS=QEXT-QSCA
            WRITE(0,FMT='(1PE9.3,1PE10.3,1PE11.3,1PE10.3,1PE11.3)')
     &         WAVE,QEXT,QABS,QSCA,GSCA
            WRITE(7,FMT='(1PE9.3,1PE10.3,1PE11.3,1PE10.3,1PE11.3)')
     &         WAVE,QEXT,QABS,QSCA,GSCA
         ENDDO ! enddo jw=1,nwave
      ENDIF
      STOP

 6010 FORMAT(
     &  1PE11.4,' = radius of core',/,
     &  1PE11.4,' = outer radius of mantle',/,
     &  1PE11.4,' = lambda',/,
     &  1PE11.4,' = x = 2*pi*rmantle/lambda',/,
     &  1P2E11.4,' = refractive index of core',/,
     &  1P2E11.4,' = refractive index of mantle',/,/,
     &  1PE11.4,' = Q_ext',/,
     &  1PE11.4,' = Q_sca',/,
     &  1PE11.4,' = Q_abs',/,
     &  1PE11.4,' = Q_bk=4*pi*dQ_sca/dOmega',/,
     &  1PE11.4,' = g = <cos(theta)>')
      END
