Program MieScat

  use vtype_module
  use mie_theory, only : CALLBH3
  use refindex, only : RefIndexInit, RefIndexFinalize, RefIndexCalc

  Implicit NONE

  ! -----CONSTANT-----
  ! mathematical constant
  real(8) :: PI

  ! -----INPUT-----
  ! (real) refractive index of surrounding medium
  real(8) :: Rmed
  ! refractive index of sphere = Rre + i*Rim
  real(8) :: Rre, Rim
  ! radius and wavelength  [m]
  real(8) :: rad, wl

  ! -----OUTPUT-----
  ! extinction efficiency
  real(8) :: QEXT
  ! single scattering albedo and co-albedo
  real(8) :: SSA
  real(8) :: SSCoA
  ! asymmetry factor
  real(8) :: AF


  ! size parameter
  real(8) :: X
  ! extinction cross section  [m2]
  real(8) :: SEXT

  real(8) :: wls
  real(8) :: wle
  real(8) :: dwl


  character(stdstr) :: NcFN = "../mkncfile_refindex/ncfiles/RefIndex_H2OLiq_Segelstein1981.nc"


  call RefIndexInit( &
    & NcFN &
    & )


  ! -----BEGIN-----
  ! set mathematical constant
  PI = acos(-1.0D0)

  ! set input parameters
  ! (real) refractive index of surrounding medium
  Rmed = 1.00045D0

  rad = 10.0d-6

  ! output parameters
  Write(*,*) '#-----input-----'
  Write(*,*) '#refractive index of surrounding medium (real)', &
    &           Rmed
  Write(*,*) '#particle radius [m]                          ', &
    &           rad
  Write(*,*)


  wls =   0.1D-6
  wle = 200.0D-6
  dwl =   1.0D-9

  wl = wls
  do while( wl <= wle )

!!$    Write(*,*) 'wavelength [m]                               ', &
!!$      &           wl

    ! Calculation of refractive index
    ! refractive index of sphere = Rre + i*Rim
    !
    call RefIndexCalc( &
      & Wl,            & ! (in)
      & Rre, Rim       & ! (out)
      & )

    ! calculation
    Call CALLBH3(                &
      & Rmed, Rre, Rim, rad, wl, &
      & QEXT, SSA, SSCoA, AF     &
      & )

    ! output results
    X = 2.0D0*PI*rad*Rmed/wl
    SEXT = PI*(rad*rad)*QEXT

!!$    Write(*,*) '-----output-----'
!!$    Write(*,*) 'size parameter               ', X
!!$    Write(*,*)
!!$    Write(*,*) 'extinction efficiency        ', QEXT
!!$    Write(*,*) 'extinction cross section [m2]', SEXT
!!$    Write(*,*)
!!$    Write(*,*) 'single scattering albedo     ', SSA
!!$    Write(*,*) 'co-single scattering albedo  ', SSCoA
!!$    Write(*,*) 'asymmetry factor             ', AF

    write( 6, * ) wl, X, QExt, SSA, SSCoA, AF


    wl = wl + dwl
  end do

  call RefIndexFinalize

  Stop
End Program MieScat
