!! Function to map projected (ray-tracing) coordinates to un-projected (clustering) coordinates in the KiDS halo mocks

function   CorrectCoordinates(X_arcmin, Y_arcmin, z_halo, z_lens)
  implicit none
  
  real(4) X_arcmin, Y_arcmin, X_arcmin_shifted, Y_arcmin_shifted
  real(4) z_halo, z_lens, delta_chi, chi_halo, chi_lens, delta_X_mpch, delta_Y_mpch, chi
  real(4), dimension(2) :: CorrectCoordinates

  real, parameter :: Omega_m = 0.2905
  real, parameter :: angle = 10.0
  real, parameter :: h = 0.6898
  real, parameter :: pi = 3.14159

  ! 1- Run basic checks
  if(X_arcmin < 0. .or. Y_arcmin < 0. .or. X_arcmin > angle*60. .or. Y_arcmin > angle*60. ) then
     write(*,*) 'Angle out of range', X_arcmin, Y_arcmin
     stop
  endif

  if(z_halo < 0.0 .or. z_halo > 3.1) then 
     write(*,*) 'Redshift out of range: z_halo = ', z_halo
     stop
  endif

  !write(*,*) 'Computing correction for ', Omega_m, h, angle, pi, X_arcmin, Y_arcmin, z_halo, z_lens

  ! 2 - Convert z to chi
  chi_lens = chi(z_lens, Omega_m, h)*h
  chi_halo = chi(z_halo, Omega_m, h)*h

  ! 3- Compute delta_chi: distance between halo and mass plane (no correction there)
  delta_chi = chi_halo - chi_lens

  ! 4- Shift angles for light cone to be centered at the center (no correction there)
   X_arcmin_shifted = X_arcmin - (angle/2.)*60.
   Y_arcmin_shifted = Y_arcmin - (angle/2.)*60.

  ! 5- Compute shift in comoving Mpc/h of the halo, caused by projection
  delta_X_mpch = tan(X_arcmin_shifted/60.*pi/180.)*delta_chi
  delta_Y_mpch = tan(Y_arcmin_shifted/60.*pi/180.)*delta_chi

  ! 6- Compute shift in apparent angle:
  CorrectCoordinates(1) = atan(delta_X_mpch/chi_halo)*180./pi*60.
  CorrectCoordinates(2) = atan(delta_Y_mpch/chi_halo)*180./pi*60.

end function CorrectCoordinates
