%
% [Sin,SNin] = Sin_vsh_ctf_fast(r_sphere,R,EX,EY,EZ,Lin)
%
% Calculate the internal SSS basis Sin for a CTF system
% using vector spherical harmonics
%
function [Sin,SNin] = Sin_vsh_ctf_fast(r_sphere,R,EX,EY,EZ,Lin)

% Copyright (c) 2016, Elekta Oy
% ---------------------------------------
% 
% Redistribution and use of the Software in source and binary forms, with or without 
% modification, are permitted for non-commercial use.
% 
% The Software is provided "as is" without warranties of any kind, either express or
% implied including, without limitation, warranties that the Software is free of defects,
% merchantable, fit for a particular purpose. Developer/user agrees to bear the entire risk 
% in connection with its use and distribution of any and all parts of the Software under this license.
% 

mu0 = 1.25664e-6; % Permeability of vacuum
%
% For numerical surface integration:
%
%baseline = 50e-3;
dx = 4.5e-3;
dy = 4.5e-3;
dz1 = 0;
dz2 = 50e-3;
D = [dx dy dz1; dx -dy dz1; -dx dy dz1; -dx -dy dz1; dx dy dz2; dx -dy dz2; -dx dy dz2; -dx -dy dz2]';
for j = 1:8
   if j <= 4
      weights(j) = 1/(4*1);
   else
      weights(j) = -1/(4*1);
   end
end
weights = weights';

for ch = 1:size(R,2)
   count = 1;
   R(:,ch) = R(:,ch) - r_sphere;
   Sin(ch,:) = -mu0*vsh_response(R(:,ch),EX(:,ch),EY(:,ch),EZ(:,ch),D,weights,Lin);
end
for j = 1:size(Sin,2)
   SNin(:,j) = Sin(:,j)/norm(Sin(:,j));
end


function Sin_elements = vsh_response(r,ex,ey,ez,D,weights,Lin)

for j = 1:length(weights)
    r_this = r + D(1,j)*ex + D(2,j)*ey;
    rn(j) = norm(r_this);
    theta(j) = acos(r_this(3)/rn(j));
    phi(j) = atan2(r_this(2),r_this(1));
    sint(j) = sin(theta(j));
    sinp(j) = sin(phi(j));
    cost(j) = cos(theta(j));
    cosp(j) = cos(phi(j));
    for l = 1:Lin
       p0{l}(:,j) = legendre(l,cos(theta(j)));
       rnv(j,l) = rn(j)^(l+2);
   end
end
Sin_elements = [];
for l = 1:Lin
  for m = -l:l
    for j = 1:length(weights)
      %vs = vsh_modified_in_fast(theta(j),phi(j),l,m,p0{l}(:,j))'/rn(j)^(l+2);
      vs = vsh_modified_in_fast(theta(j),phi(j),l,m,p0{l}(:,j))'/rnv(j,l);
      V(1,j) = vs(1)*sint(j)*cosp(j) + vs(2)*cost(j)*cosp(j) - vs(3)*sinp(j);
      V(2,j) = vs(1)*sint(j)*sinp(j) + vs(2)*cost(j)*sinp(j) + vs(3)*cosp(j);
      V(3,j) = vs(1)*cost(j) - vs(2)*sint(j);
    end
    Sin_elements = [Sin_elements dot(V*weights,ez)]; % Cartesian coordinates
                                                     %Sin_element = Sin_element/sqrt((l+1)*(2*l+1));  % Back to orthonormal presentation
  end 
end 
%Sin_element = dot(V*weights,ez); % Cartesian coordinates
%Sin_element = Sin_element/sqrt((l+1)*(2*l+1));  % Back to orthonormal presentation


