ComputeBfieldScriptIntegrate.m

%% Header
% Author: Nathan Eckert
% Date Created : 03 - 28 - 2017
% Date Edited  : 03 - 28 - 2017
% Title: MagneticFieldFunctionTest
clear all; close all; clc;
%%%%%%%%%%%% Constants / Inputs / Initial Conditions %%%%%%%%%%%%
%% Physical Parameters
% Rail Length assumed infinite
w = 0.01; %     [m] - Rail Width
h = 0.01; %     [m] - Rail Height
s = 0.012 ; %   [m] - Rail Separation
at = 0.0296 ; % [m] - Armature Thickness
ah = h; %       [m] - Armature Height
aw = s; %       [m] - Armature Width
mu0 = 4*pi*10^-7; % [N/A^2] - Magnetic Permeability of Vaccuum
I0 = 20e3; %    [A] - Constant Maximum Current
tau = 0.0007; % [Ohms-s] - Capacitor Time Constant
%% Simulation Parameters
t = 0.001; %        [s] - Simulation Time
l = 0.01; %         [m] - Simulated Armature Position 
t1 = 0.0005; %      [s] - Current Profile Times
t2 = 0.0015; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Changed For Comparison to Analytic solvers in mathematica
t3 = 0.003;
Nm = 10; %          [~] - Number of Mesh Points
Natm = 5; %        [~] - Number of Armature Thickness Values
[zp0,zpF,zpPts,dzp,zpVec] = MakeVector(0,w,'num',Nm); 
    % Flipping Direction has no impact !!!
[yp0,ypF,ypPts,dyp,ypVec] = MakeVector(-h/2,h/2,'num',Nm); 
    % Flipping Direction makes a few values non-NaN, though they compute as smaller
[s0,sF,sPts,ds,sVec] = MakeVector(0,s,'num',Nm);
[x0,xF,xPts,dx,xVec] = MakeVector(l,l+at,'num',Natm); 
    % Flipping Direction has no impact !!!
[y0,yF,yPts,dy,yVec] = MakeVector(-ah/2,ah/2,'num',Nm); 
    % Flipping Direction makes a few values non-NaN, though they compute as smaller
[z0,zF,zPts,dz,zVec] = MakeVector(aw,0,'num',Nm); 
    % Flipping direction has no impact !!!
%% Initialize Lists
B = zeros(yPts,zPts,xPts);
dB = zeros(yPts,zPts,xPts);
dBc = zeros(yPts,zPts,xPts);
%% Simulation

% Pre-Loop Computations
if t <= t1                      % Compute Current 
    I = I0 * sin(omega*t);
elseif t <= t2
    I = I0;
elseif t <= t3
    I = I0*exp(-t/tau);
else
    print('Error: The simulation time entered is outside the defined range.');
    return
end
J = I / (h*w);  % Compute current density in rail

% Integral Computation 
dI = J * abs(dzp) * abs(dyp); % Compute current element at the rail point (m)
CB = (mu0*dI) / (4*pi*h*w);

%----------- Armature Loop - Rail Values Computed with Integral2 ----------
for k = 1:zPts              % Loop through armature width -----------------
    z = zVec(k);            % Define armature width position 
    for j = 1:yPts          % Loop through armature height ----------------
        y = yVec(j);        % Define armature height position
        for i = 1:xPts      % Loop through armature thickness -------------
            x = xVec(i);    % Define armature thickness position
            
            % Compute magnetic field at each armature point (n)
            dBi = @(yp,zp) CB.*(z+zp) ./ ((y-yp).^2+(z+zp).^2).*(((x) ./ sqrt((y-yp).^2+(z+zp).^2+x.^2))-((x-l) ./ sqrt((y-yp).^2+(z+zp).^2+(x-l).^2)));
            %dBi = @(yp,zp) CB*(z+zp) / ((y-yp)^2+(z+zp)^2)*(((x) / sqrt((y-yp)^2+(z+zp)^2+x^2))-((x-l) / sqrt((y-yp)^2+(z+zp)^2+(x-l)^2)));
            
            dBc(j,k,i) = dBc(j,k,i) + integral2(dBi,yp0,ypF,zp0,zpF,'AbsTol',1e-12);
                    
        end
    end
end

dBc = dBc + fliplr(dBc);


% Simple Numeric Sum Computation 
% ------------------------- Conductor (Rail) Loop -------------------------  
for kp = 1:zpPts        % Loop through rail width -------------------------
    zp = zpVec(kp);     % Define rail width position
    for jp = 1:ypPts    % Loop through rail height ------------------------
        yp = ypVec(jp); % Define rail height position
        dI = J * abs(dzp) * abs(dyp); % Compute current element at the rail point (m)
        CB = (mu0*dI) / (4*pi*h*w);
% -------------------- Area of Interest (Armature) Loop -------------------

        for k = 1:zPts              % Loop through armature width ---------
            z = zVec(k);            % Define armature width position
            for j = 1:yPts          % Loop through armature height --------
                y = yVec(j);        % Define armature height position
                for i = 1:xPts      % Loop through armature thickness -----
                    x = xVec(i);    % Define armature thickness position
                    
                    % Compute magnetic field at each armature point (n)
                    P1 = (z+zp) / ((y-yp)^2+(z+zp)^2);
                    P2 = (x) / sqrt((y-yp)^2+(z+zp)^2+x^2);
                    P3 = (x-l) / sqrt((y-yp)^2+(z+zp)^2+(x-l)^2);
                    dB(j,k,i) = dB(j,k,i) + CB*P1*(P2-P3);
                    
                end
            end
        end

        % Add Second Rail
        dB = dB + fliplr(dB);
        % Add Contribution from each rail section
        B = B + dB;
        dB = zeros(yPts,zPts,xPts);
        
    end
end


%% Plot

% Compute average B for each armature thickness value
for i = 1:xPts
    B1avg = mean(B(1,2:end-1,i));
    BmidAvg = mean(B(end,2:end-1,i));
    BendAvg = mean(B(2:end-1,2:end-1,i));
    Bavg(i) = mean([B1avg BmidAvg BendAvg]);
end
BavgMax = max(Bavg);
disp(['The largest of the average B values is: Bavg = ', num2str(BavgMax)]);
figure()
plot(xVec,Bavg)
title('Average B Field Value at Each Armature Thickness Value')
xlabel('Armature Thickness Value (x) [m]')
ylabel('Average B Field Value [T]')

% Compute average B for each armature thickness value
for i = 1:xPts
    B1cavg = mean(dBc(1,2:end-1,i));
    BcmidAvg = mean(dBc(end,2:end-1,i));
    BcendAvg = mean(dBc(2:end-1,2:end-1,i));
    Bcavg(i) = mean([B1cavg BcmidAvg BcendAvg]);
end
BcavgMax = max(Bcavg);
disp(['The largest of the average B values is: Bavg = ', num2str(BcavgMax)]);
figure()
plot(xVec,Bcavg)
title('Average B Field Value at Each Armature Thickness Value - Integral2')
xlabel('Armature Thickness Value (x) [m]')
ylabel('Average B Field Value [T]')

figure()
hold on
% Plots -> (x,y) || Matrices -> (rows,cols) = (y,x)
for pIdx = 1:xPts
    surf(zVec,yVec,B(:,:,pIdx),'FaceColor',rand(1,3),'FaceAlpha',0.25,'EdgeColor','black','LineWidth',0.5)
    view([33 25])
end
grid on
title('B Field Value at Each Point in the Armature')
xlabel('Armature Width [m]')
ylabel('Rail Height [m]')
zlabel('Magnetic Field Strength [T]')

figure()
hold on
% Plots -> (x,y) || Matrices -> (rows,cols) = (y,x)
for pIdx = 1:xPts
    surf(zVec,yVec,dBc(:,:,pIdx),'FaceColor',rand(1,3),'FaceAlpha',0.25,'EdgeColor','black','LineWidth',0.5)
    view([33 25])
end
grid on
title('B Field Value at Each Point in the Armature - Integral2')
xlabel('Armature Width [m]')
ylabel('Rail Height [m]')
zlabel('Magnetic Field Strength [T]')
%myaa(4);

% % Plot B field Plot Between Rails
% figure((get(gcf,'Number')+1))
% hold on
% imagesc(B(:,:,1))
% contour(B(:,:,1),20)
% title('Magnetic Field Strength Between Rails')
% ylabel('Rail Height [m]')
% xlabel('Rail Separation [m]')
% colormap jet
% colorbarH = colorbar;
% xlabel(colorbarH,'B field strength [T]')
% 
% % Plot Contour of B field between rails
% figure((get(gcf,'Number')+1))
% contour(B(:,:,1),30)
% title('Magnetic Field Strength Between Rails')
% ylabel('Rail Height [m]')
% xlabel('Rail Separation [m]')

% Plot m-n
% figure((get(gcf,'Number')+1))
% plot3([0,)

% % Plot B field from top perspective
% figure((get(gcf,'Number')+1))
% %                                    Y X Z
% Btop = sum(B,1); Btop = squeeze(Btop(1,:,:));
% mesh(atIdx,sIdx,Btop)
% xlabel('Armature Thickness [m]')
% ylabel('Rail Separation [m]')
% zlabel('Magnetic Field Strength [??]')

%% GIF
gif = 0;
if gif == 1
figure((get(gcf,'Number')+1))
axis tight
% set(gca,'nextplot','replacechildren','visible','off')
% f = getframe;
% [im,map] = rgb2ind(f.cdata,256,'nodither');
% [~,map(2,:)] = rgb2ind(f.cdata,256,'nodither');
% im(1,1,1,(meshNum^2+meshNum+2^2)) = 0;
% k=0;
for kp = 1:zpPts        % Loop through rail width -------------------------
    zp = zpVec(kp);     % Define rail width position
    for jp = 1:ypPts    % Loop through rail height ------------------------
        yp = ypVec(jp); % Define rail height position
        dI = J * dzp * dyp; % Compute current element at the rail point (m)
        CB = (mu0*dI) / (4*pi*h*w);
% -------------------- Area of Interest (Armature) Loop -------------------    
        for k = 1:zPts              % Loop through armature width ---------
            z = zVec(k);            % Define armature width position
            for j = 1:yPts          % Loop through armature height --------
                y = yVec(j);        % Define armature height position
                for i = 1:xPts      % Loop through armature thickness -----
                    x = xVec(i);    % Define armature thickness position
                    k=k+1;
                    
                    plot3([-zp,z],[0,x],[yp,y]); 
                    grid on
                    %f = getframe;
                    %im(:,:,1,k) = rgb2ind(f.cdata,map,'nodither');
                    view([45 45])
                    xlim([-w w]);
                    ylim([0 0.05]);
                    zlim([-h/2 h/2]);
                    % gif utilities
                    set(gcf,'color','w'); % set figure background to white
                    drawnow;
                    frame = getframe(gcf);
                    im = frame2im(frame);
                    [imind,cm] = rgb2ind(im,256);
                    outfile = 'BfieldComputationLines.gif';

                    % On the first loop, create the file. In subsequent loops, append.
                    if kp == 1 && i == 1
                        imwrite(imind,cm,outfile,'gif','DelayTime',0,'loopcount',inf);
                    else
                        imwrite(imind,cm,outfile,'gif','DelayTime',0,'writemode','append');
                    end
                end
            end
        end
    end
end
xlabel('Z')
ylabel('X')
zlabel('Y')
%imwrite(im,map,'DancingPeaks.gif','DelayTime',0,'LoopCount',inf)
end