!https://www.science-emergence.com/Codes/Lorentz-Mie-scattering-calculation-fortran-code/
program mie_code

  ! reference: Hong Du: Mie-Scattering Calculation, Applied Optics,
  !                     Vol. 43, Issue 9, pp. 1951-1956 (2004)

  implicit none

  integer :: i,j,k

  integer, parameter :: n = 100  ! number of coefficients

!!$  double precision :: x = 5.0  ! Mie Parameter (2*pi*r/wl)
!!$
!!$  double complex :: ri = (1.33d0,-0.0d0)  ! refractive index

  double precision :: x = 5.6  ! Mie Parameter (2*pi*r/wl)
  double complex :: ri = (1.43d0,-1.85d-10)  ! refractive index

  double complex, dimension(n) :: ma , mb

  double precision, dimension(:,:), allocatable :: pi , tau

  double precision, dimension(0:n) :: phi , chi
  double complex, dimension(0:n) :: zeta

  double complex, dimension(0:n) :: rn

  double complex, dimension(:), allocatable :: sa1 , sa2

  double precision, dimension(:), allocatable :: p11,p12,p33,p34

  double precision :: qe

  double precision :: albedo , gfact

  integer :: nb_angles

  double precision :: teta
  double precision :: teta_resol = 1.0 ! scattering angle resolution

  double precision :: s , t

  double precision :: temp1 , temp2 , temp3
  double complex :: ctemp1 , ctemp2 ,ctemp3

  double precision :: coef1 , coef2 , coef3 , coef4

  double precision :: norm

  double precision :: cdr = 0.01745329252

  double complex :: ci = (0.0d0,1.0d0)

  !----------------------------------------------------------------------------------------!

  if (n < int(x)) write(6,*) 'Warning: number of coeffs n is not enough !'

  nb_angles = int(180/teta_resol) + 1

  allocate( pi(nb_angles,n) )
  allocate( tau(nb_angles,n) )

  allocate( sa1(nb_angles) )
  allocate( sa2(nb_angles) )

  allocate( p11(nb_angles) )
  allocate( p12(nb_angles) )
  allocate( p33(nb_angles) )
  allocate( p34(nb_angles) )

  do i = 1 , n
    ma(i)  = (0.0d0,0.0d0)
    mb(i)  = (0.0d0,0.0d0)
  end do

  do i = 1 , nb_angles
    do j = 1 , n
      pi(i,j)  = 0.0d0
      tau(i,j) = 0.0d0
    end do
  end do

  do i = 1 ,nb_angles
    sa1(i) = 0.0d0
    sa2(i) = 0.0d0
  end do

  !----------------------------------------------------------------------------------------!
  ! Step 1: Calculations of the Mie angular functions.

  teta = 0.0d0

  do i = 1 , nb_angles

    !---- pi function

    temp1 = dcos( teta * cdr )

    temp2 = 0.0d0
    temp3 = 1.0d0
    pi(i,1) = temp3
    do j = 2 , n
      s = temp1 * temp3
      t = s - temp2
      pi(i,j) = s + t + t / ( j - 1 )
      temp2 = temp3
      temp3 = pi(i,j)
    end do

    !----- tau function

    tau(i,1) = temp1
    do j = 2 , n
      t = temp1 * pi(i,j) - pi(i,j-1)
      tau(i,j) = j * t - pi(i,j-1)
    end do

    teta = teta + teta_resol
  end do

  !----------------------------------------------------------------------------------------!
  ! Step 2: Calculations of the Riccati-Bessel functions by upward recurrence.

  !----- phi function

  temp1 = dcos( x )
  phi(0) = dsin( x )
  temp2 = phi(0)
  do i = 1 , n
    phi(i) = ( ( 2 * i - 1 ) / x ) * temp2  - temp1
    temp1 = temp2
    temp2 = phi(i)
  end do

  !----- chi function

  chi(0) = dcos(x)
  temp1  = chi(0)
  chi(1) = dcos(x) / x + dsin(x)
  temp2  = chi(1)
  do i = 2 , n
    chi(i) = ( ( 2 * i - 1 ) / x ) * temp2  - temp1
    temp1 = temp2
    temp2 = chi(i)
  end do

  !----- zeta fcuntion

  do i = 0 , n
    zeta(i) = phi(i) + ci * chi(i)
  end do

  !----------------------------------------------------------------------------------------!
  ! Step 3: Calculations of the complex functions rn(mx) by upward recurrence.

  temp1 = dble( ri ) * x
  temp2 = aimag( ri ) * x

  coef1 = dtan( temp1 ) - dexp( - 2.0 * temp2 ) * dtan( temp1 )
  coef2 = 1.0d0 + dexp( - 2.0 * temp2 )
  coef3 = - 1.0d0 + dexp( - 2.0 * temp2 )
  coef4 = dtan( temp1 ) + dexp( - 2.0 * temp2 ) * dtan( temp1 )

  temp1 = coef1 * coef3 + coef2 * coef4
  temp2 = coef2 * coef3 - coef1 * coef4
  temp3 = coef3 * coef3 + coef4 * coef4

  rn(0) = ( temp1 / temp3 ) + ci * ( temp2 / temp3 )

  do i = 1 , n
    ctemp1 = ( 2 * i - 1 ) / ( ri * x ) - rn(i-1)
    rn(i) = 1.0d0 / ctemp1
  end do

  !----------------------------------------------------------------------------------------!
  ! Step 4: Calculations of the Mie coefficients an et bn.

  do i = 1 , n
    ctemp1 = rn(i) / ri + i * ( 1.0d0 - 1.0d0 / ( ri * ri ) ) / x
    ma(i) = ( ctemp1 * phi(i) - phi(i-1) ) / ( ctemp1 * zeta(i) - zeta(i-1) )
    mb(i) = ( rn(i) * ri  * phi(i) - phi(i-1) ) / ( rn(i) * ri  * zeta(i) - zeta(i-1) )
  end do

  !----------------------------------------------------------------------------------------!
  ! Step 5: scattering amplitudes

  do i = 1 , nb_angles
    ctemp2 = (0.0d0,0.0d0)
    ctemp3 = (0.0d0,0.0d0)
    do j = 1 , n
      temp1 = DBLE( 2 * j + 1 ) / DBLE( j * ( j + 1 ) )
      ctemp2 = ctemp2 + temp1 * ( ma(j) * pi(i,j) + mb(j) * tau(i,j) )
      ctemp3 = ctemp3 + temp1 * ( mb(j) * pi(i,j) + ma(j) * tau(i,j) )
    end do
    sa1(i) = ctemp2
    sa2(i) = ctemp3
  end do

  !----------------------------------------------------------------------------------------!
  ! Step 6: phase function.

  !----- norm (need to be improved <- integration too simple)

  norm = 0.0
  do i = 1 , n
    temp1 = abs( ma(i) ) * abs( ma(i) ) + abs( mb(i) ) * abs( mb(i) )
    norm = norm + ( 2.0d0 * i + 1.0d0 ) * temp1
  end do

  do i = 1 , nb_angles
    p11(i) = ( sa1(i) * conjg( sa1(i) ) + sa2(i) * conjg( sa2(i) ) ) / norm
    p12(i) = ( sa2(i) * conjg( sa2(i) ) - sa1(i) * conjg( sa1(i) ) ) / norm
    p33(i) = ( sa2(i) * conjg( sa1(i) ) + sa1(i) * conjg( sa2(i) ) ) / norm
    p34(i) = ci * ( sa1(i) * conjg( sa2(i) ) - sa2(i) * conjg( sa1(i) ) ) / norm
  end do

  !----------------------------------------------------------------------------------------!
  ! extinction efficiency.

  temp1 = 0.0d0
  do i = 1 , n
    temp1 = temp1 + ( 2 * i + 1 ) * dble( ma(i) + mb(i) )
  end do

  qe = ( 2.0d0 / ( x * x ) ) * temp1

  !----------------------------------------------------------------------------------------!
  ! asymmetry factor.

  temp1 = 0.0d0
  do i = 1 , n - 1
    ctemp1 = ma(i) * conjg( ma(i+1) ) + mb(i) * conjg( mb(i+1) )
    ctemp2 = ma(i) * conjg( mb(i) )
    temp2 = dble( i * ( i + 2.0d0 ) ) / dble( i + 1.0d0 )
    temp3 = dble( 2.0d0 * i + 1.0d0 ) / dble( i * ( i + 1.0d0 ) )
    temp1 = temp1 + ( temp2 * dble( ctemp1 ) + temp3 * dble( ctemp2 ) )
  end do

  gfact = ( 2.0d0 / norm ) * temp1

  !----------------------------------------------------------------------------------------!
  ! single scattering albedo.

  albedo = ( 2.0d0 / ( x * x ) ) * norm / qe

  !----------------------------------------------------------------------------------------!
  ! save data in file

  teta = 0.0d0
  open(1,file='mie_code_output.txt')
  write(1,*) x, real(real(ri)), real(AIMAG(ri)) !, wl
  write(1,*) qe,albedo,gfact
!!$  do i = 1 , nb_angles
!!$    write(1,*) teta,p11(i),p12(i),p33(i),p34(i)
!!$    teta = teta + teta_resol
!!$  end do
  close(1)

end program mie_code
