!+
! Program to produce a spline table of the integrated energy probability function.
!

program create_splines

use spline_mod
use sim_utils

implicit none

type (spline_struct), allocatable, target :: prob_spline1(:), pl_spline1(:), pc_spline1(:), pl45_spline1(:)
type (spline_struct), pointer :: es

real(rp) dE_spline_max, dP_spline_max, dum
real(rp), allocatable :: x(:), x2(:,:), y2(:, :)

integer i, j, n, nx, ny, ix, num_rows_energy, num_rows_y_angle, num_rows_x_angle, ios
integer spline_space_dimensions

logical ok

character(20) source_type
character(100) spectra_energy_file, energy_spline_file, master_dir, file_name
character(100) spectra_angle_file, path, angle_spline_file

namelist / master_params / source_type, dE_spline_max, dP_spline_max, num_rows_energy, &
                num_rows_x_angle, num_rows_y_angle, spline_space_dimensions

! Read parameters

call cesr_getarg (1, master_dir)
n = len_trim(master_dir)
if (master_dir(n:n) /= '/') master_dir = trim(master_dir) // '/'
file_name = trim(master_dir) // 'spline.params'
open (1, file = file_name, status = 'old')
read (1, nml = master_params)
close (1)

!--------------------------------------------------
! Energy spline.
! Read in data points. First 10 lines are header.

spectra_energy_file = trim(master_dir) // 'spectra_output/energy.dc0'
open (1, file = spectra_energy_file, status = 'old')
do i = 1, 10
  read (1, *)
enddo

allocate (prob_spline1(num_rows_energy))

n = 0
do i = 1, num_rows_energy
  es => prob_spline1(n+1)
  read (1, *, iostat = ios) es%x, dum, dum, dum, dum, es%y
  if (ios /= 0) exit
  if (n > 0) then
    if ((es%x - prob_spline1(n)%x) < dE_spline_max .and. (prob_spline1(n)%y - es%y) < dP_spline_max * prob_spline1(1)%y) cycle
  endif
  n = n + 1
enddo

close (1)

! Write spline

prob_spline1%y = (prob_spline1(1)%y - prob_spline1%y) / (prob_spline1(1)%y - prob_spline1(n)%y)
call spline_akima (prob_spline1(1:n), ok)

energy_spline_file = trim(master_dir) // 'spline/energy.spline'

open (1, file = energy_spline_file, recl = 128)
write (1, '(a)') '&spline'
do j = 1, n
  write (1, '(2x, a, i0, a, 6es14.6)') 'prob_spline1(', j, ') =', prob_spline1(j)
enddo
write (1, '(a)') '/'

close (1)


deallocate(prob_spline)

!--------------------------------------------------
! Bend vertical angle spline

if (spline_space_dimensions == 2) then

  n = num_rows_y_angle
  allocate (prob_spline1(n), x(n))
  allocate (pl_spline1(n), pc_spline1(n), pl45_spline1(n))

  do i = 1, num_rows_energy  ! Number of angle files = number of energy points
    write (spectra_angle_file, '(2a, i0, a)') = trim(master_dir), 'spectra_output/angle', i, '.dty'

    open (1, file = spectra_angle_file, status = 'old')
    do j = 1, 10   ! Skip header
      read (1, *)
    enddo

    do j = 1, n
      read (1, *) x(j), prob_spline1(j)%y, pl_spline1(j)%y, pc_spline1(j)%y, pl45_spline1(j)%y
    enddo

    close (1)

    x = x * 1d-3 ! Convert mm -> rad (observation distance is 1 meter).
    prob_spline1%x = x
    pl_spline1%x   = x
    pc_spline1%x   = x
    pl45_spline1%x = x

    ! Integrate and normalize probability

    call integrate_and_normalize_prob (prob_spline1)

    call spline_akima (prob_spline1(1:n), ok)
    call spline_akima (pl_spline1(1:n), ok)
    call spline_akima (pc_spline1(1:n), ok)
    call spline_akima (pl45_spline1(1:n), ok)

    write (angle_spline_file, '(2a, i0, a)') trim(master_dir), 'spline/e', i, '_y.spline'
    open (1, file = angle_spline_file, recl = 128)

    write (1, '(a)') '&spline'
    do j = 1, n
      write (1, '(2x, a, i0, a, 6es14.6)') 'prob_spline(', j, ') =', prob_spline1(j)
    enddo

    do j = 1, n
      write (1, '(2x, a, i0, a, 6es14.6)') 'pl_spline(', j, ') =', pl_spline1(j)
    enddo

    do j = 1, n
      write (1, '(2x, a, i0, a, 6es14.6)') 'pc_spline(', j, ') =', pc_spline1(j)
    enddo

    do j = 1, n
      write (1, '(2x, a, i0, a, 6es14.6)') 'pl45_spline(', j, ') =', pl45_spline1(j)
    enddo

    write (1, '(a)') '/'

    close (1)

  enddo

  stop

endif

!--------------------------------------------------
! Non-bend far field angle spline 

if (spline_space_dimensions == 3) then

  nv = num_rows_y_angle
  nx = num_rows_h_angle

  allocate (prob_spline2(nx,ny), x2(nx,ny), y2(nx,ny))
  allocate (pl_spline2(nx, ny), pc_spline2(nx, ny), pl45_spline2(nx, ny))

  do i = 1, num_rows_energy  ! Number of angle files = number of energy points
    write (spectra_angle_file, '(2a, i0, a)') trim(master_dir), 'spectra_output/angle', i, '.dta'

    open (1, file = spectra_angle_file, status = 'old')
    do j = 1, 10   ! Skip header
      read (1, *)
    enddo

    do k = 1, ny
      do j = 1, nx
        read (1, *) x2(j,k), y2(j,k), prob_spline2(j,k)%y, pl_spline2(j,k)%y, pc_spline2(j,k)%y, pl45_spline2(j,k)%y
      enddo
    enddo

    allocate (prob_spline1(nx))

    do j = 1, nx
      prob_spline1(j)%y = sum(prob_spline2(j,:))%y
      prob_spline1%x = x2(j,1) * 1e3
    enddo

    ! Integrate and normalize probability

    call integrate_and_normalize_prob (prob_spline1)

    call spline_akima (prob_spline1(1:nx), ok)     

    write (angle_spline_file, '(2a, i0, a)') trim(master_dir), 'spline/e', i, '_x.spline'
    open (1, file = angle_spline_file, recl = 128)

    write (1, '(a)') '&spline'
    do j = 1, n
      write (1, '(2x, a, i0, a, 6es14.6)') 'prob_spline(', j, ') =', prob_spline1(j)
    enddo
    write (1, '(a)') '/'

    close (1)

    deallocate (prob_spline1)

    !

    allocate (prob_spline1(ny), pl_spline1(ny), pc_spline1(ny), pl45_spline1(ny))

    prob_spline1%x = y(1,:)
    pl_spline1%x   = y(1,:)
    pc_spline1%x   = y(1,:)
    pl45_spline1%x = y(1,:)

    do j = 1, nx
      prob_spline1%y = prob_spline2(j,:)%y
      pl_spline1%y   = pl_spline2(j,:)%y
      pc_spline1%y   = pc_spline2(j,:)%y
      pl45_spline1%y = pl45_spline2(j,:)%y

      call integrate_and_normalize_prob(prob_spline1)

      call spline_akima (prob_spline1(1:n), ok)
      call spline_akima (pl_spline1(1:n), ok)
      call spline_akima (pc_spline1(1:n), ok)
      call spline_akima (pl45_spline1(1:n), ok)

      write (angle_spline_file, '(2a, i0, a, i0, a)') trim(master_dir), 'spline/e', i, '_x', j, '_y.spline'
      open (1, file = angle_spline_file, recl = 128)

      write (1, '(a)') '&spline'
      do j = 1, n
        write (1, '(2x, a, i0, a, 6es14.6)') 'prob_spline(', j, ') =', prob_spline1(j)
      enddo

      do j = 1, n
        write (1, '(2x, a, i0, a, 6es14.6)') 'pl_spline(', j, ') =', pl_spline1(j)
      enddo

      do j = 1, n
        write (1, '(2x, a, i0, a, 6es14.6)') 'pc_spline(', j, ') =', pc_spline1(j)
      enddo

      do j = 1, n
        write (1, '(2x, a, i0, a, 6es14.6)') 'pl45_spline(', j, ') =', pl45_spline1(j)
      enddo

      write (1, '(a)') '/'

      close (1)

    enddo

  enddo

  stop

endif

!------------------------------------------------------------------
contains

subroutine integrate_and_normalize_prob (spline)

type (spline_struct) spline(:)
real(rp), allocatable :: dy(:)
integer n, j

n = size(spline)
allocate (dy(n+1))
dy = [0.0_rp, spline%y] + [spline%y, 0.0_rp]

do j = 2, n
  spline(j)%y = spline(j-1)%y + dy(j)
enddo
spline%y = (spline(1:n)%y - spline(1)%y) / (spline(n)%y - spline(1)%y)

end subroutine integrate_and_normalize_prob

end program