dirichlet_neumann_difference

 TODO:
  - consider to implement Harach https://iopscience.iop.org/article/10.1088/0266-5611/31/11/115008/pdf
  - maybe other approaches 

 (C) 2023 Joeran Rixen. Licensed under GPL version 2 or 3

0001 function del_L = dirichlet_neumann_difference(fmdl,data1,data2);
0002 % TODO:
0003 %  - consider to implement Harach https://iopscience.iop.org/article/10.1088/0266-5611/31/11/115008/pdf
0004 %  - maybe other approaches
0005 %
0006 % (C) 2023 Joeran Rixen. Licensed under GPL version 2 or 3
0007 
0008 % TODO: allow choosing the approach, only have one for now
0009 del_L =  d2n_from_full_data(fmdl,data1,data2);
0010 
0011 end
0012 
0013 function del_L = d2n_from_full_data(fmdl,data1,data2);
0014     %% Calculations regarding the Dirichlet to Neumann map
0015     stim_pattern_mat = create_stim_pattern_matrix(fmdl);
0016     stim_pattern_mat = stim_pattern_mat./sqrt(sum(stim_pattern_mat.^2,1));
0017 
0018     n_elec = num_elecs(fmdl);
0019 
0020     data1 = calc_difference_data(0,data1,fmdl);
0021     data2 = calc_difference_data(0,data2,fmdl);
0022 
0023     % reshaping voltages for easy BN matrix assembly
0024     data1 = reshape(data1, [n_elec, n_elec-1]);
0025     data2 = reshape(data2, [n_elec, n_elec-1]);
0026     
0027     data1 = make_voltage_mean_free(data1);
0028     data2 = make_voltage_mean_free(data2);
0029     
0030     data1 = data1./sqrt(sum(stim_pattern_mat.^2,1));
0031     data2 = data2./sqrt(sum(stim_pattern_mat.^2,1));
0032     
0033     %either calcualte or give, as this is rather imporant for scaling of
0034     %the output of the d-bar algorithm.
0035     elec_area = pi*0.25^2;
0036 
0037     NtoD1 = (1/elec_area)*(full(stim_pattern_mat)'*data1);
0038     NtoD2 = (1/elec_area)*(full(stim_pattern_mat)'*data2);
0039     
0040     DN1 = inv(NtoD1);
0041     DN2 = inv(NtoD2);
0042 
0043     del_L = (DN1 - DN2);
0044 end
0045     
0046 
0047 function stim_pattern_mat = create_stim_pattern_matrix(fmdl)
0048     % converts the stimpattern in the model into matrix form for ease of
0049     % computation
0050     n_elec = length(fmdl.electrode);
0051     stim_pattern_mat = zeros(n_elec, n_elec-1);
0052     
0053     for i = 1:n_elec-1
0054         stim_pattern_mat(:, i) = full(fmdl.stimulation(i).stim_pattern);
0055     end
0056 end
0057 
0058 function data_mat_free = make_voltage_mean_free(data_mat)
0059     % makes each stimulation voltage vector mean free
0060     U_mean = mean(data_mat, 1);
0061     data_mat_free = data_mat - U_mean;
0062 end
0063