function [Jx,Jy,Jz] = J_coeff(x,y,z,N_wire,I_mag,mu_o)
%Physically, this function assigns the current density to nodes inside each
%specific coil - specifically coils with a rectangular, hollow-cored cross
%section.
%As a physical model, each point in the coil has a constant magnitude of
%the current density, with associated vector components used to create a
%continous current loop around the coil in the XY plane.
%N_wire is number or turns, Imag is the current for that specific coil,
%mu_o is the standard variable but useful to = 1 when testing code.
%L_rat is used to calculate the inner length of the coil regardless of
%specific scale of the mesh.

[Corner,~,~] = cubecorner(x,y,z);%Find Coil Corners
Z_height = range(Corner(:,3));       %Coil Variables
L_outer = range(Corner(:,1));
L_rat = 1.46257197696737;            %Fixed for this coil
L_inner = L_outer ./ L_rat ;
CA = Z_height.*(L_outer-L_inner)/2;  %Cross Sectional Area of Coil
jo = -mu_o*N_wire*I_mag ./ CA;   %Current Density Magnitude

f_x = (abs(x) <= abs(y)).*sign(y);
f_y = (abs(y) <= abs(x)).*sign(x);
J(1,:) = jo.*f_x;   %Jx
J(2,:) = -jo*f_y;   %Jy
J(3,:) = 0;         %Jz
clearvars -except J

%Ideally points on the perimeter at the corners of the coil should be
%reassigneed so they have a single vector component directing flow around 
%the perimeter but I can't find a nice way to do this.


%For mag(J) to be constant, this scaling is required where both Jx and Jy
%components are defined:
if ~isempty((J(1,:) ~= 0 & J(2,:)~=0))
    doub_vect = (J(1,:) ~= 0 & J(2,:)~=0);
    J(1:2,doub_vect) = (1/sqrt(2)) .* J(1:2,doub_vect);
end

Jx = J(1,:); Jy = J(2,:); Jz = J(3,:);
%
% % %Quiver Plot of Current Assignment:
% figure; quiver3(x,y,z,J(1,:),J(2,:),J(3,:)); hold on;
% scatter3(x,y,z,'kx');
% xlabel('x'); ylabel('y'); zlabel('z');

end