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/>.
!
 
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
 
    subroutine print_help
 
        print '(A)', 'Calculate discrete RABITT time delays using the two-photon raw dipole elements'
        print '( )'
        print '(A)', 'Usage:'
        print '( )'
        print '(A)', '    multidip_rabitt_delays [OPTIONS] --xuv-plus-ir [FILES1] --xuv-minus-ir [FILES2]'
        print '( )'
        print '(A)', 'Options:'
        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)', '    --spherical             Assume matrix elements in the spherical basis'
        print '(A)', '    --cartesian             Assume matrix elements in the Cartesian basis'
        print '(A)', '    --xuv-plus-ir FILE      File with paths to files containing raw matrix elements for the XUV+IR pathway'
        print '(A)', '    --xuv-minus-ir FILE     File with paths to files containing raw matrix elements for the XUV-IR pathway'
        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 '(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 --xuv-plus-ir <(ls XUV+IR/rawdip-*) --xuv-minus-ir <(ls XUV-IR/rawdip-*)'
        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
 
 
    subroutine rabitt_main
 
        use iso_fortran_env,   only: output_unit
        use mpi_gbl,           only: mpi_mod_start, mpi_mod_finalize
        use multidip_routines, only: calculate_quadratic_dipole_sph, calculate_quadratic_dipole_xyz, convert_xyz_to_sph
 
        type(string), allocatable :: XUV_plus_IR_filenames(:), XUV_minus_IR_filenames(:)
 
        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(:)
        real(real64)              :: polar(2), polarIR(2), axis(3), axisIR(3)
 
        integer, allocatable :: target_states(:), chains1(:, :), chains2(:, :)
 
        integer :: i, L, order1, order2, maxl, ntarg, nesc, maxlab
        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'
 
        if (.not. process_command_line(order1, order2, sph, maxlab, polar, polarir, theta, phi, &
                                       xuv_plus_ir_filenames, xuv_minus_ir_filenames)) 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
            m_xuv_plus_ir_sph = m_xuv_plus_ir_xyz
            m_xuv_minus_ir_sph = m_xuv_minus_ir_xyz
        else
            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 + order2)
 
            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)
            call write_quadratic_dipoles(l, mm, str, 'sph')
 
            if (l == 0) then
               mm0 = mm
            else
               do i = 3, 2 + ntarg
                  mm(i,:) = mm(i,:) / mm0(i,:)
               enddo
 
               str = 'rabitt-beta-L'
               call write_quadratic_dipoles(l, mm, str, '')
            endif
 
        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
 
 
    logical function process_command_line (ord1, ord2, sph, maxlab, polar, polarIR, theta, phi, &
                                           XUV_plus_IR_filenames, XUV_minus_IR_filenames)
 
        integer,      intent(out) :: ord1, ord2, maxlab
        real(real64), intent(out) :: polar(2), polarir(2)
        logical,      intent(out) :: sph
 
        real(real64), allocatable, intent(inout) :: theta(:), phi(:)
        real(real64), allocatable                :: real_array(:)
 
        type(string), allocatable, intent(inout) :: xuv_plus_ir_filenames(:), xuv_minus_ir_filenames(:)
        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
        basis = -1
        maxlab = 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 ('--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(theta)) real_array(1:n) = phi(1:n)
                        call move_alloc(real_array, phi)
                        read (arg, *) phi(n + 1)
                        iarg = iarg + 1
                    end do
 
                case ('--xuv-plus-ir')
                    iarg = iarg + 1
                    call get_command_argument(iarg, filenames(1))
 
                case ('--xuv-minus-ir')
                    iarg = iarg + 1
                    call get_command_argument(iarg, filenames(2))
 
                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/--spehrical'
            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,2x,i0,2x,i0)', 'Photon orders:', ord1, ord2
        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
 
 
    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
 
 
    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
 
 
    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
 
 
    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
 
 
    subroutine write_quadratic_dipoles (L, MM, prefix, stem)
 
        complex(real64), allocatable, intent(in) :: MM(:, :)
        character(len=*),             intent(in) :: prefix
        character(len=*),             intent(in) :: stem
        integer,                      intent(in) :: L
 
        character(len=200) :: filename
        integer            :: u, ie
 
        write (filename, '(A,I0,3A)') trim(prefix), l, '-', trim(stem), '.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
 
 
    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
 
 
    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