params_C.m

%------------------------------------------------------------------------%  
%   [K_1, Lc_1, Rc_1] = params_C(J_sol, L_sol, R1_sol, R2_sol)
%   solves for the cylinder parameters for 1 cylinder using the equations
%   defined by case C
%   case C equations [b1; b3; Bz];
%
%   INPUTS
%   J_sol: J value for solenoid
%   L_sol: length of solenoid
%   R_1: inner radius of solenoid
%   R_2: outer radius of solenoid
%
%   OUTPUTS
%   K_1: current for cylinder 1
%   Lc_1: length of cylinder 1
%   Rc_1: radius of cylinder 1
%
%   REMARKS
%   case C = [b1; b3; Bz];
%   The optimal paramters of solved by using an LM algorithum
%
%   AUTHOR(S): Paige Husa
%
%                  v1.0 4/29/2019
% ----------------------------------------------------------------------- %

function [K_1, Lc_1, Rc_1, soln_found] = params_C(J_sol, L_sol, R1_sol, R2_sol)
format long;
%% solenoid parameters
realS_J = J_sol; 
realS_L = L_sol; 
realS_R1 = R1_sol;
realS_R2 = R2_sol; 

delta_R = realS_R2 - realS_R1;
real_z = 0; 
real_mu = 1;%(4*pi)*1e-7;

%% solenoid side of equation
c = c_vector_C(realS_J, realS_L, realS_R1, realS_R2, real_mu, real_z);

%% inputs to LM function
ravg = (realS_R2+realS_R1)/2;
rdiff = (realS_R2-realS_R1)/2;

%initial guess the optimal paramters
cur_scale =  c(1)/(pi*realS_L*(ravg)^2) ;
gammaGuess = [1; realS_L; ravg];
% parameters K_1,Lc_1,Rc_1

tol = 1e-5;

%is basically going to run forever unless you stop the algorithum or
%solution is found
maxIterations = 500;
maxTime = 60*5;
dispVals = 0;
maxChange = [];

%% LM method
[ currentEstimate, cost, soln_found ] = LM_LeastSquares( gammaGuess, @(x)costF(x,c,cur_scale), tol, maxIterations, maxTime, dispVals, maxChange );

if( any(currentEstimate <0 ))
    warning('Solution Found is non-physical. Resolving...');
    [ currentEstimate, cost, soln_found ] = LM_LeastSquares( abs(currentEstimate), @(x)costF(x,c,cur_scale), tol, maxIterations, maxTime, dispVals, maxChange );
    currentEstimate = abs(currentEstimate);
  
end
% currentEstimate
% cost
soln_found = cost < 1e-3;

%% final optimal values
K_1 = currentEstimate(1)*cur_scale;
Lc_1 = currentEstimate(2);
Rc_1 = currentEstimate(3);
end

%% cost function
function [error, Jacobian] = costF(gammaGuess,c,cur_scale)
real_z = 0; 
real_mu = 1;%(4*pi)*1e-7;

f =              f_vector_C(gammaGuess(1)*cur_scale, gammaGuess(2), gammaGuess(3), real_mu, real_z);

J_gamma = jacobian_matrix_C(gammaGuess(1)*cur_scale, gammaGuess(2), gammaGuess(3), real_mu, real_z);

J_gamma(:,1) = J_gamma(:,1)*cur_scale;
Jacobian = J_gamma;

% error  = Jacobian*gammaGuess;
error = f -c;

error = error/norm(c);
Jacobian = Jacobian/norm(c);

%scale the jabocian appropriately so don't get a singularity
for i = 1:length(error)
   sc = norm(Jacobian(i,:));
   error(i) = error(i)/sc;
   Jacobian(i,:) = Jacobian(i,:)/sc;
end

end

%------------------------------------------------------------------------%  
%   c_1 = c_vector_C(J,Lsol,R1,R2,mu,z)
%   constants in equation (solenoid side of system of equations) for case (C) where the
%   optimal solution is 1 cylindrical shell
%
%   INPUTS
%   J: J value for solenoid
%   Lsol: length of solenoid
%   R1: inner radius of solenoid
%   R2: outer radius of solenoid
%   mu: magnetic permeability 
%   z: z location in solenoid (0 if at center - origin is at center)
%
%   OUTPUTS
%   c_1 = 3x1 vector of constants based on solenoid paramters for the 3
%   equations used in case C
%
%   REMARKS
%   equation solving for C case [b1; b3; Bz;];
%
%   AUTHOR(S): Paige Husa
%
%   MODIFICATIONS: 
%                  v1.0 4/29/2019
% ----------------------------------------------------------------------- %

function c_1 = c_vector_C(J,Lsol,R1,R2,mu,z)
%C_VECTOR_C
%    C_1 = C_VECTOR_C(J,LSOL,R1,R2,MU,Z)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    29-Apr-2019 19:53:38

t2 = R1.^2;
t3 = R2.^2;
t4 = Lsol.^2;
t5 = t2.^2;
t6 = t3.^2;
t7 = Lsol.*z;
t8 = t4./4.0;
t9 = z.^2;
t10 = t2+t7+t8+t9;
t11 = sqrt(t10);
t12 = R1+t11;
t13 = log(t12);
t14 = t3+t7+t8+t9;
t15 = sqrt(t14);
t16 = R2+t15;
t17 = log(t16);
t18 = t13-t17;
t19 = t2-t7+t8+t9;
t20 = sqrt(t19);
t21 = R1+t20;
t22 = log(t21);
t23 = t3-t7+t8+t9;
t24 = sqrt(t23);
t25 = R2+t24;
t26 = log(t25);
t27 = t22-t26;
c_1 = [J.*Lsol.*pi.*(R1.*t2-R2.*t3).*(-1.0./3.0);J.*Lsol.*R1.*t5.*pi.*(-3.0./4.0e1)+J.*Lsol.*R2.*t6.*pi.*(3.0./4.0e1)+(J.*Lsol.*R1.*t2.*t4.*pi)./2.4e1-(J.*Lsol.*R2.*t3.*t4.*pi)./2.4e1;(mu.*(J.*Lsol.*t18.*pi+J.*Lsol.*t27.*pi+J.*t18.*z.*pi.*2.0-J.*t27.*z.*pi.*2.0).*(-1.0./4.0))./pi];
end

%------------------------------------------------------------------------%  
%   f_1 = f_vector_C(K_1,Lc_1,Rc_1,mu,z)
%   approx geometry side of system of equations with paramters for 1 cyclinder (case C)
%
%   INPUTS
%   K_1: K value for cylinder 
%   Lc_1: length of cylinder 
%   Rc_1: radius of cylinder 
%   mu: magnetic permeability 
%   z: z location in cylinders (0 if at center - origin is at center)
%
%   OUTPUTS
%   f1 = 3x1 vector of equations based on cylinder paramters for the 3
%   equations used in case C
%
%   REMARKS
%   equations solving for C case [b1; b3; Bz];
%
%   AUTHOR(S): Paige Husa
%
%   MODIFICATIONS: 
%                  v1.0 4/29/2019
% ----------------------------------------------------------------------- %
function f_1 = f_vector_C(K_1,Lc_1,Rc_1,mu,z)
%F_VECTOR_C
%    F_1 = F_VECTOR_C(K_1,LC_1,RC_1,MU,Z)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    29-Apr-2019 19:53:38

t2 = Rc_1.^2;
t3 = Lc_1.^2;
t4 = Lc_1./2.0;
t5 = t3./4.0;
t6 = z.^2;
f_1 = [K_1.*Lc_1.*t2.*pi;(K_1.*Lc_1.*t2.*pi.*(t2.*3.0-t3))./8.0;(mu.*(K_1.*pi.*(t4+z).*1.0./sqrt(t2+t5+t6+Lc_1.*z).*2.0+K_1.*pi.*(t4-z).*1.0./sqrt(t2+t5+t6-Lc_1.*z).*2.0))./(pi.*4.0)];
end

%------------------------------------------------------------------------%  
%   Jacob1 = jacobian_matrix_C(K_1,Lc_1,Rc_1,mu,z)
%   equations for jacobian martix for case C
%
%   INPUTS
%   K_1: K value for cylinder 
%   Lc_1: length of cylinder 
%   Rc_1: radius of cylinder 
%   mu: magnetic permeability 
%   z: z location in cylinders (0 if at center - origin is at center)
%
%   OUTPUTS
%   Jacob1 = 3x3 matrix of equations based on cylinder paramters for the 3
%   equations used in case C in f_vector and derivatives with respect to the
%   3 parameters
%
%   REMARKS
%   case C equations [b1; b3; Bz];
%   parameter order [K_1, Lc_1, Rc_1];
%
%   AUTHOR(S): Paige Husa
%
%   MODIFICATIONS: 
%                  v1.0 4/29/2019
% ----------------------------------------------------------------------- %

function Jacob1 = jacobian_matrix_C(K_1,Lc_1,Rc_1,mu,z)
%JACOBIAN_MATRIX_C
%    JACOB1 = JACOBIAN_MATRIX_C(K_1,LC_1,RC_1,MU,Z)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    29-Apr-2019 19:53:38

t2 = Rc_1.^2;
t3 = Lc_1.^2;
t4 = t2.*3.0;
t5 = t3-t4;
t6 = Lc_1./2.0;
t7 = t3./4.0;
t8 = z.^2;
t9 = 1.0./pi;
t10 = Lc_1.*z;
t11 = t2+t7+t8+t10;
t12 = 1.0./sqrt(t11);
t13 = t6+z;
t14 = t6-z;
t15 = t2+t7+t8-t10;
t16 = 1.0./t11.^(3.0./2.0);
t17 = 1.0./t15.^(3.0./2.0);
Jacob1 = reshape([Lc_1.*t2.*pi,Lc_1.*t2.*t5.*pi.*(-1.0./8.0),(mu.*t9.*(t14.*pi.*1.0./sqrt(t2+t7+t8-Lc_1.*z).*2.0+t12.*t13.*pi.*2.0))./4.0,K_1.*t2.*pi,K_1.*t2.*t3.*pi.*(-1.0./4.0)-(K_1.*t2.*t5.*pi)./8.0,(mu.*t9.*(K_1.*t12.*pi+K_1.*1.0./sqrt(t15).*pi-K_1.*t13.^2.*t16.*pi-K_1.*t14.^2.*t17.*pi))./4.0,K_1.*Lc_1.*Rc_1.*pi.*2.0,K_1.*Lc_1.*Rc_1.*t2.*pi.*(3.0./4.0)-(K_1.*Lc_1.*Rc_1.*t5.*pi)./4.0,mu.*t9.*(K_1.*Rc_1.*t13.*t16.*pi.*2.0+K_1.*Rc_1.*t14.*t17.*pi.*2.0).*(-1.0./4.0)],[3,3]);
end