function M = green3d(r1,z1,t1,r2,z2,t2)
% GREEN3D  Mutual inductance of 3D loops
%  GREEN3D(R1,Z1,PHI1,R2,Z2,PHI2) returns the MxN mutual inductances of the 3D
%  loops (R1(1:M,:),Z1,PHI1) and (R2(1:N,:),Z2,PHI2) where rows are paded with
%  NaNs, using mu0/(4*pi)*int(dl1.dl2/|x2-x1|).
%  If PHI1 (or PHI2) is empty, the loops are taken as circles at (R1,Z1) and
%  int(dphi2/(2*pi)*M(r1,z1,r2(phi2),z2(phi2))) is used.
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
 
 if isempty(t1) && isempty(t2)
  M = greenem('mut',r1,z1,r2,z2);
  
 elseif isempty(t1)
  M = NaN(numel(r1),size(r2,2));
  r0 = r1(:);
  z0 = z1(:);
  parfor l = 1:size(M,2)
   M(:,l) = green2d5(r0,z0,r2(:,l),z2(:,l),t2(:,l));
  end
  
 elseif isempty(t2)
  M = NaN(size(r1,2),numel(r2));
  r0 = r2(:);
  z0 = z2(:);
  parfor k = 1:size(M,1)
   M(k,:) = green2d5(r0,z0,r1(:,k),z1(:,k),t1(:,k));
  end
  
 else

  assert(~verLessThan('MATLAB','9.1'),'green3d with only 3d loops needs MATLAB 9.1 or later');
  
  x1 = r1 .* cos(t1);
  y1 = r1 .* sin(t1);
  x2 = r2 .* cos(t2);
  y2 = r2 .* sin(t2);
  M = NaN(size(r1,2),size(r2,2));
  for k = 1:size(M,1)
   k1 = ~isnan(x1(:,k));
   x0 = x1(k1,k);
   y0 = y1(k1,k);
   z0 = z1(k1,k);
   parfor l = 1:size(M,2)
    M(k,l) = greenxyz(x0,y0,z0,x2(:,l),y2(:,l),z2(:,l));
   end
  end
 end
 
end

function M = greenxyz(x1,y1,z1,x2,y2,z2)
 dx1 = diff(x1);
 dy1 = diff(y1);
 dz1 = diff(z1);
 dx2 = diff(x2');
 dy2 = diff(y2');
 dz2 = diff(z2');
 M = sum(sum((dx1*dx2 + dy1*dy2 + dz1*dz2)./ ...
  sqrt((x1(1:end-1)-x2(1:end-1)').^2 + ...
  (y1(1:end-1)-y2(1:end-1)').^2 + ...
  (z1(1:end-1)-z2(1:end-1)').^2),'omitnan'),'omitnan')*1e-7;
end

function M = green2d5(r1,z1,r2,z2,t2)
 l = find(~isnan(r2),1,'last');
 dt2 = [0; diff(t2(1:l))];
 dt2 = 0.5 * (dt2 + [dt2(2:end); 0]);
 M = greenem('mut',r1,z1,r2(1:l),z2(1:l))*dt2/(2*pi);
end