multidip_rabitt_delays.F90

Go to the documentation of this file.

! Copyright 2021
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-out (UKRmol+ suite).
!
!     UKRmol-out 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.
!
!     UKRmol-out 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  UKRmol-out (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
!
 
!> \brief   Calculate RABITT delays from raw multipole elements
!> \author  J Benda
!> \date    2021 - 2024
!>
!> Having (spherical or Cartesian) raw multipole elements calculated by MULTIDIP for the two different absorption/emission pathways,
!> evaluate the corresponding discrete time delays, unresolved (i.e. averaged) in emission direction as well as in molecular
!> orientation. The raw multipole elements filenames are provided on the command line and they are expected to bear either the
!> original file name assigned by MULTIDIP or additional stem suffixes, for instance "-smooth" as added by the utility program
!> *multidip_smooth_rawdips*. Example usage:
!> ```
!>   multidip_rabitt_delays --cartesian \
!>                          --two-photon \
!>                          --xuv-plus-ir <(ls multidip.2.avg.XUV+IR.r100/rawdip-[xyz]*-*) \
!>                          --xuv-minus-ir <(ls multidip.2.avg.XUV-IR.r100/rawdip-[xyz]*-*)
!> ```
!> See the program help for details. Note that the program can be used to calculate one-photon delays, too. For this to work,
!> one needs to produce two sets of the one-photon raw transition dipoles that are offset in energies by 2*omega_IR. These
!> two sets are then used with `--xuv-plus-ir` and `--xuv-minus-ir` (with smaller and larger energies, respectively). The results
!> are then for energies in between.
!>
!> When the parameters `--polar`, `--theta` and `--phi` are specified, molecular frame delays will be calculated as well. The
!> parameter `--polar` specifies the direction of the polarization axis the polar and azimuthal angle in degrees. The parameters
!> `--theta` and `--phi` specify the emission directions in the molecular frame; arbitrary number of angles in degrees can be given.
!> If one of these is omitted, it is assumed to be zero. If the two arrays of angles have unequal lengths, the shorter one is
!> cycled over as long as necessary. When `--theta` and/or `--phi` are given, the program will produce filws with the RABITT
!> interference term called "rabitt-signal-mol-?.txt", where "?" is the index of the emission direction. Then there weill be
!> two more sets
!> of results:
!>    - "rabitt-signal-axs-?.txt" - integrated over orientations of the polarization, preserving the polar angle of the polarization
!>      vector
!>    - "rabitt-signal-rec-?.txt" - integrated over orientations of the polarizatoin (both angles)
!> These two semi-averaged results also include integration over azimuthal emission angle. Molecular frame output is only
!> computed with Cartesian angular basis.
!>
program multidip_rabitt_delays
 
    use iso_fortran_env, only: real64
#ifdef HAVE_NO_FINDLOC
    use algorithms,      only: findloc
#endif
 
    implicit none
 
    type string
        character(:), allocatable :: str
    end type string
 
    call rabitt_main
 
contains
 
    !> \brief   Print program usage information
    !> \author  J Benda
    !> \date    2021 - 2024
    !>
    subroutine print_help
 
        print '(A)', 'Calculate discrete RABITT time delays using the two- (or multi-)photon raw dipole elements'
        print '( )'
        print '(A)', 'Usage:'
        print '( )'
        print '(A)', '    multidip_rabitt_delays [OPTIONS] --xuv+ir <FILE1> --xuv-ir <FILE2>'
        print '( )'
        print '(A)', 'General options:'
        print '( )'
        print '(A)', '    --help                  Display this help'
        print '(A)', '    --one-photon            Calculate one-photon delays (equivalent to --orders 1 1)'
        print '(A)', '    --two-photon            Calculate two-photon delays (equivalent to --orders 2 2)'
        print '(A)', '    --orders M N            Calculate mixed-photon delays'
        print '(A)', '    --extend P Q            Extend the orders to M+|P| and N+|Q| using per-pw asymptotic approximation'
        print '(A)', '    --spherical             Assume matrix elements in the spherical basis'
        print '(A)', '    --cartesian             Assume matrix elements in the Cartesian basis'
        print '(A)', '    --rawdips FILE1 FILE2   File with paths to files containing raw matrix elements for the XUV+/-IR pathways'
        print '( )'
        print '(A)', 'Higher-order extension control options:'
        print '( )'
        print '(A)', '    --ir-energy EIR         Energy of IR photon (in au, needed for --extend)'
        print '(A)', '    --properties FILE       Path to the residual ion properties file (needed'
        print '(A)', '    --energies FILE1 FILE2  Files with photoelectron energies (in au) in the 1st channel (for --extend)'
        print '(A)', '    --out-dir DIR1 DIR2     Directory for output of XUV+/-IR matrix elements calculated via --extend'
        print '( )'
        print '(A)', 'Angular distribution control options:'
        print '( )'
        print '(A)', '    --max-lab MAXL          Maximal angular momentum for asymmetry parameer (default: 0)'
        print '(A)', '    --polar THETA PHI       Polarization direction of XUV in molecular frame (degrees, default: 0 0, i.e. z)'
        print '(A)', '    --polar-ir THETA PHI    As above but for IR, by default oriented in the same direction as XUV'
        print '(A)', '    --theta ANGLE(s)        Polar emission angles (degrees, default: 0)'
        print '(A)', '    --phi ANGLE(s)          Azimuthal emission angles (degrees, default: 0)'
        print '( )'
        print '(A)', 'The files containing the list of paths to raw dipole data are assumed to contain one path on each line.&
                     & A convenient way how to pass all files of a given mask to this program in some Unix-like shells is using&
                     & the so-called process substitution like this:'
        print '( )'
        print '(A)', '    multidip_rabitt_delays --rawdips <(ls XUV+IR/rawdip-[xyz]-*) <(ls XUV-IR/rawdip-[xyz]-*)'
        print '( )'
        print '(A)', 'Note that the file names provided on the command line need to satisfy the original MULTIDIP naming scheme,&
                     & so that it is possible to extract the partial wave quantum numbers.'
        print '( )'
        print '(A)', 'The resulting RABITT cross section interference term will be written to file "rabitt-signal-L0-xyz.txt" or&
                     & "rabitt-signal-L0-sph.txt" depending on the selected mode. To obtain the resulting&
                     & molecular orientation averaged photoionization delays, one needs to take the complex argument and divide&
                     & it by 2*omega_IR.'
 
    end subroutine print_help
 
 
    !> \brief   Main program
    !> \author  J Benda
    !> \date    2021 - 2024
    !>
    !> Read the intermediate data from disk and produces the (orientation averaged) delays.
    !>
    subroutine rabitt_main
 
        use iso_fortran_env,   only: output_unit
        use mpi_gbl,           only: mpi_mod_start, mpi_mod_finalize
        use multidip_io,       only: read_target_properties
        use multidip_routines, only: calculate_quadratic_dipole_sph, calculate_quadratic_dipole_xyz, convert_xyz_to_sph
        use multidip_params,   only: nmaxphotons
 
        type(string), allocatable :: XUV_plus_IR_filenames(:), XUV_minus_IR_filenames(:)
        type(string) :: XUV_plus_IR_outdir, XUV_minus_IR_outdir, properties
 
        complex(real64), allocatable :: MM(:, :), MM0(:,:)
        complex(real64), allocatable :: M_XUV_plus_IR_sph(:, :, :, :), M_XUV_minus_IR_sph(:, :, :, :)
        complex(real64), allocatable :: M_XUV_plus_IR_xyz(:, :, :, :), M_XUV_minus_IR_xyz(:, :, :, :)
 
        real(real64), parameter   :: rone = 1, pi = 4*atan(rone)
        real(real64), allocatable :: theta(:), phi(:), Ek1(:), Ek2(:), prop(:, :, :), echl(:)
        real(real64)              :: polar(2), polarIR(2), axis(3), axisIR(3), wIR
 
        integer, allocatable :: target_states(:), chains1(:, :), chains2(:, :)
 
        integer :: i, L, order1, order2, extend1, extend2, maxl, ntarg, nesc, maxlab, lutarg, p(nMaxPhotons)
        logical :: sph
        character(len=200) :: str
 
        if (command_argument_count() == 0) then
            call print_help
            return
        end if
 
        call print_ukrmol_header(output_unit)
        call mpi_mod_start
 
        print '(/,A,/)', 'Program RABITT: Calculation of two-photon interference amplitudes'
 
        ! set lab-frame polarisations to zero
        p = 0
 
        if (.not. process_command_line(order1, order2, extend1, extend2, sph, maxlab, polar, polarir, &
                                       theta, phi, ek1, ek2, wir, properties, &
                                       xuv_plus_ir_filenames, xuv_minus_ir_filenames, &
                                       xuv_plus_ir_outdir, xuv_minus_ir_outdir)) then
            return
        end if
 
        ! read raw multipole elements from disk
        call read_raw_multipoles(order1, order2, sph, ntarg, maxl, nesc, target_states, chains1, chains2, &
                                 xuv_plus_ir_filenames, xuv_minus_ir_filenames, &
                                 m_xuv_plus_ir_xyz, m_xuv_minus_ir_xyz)
        if (nesc == 0) return
 
        ! convert Cartesian matrix elements to spherical basis
        if (sph) then
            if (extend1 /= 0 .or. extend2 /= 0) then
                print '(/,a)', 'Error: Asymptotic extension not implemented for --spherical'
                stop 1
            end if
            m_xuv_plus_ir_sph = m_xuv_plus_ir_xyz
            m_xuv_minus_ir_sph = m_xuv_minus_ir_xyz
        else
            if (extend1 /= 0 .or. extend2 /= 0) then
                if (.not. allocated(ek1) .or. .not. allocated(ek2)) then
                    print '(/,a)', 'Error: Missing photoelectron energies on the command line (--energies) for extending'
                    stop 1
                end if
                if (wir <= 0) then
                    print '(/,a)', 'Error: Need a positive IR energy (--ir-energy) on the command line for extending'
                    stop 1
                end if
                if (len_trim(properties % str) == 0) then
                    print '(/,a)', 'Error: Need the residual ion properties file (--properties) for extending'
                    stop 1
                end if
                open (newunit = lutarg, file = properties % str, action = 'read', form = 'formatted')
                call read_target_properties(lutarg, prop, echl)
                close (lutarg)
                echl = echl - echl(1)
                call extend_order(order1, m_xuv_plus_ir_xyz, chains1, extend1, ek1, wir, &
                                  target_states, echl, prop, maxl, xuv_plus_ir_outdir)
                call extend_order(order2, m_xuv_minus_ir_xyz, chains2, extend2, ek2, wir, &
                                  target_states, echl, prop, maxl, xuv_minus_ir_outdir)
            end if
            allocate (m_xuv_plus_ir_sph, mold = m_xuv_plus_ir_xyz)
            allocate (m_xuv_minus_ir_sph, mold = m_xuv_minus_ir_xyz)
            call convert_xyz_to_sph(m_xuv_plus_ir_xyz, m_xuv_plus_ir_sph, maxl, chains1)
            call convert_xyz_to_sph(m_xuv_minus_ir_xyz, m_xuv_minus_ir_sph, maxl, chains2)
        end if
 
        ! allocate final storage for the quadratic dipoles
        allocate (mm(2 + ntarg, nesc), mm0(2 + ntarg, nesc))
 
        do l = 0, min(maxlab, order1 + abs(extend1) + order2 + abs(extend2))
 
            str = 'rabitt-signal-L'
 
            ! for checking purposes, evaluate the isotropic parameter in the Cartesian basis
            if (l == 0 .and. .not. sph) then
                call calculate_quadratic_dipole_xyz(mm, l, maxl, chains1, chains2, ntarg, nesc, &
                                                    m_xuv_plus_ir_xyz, m_xuv_minus_ir_xyz)
                call write_quadratic_dipoles(l, mm, str, '-xyz')
            end if
 
            ! evaluate the asymmetry parameter in the spherical basis
            call calculate_quadratic_dipole_sph(mm, l, maxl, chains1, chains2, ntarg, nesc, &
                                                m_xuv_plus_ir_sph, m_xuv_minus_ir_sph, p)
            call write_quadratic_dipoles(l, mm, str, '-sph')
 
            if (l == 0) then
                mm0 = mm
            else
                do i = 3, 2 + ntarg
                    where (mm0(i, :) /= 0)
                        mm(i,:) = mm(i,:) / mm0(i,:)
                    elsewhere
                        mm(i,:) = 0
                    end where
                end do
 
                str = 'rabitt-beta-L'
                call write_quadratic_dipoles(l, mm, str, '')
            end if
 
        end do
 
        if (.not. sph) then
            polar = polar * pi / 180
            polarir = polarir * pi / 180
            axis  = [ sin(polar(1))*cos(polar(2)), sin(polar(1))*sin(polar(2)), cos(polar(1)) ]
            axisir = [ sin(polarir(1))*cos(polarir(2)), sin(polarir(1))*sin(polarir(2)), cos(polarir(1)) ]
            call calculate_oriented_dipoles(maxl, chains1, chains2, ntarg, nesc, &
                                            m_xuv_plus_ir_xyz, m_xuv_minus_ir_xyz, axis, axisir, theta, phi)
        end if
 
        print '(/,A)', 'Done.'
 
        call mpi_mod_finalize
 
    end subroutine rabitt_main
 
 
    !> \brief  Read options from the command line
    !> \author J Benda
    !> \date   2022 - 2024
    !>
    !> Obtain command line options, in particular the list of transition multipole files written by MULTIDIP.
    !>
    !> \param[out] ord1    Photon absorption order for conjugated pathway
    !> \param[out] ord2    Photon absorption order for non-conjugated pathway
    !> \param[out] ext1    Requested asymptotic extension of the absorption order in conjugated pathway
    !> \param[out] ext2    Requested asymptotic extension of the absorption order in non-conjugated pathway
    !> \param[out] sph     Whether the spherical basis is used
    !> \param[out] maxlab  Highest angular momentum to consider in the laboratory frame
    !> \param[out] polar   Polarization vector of the XUV field in the molecular frame
    !> \param[out] polarIR Polarization vector of the IR field in the molecular frame
    !> \param[inout] theta List of polar emission angles
    !> \param[inout] phi   List of azimuthal emission angles
    !> \param[inout] Ek1   List of photoelectron kinetic energies on input (needed for extension to higher orders)
    !> \param[inout] Ek2   List of photoelectron kinetic energies on input (needed for extension to higher orders)
    !> \param[inout] properties Target properties file path
    !> \param[out]   wIR   Energy of the IR quantum (needed for extension to higher orders)
    !> \param[inout] XUV_plus_IR_filenames  List of files with transition multipoles for XUV+IR pathway
    !> \param[inout] XUV_minus_IR_filenames List of files with transition multipoles for XUV-IR pathway
    !> \param[inout] XUV_plus_IR_outdir     Output directory for extended XUV+IR amplitudes
    !> \param[inout] XUV_minus_IR_outdir    Output directory for extended XUV-IR amplitudes
    !>
    logical function process_command_line (ord1, ord2, ext1, ext2, sph, maxlab, polar, polarIR, theta, phi, Ek1, Ek2, wIR, &
                                           properties, XUV_plus_IR_filenames, XUV_minus_IR_filenames, XUV_plus_IR_outdir, &
                                           XUV_minus_IR_outdir)
 
        integer,      intent(out) :: ord1, ord2, ext1, ext2, maxlab
        real(real64), intent(out) :: polar(2), polarir(2), wir
        logical,      intent(out) :: sph
 
        real(real64), allocatable, intent(inout) :: theta(:), phi(:), ek1(:), ek2(:)
        real(real64), allocatable                :: real_array(:)
 
        type(string), allocatable, intent(inout) :: xuv_plus_ir_filenames(:), xuv_minus_ir_filenames(:)
        type(string),              intent(inout) :: xuv_plus_ir_outdir, xuv_minus_ir_outdir, properties
        type(string), allocatable                :: string_array(:)
 
        character(len=8192) :: arg, filenames(2)
 
        integer :: iarg, narg, basis, n, ifile, stat, u
 
        narg = command_argument_count()
        iarg = 1
        polar = 0
        polarir = 0
        ord1 = 2
        ord2 = 2
        ext1 = 0
        ext2 = 0
        basis = -1
        maxlab = 0
        wir = 0
        process_command_line = .true.
 
        main_loop: do while (iarg <= narg)
 
            call get_command_argument(iarg, arg)
 
            select case (trim(arg))
 
                case ('--help')
                    call print_help
                    process_command_line = .false.
                    return
 
                case ('--one-photon')
                    ord1 = 1
                    ord2 = 1
 
                case ('--two-photon')
                    ord1 = 2
                    ord2 = 2
 
                case ('--orders')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) ord1
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) ord2
 
                case ('--extend')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) ext1
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) ext2
 
                case ('--cartesian')
                    basis = 1
 
                case ('--spherical')
                    basis = 2
 
                case ('--max-lab')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) maxlab
 
                case ('--polar')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) polar(1)
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) polar(2)
                    polarir = polar
 
                case ('--polar-ir')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) polarir(1)
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) polarir(2)
 
                case ('--theta')
                    iarg = iarg + 1
                    do while (iarg <= narg)
                        call get_command_argument(iarg, arg)
                        if (arg(1:2) == '--') cycle main_loop
                        n = 0
                        if (allocated(theta)) n = size(theta)
                        allocate (real_array(n + 1))
                        if (allocated(theta)) real_array(1:n) = theta(1:n)
                        call move_alloc(real_array, theta)
                        read (arg, *) theta(n + 1)
                        iarg = iarg + 1
                    end do
 
                case ('--phi')
                    iarg = iarg + 1
                    do while (iarg <= narg)
                        call get_command_argument(iarg, arg)
                        if (arg(1:2) == '--') cycle main_loop
                        n = 0
                        if (allocated(phi)) n = size(phi)
                        allocate (real_array(n + 1))
                        if (allocated(phi)) real_array(1:n) = phi(1:n)
                        call move_alloc(real_array, phi)
                        read (arg, *) phi(n + 1)
                        iarg = iarg + 1
                    end do
 
                case ('--energies')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    call read_real_array(arg, ek1)
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    call read_real_array(arg, ek2)
 
                case ('--ir-energy')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    read (arg, *) wir
 
                case ('--xuv-plus-ir', '--xuv+ir')
                    iarg = iarg + 1
                    call get_command_argument(iarg, filenames(1))
 
                case ('--xuv-minus-ir', '--xuv-ir')
                    iarg = iarg + 1
                    call get_command_argument(iarg, filenames(2))
 
                case ('--rawdips')
                    iarg = iarg + 1
                    call get_command_argument(iarg, filenames(1))
                    iarg = iarg + 1
                    call get_command_argument(iarg, filenames(2))
 
                case ('--out-dir')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    xuv_plus_ir_outdir % str = trim(arg)
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    xuv_minus_ir_outdir % str = trim(arg)
 
                case ('--properties')
                    iarg = iarg + 1
                    call get_command_argument(iarg, arg)
                    properties % str = trim(arg)
 
                case default
                    print '(2a)', 'Unknown command line argument ', trim(arg)
                    stop 1
 
            end select
 
            iarg = iarg + 1
 
        end do main_loop
 
        if (basis == -1) then
            print '(a)', 'Missing one of --cartesian/--spherical'
            stop 1
        end if
 
        sph = basis == 2
 
        ! read filenames from the XUV+IR and XUV-IR files
        do ifile = 1, 2
            open (newunit = u, file = trim(filenames(ifile)), action = 'read', iostat = stat)
            if (stat /= 0) then
                print '(a)', 'Failed to open file "', trim(filenames(ifile)), '"'
                stop 1
            end if
            if (ifile == 0) allocate (xuv_plus_ir_filenames(0))
            if (ifile == 1) allocate (xuv_minus_ir_filenames(0))
            n = 0
            do
                read (u, '(a)', iostat = stat) arg
                if (is_iostat_end(stat)) exit
                if (stat /= 0) then
                    print '(a,i0,3a)', 'Error ', stat, ' while reading from file "', filenames(ifile), '"'
                    stop 1
                end if
                allocate (string_array(n + 1))
                if (ifile == 1 .and. n >= 1) string_array(1:n) = xuv_plus_ir_filenames
                if (ifile == 2 .and. n >= 1) string_array(1:n) = xuv_minus_ir_filenames
                n = n + 1
                string_array(n) % str = trim(arg)
                if (ifile == 1) call move_alloc(string_array, xuv_plus_ir_filenames)
                if (ifile == 2) call move_alloc(string_array, xuv_minus_ir_filenames)
            end do
            close (u)
        end do
 
        print '(a,*(1x,i0,a))', 'Photon orders:', ord1, ' (extend to', ord1+abs(ext1), '),', ord2, ' (extend to', ord2+abs(ext2),')'
        print '(a,1x,a)',  'Angular basis:', merge('spherical', 'cartesian', sph)
        print '(a,1x,2(1x,f0.3))', 'XUV polarization:', polar
        print '(a,1x,2(1x,f0.3))', 'IR polarization:', polarir
        print '(a,5x,i0)', 'Max lab L:', maxlab
        if (allocated(theta)) print '(a,1x,*(1x,f0.1))', 'Polar angles:', theta
        if (allocated(phi)) print '(a,1x,*(1x,f0.1))', 'Azimut angles:', phi
        print '()'
 
    end function process_command_line
 
 
    !> \brief   Read multipoles from disk
    !> \author  J Benda
    !> \date    2021 - 2022
    !>
    !> Read all raw multipole files specified on the command line and store the data in `M_XUV_plus_IR` and `M_XUV_minus_IR`.
    !>
    subroutine read_raw_multipoles (order1, order2, sph, ntarg, maxl, nesc, target_states, chains1, chains2, &
                                    XUV_plus_IR_filenames, XUV_minus_IR_filenames, M_XUV_plus_IR, M_XUV_minus_IR)
 
        type(string), allocatable, intent(in) :: XUV_plus_IR_filenames(:), XUV_minus_IR_filenames(:)
 
        complex(real64), allocatable, intent(inout) :: M_XUV_plus_IR(:, :, :, :), M_XUV_minus_IR(:, :, :, :)
        integer,         allocatable, intent(inout) :: target_states(:), chains1(:, :), chains2(:, :)
 
        logical, intent(in)  :: sph
        integer, intent(in)  :: order1, order2
        integer, intent(out) :: ntarg, maxl, nesc
 
        character(len=:),allocatable :: filename
        complex(real64), allocatable :: buffer(:)
        integer,         allocatable :: new_target_states(:)
        integer                      :: i, idx, ichain, nchain1, nchain2, l, m, nene, ifile, nplus, nminus
 
        ntarg = 0
        maxl = 0
        nesc = 0
 
        allocate (target_states(ntarg))
 
        nplus = size(xuv_plus_ir_filenames)
        nminus = size(xuv_minus_ir_filenames)
 
        call setup_chains(order1, nchain1, chains1)
        call setup_chains(order2, nchain2, chains2)
 
        if (nplus /= nminus) then
            print '(a,i0,a,i0,a,/)', 'Warning: Number of XUV-IR and XUV+IR files is different (', nminus, ' vs ', nplus, ').'
        end if
 
        ! first pass over all files to get limiting values
        do ifile = 1, nplus + nminus
            if (ifile <= nplus) filename = xuv_plus_ir_filenames(ifile) % str
            if (ifile > nplus) filename = xuv_minus_ir_filenames(ifile - nplus) % str
            print '(3A)', 'Reading file "', filename, '"...'
            call read_file_data(filename, nene)
            if (ifile <= nplus) call parse_file_name(filename, order1, sph, ichain, i, l, m)
            if (ifile > nplus) call parse_file_name(filename, order2, sph, ichain, i, l, m)
            print '(2(A,I0),A,SP,I0)', '  - target state: ', i, ', partial wave: ', l, ', ', m
            print '(A,I0)', '  - number of energies: ', nene
            if (ifile <= nplus) print '(A,*(1x,I0))', '  - absorption chain:', chains1(:, ichain)
            if (ifile > nplus) print '(A,*(1x,I0))', '  - absorption chain:', chains2(:, ichain)
            maxl = max(maxl, l)
            nesc = max(nesc, nene)
            if (findloc(target_states, i, 1) == 0) then
                allocate (new_target_states(ntarg + 1))
                new_target_states(1 : ntarg) = target_states(1:ntarg)
                new_target_states(1 + ntarg) = i
                call move_alloc(new_target_states, target_states)
                ntarg = ntarg + 1
            end if
            deallocate (filename)
        end do
 
        ! resize the multipole storage
        allocate (m_xuv_plus_ir(nesc, (maxl + 1)**2, nchain1, ntarg), m_xuv_minus_ir(nesc, (maxl + 1)**2, nchain2, ntarg))
 
        m_xuv_plus_ir = 0
        m_xuv_minus_ir = 0
 
        ! actually read the dipoles now
        do ifile = 1, nplus + nminus
            if (ifile <= nplus) then
                filename = xuv_plus_ir_filenames(ifile) % str
                call parse_file_name(filename, order1, sph, ichain, i, l, m)
            else
                filename = xuv_minus_ir_filenames(ifile - nplus) % str
                call parse_file_name(filename, order2, sph, ichain, i, l, m)
            end if
            call read_file_data(filename, nene, buffer)
            idx = findloc(target_states, i, 1)
            if (ifile <= nplus) m_xuv_plus_ir(1:nene, l*l + l + m + 1, ichain, idx) = buffer(1:nene)
            if (ifile > nplus) m_xuv_minus_ir(1:nene, l*l + l + m + 1, ichain, idx) = buffer(1:nene)
            deallocate (filename)
        end do
 
    end subroutine read_raw_multipoles
 
 
    !> \brief   Assemble absorption chains
    !> \author  J Benda
    !> \date    2021 - 2023
    !>
    !> Construct database of all possible dipole component absorption chains. This is a two-dimensional
    !> array, where the second index numbers the individual pathway, while the first index numbers the
    !> absorbed photons. Polarizations of the photons are denoted by labels -1, 0, +1. The pathways are
    !> organized in lexicographical order. For example, the two-photon pathways look like this:
    !>
    !> \verbatim
    !>    1:  -1, -1
    !>    2:  -1,  0
    !>    3:  -1, +1
    !>    4:   0, -1
    !>    5:   0,  0
    !>    6:   0, +1
    !>    7:  +1, -1
    !>    8:  +1,  0
    !>    9:  +1, +1
    !> \endverbatim
    !>
    subroutine setup_chains (order, nchain, chains)
 
        integer,              intent(in)    :: order
        integer,              intent(out)   :: nchain
        integer, allocatable, intent(inout) :: chains(:, :)
 
        integer :: i, idx
 
        nchain = 3**order
 
        allocate (chains(order, nchain))
 
        ! the first pathway uses the lowest labels for all photons
        chains(:, 1) = -1
 
        ! build other pathways one after another
        do idx = 2, 3**order
 
            ! use the previous pathway as starting point
            chains(:, idx) = chains(:, idx - 1)
 
            ! atempt to increment the sequence (least significant digit last)
            increment_chain: do i = order, 1, -1
                if (chains(i, idx) < 1) then
                    ! if possible, increment the right-most value and exit
                    chains(i, idx) = chains(i, idx) + 1
                    exit increment_chain
                else
                    ! otherwise wrap around to the lowest index and try next
                    chains(i, idx) = -1
                end if
            end do increment_chain
 
        end do
 
    end subroutine setup_chains
 
 
    !> \brief   Extract quantum numbers from raw multipole file name
    !> \author  J Benda
    !> \date    2021
    !>
    !> Expect the multipole to be named "XXX-ab-(u,v,w)YYY.txt" and extract the characters 'a' and 'b' as well as the numbers
    !> 'u', 'v' and 'w'.
    !>
    subroutine parse_file_name (filename, order, sph, ichain, i, l, m)
 
        character(len=*), intent(in)  :: filename
        logical,          intent(in)  :: sph
        integer,          intent(in)  :: order
        integer,          intent(out) :: ichain, i, l, m
 
        character(len=1) :: c(order)
        integer          :: j, k, q(order), u, v
 
        k = scan(filename, '(', back = .true.)
 
        v = scan(filename(1:k-1), '-', back = .true.)
        u = scan(filename(1:v-1), '-', back = .true.)
 
        if (v - u /= order + 1) then
            print '(3A,I0)', 'Name of the file "', filename, '" is not consistent with the given order ', order
            stop 1
        end if
 
        do j = 1, order
            c(j) = filename(k - 2 - order + j : k - 2 - order + j)
        end do
 
        do j = 1, order
            if (sph) then
                select case (c(j))
                    case ('m'); q(j) = -1
                    case ('0'); q(j) =  0
                    case ('p'); q(j) = +1
                    case default
                        print '(3A)', 'Dipole component "', c(j), '" is not valid in spherical basis.'
                        stop 1
                end select
            else
                select case (c(j))
                    case ('y'); q(j) = -1
                    case ('z'); q(j) =  0
                    case ('x'); q(j) = +1
                    case default
                        print '(3A)', 'Dipole component "', c(j), '" is not valid in Cartesian basis.'
                        stop 1
                end select
            end if
        end do
 
        ichain = 1
        do j = 1, order
            ichain = 3*(ichain - 1) + q(j) + 2
        end do
 
        j = k + 1
        k = j + scan(filename(j:), ',') - 1
        read (filename(j : k - 1), *) i
 
        j = k + 1
        k = j + scan(filename(j:), ',') - 1
        read (filename(j : k - 1), *) l
 
        j = k + 1
        k = j + scan(filename(j:), ')') - 1
        read (filename(j : k - 1), *) m
 
    end subroutine parse_file_name
 
 
    !> \brief   Read raw multipole from file
    !> \author  J Benda
    !> \date    2021
    !>
    !> Read real and imaginary part of the raw multipole for all energies from the given file. The number of read elements
    !> is returned via the parameter "n", the values themselves via "buffer". If "buffer" is not allocated or is too short,
    !> it will be reallocated by this subroutine.
    !>
    subroutine read_file_data (filename, n, buffer)
 
        character(len=*),             intent(in)              :: filename
        integer,                      intent(out)             :: n
        complex(real64), allocatable, intent(inout), optional :: buffer(:)
 
        integer      :: ierr, u, i
        real(real64) :: re_d, im_d
 
        open (newunit = u, file = filename, action = 'read', form = 'formatted', iostat = ierr)
 
        if (ierr /= 0) then
            print '(3A)', 'Failed to open file "', filename, '"'
            stop 1
        end if
 
        ! first pass: count valid lines only
        ierr = 0
        n = 0
        do while (ierr == 0)
            read (u, *, iostat = ierr) re_d, im_d
            if (ierr == 0) n = n + 1
        end do
 
        ! if no buffer was given, do not read the data
        if (.not. present(buffer)) then
            close (u)
            return
        end if
 
        ! resize the buffer
        if (allocated(buffer)) then
            if (size(buffer) < n) then
                deallocate (buffer)
            end if
        end if
        if (.not. allocated(buffer)) then
            allocate (buffer(n))
        end if
 
        rewind(u)
 
        ! actually read the data now
        do i = 1, n
            read (u, *) re_d, im_d
            buffer(i) = cmplx(re_d, im_d, real64)
        end do
 
        close (u)
 
    end subroutine read_file_data
 
 
    !> \brief   Calculate multi-photon matrix elements by asymptotic approximation
    !> \authors J Benda
    !> \date    2024
    !>
    !> Use the provided matrix elements of order `ord` to calculate matrix elements of order `ord`+`ext`. The results are
    !> returned in the same (reallocated) array. If `outdir` is given, write out the resulting matrix elements to that
    !> directory.
    !>
    subroutine extend_order (ord, MM, chains, ext, Ek, wIR, target_states, echl, prop, maxl, outdir)
 
        use coupling_obj_gbl, only: couplings_type
        use multidip_asy,     only: akkl
        use multidip_params,  only: pi
 
        complex(real64), allocatable, intent(inout) :: MM(:, :, :, :)
        integer,                      intent(in)    :: ord, ext, target_states(:)
        integer,                      intent(in)    :: maxl
        real(real64),                 intent(in)    :: wIR
        real(real64),    allocatable, intent(in)    :: Ek(:), echl(:), prop(:, :, :)
        type(string),                 intent(in)    :: outdir
        integer,         allocatable, intent(inout) :: chains(:, :)
 
        complex(real64), allocatable :: M0(:, :, :, :), Apws(:, :), Aion(:, :)
        integer,         allocatable :: chains0(:, :)
 
        type(couplings_type) :: cpl
        character(len=1)     :: cmptname(3) = ['y', 'z', 'x']
        character(len=1024)  :: filename, strchain
 
        integer      :: order, nesc, nchain, nchain0, ntarg, l, m, n, u, lambda, mu, ipw, jpw, itarg, jtarg, &
                        ichain, jchain, ie, iex, compt
        real(real64) :: gaunt, dEIR, Qion, Qpws, kappa, k
 
        nesc = size(mm, 1)
        ntarg = size(mm, 4)
        order = ord
        deir = merge(wir, -wir, ext > 0)
 
        call cpl % prec_cgaunt(maxl)
 
        allocate (apws(nesc, 0:maxl), aion(nesc, ntarg))
 
        apws = 0
        aion = 0
 
        do iex = 1, abs(ext)
 
            nchain0 = size(chains, 2)
            order = order + 1
 
            print '(/,a,i0,/)', 'Extending pathway to order ', order
 
            call move_alloc(mm, m0)
            call move_alloc(chains, chains0)
            call setup_chains(order, nchain, chains)
 
            allocate (mm(nesc, (maxl + 1)**2, nchain, ntarg))
 
            mm = 0
 
            do itarg = 1, ntarg
                do l = 0, maxl
                    print '(2x,2(a,i0))', 'Precomputing Akkl factors for target ', itarg, ', and l = ', l
                    ! precompute asymptotic factors
                    !$omp parallel do schedule(dynamic) private(kappa, k, lambda, jtarg)
                    do ie = 1, nesc
                        if (ek(ie) - echl(itarg) + deir*iex <= 0) cycle
                        k = sqrt(2*(ek(ie) - echl(itarg) + deir*iex))
                        do lambda = abs(l - 1), min(maxl, l + 1)
                            if (ek(ie) - echl(itarg) + deir*(iex - 1) <= 0) cycle
                            kappa = sqrt(2*(ek(ie) - echl(itarg) + deir*(iex - 1)))
                            apws(ie, lambda) = akkl(kappa, lambda, k, l, 1)
                        end do
                        do jtarg = 1, ntarg
                            if (ek(ie) - echl(jtarg) + deir*(iex - 1) <= 0) cycle
                            kappa = sqrt(2*(ek(ie) - echl(jtarg) + deir*(iex - 1)))
                            aion(ie, jtarg) = akkl(kappa, l, k, l, 0)
                        end do
                    end do
                    do m = -l, l
                        ipw = l*l + l + m + 1
                        ! continuum-continuum transition
                        do lambda = abs(l - 1), min(maxl, l + 1)
                            do mu = -lambda, lambda
                                jpw = lambda*lambda + lambda + mu + 1
                                do compt = -1, +1
                                    gaunt = cpl % rgaunt(l, lambda, 1, m, mu, compt)
                                    if (gaunt /= 0) then
                                        qpws = sqrt(4*pi/3) * gaunt
                                        do jchain = 1, nchain0
                                            ichain = 3*(jchain - 1) + compt + 2
                                            mm(:, ipw, ichain, itarg) = mm(:, ipw, ichain, itarg) &
                                                + qpws * apws(:, lambda) * m0(:, jpw, jchain, itarg)
                                        end do
                                    end if
                                end do
                            end do
                        end do
                        ! ion-ion transition
                        do jtarg = 1, ntarg
                            do compt = -1, +1
                                qion = prop(target_states(itarg), target_states(jtarg), compt + 2)
                                if (qion /= 0) then
                                    do jchain = 1, nchain0
                                        ichain = 3*(jchain - 1) + compt + 2
                                        mm(:, ipw, ichain, itarg) = mm(:, ipw, ichain, itarg) &
                                            + qion * aion(:, jtarg) * m0(:, ipw, jchain, jtarg)
                                    end do
                                end if
                            end do
                        end do
                    end do
                end do
            end do
 
        end do
 
        if (len_trim(outdir % str) /= 0) then
            print '(/,3a)', 'Writing extended partial-wave amplitudes to "', outdir % str, '"'
            do itarg = 1, ntarg
                do ichain = 1, size(mm, 3)
                    do l = 0, maxl
                        do m = -l, l
                            ipw = l*l + l + m + 1
                            if (all(mm(:, ipw, ichain, itarg) == 0)) cycle
                            write (strchain, '(*(a))') [(cmptname(mod(ichain - 1, 3**n)/3**(n - 1) + 1), n = 1, order)]
                            write (filename, '(3a,2(i0,a),sp,i0,a)') 'rawdip-', trim(strchain), '-(', itarg, ',', l, ',', m, ').txt'
                            open (newunit = u, file = outdir % str // '/' // filename, action = 'write', form = 'formatted')
                            write (u, '(2e25.15)') mm(:, ipw, ichain, itarg)
                            close (u)
                        end do
                    end do
                end do
            end do
            open (newunit = u, file = outdir % str // '/' // 'pe_energies.txt', action = 'write', form = 'formatted')
            write (u, '(e25.15)') [(ek(ie) + ext*wir, ie = 1, nesc)]
            close (u)
        end if
 
    end subroutine extend_order
 
 
    !> \brief   Write results to file
    !> \author  J Benda
    !> \date    2021 - 2024
    !>
    subroutine write_quadratic_dipoles (L, MM, prefix, suffix)
 
        complex(real64), allocatable, intent(in) :: MM(:, :)
        character(len=*),             intent(in) :: prefix, suffix
        integer,                      intent(in) :: L
 
        character(len=200) :: filename
        integer            :: u, ie
 
        write (filename, '(A,I0,3A)') trim(prefix), l, trim(suffix), '.txt'
 
        print '(/,3A)', 'Writing M[XUV+IR]* M[XUV-IR] to file "', trim(filename), '"'
 
        open (newunit = u, file = trim(filename), action = 'write', form = 'formatted')
 
        do ie = 1, size(mm, 2)
            write (u, '(2E25.15)') sum(mm(:, ie))
        end do
 
        close (u)
 
    end subroutine write_quadratic_dipoles
 
 
    !> \brief   Molecular- and recoil-frame RABITT delays
    !> \author  J Benda
    !> \date    2022 - 2024
    !>
    !> Calculate the interference term of the photoelectron ionization probability for the following configurations:
    !>   1. If the polarization axis in the molecular frame is given, evalutes the distribution for all emission
    !>      directions provided by the arguments `theta` and `phi` (for molecular-frame observables). These are
    !>      written to files "rabitt-signal-mol-?.txt", one file per a given angle, one line per energy.
    !>   2. If the polarization axis in the molecular frame is given, evalutes the emission angle-integrated RABITT
    !>      interference term. The emission angle input from the command line is not used. The results are
    !>      written to file "rabitt-signal-int-1.txt".
    !>   3. If the polarization axis in the molecular frame is given, evaluates the distribution integrated around
    !>      the axis, for values of the polar angle `theta`. The values of `phi` are ignored. This is useful when
    !>      partially resolved recoil frame, where the recoil axis coincides with the main quantization axis (z).
    !>      The outputs are written to files "rabitt-signal-axs-?.txt", one file per a given angle, one line per energy.
    !>   4. Finally, regardless of the specification of the polarization axis, evaluate the quantity integrated over
    !>      all possible linear polarization as well as over all azimuthal emission angles. This is useful for
    !>      recoil-frame quantities. The outputs are written to files "rabitt-signal-rec-?.txt", one file per a given angle,
    !>      one line per energy.
    !>
    !> Colinear polarization of XUV and IR fields is assumed.
    !>
    !> \param[in] maxl    Highest angular momentum
    !> \param[in] chains1 List of available absorption paths for first (conjugated) pathway
    !> \param[in] chains2 List of available absorption paths for second (non-conugated) pathway
    !> \param[in] ntarg   Number of targets to consider
    !> \param[in] nesc    Number of scattering energies
    !> \param[in] Mp      Raw ionization matrix elements for the XUV+IR pathway (Cartesian basis)
    !> \param[in] Mm      Raw ionization matrix elements for the XUV-IR pathway (Cartesian basis)
    !> \param[in] axis    Unit vector in direction of the XUV polarization
    !> \param[in] axisIR  Unit vector in direction of the IR polarization
    !> \param[in] theta   List of photoelectron emission polar angles (in degrees)
    !> \param[in] phi     List of photoelectron emission azimuthal angles (in degrees)
    !>
    subroutine calculate_oriented_dipoles (maxl, chains1, chains2, ntarg, nesc, Mp, Mm, axis, axisIR, theta, phi)
 
        use dipelm_defs,              only: pi
        use dipelm_special_functions, only: a_re_sp_harm, a_sp_harm
        use multidip_special,         only: cartesian_vector_component_product_average
 
        integer,                      intent(in) :: maxl, ntarg, nesc
        integer,         allocatable, intent(in) :: chains1(:, :), chains2(:, :)
        real(real64),                 intent(in) :: axis(3), axisIR(3)
        real(real64),    allocatable, intent(in) :: theta(:), phi(:)
        complex(real64), allocatable, intent(in) :: Mp(:, :, :, :), Mm(:, :, :, :)
 
        real(real64),    allocatable :: Xlm(:, :), Yl0(:, :), sins(:), coss(:)
        complex(real64), allocatable :: Qmol(:, :), Qint(:, :), Qaxs(:, :), Qrec(:, :), QM(:), Xvalues(:)
 
        integer            :: l1, l2, m1, m2, pw1, pw2, pw1a, pw2a, c1, c2, i1, i2, icomp, iang, nang, ntheta, nphi, itg, ncompt(3)
        real(real64)       :: cav, cai, cax, rtheta, rphi, zero = 0
 
        ntheta = 0
        nphi = 0
 
        if (allocated(theta)) ntheta = size(theta)
        if (allocated(phi)) nphi = size(phi)
 
        nang = max(ntheta, nphi)
 
        if (nang == 0) return
 
        allocate (xlm((maxl + 1)**2, nang), yl0((maxl + 1)**2, nang), qmol(nesc, nang), qaxs(nesc, nang), qrec(nesc, nang), &
                  qint(nesc, 1), sins(nang), coss(nang), qm(nesc))
 
        qmol = 0
        qint = 0
        qaxs = 0
        qrec = 0
 
        ! precalculate angular functions
        do iang = 1, nang
            rtheta = 0
            rphi = 0
 
            if (ntheta > 0) rtheta = theta(1 + mod(iang - 1, ntheta)) * pi / 180
            if (nphi > 0) rphi = phi(1 + mod(iang - 1, nphi)) * pi / 180
 
            call a_re_sp_harm(maxl, rtheta, rphi, xvalues)
            xlm(:, iang) = real(xvalues, real64)    ! X_{lm}(theta, phi)
            deallocate (xvalues)
 
            call a_sp_harm(maxl, rtheta, zero, xvalues)
            yl0(:, iang) = real(xvalues, real64)    ! Y_{lm}(theta, 0)
            deallocate (xvalues)
 
            sins(iang) = sin(rtheta)
            coss(iang) = cos(rtheta)
        end do
 
        ! evaluate the interference terms
        do itg = 1, ntarg
            do l1 = 0, maxl
                do m1 = -l1, l1
                    pw1 = l1*l1 + l1 + m1 + 1
                    pw1a = l1*l1 + l1 + abs(m1) + 1
                    do l2 = 0, maxl
                        do m2 = -l2, l2
                            pw2 = l2*l2 + l2 + m2 + 1
                            pw2a = l2*l2 + l2 + abs(m2) + 1
                            do c1 = 1, size(chains1, 2)
                                do c2 = 1, size(chains2, 2)
                                    cax = 0
                                    cai = 0
                                    ncompt = 0
 
                                    if (sum(axis**2) > 0 .and. sum(axisir**2) > 0) then
                                        ! product of oriented polarization unit vectors
                                        cax = 1
                                        do icomp = 1, size(chains1, 1)
                                            i1 = 1 + mod(chains1(icomp, c1) + 2, 3)  ! m projection -> Cartesian index
                                            ncompt(i1) = ncompt(i1) + 1
                                            cax = cax * merge(axis(i1), axisir(i1), icomp == 1)
                                        end do
                                        do icomp = 1, size(chains2, 1)
                                            i2 = 1 + mod(chains2(icomp, c2) + 2, 3)  ! m projection -> Cartesian index
                                            ncompt(i2) = ncompt(i2) + 1
                                            cax = cax * merge(axis(i2), axisir(i2), icomp == 1)
                                        end do
 
                                        ! azimuthal integral over the product of polarization vectors (assume axis == axisIR)
                                        if (ncompt(1) == 0 .and.  ncompt(2) == 0) cai = 2*pi     ! zz, zzzz
                                        if (ncompt(1) == 2 .neqv. ncompt(2) == 2) cai = pi       ! xx, yy, xxzz, yyzz
                                        if (ncompt(1) == 2 .and.  ncompt(2) == 2) cai = pi/4     ! xxyy
                                        if (ncompt(1) == 4 .or.   ncompt(2) == 4) cai = 3*pi/4   ! xxxx, yyyy
                                        if (any(mod(ncompt, 2) /= 0)) cai = 0
                                        cai = cai * axis(3)**ncompt(3) * hypot(axis(1), axis(2))**(ncompt(1) + ncompt(2))
                                    end if
 
                                    ! angular average of product of polarization unit vectors
                                    cav = cartesian_vector_component_product_average([chains1(:, c1), chains2(:, c2)])
 
                                    ! update the value of the RABITT interference term
                                    if (cax /= 0 .or. cai /= 0 .or. cav /= 0) then
                                        qm = conjg(mp(:, pw1, c1, itg)) * mm(:, pw2, c2, itg)
                                        if (all(qm == 0)) cycle
                                        !$omp parallel do
                                        do iang = 1, nang
                                            if (cax /= 0) then
                                                qmol(:, iang) = qmol(:, iang) + qm * cax * xlm(pw1, iang) * xlm(pw2, iang)
                                            end if
                                            if (cax /= 0 .and. l1 == l2 .and. m1 == m2 .and. iang == 1) then
                                                qint(:, 1) = qint(:, 1) + qm * cax
                                            end if
                                            if (m1 == m2 .and. cai /= 0) then
                                                qaxs(:, iang) = qaxs(:, iang) + qm * cai * yl0(pw1a, iang) * yl0(pw2a, iang)
                                            end if
                                            if (m1 == m2 .and. cav /= 0) then
                                                qrec(:, iang) = qrec(:, iang) + qm * cav * yl0(pw1a, iang) * yl0(pw2a, iang)
                                            end if
                                        end do
                                    end if
                                end do  ! c2
                            end do  ! c1
                        end do  ! m2
                    end do  ! l2
                end do  ! m1
            end do  ! l1
        end do  ! itg
 
        ! write fixed-polarization results to files (one per emission direction)
        call write_to_file('rabitt-signal-mol-', qmol)
        call write_to_file('rabitt-signal-int-', qint)
        call write_to_file('rabitt-signal-axs-', qaxs)
 
        ! write fully polarization-averaged results to files (one per emission direction)
        call write_to_file('rabitt-signal-rec-', qrec)
 
    end subroutine calculate_oriented_dipoles
 
 
    !> \brief   Read a column from the text file
    !> \authors J Benda
    !> \date    2024
    !>
    subroutine read_real_array (filename, array)
 
        character(len=*), intent(in) :: filename
        real(real64), allocatable, intent(inout) :: array(:)
 
        integer :: u, i, stat, nlines
 
        nlines = 0
 
        open (newunit = u, file = filename, action = 'read', form = 'formatted')
 
        ! count lines first
        do
            read (u, *, iostat = stat)
            if (stat /= 0) exit
            nlines = nlines + 1
        end do
 
        allocate (array(nlines))
 
        rewind(u)
 
        ! now read the data
        do i = 1, nlines
            read (u, *) array(i)
        end do
 
        close (u)
 
    end subroutine read_real_array
 
 
    !> \brief   Write the interference term to file
    !> \author  J Benda
    !> \date    2022
    !>
    !> Produce one text file per angular combination. Write the real and imaginary part of the interference term, one line
    !> per scattering energy.
    !>
    !> \param[in]  stem  Initial string to be used in the name of the file to write
    !> \param[in]  Q     Interference terms (number of energies times number fo angles)
    !>
    subroutine write_to_file (stem, Q)
 
        character(len=*),             intent(in) :: stem
        complex(real64), allocatable, intent(in) :: Q(:, :)
 
        character(len=256) :: filename
        integer            :: iang, iene, u
 
        do iang = 1, size(q, 2)
            write (filename, '(A,I0,A)') stem, iang, '.txt'
            open (newunit = u, file = trim(filename), form = 'formatted')
            do iene = 1, size(q, 1)
                write (u, '(2E25.15)') q(iene, iang)
            end do
            close (u)
        end do
 
    end subroutine write_to_file
 
end program multidip_rabitt_delays