function [rfjac_lnA, rfjac_phase] = mcxrfreplay(cfg, f_mod, detp, seeds, detnums)
%
% [rfjac_lnA, rfjac_phase]=mcxrfreplay(cfg,f_mod,detp,seeds,detnums)
%
% Compute the frequency domain (FD) log-amplitude and phase shift Jacobians
% with respect to voxel-wise absorption coefficients using the radio
% frequency (RF) replay algorithm.
%
% Author: Pauliina Hirvi (pauliina.hirvi <at> aalto.fi)
%         Qianqian Fang (q.fang <at> neu.edu)
%
% Ref.: Hirvi et al. (2023). Effects of atlas-based anatomy on modelled
% light transport in the neonatal head. Phys. Med. Biol.
% https://doi.org/10.1088/1361-6560/acd48c
%
% Input:
%     cfg: struct used in the main forward simulation
%     f_mod: RF modulation frequency
%     detp: the 2nd output from mcxlab, must be a struct
%     seeds: the 4th output from mcxlab
%     detnums: row array with the indices of the detectors to replay and obtain Jacobians
%
% Output:
%     rfjac_lnA: a 4D array with dimensions specified by [size(vol) numb-of-detectors];
%                each 3D array contains the Jacobians for log-amplitude measurements
%     rfjac_phase: a 4D array with dimensions specified by [size(vol) numb-of-detectors];
%                  each 3D array contains the Jacobians for phase shift measurements
%
% License: GPLv3, see http://mcx.space/ for details
%

if (nargin < 5)
    error('you must provide all 5 required input parameters');
end

if (~isfield(cfg, 'unitinmm'))
    cfg.unitinmm = 1;
end

% Initialize the 4D arrays for collecting the Jacobians. The 4th dimension
% corresponds to the detector index.
rfjac_lnA = zeros([size(cfg.vol), length(detnums)]);
rfjac_phase = zeros([size(cfg.vol), length(detnums)]);

% Collect Jacobians one detector index at a time.
for d = detnums
    % MCXLAB REPLAY SETTINGS
    clear cfg_jac;
    cfg_jac = cfg;
    cfg_jac.seed = seeds.data;
    cfg_jac.detphotons = detp.data;
    cfg_jac.replaydet = d;
    cfg_jac.outputtype = 'rf';
    cfg_jac.omega = 2 * pi * f_mod; % 100 MHz RF modulation
    cfg_jac.isnormalized = 0; % !

    % REPLAY SIMULATION
    [rfjac_d, detp_d, vol_d, seeds_d] = mcxlab(cfg_jac);
    detw = mcxdetweight(detp_d, cfg_jac.prop, cfg_jac.unitinmm); % array with detected photon weights
    dett = mcxdettime(detp_d, cfg_jac.prop, cfg_jac.unitinmm); % array with photon time-of-flights

    % FD MEASUREMENT ESTIMATES
    X = dot(detw, cos((2 * f_mod) .* dett * pi));
    Y = dot(detw, sin((2 * f_mod) .* dett * pi));
    A = sqrt(X^2 + Y^2); % amplitude [a.u.]
    phase = atan2(Y, X) + (double(atan2(Y, X) < 0)) * 2 * pi; % phase shift in [0,2*pi] [rad]
    if A == 0
        fprintf('MCX WARNING: No detected photons for detector %d.\n', d);
        continue
    end

    % FD JACOBIANS
    % Compute the Jacobians with the rf replay feature.
    rfjac_d = sum(rfjac_d.data, 4); % sum over time instances
    if cfg_jac.isnormalized == 0
        rfjac_d = (cfg_jac.unitinmm) .* rfjac_d; % correct units to [mm]
    end
    % Jacobians for X and Y wrt mua:
    rfjac_X = rfjac_d(:, :, :, :, :, 1);
    rfjac_Y = rfjac_d(:, :, :, :, :, 2);
    % Jacobians for log-amplitude and phase shift wrt mua:
    rfjac_lnA(:, :, :, d) = (1 / (A^2)) .* (X .* rfjac_X + Y .* rfjac_Y);
    rfjac_phase(:, :, :, d) = (1 / (A^2)) .* (X .* rfjac_Y - Y .* rfjac_X);
end