update mcxrfreplay, update mcxlab jupyter notebook, add PKG_ADD

mcxlab/PKG_ADD

@@ -0,0 +1,5 @@
+if(exist(file_in_loadpath('mcxcl.mex')))
+    autoload('mcxcl',file_in_loadpath('mcxcl.mex'))
+else
+    autoload('mcxcl',file_in_loadpath(['octave' filesep regexprep(computer('arch'), 'darwin[0-9.]+-', 'darwin-') filesep 'mcxcl.mex']))
+end

mcxlab/mcxlab.m

@@ -224,7 +224,7 @@
 %                   light (-1 <= Q <= 1)
 %      cfg.lambda: source light wavelength (nm) for polarized MC
 %      cfg.issrcfrom0: 1-first voxel is [0 0 0], [0]- first voxel is [1 1 1]
-%      cfg.replaydet:  only works when cfg.outputtype is 'jacobian', 'wl', 'nscat', or 'wp' and cfg.seed is an array
+%      cfg.replaydet:  only works when cfg.outputtype is 'jacobian', 'wl', 'nscat', 'wp' or 'rf' and cfg.seed is an array
 %                      -1 replay all detectors and save in separate volumes (output has 5 dimensions)
 %                       0 replay all detectors and sum all Jacobians into one volume
 %                       a positive number: the index of the detector to replay and obtain Jacobians

utils/mcxrfreplay.m

@@ -1,8 +1,7 @@
-
 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.
@@ -20,7 +19,7 @@
 %     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
@@ -30,6 +29,14 @@
 % 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)]);
@@ -40,9 +47,9 @@
     % 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.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; % !
@@ -53,8 +60,8 @@
     dett=mcxdettime(detp_d,cfg_jac.prop,cfg_jac.unitinmm); % array with photon time-of-flights
 
     % FD MEASUREMENT ESTIMATES
-    X=dot(detw,cospi((2*f_mod).*dett));
-    Y=dot(detw,sinpi((2*f_mod).*dett));
+    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