function source = FldT_source_coherence(data, freq, fd, fmin, fmax, fmid, ref_chan, hdm_file, range, step, trial_ind),

% FieldTrip function to generate source-EMG coherence estimate
%
%   Input parameters
%       data = output structure from FldT_read_MEGEMG.m
%       freq = output structure from FldT_MEGEMG_coherence.m
%       fd = output structure from FldT_MEGEMG_coherence.m
%       fmin/fmax = minimum/maximum frequency for source estimate
%       fmid = frequency of coherence calculation
%       ref_chan = reference channel name (string)
%       hdm_file = head model (as generated by localSpheres) filename (string)
%       range = [3 x 2] array of x/y/z minimum/maximum for ROI
%       step = grid spacing
%	trial_ind = trials to be included in the CMC calculation
%
%   Output parameters
%       source = output structure includes coherence values on ROI
%

    cfg = [];
    cfg.method = 'fft';
    cfg.output = 'powandcsd';
    cfg.tapsmofrq = 5;
    if exist('trial_ind', 'var'),
    	cfg.trials = trial_ind;
    else,
    	cfg.trials = 'all';
    end,
    cfg.foilim = [fmin fmax]; %cfg.foilim = [2 30];
    cfg.keeptrials = 'yes';
    cfg.sgncmb = {'MEG' 'MEG';'MEG' ref_chan};
    FDIfreq = freqanalysis(cfg,data);

    r=step;
    cfg = [];
    cfg.xgrid = range(1,1):r:range(1,2);
    cfg.ygrid = range(2,1):r:range(2,2);
    cfg.zgrid = range(3,1):r:range(3,2);
    cfg.hdmfile = hdm_file;
    [grid] = prepare_leadfield(cfg, FDIfreq);

    cfg = [];
    cfg.method = 'dics';
    cfg.refchannel = ref_chan;
    cfg.projectnoise ='yes';
    cfg.hdmfile = hdm_file;
    cfg.grid = grid;
    cfg.frequency = fmid;
    cfg.keepleadfield = 'no';
    cfg.feedback ='none';
    [source]= sourceanalysis(cfg,FDIfreq);

    [source]= sourcedescriptives([],source);

return,