Tank/Field Simulations

calc_tankfield.m

function calc_tankfield

global EROOT;

if ~exist('EROOT')
    EROOT = alabstartup('snyd07a');
end

% Enable/Disable Features
tank = 1;

% These are the locations of the simulated tank edges
tanklims_x = [-0.380 0.380];
tanklims_y = [-0.370 0.370];
tanklims_z = [-0.060 0.060];

% Plotting Range/Scale
view_scale = 1;
x_min = tanklims_x(1)*view_scale;
x_max = tanklims_x(2)*view_scale;
y_min = tanklims_y(1)*view_scale;
y_max = tanklims_y(2)*view_scale;

 
pscale = 1.2;

% Electrode Configuration
elec_sep = 0.0254;
elec_lat_dist = 0.035;

% Object Position Parameters
nsteps = 100;
obj_lat_dist = 0.350; % m
obj_height = 0; % m
obj_radius = 0.010; % m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pole, Electode & Object Parameters/Positions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Positions of Poles (mm) in x,y,z
poles = [0  0.090 0; 
         0 -0.090 0];
     
q = 0.012; % This is somewhat of a fudge factor

q_poles = [-q ; q];


% Create Electrode Arrangement
elec_pos = [elec_lat_dist -elec_sep*2 0; -elec_lat_dist -elec_sep*2 0;
            elec_lat_dist -elec_sep   0; -elec_lat_dist -elec_sep   0; 
            elec_lat_dist 0           0; -elec_lat_dist 0           0;
            elec_lat_dist elec_sep    0; -elec_lat_dist elec_sep    0;
            elec_lat_dist elec_sep*2  0; -elec_lat_dist elec_sep*2  0];


% Generate Object Positions
objpos = [linspace(0,0,nsteps)'+obj_lat_dist linspace(-0.130,0.130,nsteps)' linspace(0,0,nsteps)'+obj_height];
xyz = objpos;

% Compute primary electrostatic field from poles
exyz  =  compute_efield_general(xyz, q_poles, poles);

% Compute Tank/Plane Insulator Effects
plane_dists = [tanklims_x' zeros(2,2); zeros(2,1) tanklims_y' zeros(2,1); zeros(2,2) tanklims_z'];
plane_norms = [1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1];

if tank == 1
    [exyz, tank_poles_all] = apply_tank_effect(xyz, exyz, poles, q_poles, plane_dists, plane_norms);
end

[all_dphi_diff, all_dphi] = compute_electrode_potentials(exyz,elec_pos,objpos,obj_radius);

% Generate Electrode Potentials Plot
figure(1);
plot(objpos(:,2),all_dphi_diff);
xlabel('Position (mm)');
ylabel('Voltage');
if tank == 1
    title(['Chen Model Object Pass, dist: ' sprintf('%2.1f',(obj_lat_dist - elec_lat_dist)*1000) ' mm (center from elec)'])
else
    title(['Chen Model Object Pass (notank), dist: ' sprintf('%2.1f',(obj_lat_dist - elec_lat_dist)*1000) ' mm (center from elec)'])
end
drawnow;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate Field Slice
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generate Second Plot of Field Mag
[x, y, z] = meshgrid(linspace(x_min,x_max,50),linspace(y_min,y_max,50),0);
[nx,ny,nz] = size(x);

xyz = [x(:) y(:) z(:)];

% Compute primary electrostatic field from poles
exyz  =  compute_efield_general(xyz, q_poles, poles);
if tank == 1
    [exyz, tank_poles_all] = apply_tank_effect(xyz, exyz, poles, q_poles, plane_dists, plane_norms);
end

mag_exyz = sqrt(sum(exyz.^2,2)).*sign(exyz(:,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate Field Slice Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f2 = figure(2); clf;
f2pos = get(f2,'Position');
set(f2, 'Name','Electric Field Slice: X-Y Plane','NumberTitle','off', ...
    'Units', 'pixels', ...
    'Position', [f2pos(1:2) 650 630]);
h = subplot(2,2,1);
contourf(squeeze(x),squeeze(y),log(real(squeeze(reshape(mag_exyz,nx,ny,nz)))),20,'LineStyle','--','LineWidth',0.2,'LineColor',[0.2 0.2 0.2]);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
if tank == 1
    scatter(tank_poles_all(:,1),tank_poles_all(:,2),'wo')
end
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
title('Field Magnitude')
fix_subplot(h,pscale);

% Generate Field Angles Plot
h = subplot(2,2,2);
contourf(squeeze(x),squeeze(y),squeeze(reshape(atan2(exyz(:,2),exyz(:,1)),nx,ny,nz)).*(360/(2*pi)),[-pi:pi/8:pi].*(360/(2*pi)),'LineStyle','--','LineWidth',0.2,'LineColor',[0.2 0.2 0.2]);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
if tank == 1
    scatter(tank_poles_all(:,1),tank_poles_all(:,2),'wo')
end
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
title('Field Angles')
fix_subplot(h,pscale);


% X Vectors
h = subplot(2,2,3);
contourf(squeeze(x),squeeze(y),log(real(squeeze(reshape(exyz(:,1),nx,ny,nz)))),20,'LineStyle','--','LineWidth',0.2,'LineColor',[0.2 0.2 0.2]);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
if tank == 1
    scatter(tank_poles_all(:,1),tank_poles_all(:,2),'wo')
end
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
title('X Field Vectors')
fix_subplot(h,pscale);


% Y Vectors
h = subplot(2,2,4);
contourf(squeeze(x),squeeze(y),log(real(squeeze(reshape(exyz(:,2),nx,ny,nz)))),20,'LineStyle','--','LineWidth',0.2,'LineColor',[0.2 0.2 0.2]);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
if tank == 1
    scatter(tank_poles_all(:,1),tank_poles_all(:,2),'wo')
end
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
title('Y Field Vectors')
fix_subplot(h,pscale);


function fix_subplot(h,pscale)
p = get(h, 'pos');
set(h, 'pos', [p(1:2)-(p(3:4).*(pscale-1))./2 p(3:4).*pscale]);

set(h,'FontName', 'Helvetica' );
set(h,'FontSize', 8 );

axis equal;




function [exyz, tank_poles_all] = apply_tank_effect(xyz, exyz, poles, q_poles, plane_dists, plane_norms)
% Apply tank effects to field

tank_poles_all = [];

for i = 1:length(plane_dists)
    plane_col = plane_dists(i,:)~=0;
    tank_poles = poles - 2.*repmat((poles*plane_norms(i,:)' + abs(plane_dists(i,plane_col))),1,3).*repmat(plane_norms(i,:),size(poles,1),1);
    exyz = exyz + compute_efield_general(xyz, -q_poles, tank_poles);

    tank_poles_all = [tank_poles_all; tank_poles];
end

function [all_dphi_diff, all_dphi] = compute_electrode_potentials(exyz,elec_pos,objpos,obj_radius)

n_elec = size(elec_pos,1);

% Iteratively Compute potentials at Electrodes
all_dphi = [];
all_dphi_diff = [];

for i = 1:size(objpos,1)
    % Compute potential with dipole
    r       =  repmat(objpos(i,:),n_elec,1) - elec_pos;
    rnorm   =  sqrt(sum(r.^2,2));
    rht     =  (obj_radius./rnorm).^3;
    dphi    =  -sum(repmat(exyz(i,:), n_elec, 1).*r,2) .* rht;
    
    dphi_diff = dphi(1:2:end)-dphi(2:2:end);
    all_dphi  = [all_dphi; dphi'];
    all_dphi_diff = [all_dphi_diff; dphi_diff'];
end

calc_tankfield_evenfield.m

% PARAMS
%
% Electrode Properties
elec_sep      = 0.04; % m
elec_lat_dist = 0.02; % m
elec_rows     = 5; % m

% Object Properties
obj_lat_dist = 0.04; % m
obj_height   = 0;     % m
obj_radius   = 0.010; % m

% Field Limits
x_min = -0.1;
x_max = 0.1;
y_min = -0.1;
y_max = 0.1;

field_max = 1;
field_min = -1;



% Create Electrode Arrangement
if mod(elec_rows,2)
    elec_y = repmat(((-(elec_rows-1)/2):((elec_rows-1)/2))*elec_sep,2,1);
else
    elec_y = repmat((((-(elec_rows)/2):((elec_rows-1)/2))+0.5)*elec_sep,2,1);
end

elec_pos = [repmat([elec_lat_dist; -elec_lat_dist],elec_rows,1) elec_y(:) repmat(0,elec_rows*2,1)];
n_elec = size(elec_pos,1);

% Create Object Trajectory
nsteps = 100;
objpos = [linspace(0,0,nsteps)'+obj_lat_dist linspace(-0.1,0.1,nsteps)' linspace(0,0,nsteps)'+obj_height];

% Make grid for plotting and create "field"
[x, y, z] = meshgrid(linspace(x_min,x_max,100),linspace(y_min,y_max,100),0);
[nx,ny,nz] = size(x);
[exyz_grid_x, exyz_grid_y, exyz_grid_z] = meshgrid(linspace(field_max,field_min,nx),linspace(field_min,field_max,ny),0);

exyz = [interp2(x,y,exyz_grid_x,objpos(:,1),objpos(:,2)) ... % x component
        interp2(x,y,exyz_grid_y,objpos(:,1),objpos(:,2)) ... % y component
        interp2(x,y,exyz_grid_z,objpos(:,1),objpos(:,2))];   % z component

% Iteratively Compute potentials at Electrodes Using Dipoles
all_dphi = [];
all_dphi_diff = [];

for i = 1:size(objpos,1)
    % Compute potential with dipole
    r       =  repmat(objpos(i,:),n_elec,1) - elec_pos;
    rnorm   =  sqrt(sum(r.^2,2));
    rht     =  (obj_radius./rnorm).^3;
    dphi    =  -sum(repmat(exyz(i,:), n_elec, 1).*r,2) .* rht;
    
    dphi_diff = dphi(1:2:end)-dphi(2:2:end);
    all_dphi  = [all_dphi; dphi'];
    all_dphi_diff = [all_dphi_diff; dphi_diff'];
end


%% Generate Differential Electrode Output Plot
figure(1);
plot(objpos(:,2),all_dphi_diff);
xlabel('Position (mm)');
ylabel('Voltage');
title(['Even field Object Pass, dist: ' num2str(obj_lat_dist - elec_lat_dist)*1000 ' mm (center from elec)'])

% Generate Second Plot of Field Mag
xyz = [x(:) y(:) z(:)];
exyz = [exyz_grid_x(:) exyz_grid_y(:) exyz_grid_z(:)];
mag_exyz = sqrt(sum(exyz.^2,2));


%% Generate Field Slice Plots
figure(2); clf;
subplot(2,2,1)
contourf(squeeze(x),squeeze(y),real(squeeze(reshape(mag_exyz,nx,ny,nz))),20);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
ylabel('Y position (m)')
title('Field Magnitude')

% Generate Field Angles Plot
subplot(2,2,2)
contourf(squeeze(x),squeeze(y),squeeze(reshape(atan2(exyz(:,2),exyz(:,1)),nx,ny,nz)).*(360/(2*pi)),[-pi:pi/8:pi].*(360/(2*pi)));
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
title('Field Angles')

% X Vectors
subplot(2,2,3)
contourf(squeeze(x),squeeze(y),real(squeeze(reshape(exyz(:,1),nx,ny,nz))),20);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
ylabel('Y position (m)')
title('X Field Vectors')

% Y Vectors
subplot(2,2,4)
contourf(squeeze(x),squeeze(y),real(squeeze(reshape(exyz(:,2),nx,ny,nz))),20);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'w.')
plot(objpos(:,1),objpos(:,2),'k','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
ylabel('Y position (m)')
title('Y Field Vectors')

calc_tankfield_uniform.m

% PARAMS
%
% Electrode Properties
elec_sep      = 0.025; % m
elec_lat_dist = 0.035; % m
elec_rows     = 5; % m

% Object Properties
obj_lat_dist = 0.065; % m
obj_height   = 0;     % m
obj_radius   = 0.010; % m
obj_excursion= 0.1;   % m, obj travels from -obj_excursion in y to obj_excursion


% Field Limits
x_min = -0.1;
x_max = 0.1;
y_min = -0.1;
y_max = 0.1;

field_mag = 1;



% Create Electrode Arrangement
if mod(elec_rows,2)
    elec_y = repmat(((-(elec_rows-1)/2):((elec_rows-1)/2))*elec_sep,2,1);
else
    elec_y = repmat((((-(elec_rows)/2):((elec_rows-1)/2))+0.5)*elec_sep,2,1);
end

elec_pos = [repmat([elec_lat_dist; -elec_lat_dist],elec_rows,1) elec_y(:) repmat(0,elec_rows*2,1)];
n_elec = size(elec_pos,1);

% Create Object Trajectory
nsteps = 100;
objpos = [linspace(0,0,nsteps)'+obj_lat_dist linspace(-obj_excursion,obj_excursion,nsteps)' linspace(0,0,nsteps)'+obj_height];

% Make grid for plotting and create "field"
[x, y, z] = meshgrid(linspace(x_min,x_max,100),linspace(y_min,y_max,100),0);
[nx,ny,nz] = size(x);
[exyz_grid_x, exyz_grid_y, exyz_grid_z] = meshgrid(linspace(0,0,nx),linspace(field_mag,field_mag,ny),0);

exyz = [interp2(x,y,exyz_grid_x,objpos(:,1),objpos(:,2)) ... % x component
        interp2(x,y,exyz_grid_y,objpos(:,1),objpos(:,2)) ... % y component
        interp2(x,y,exyz_grid_z,objpos(:,1),objpos(:,2))];   % z component

% Iteratively Compute potentials at Electrodes Using Dipoles
all_dphi = [];
all_dphi_diff = [];

for i = 1:size(objpos,1)
    % Compute potential with dipole
    r       =  repmat(objpos(i,:),n_elec,1) - elec_pos;
    rnorm   =  sqrt(sum(r.^2,2));
    rht     =  (obj_radius./rnorm).^3;
    dphi    =  -sum(repmat(exyz(i,:), n_elec, 1).*r,2) .* rht;
    
    dphi_diff = dphi(1:2:end)-dphi(2:2:end);
    all_dphi  = [all_dphi; dphi'];
    all_dphi_diff = [all_dphi_diff; dphi_diff'];
end


%% Generate Differential Electrode Output Plot
figure(1);
subplot(3,1,1)
plot(objpos(:,2),all_dphi_diff);
xlabel('Position (mm)');
ylabel('Voltage (diff)');
title(['Uniform Field Object Pass, dist: ' num2str(obj_lat_dist - elec_lat_dist)*1000 ' mm (center from elec)'])

subplot(3,1,2);
plot(objpos(:,2),all_dphi(:,2:2:end));
xlabel('Position (mm)');
ylabel('Voltage (left)');
title(['Uniform Field Object Pass, dist: ' num2str(obj_lat_dist - elec_lat_dist)*1000 ' mm (center from elec)'])

subplot(3,1,3);
plot(objpos(:,2),all_dphi(:,1:2:end));
xlabel('Position (mm)');
ylabel('Voltage (right)');
title(['Uniform Field Object Pass, dist: ' num2str(obj_lat_dist - elec_lat_dist)*1000 ' mm (center from elec)'])


% Generate Second Plot of Field Mag
xyz = [x(:) y(:) z(:)];
exyz = [exyz_grid_x(:) exyz_grid_y(:) exyz_grid_z(:)];
mag_exyz = sqrt(sum(exyz.^2,2));


%% Generate Field Slice Plots
figure(2); clf;
subplot(2,2,1)
contourf(squeeze(x),squeeze(y),real(squeeze(reshape(mag_exyz,nx,ny,nz))),20);
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'k.')
plot(objpos(:,1),objpos(:,2),'b','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
ylabel('Y position (m)')
title('Field Magnitude')

% Generate Field Angles Plot
subplot(2,2,2)
contourf(squeeze(x),squeeze(y),squeeze(reshape(atan2(exyz(:,2),exyz(:,1)),nx,ny,nz)).*(360/(2*pi)),[-pi:pi/8:pi].*(360/(2*pi)));
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'k.')
plot(objpos(:,1),objpos(:,2),'b','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
title('Field Angles')

% X Vectors
subplot(2,2,3)
%quiver(squeeze(x),squeeze(y),real(squeeze(reshape(exyz(:,1),nx,ny,nz))),20);
quiver(x(:),y(:),exyz(:,1),exyz(:,2).*0,0.5)
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'k.')
plot(objpos(:,1),objpos(:,2),'b','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
ylabel('Y position (m)')
title('X Field Vectors')

% Y Vectors
subplot(2,2,4)
quiver(x(:),y(:),exyz(:,1).*0,exyz(:,2),0.5)
hold on;
scatter(elec_pos(:,1),elec_pos(:,2),'k.')
plot(objpos(:,1),objpos(:,2),'b','LineWidth',2);
axis equal;
set(gca,'xlim',[x_min x_max],'ylim',[y_min y_max])
xlabel('X Position (m)')
ylabel('Y position (m)')
title('Y Field Vectors')