! Module with some useful functions and global variables.
! This is used to automatically create a function interface.
module utils
  implicit none

  real(8) :: r_max  ! radius of the source
  real(8) :: pi     
  real(8) :: s      ! source-slits distance
  real(8) :: D      ! slit-screen distance
  real(8) :: l      ! slit-slit distance
  real(8) :: lambda ! wavelength

  real(8) :: r_min  ! only used for ring source

  real(8), dimension(:,:), allocatable :: source  ! custom source
  integer :: SX, SY  ! custom source dimensions
  real(8) :: s_scale ! custom source scale

contains

  ! Loads a matrix from file f in the array A
  ! Writes the size of A in x and y
  subroutine load_file(f, A, x, y)
    implicit none

    real(8), dimension(:,:), allocatable :: A
    character(len=*) :: f
    integer :: x,y,i,j,xi,yj

    open(1, file=f)

    read(1,*) x, y

    allocate(A(0:x-1,0:y-1))

    do i = 0, x - 1
      do j = 0, y - 1
        read(1,*) xi,yj,A(xi,yj)
      enddo
    enddo

    close(1)

  end

end

! Computes an integral using the trapezoidal rule
real(8) function trapeze(f,a,b,n)

  integer :: n,i
  real(8) :: a,b,res
  real(8),external :: f

  res = (f(a) + f(b))/2.0

  do i = 1, n - 1
    res = res + f(a + real(i) * (b - a) / real(n))
  enddo

  trapeze = res * (b - a) / real(n)
end


! Computes a 2D integral using the trapezoidal rule
real(8) function trapeze2(f,a,b,c,d,n,m)

  integer :: n,m,i,j
  real(8) :: a,b,c,d,res,dx,dy
  real(8),external :: f

  res = (f(a,c) + f(a,d) + f(b,c) + f(b,d)) / 4.0
  
  dx = (b - a) / real(n)
  dy = (d - c) / real(m)

  do i = 1, n - 1
    res = res + (f(a + real(i)*dx, c) + f(a + real(i)*dx, d)) / 2.0
  enddo

  do j = 1, m - 1
    res = res + (f(a, c + real(j)*dy) + f(b, c + real(j)*dy)) / 2.0
  enddo

  do i = 1, n - 1
    do j = 1, m - 1
      res = res + f(a + real(i)*dx, c + real(j)*dy)
    enddo
  enddo

  trapeze2 = res * dx * dy
end


! Computes a 2D complex integral using the trapezoidal rule
complex(8) function trapeze2C(f,a,b,c,d,n,m)

  integer :: n,m,i,j
  real(8) :: a,b,c,d,dx,dy
  complex(8)::res
  complex(8),external :: f

  res = (f(a,c) + f(a,d) + f(b,c) + f(b,d)) / 4.0
 
  dx = (b - a) / real(n)
  dy = (d - c) / real(m)

  do i = 1, n - 1
    res = res + (f(a + real(i)*dx, c) + f(a + real(i)*dx, d)) / 2.0
  enddo

  do j = 1, m - 1
    res = res + (f(a, c + real(j)*dy) + f(b, c + real(j)*dy)) / 2.0
  enddo

  do i = 1, n - 1
    do j = 1, m - 1
      res = res + f(a + real(i)*dx, c + real(j)*dy)
    enddo
  enddo

  trapeze2C = res * dx * dy
end


! Some example functions
real(8) function my_sine(x)
  real(8) :: x
  my_sine = sin(x)
end

real(8) function expx2(x)
  real(8) :: x
  expx2 = exp(-x*x)
end

real(8) function my_cosine(x)
  real(8) :: x
  my_cosine = cos(x)
end

real(8) function my_sine2(x,y)
  real(8) :: x,y
  my_sine2 = sin(x*y)
end


! Optical path difference
real(8) function optical_path(r, phi, X, Y)
  use utils
  real(8) :: r, phi, X, Y, delta, u, v

  u = r * cos(phi)
  v = r * sin(phi)

  delta = sqrt(X*X + (Y + l/2.0)**2.0 + D*D) - sqrt(X*X + (Y - l/2.0)**2.0 + D*D)
  delta = sqrt(u*u + (v + l/2.0)**2.0 + s*s) - sqrt(u*u + (v - l/2.0)**2.0 + s*s) + delta

  optical_path = delta
end


! Circular uniform intensity
real(8) function uniform_intensity(r, phi)
  use utils
  real(8) :: r,phi

  if (r < r_max) then
    uniform_intensity = 1.0 / (pi * r_max * r_max)
  else
    uniform_intensity = 0.0
  endif
end


! Ring uniform intensity
real(8) function ring_intensity(r, phi)
  use utils
  implicit none
  real(8) :: r,phi

  if ((r < r_max) .and. (r > r_min)) then
    ring_intensity = 1.0 / ((pi * r_max * r_max) - (pi * r_min * r_min))
  else
    ring_intensity = 0.0
  endif
end


! Custom source intensity from the source array
real(8) function custom_intensity(r, phi)
  use utils
  implicit none

  real(8) :: r,phi,x,y
  integer :: xi, yi

  x = r * cos(phi) / s_scale
  y = r * sin(phi) / s_scale

  xi = int(x) + SX / 2
  yi = int(y) + SY / 2

  if((xi < SX) .and. (xi >= 0) .and. (yi < SY) .and. (yi >= 0)) then
    custom_intensity = source(xi, yi)
  else
    custom_intensity = 0
  endif
end


! Computes the intensity at a point of the screen
! given an intensity function
real(8) function screen_intensity(I,X,Y)
  use utils
  implicit none
  real(8) :: X,Y
  real(8), external :: I, trapeze2, optical_path

  screen_intensity = 2.0 * trapeze2(pointwise_intensity, real(0.,8), r_max, real(0.,8), 2.*pi, 40, 40)

  contains 

  real(8) function pointwise_intensity(r, phi)
    real(8) :: r,phi

    pointwise_intensity = r * I(r,phi) * (1.0 + cos(2.0 * pi * optical_path(r,phi,X,Y) / lambda))

  end

end


! Computes the amplitude A received by the screen given an intensity function
complex(8) function screen_amplitude(I)
  use utils
  implicit none
  real(8), external :: I, optical_path
  complex(8), external :: trapeze2C

  screen_amplitude = 2.0 * trapeze2C(pointwise_amplitude, real(0.,8), r_max, real(0.,8), 2.*pi, 200, 200)

  contains 

  complex(8) function pointwise_amplitude(r, phi)
    real(8) :: r,phi

    pointwise_amplitude = r * I(r,phi) * exp(complex(real(0.0,8), 2.0 * pi *& 
                          optical_path(r,phi,real(0.,8),real(0.,8)) / lambda))

  end

end


! Computes the contrast at a point of the screen 
! given an intensity function
real(8) function contrast(I)
  use utils
  implicit none
  real(8) :: total
  real(8), external :: I, trapeze2
  complex(8), external :: screen_amplitude
  complex(8) :: amp

  amp = screen_amplitude(I)
  total = trapeze2(pointwise_total, real(0.,8), r_max, real(0.,8), 2.*pi, 40, 40)

  contrast = abs(amp)/total

  contains 

  real(8) function pointwise_total(r, phi)
    real(8) :: r,phi

    pointwise_total = r * I(r,phi) 

  end

end


! Main program
program main
  use utils
  implicit none

  integer :: x,y,n,m
  real(8),external :: trapeze, trapeze2, my_sine, my_cosine, expx2, my_sine2
  real(8),external :: screen_intensity, uniform_intensity, optical_path 
  real(8),external :: custom_intensity, contrast, ring_intensity
  real(8) :: last_val, val

  pi     = acos(-1.0)
  r_max  = 1e-3
  r_min  = 5e-4
  s      = 1000.
  D      = 10.
  l      = 1.
  lambda = 5e-4

  call load_file('resources/elliptical.data', source, SX, SY)
  r_max   = 0.002
  s_scale = 1e-4


  ! We compute the value of l such that the contrast is near zero
  last_val = 1000.
  do x = 1000, 4000
    l = real(x,8) / 10.
    val = contrast(custom_intensity)
    write(*,*) val
    if (last_val < val) then
      exit
    else
      last_val = val
    endif
  enddo

  write(*,*) last_val, l

end
