SOS2DFTN

! This file is part of SOS2DFTN.
! Copyright (C) 2014 Adam Hirst <[email protected]>
!
! SOS2DFTN is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! SOS2DFTN is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with SOS2DFTN. If not, see <http://www.gnu.org/licenses/>.

module Constants
  implicit none

  real, parameter :: pi = 4 * atan(1.0)

end module Constants

module Types
  use Constants
  implicit none

  type Emitter
    ! The extended emitter lies on the x-axis
    ! edges of element are `-R` and `+R`, units in mm
    ! `-R` also referred to as `R'` in literature
    real :: R
    ! prescribed response at the z-axis
    real :: L_z
    ! prescribed response is defined from the z-axis down to some minimum angle
    real :: theta_min
  contains
    ! prescribed response in terms of the emitter's angular projected width
    procedure :: L => Response
  end type Emitter

  type SurfacePoint
    ! points exist in x-z plane, units in mm
    real :: P(2)
    ! we also store the normal unit-vector
    real :: N(2)
  end type SurfacePoint

  type Surface
    ! the surface points themselves
    type(SurfacePoint), allocatable :: A(:)
    ! total number of points comprising the surface
    integer :: J
  end type Surface

  type, extends(Surface) :: Lens
    ! Refractive index of the lens.
    real    :: n_1
    ! Number of oval-segments for the initial lens portion, where # oval points = # oval segments + 1
    integer :: N
  contains
    ! The surface is "seeded" with a Cartesian oval from the x-axis up to the emitter's minimum prescription angle.
    procedure :: Seed => SeedLens
    procedure :: Generate
  end type Lens

contains

  pure function Reflect(i, n) result (rfx)
    ! Returns the unit vector of a ray after reflection, given unit vectors of the `i`ncident ray and the surface `n`ormal.
    ! "Introduction to Nonimaging Optics" (Chaves, 2008), Section 17.15
    real, intent(in) :: i(2), n(2)
    real             :: rfx(2)

    rfx = i - 2*dot_product(i,n)*n
  end function Reflect

  pure function Refract(i, n_S, n_1) result (rfr)
    ! Returns the unit-vector of a ray after refraction, given unit vectors of the `i`ncident ray and surface normal `n_S`, and the
    ! refractive index of the incident medium `n_1`.
    ! We assume that the second medium is air, i.e. that `n_2 = 1`.
    ! "Introduction to Nonimaging Optics" (Chaves, 2008), Section 17.15
    real, intent(in) :: i(2), n_S(2), n_1
    real             :: rfr(2), n(2), delta

    if (dot_product(i,n_S) >= 0) then
      n = n_S
    else
      n = -n_S
    end if
    delta = 1 - (n_1**2) * (1 - dot_product(i,n)**2)
    if (delta > 0) then
      rfr = i*n_1 + n*(sqrt(delta) - n_1*dot_product(i,n))
    else
      rfr = Reflect(i,n_S)
    end if
  end function Refract

  pure function IntersectLines(P, v, Q, u) result (isl)
    ! Returns the point of intersection of two lines, each defined as a point (`P` and `Q`) and a direction (`v` and `u`).
    ! "Introduction to Nonimaging Optics" (Chaves, 2008), Section 17.15
    real, intent(in) :: P(2), v(2), Q(2), u(2)
    real             :: isl(2), u_P(2)

    u_P = [-u(2), u(1)]
    isl = P + dot_product((Q-P),u_P)/dot_product(v,u_P) * v
  end function IntersectLines

  pure function RefractionNormal(i, r, n_1) result (rfrnrm)
    ! Returns the surface normal unit-vector given that an `i`ncident ray is `r`efracted from a medium of index `n_1` into air.
    ! "Introduction to Nonimaging Optics" (Chaves, 2008), Section 17.15
    real, intent(in) :: i(2), r(2), n_1
    real             :: rfrnrm(2)

    rfrnrm = n_1*i - r; rfrnrm = rfrnrm / norm2(rfrnrm)
  end function RefractionNormal

  pure function CartesianOval(F, n_1, P, alpha, phi) result (cop)
    ! Returns the surface point of the Cartesian oval which collimates rays from `F` in a medium of refractive index `n_1` into air
    ! at an angle `alpha`, which intersects the point `P`, in terms of `phi` (the angle about `F` above `alpha`).
    ! "Introduction to Nonimaging Optics" (Chaves, 2008), Section 17.15
    real, intent(in) :: F, n_1, P, alpha, phi
    ! Note that `F` and `P` are scalar, as are defined to lie on the x-axis.
    ! We additionally assume that `F` lies on `-x` and that P lies on `+x`.
    real             :: cop(2)

    cop = (P - F)*(n_1 - cos(alpha))*[cos(phi+alpha), sin(phi+alpha)] / (n_1 - cos(phi)) + [F,0.0]
  end function CartesianOval

  pure function CartesianOvalNormal(F, n_1, P, alpha, phi) result (copn)
    ! Returns the unit-vector of the normal to the surface of the above-described Cartesian oval at a given parameter `phi`.
    real, intent(in) :: F, n_1, P, alpha, phi
    real             :: copn(2)

    copn = (P - F)*(n_1 - cos(alpha))*[n_1*cos(phi+alpha) - cos(alpha), n_1*sin(phi+alpha) - sin(alpha)] / (n_1 - cos(phi))**2
    copn = copn / norm2(copn)
  end function CartesianOvalNormal

  function CartesianOvalProjection(F, n_1, P, alpha, L) result (phi)
    ! Returns the parameter `phi` corresponding to the end of the Cartesian oval with projected width `L` in the direction `alpha`.
    real, intent(in) :: F, n_1, P, alpha, L
    real             :: phi, C, phis(2), phi_C

    C = (P - F)*(n_1 - cos(alpha)) / ((P - F)*sin(alpha) - L)
    if (C**2 < n_1**2 - 1) then
      stop 'Projection and oval do not intersect.'
    end if
    ! `[1,-1]` indicates `\pm`, as lines intersect closed curves twice.
    phis(:) = 2 * atan( ([1, -1]*sqrt(C**2 - n_1**2 + 1.0) - C) / (n_1+1.0) )
    ! The intersecting line necessarily has positive slope (`alpha` is always in the positive x-z quadrant), so take the
    ! intersection closest to `phi = 0`.
    phi = phis( minloc(abs(phis),1) )
    ! `phi` must also correspond to a point on the oval within a valid range, namely that which permits refraction onto `F`.
    phi_C = pi/2 - asin(min(1.0,n_1)/max(1.0,n_1))
    if (abs(phi) > phi_C) then
      stop 'Projection does not intersect with refractive portion of oval.'
    else if (phi < -alpha) then
      stop 'Projection does not intersect above the x-axis.'
    end if
  end function CartesianOvalProjection

  pure real function Response(this, theta) result (L)
    ! Returns the prescribed angular response of the emitter at the given `theta`.
    ! "Nonimaging Optics" (Winston et al., 2005), Section 7.4.2
    class(Emitter), intent(in) :: this
    real,           intent(in) :: theta

    if (theta >= this%theta_min) then
      L = this%L_z * (1 + (2*theta - pi)/(2*this%theta_min - pi))
    else
      L = 0.0
    end if
  end function Response

  subroutine SeedLens(this, my_emitter, A_0)
    class(Lens),   intent(in out) :: this
    type(Emitter), intent(in)     :: my_emitter
    real,          intent(in)     :: A_0
    real                          :: phi_max
    integer :: i

    associate (R => my_emitter%R, n_1 => this%n_1, alpha => my_emitter%theta_min, N => this%N)
      phi_max = CartesianOvalProjection(-R, n_1, A_0, alpha, my_emitter%L(alpha))
      allocate(this%A(0:N))
      this%A(0)%P = [A_0, 0.0]
      do concurrent (i = 1:N)
        this%A(i)%P = CartesianOval(-R, n_1, A_0, alpha, i*(phi_max+alpha)/N - alpha)
      end do
      do concurrent (i = 0:N)
        this%A(i)%N = CartesianOvalNormal(-R, n_1, A_0, alpha, i*(phi_max+alpha)/N - alpha)
      end do
      this%J = N
    end associate
  end subroutine SeedLens

  pure subroutine Generate(this, my_emitter)
    class(Lens),         intent(in out)              :: this
    type(Emitter),       intent(in)                  :: my_emitter
    integer                                          :: i
    real                                             :: i_i(2), r_i(2), theta_i
    type(SurfacePoint)                               :: B
    type(SurfacePoint),                  allocatable :: A_new(:)

    ! for `i = 1` upwards
    i = 1
    do
      ! Refract from `R` through `A(i)` to get projection angle `theta_i` (and `r_i`, the unit-vector in that direction).
      ! `i_i` and `r_i` are the `i`th `i`ncident and `r`efracted ray directions respectively
      i_i = (this%A(i)%P - [my_emitter%R, 0.0]); i_i = i_i / norm2(i_i)
      r_i = Refract(i_i, this%A(i)%N, this%n_1)
      ! `r_i` is a unit vector, and unit vectors are equal to `[cos(theta), sin(theta)]`
      theta_i = acos(r_i(1))
      ! Intersect the line parallel to `r_i` perpendicularly separated from `A(i)` by `L(theta_i)` with the tangent to `A(N+i-1)` to
      ! get the next point in the surface, `A(N+i)`.
      associate (L => my_emitter%L(theta_i), P => this%A(i)%P, Q => this%A(this%N+i-1)%P, N => this%A(this%N+i-1)%N)
        ! TODO: Use a different offset for `Q` to avoid the need to compute `sin(theta_i)`
        B%P = IntersectLines(P - [L/sin(theta_i), 0.0], r_i, Q, [-N(2), N(1)])
      end associate
      ! Compute the normal to `A(N+i)`, from the fact that light from `R'` through `A(N+i)` must have resulting direction `theta_L`.
      i_i = B%P + [my_emitter%R, 0.0]; i_i = i_i / norm2(i_i)
      B%N = RefractionNormal(i_i, r_i, this%n_1)
      ! if the new point has not crossed the z-axis, simply store it
      if (B%P(1) >= 0.0) then
        ! TODO: implement chunked reallocation of `type(SurfacePoint)` array
        this%J = this%J + 1
        allocate(A_new(0:this%J))
        A_new(0:) = [this%A(0:), B]
        call move_alloc(A_new, this%A)
      else
        ! TODO: store the intersection of the new surface segment with the z-axis, and then terminate
        return
      end if
      ! keep increasing `i` until done
      i = i + 1
    end do
  end subroutine Generate

end module Types

program SOS2DFTN
  use Types
  implicit none

  type(Emitter) :: LED = Emitter(1.0, 6.0, pi/6)
  type(Lens)    :: glass
  integer       :: i

  ! TODO: %Seed should also store `N` and `n_1` in the `type(Lens)` object
  glass%N = 50; glass%n_1 = 1.49
  call glass%Seed(LED, 15.0)
  call glass%Generate(LED)
  do i = 0, glass%J
    print *, 'Point: ', glass%A(i)%P
    if (i == glass%N) print *, '##########'
  end do

end program SOS2DFTN