function varargout = greenem(mode,r1,z1,r2,z2,T)

%GREENEM  Electromagnetic Green's functions
% [Y1,...] = GREENEM({G1,...},R1,Z1[,R2,Z2]) computes according to the G
% selectors:
% 'm[ut]' the mutual inductance between the turn set R1,Z1 and R2,Z2 (R1,Z1
% if omitted)
% 'br','bz' the field at R1,Z1 producted by unit currents at R2,Z2
% 'dr1m','dr1br','dr1bz','dz1..','dr2..','dz2..' the partial derivates of
%  the mutual inductance or magnetic field
% 'dr1r1m','dr1r2m',... the second derivatives of the mutual inductance
% GREENEM({G1,...},R1,Z1,R2,Z2,T) post multiplies the output with T.
% GREENEM('m[ut]',R,Z,A|[A B]) computes the mutual between loops at R,Z
%  and their self with a circular cross section of radius A or a
%  rectangular cross section AxB
% GREENEM('s[elf]',R,~,A|[A B][,B]) returns only the self
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if ischar(mode), mode = {mode}; end
m = zeros(1,length(mode),'int32');
for k = 1:length(m)
 switch mode{k}
  case {'self' 's'}              , m(k) = -1;
  case {'mut' 'm'}               , m(k) =  0;
  case {'dr1m' 'dr1mut' 'dmdr1'} , m(k) =  1;
  case {'dz1m' 'dz1mut' 'dmdz1'} , m(k) =  2;
  case {'dr2m' 'dr2mut' 'dmdr2'} , m(k) =  3;
  case {'dz2m' 'dz2mut' 'dmdz2'} , m(k) =  4;
  case  'br'                     , m(k) =  5;
  case  'bz'                     , m(k) =  6;
  case {'dr1r1m' 'dr1r1mut'}     , m(k) = 11;
  case {'dr1z1m' 'dr1z1mut'}     , m(k) = 12;
  case {'dr1r2m' 'dr1r2mut'}     , m(k) = 13;
  case {'dr1z2m' 'dr1z2mut'}     , m(k) = 14;
  case {'dz1r1m' 'dz1r1mut'}     , m(k) = 12; % 21
  case {'dz1z1m' 'dz1z1mut'}     , m(k) = 22;
  case {'dz1r2m' 'dz1r2mut'}     , m(k) = 23;
  case {'dz1z2m' 'dz1z2mut'}     , m(k) = 24;
  case {'dr2r1m' 'dr2r1mut'}     , m(k) = 13; % 31
  case {'dr2z1m' 'dr2z1mut'}     , m(k) = 23; % 32
  case {'dr2r2m' 'dr2r2mut'}     , m(k) = 33;
  case {'dr2z2m' 'dr2z2mut'}     , m(k) = 34;
  case {'dz2r1m' 'dz2r1mut'}     , m(k) = 14; % 41
  case {'dz2z1m' 'dz2z1mut'}     , m(k) = 24; % 42
  case {'dz2r2m' 'dz2r2mut'}     , m(k) = 34; % 43
  case {'dz2z2m' 'dz2z2mut'}     , m(k) = 22; % 44
  case {'dr1br' 'dbrdr' 'dbrdr1'}, m(k) = 51;
  case {'dz1br' 'dbrdz' 'dbrdz1'}, m(k) = 52;
  case {'dr2br' 'dbrdr2'}        , m(k) = 53;
  case {'dz2br' 'dbrdz2'}        , m(k) = 54;
  case {'dr1bz' 'dbzdr' 'dbzdr1'}, m(k) = 61;
  case {'dz1bz' 'dbzdz' 'dbzdz1'}, m(k) = 62;
  case {'dr2bz' 'dbzdr2'}        , m(k) = 63;
  case {'dz2bz' 'dbzdz2'}        , m(k) = 64;
  case {'0' 'zero'}              , m(k) = 99;
  otherwise,	error('Invalid G specifier: %s',mode{k})
 end
end

switch nargin
 case 3
  r2 = r1; z2 = z1; 
 case 4
  c = r2;
  r2 = r1; z2 = z1;
  nr = size(r1); nc = size(c);
  if isequal(nr,nc)
   a = c; b = [];
  elseif nr(1) == 1 && nc(1) == 2 && nr(2) == nc(2)
   a = c(1,:); b = c(2,:);
  elseif nr(2) == 1 && nc(2) == 2 && nr(1) == nc(1)
   a = c(:,1); b = c(:,2);
  else
   error('R,Z A|[A B] size mismatch')
  end
 case 5
  if m == -1
   a = r1; b = z1;
  end
 case 6
 otherwise
  error('Invalid arguments')
end

if any(m(1) >= 0)
 if ~isequal(size(r1),size(z1)) || ~isequal(size(r2),size(z2)), error('R1,Z1 R2,Z2 size mismatch'), end
 varargout = cell(1,length(m));
 % greenemmex assumes doubles as inputs
 [varargout{:}] = greenemmex(double(r1),double(z1),...
                             double(r2),double(z2),...
                             m);
end

k0 = find(m ==  0);
k1 = find(m == -1);
if ~isempty(k1) || ~isempty(k0) && nargin == 4
 if isempty(b)
  s = 4e-7*pi * r1 .* (log(8*r1./a    ) - 1.75);
 else
  s = 4e-7*pi * r1 .* (log(8*r1./(a+b)) - 0.50);
 end
 if k0
  n = length(r1);
  varargout{k0}(sub2ind([n n],1:n,1:n)) = s;
 else
  varargout{k1} = s;
 end
end

if nargin == 6
 for k = 1:nargout
  varargout{k} = varargout{k} * T;
 end
end
