function [tau, g1] = mcxdcsg1(detps, tau, disp_model, DV, lambda, format, varargin)
%
%   [tau,g1]=mcxdcsg1(detps,tau, disp_model, DV, lambda, format)
%
%   Compute simulated electric-field auto-correlation function using
%   simulated photon pathlengths and scattering momentum transfer
%
%   author: Stefan Carp (carp <at> nmr.mgh.harvard.edu)
%
%   input:
%       detps:      the file name of the output .mch file or the 2nd output from mcxlab
%       tau:        correlation times at which to compute g1
%                   (default: 1e-7 to 1e-1 seconds, log equidistant)
%       disp_model: displacement model ('brownian', 'random_flow', <custom>)
%                   (default: brownian, see further explanation below)
%       disp_var:   value of displacement variable using mm as unit of
%                   length and s as unit of time
%                   (default: 1e-7 mm^2/s, see further explanation below)
%       lambda:     wavelenght of light used in nm
%                   (default: 785)
%       format:     the format used to save the .mch file
%                   (default: 'float')
%
%   output:
%
%       tau:        correlation times at which g1 was computed provided for
%                   convenience (copied from input if set, otherwise
%                   outputs default)
%       g1:         field auto-correlation curves, one for each detector
%
%   The displacement model indicates the formula used to compute the root
%   mean square displacement of scattering particles during a given delay
%
%   brownian:       RMS= 6 * DV * tau;
%                   DV(displacement variable)=Db (brownian diffusion coeff)
%   random_flow:    RMS= DV^2 * tau^2;
%                   DV = V (first moment of velocity distribution)
%   <custom>:       any string other than 'brownian' or 'random_flow' will
%                   be evaluate as is using Matlab evalf, make sure it uses
%                   'DV' as the flow related independent variable, tau is
%                   indexed as tau(J). Any additional parameters can be
%                   sent via "varargin"
%
%   This file is part of Mesh-Based Monte Carlo
%   License: GPLv3, see http://mcx.space for details
%

if nargin < 6
    format = 'float';
end
if nargin < 5
    lambda = 785;
end
if nargin < 4
    DV = 1e-7;
end
if nargin < 3
    disp_model = 'brownian';
end
if nargin < 2
    tau = logspace(-7, -1, 200);
end

if (ischar(detps))
    [mch_data, mch_header] = loadmch(detps, format);
    [fpath, fname, fext] = fileparts(detps);

    cfg = loadjson([fpath filesep fname '.json']);
    prop = cell2mat(cfg.Domain.Media);
    mua = cell2mat({prop.mua});
    mua = mua(2:end);
    n = cell2mat({prop.n});
    n = n(2:end);
    medianum = length(mua) - 1;
    detps = struct('detid', mch_data(1, :)');
    detps.ppath = mch_data(3:2 + medianum, :)';
    if (size(mch_data, 1) >= 2 * medianum + 2)
        detps.mom = mch_data(medianum + 3:2 * medianum + 2, :)';
    end
else
    mua = detps.prop(2:end, 1)';
    n = detps.prop(2:end, 4)';
end

if (~isfield(detps, 'mom'))
    error('No momentum transfer data are found, please rerun your simulation and set cfg.ismomentum=1.');
end

if strcmp(disp_model, 'brownian')
    disp_str = 'rmsdisp=6*DV.*tau(J);';
elseif strcmp(disp_model, 'random_flow')
    disp_str = 'rmsdisp=DV.^2.*tau(J).^2;';
else
    disp_str = ['rmsdisp=' disp_model ';'];
end

k0 = 2 * pi * n / (lambda * 1e-6);

detlist = sort(unique(detps.detid));
g1 = zeros(max(detlist), length(tau));

for detid = 1:length(detlist)
    I = detlist(detid);
    idx = find(detps.detid == I);
    fprintf('Processing detector %.0f: %.0f photons\n', I, length(idx));

    for J = 1:length(tau)
        eval(disp_str);
        g1(I, J) = sum(exp(-(k0.^2 .* rmsdisp / 3) * detps.mom' - mua * detps.ppath'));
    end
    g1_norm = sum(exp(-mua * detps.ppath'));
    g1(I, :) = g1(I, :) ./ g1_norm;
end