function snubani(Rs)
%SNUBANI(Rs)
%Runs an animation of the snubber circuit from ENGS22 lab
%Rs is the value of the snubber resistor.
%Scale is 10:1; i.e., max voltage is really 10 V even though it says 1.

%Copyright 2001, Charles R. Sullivan
%Version 1.0, Aug 6, 2001

% Components may be placed between any of the numbered nodes
% in this pre-defined grid.  Wires must also be placed
% between numbered nodes to complete the circuit.
% 7 8 9
% 4 5 6
% 1 2 3

%%%%%%%%Set up and clear figure%%%%%%%%%%%
figure(1)
clf
reset(gcf)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
ca = cell(0);
ca = addV(ca,1,4);
ca = addwire(ca,4,7);
ca = addR(ca,7,8);
ca = addL(ca,8,5);
ca = addswitch(ca,5,2,1);
ca = addR(ca,5,6);
ca = addC(ca,6,3);
ca = addwire(ca,3,1);

ia = cell(0);
ia = currentarrow(ia,7,8,1,2);
%load current, numbered current 1
ia = currentarrow(ia,8,5,1,2);
%same
ia = currentarrow(ia,5,6,2,2);
%snubber current, 2
ia = currentarrow(ia,5,2,3,2);
%switch current, 3

% calculate voltage(time)
t = linspace(0,200e-6,400);
t = t'; % make column, consistent with ode45 output
L = 575e-6;
C = 0.1366e-6;
R = 8 + Rs;
Vs = 1.2;
IL0 = Vs/8;

a = R/2/L;
b = sqrt(1/L/C - R^2/4/L/L);
alpha = Vs/L/IL0;
phi = atan2(b,alpha-a)-atan2(b,-a);
phiI = atan2(b,alpha-a);

vc =  IL0/C*(alpha/(a^2 + b^2) + 1/b*sqrt((alpha-a)^2 +b^2)/sqrt(a^2+b^2)...
   *exp(-a*t).*sin(b*t+phi));

ic = IL0*sqrt((alpha-a)^2 +b^2)/b*exp(-a*t).*sin(b*t+phiI);
vc = [0;vc];
ic = [0;ic];
t = [-eps;t];
is = zeros(size(ic));
is(1) = IL0;
iL = ic;
iL(1) = IL0;
vs = vc + ic*Rs;
vn8 = Vs - iL*8;
Vs = Vs * ones(size(t));
z = zeros(size(t));
%keyboard;
vm = 0.1*[z z z Vs vs vc Vs vn8 z];
%keyboard
% ones indicate nodes that get that voltage as a function of time.  
% in general create matrix with trajectories of all node voltages.

im = 1/max(iL)*[iL ic is];
%create a current matrix; sort of like voltage matrix but with 
%columns only for branches for which current arrows were defined.
switchm = zeros(size(t));
switchm(1) = 1;

rate = 20;

%%%%%%%%%SET UP 
zmax=1;
zmin=0;
coords=axis;
coords(5)=zmin;
coords(6)=zmax;
axis(coords);
axis equal;
set(gcf,'DoubleBuffer','on')
cvm=vm/max(max(abs(vm)));
zlabel('voltage')
view([ 0.8801 2.1249  0 -1.5025; -0.8643 0.3580 0.9135 -0.2036;   -1.9412    0.8040   -0.4067   17.7867; 0 0 0  1]);

%%%%%%%%%%%%%RUN
snubanimate(t,vm,cvm,im,ca,ia,rate,switchm);


%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%  END OF MAIN, Start of Functions

%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function snubanimate(t,vm,cvm,im,ca,ia,rate,switches);
for ti=1:length(t)
	[junk,nc]=size(ca);
	for i = 1:nc
   	updatec(ca(:,i),vm(ti,:),cvm(ti,:),switches(ti,:))
	end
	[junk,nc]=size(ia);
   for i = 1:nc
      updatei(ia(:,i),vm(ti,:),im(ti,:))
   end
   pause(1/rate)
   if ti==1
      title('t = 0- ; Hit any key (in command window) to go to 0+')
      %waitforbuttonpress
      %ginput(1);
      pause
      title('')
   
   elseif ti==2
      title('t = 0+ ;  Hit any key (in command window) to run')
      %waitforbuttonpress
      %ginput(1);
      pause
      title('')
   else
      title('Click and hold on titlebar to pause; control-C in command window to stop')
   end
end

%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
function ca = addC(ca,n1,n2)

xes = [ 0 1 2 0 1 2 0 1 2];
ys =  [ 0 0 0 1 1 1 2 2 2];

%plate1 and wire1

x1 = [0 0 -0.1 0.1];
y1 = [-0.5 -0.05 -0.05 -0.05];
z1 = zeros(size(x1));

x=[xes(n1),xes(n2)];
y=[ys(n1),ys(n2)];

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'k','linewidth',2);
hold on

zvm = zeros(9,length(z1));
zvm(n1,:)=1;
colselect=n1;

aux = zeros(0); %placeholder for future features, e.g. switches and exploding wires.
if length(ca)>0,
   ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
else
   ca = {h; [n1 n2]; zvm; colselect; aux};
end

%%%%%%%%%%Done plate1 and wire1
%%%%%%%%%%plate2 and wire2

x2 = [0 0 -0.1 0.1];
y2 = [0.5 0.05 0.05 0.05];
z2 = zeros(size(x2));

[x2,y2] = move(x,y,x2,y2);
   
h = plot3(x2,y2,z2,'k','linewidth',2);

zvm = zeros(9,length(z2));
zvm(n2,:)=1;
colselect=n2;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
function ca = addL(ca,n1,n2)

xes = [ 0 1 2 0 1 2 0 1 2];
ys =  [ 0 0 0 1 1 1 2 2 2];

%wire1

x1 = [0 0];
y1 = [-0.5 -0.35];
z1 = zeros(size(x1));

x=[xes(n1),xes(n2)];
y=[ys(n1),ys(n2)];

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'k','linewidth',2);
hold on

zvm = zeros(9,length(z1));
zvm(n1,:)=1;
colselect=n1;

aux = zeros(0); %placeholder for future features, e.g. switches and exploding wires.
if length(ca)>0,
   ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
else
   ca = {h; [n1 n2]; zvm; colselect; aux};
end

%%%%%%%%%%Done wire1

%%%%%%%%%%Now wire2

x2 = [0 0];
y2 = [0.5 0.35];
z2 = zeros(size(x2));

[x2,y2] = move(x,y,x2,y2);
   
h = plot3(x2,y2,z2,'k','linewidth',2);

zvm = zeros(9,length(z2));
zvm(n2,:)=1;
colselect=n2;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%Done wire2

%%%%%%%%%%Now Inductor!
x3 = [0 1 1.2 1 -0.5 1 1.2 1 -0.5 1 1.2 1 0]*0.1;
y3 = linspace(-0.35,0.35,13);
z3 = zeros(size(x3));

[x3,y3] = move(x,y,x3,y3);
   
h = plot3(x3,y3,z3,'color',[0.6,0.4,0],'linewidth',2);

zvm = zeros(9,length(z3));
zvm(n2,:)= linspace(0,1,13);
zvm(n1,:)=fliplr(zvm(n2,:));
colselect=0;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%


function ca = addR(ca,n1,n2)

xes = [ 0 1 2 0 1 2 0 1 2];
ys =  [ 0 0 0 1 1 1 2 2 2];

%wire1

x1 = [0 0];
y1 = [-0.5 -0.35];
z1 = zeros(size(x1));

x=[xes(n1),xes(n2)];
y=[ys(n1),ys(n2)];

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'k','linewidth',2);
hold on

zvm = zeros(9,length(z1));
zvm(n1,:)=1;
colselect=n1;

aux = zeros(0); %placeholder for future features, e.g. switches and exploding wires.
if length(ca)>0,
   ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
else
   ca = {h; [n1 n2]; zvm; colselect; aux};
end

%%%%%%%%%%Done wire1

%%%%%%%%%%Now wire2

x2 = [0 0];
y2 = [0.5 0.35];
z2 = zeros(size(x2));

[x2,y2] = move(x,y,x2,y2);
   
h = plot3(x2,y2,z2,'k','linewidth',2);

zvm = zeros(9,length(z2));
zvm(n2,:)=1;
colselect=n2;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%Done wire2

%%%%%%%%%%Now Resistor!
x3 = [0 -1 1 -1 1 -1 1 -1 0]*0.1;
y3 = [-0.35 -0.3:0.1:0.3 0.35];
z3 = zeros(size(x3));

[x3,y3] = move(x,y,x3,y3);
   
h = plot3(x3,y3,z3,'b','linewidth',2);

zvm = zeros(9,length(z3));
zvm(n2,:)= [0 0.0714 0.2143 0.3571 0.5000 0.6429 0.7857  0.9286 1.0000];
zvm(n1,:)=fliplr(zvm(n2,:));
colselect=0;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ca = addswitch(ca,n1,n2,swn)

xes = [ 0 1 2 0 1 2 0 1 2];
ys =  [ 0 0 0 1 1 1 2 2 2];

%ls = length of switch/2
ls = 0.2;
%ws = width of switch
ws = 0.25;
%wire1

x1 = [0 0];
y1 = [-0.5 -ls];
z1 = zeros(size(x1));

x=[xes(n1),xes(n2)];
y=[ys(n1),ys(n2)];

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'k','linewidth',2);
hold on

zvm = zeros(9,length(z1));
zvm(n1,:)=1;
colselect=n1;

aux = zeros(0); %placeholder for future features, e.g. switches and exploding wires.
if length(ca)>0,
   ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
else
   ca = {h; [n1 n2]; zvm; colselect; aux};
end

%%%%%%%%%%Done wire1
%%%%%%%%%%wire2

x1 = [0 0];
y1 = [0.5 ls];
z1 = zeros(size(x1));

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'k','linewidth',2);

zvm = zeros(9,length(z1));
zvm(n2,:)=1;
colselect=n2;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%wire2 terminal

x1 = [0];
y1 = [ls];
z1 = zeros(size(x1));

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'bo','linewidth',2);

zvm = zeros(9,length(z1));
zvm(n2,:)=1;
colselect=n2;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
%%%%%%%%%wire 1 terminal

x1 = [0];
y1 = [-ls];
z1 = zeros(size(x1));

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'bo','linewidth',2);

zvm = zeros(9,length(z1));
zvm(n1,:)=1;
colselect=n1;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%%%%%%%%%%%%%%%%%% switchleg
x1 = [ 0 -ws];
y1 = [-ls -ls+sqrt(4*ls^2-ws^2)];
z1 = zeros(size(x1));

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'b-','linewidth',2);

zvm = zeros(9,length(z1));
zvm(n1,:)=1;
colselect=0;
aux(1,1) = swn;
aux(1,2) = 0;
aux(2,:) = x1;  %x1 with switch off
aux(3,:) = [x1(1) x1(1)]; % x1 with swtich on
aux(4,:) = y1;
aux(5,:) = [y1(1), y1(1) + sign(y1(2)-y1(1))*2*ls ];
ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%
function ca = addV(ca,n1,n2)

xes = [ 0 1 2 0 1 2 0 1 2];
ys =  [ 0 0 0 1 1 1 2 2 2];

% do three plates
%wire1 & plate 1

x1 = [0 0 -0.05 0.05]*2;
y1 = [-0.5 -0.325 -0.325 -0.325];
z1 = zeros(size(x1));

x=[xes(n1),xes(n2)];
y=[ys(n1),ys(n2)];

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'k','linewidth',2);
hold on

zvm = zeros(9,length(z1));
zvm(n1,:)=1;
colselect=n1;

aux = zeros(0); %placeholder for future features, e.g. switches and exploding wires.
if length(ca)>0,
   ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
else
   ca = {h; [n1 n2]; zvm; colselect; aux};
end

%%%%%%%%%%Done wire1

%%%%%%%%%%Now wire2, top plate

x1 = [0 0 -0.1 0.1]*2;
y1 = [0.5 0.325 0.325 0.325];
z1 = zeros(size(x1));

[x1,y1] = move(x,y,x1,y1);
   
h = plot3(x1,y1,z1,'k','linewidth',2);
hold on

zvm = zeros(9,length(z1));

zvm(n2,:)=1;
colselect=n2;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%Done wire2

%%%%%%%%%%Now mid plates 1
x3 = [-0.1 0.1 0 0 -0.05 0.05]*2;
y3 = [-0.3 -0.3 -0.3 -0.2 -0.2 -0.2]+.025;
z3 = zeros(size(x3));

[x3,y3] = move(x,y,x3,y3);
   
h = plot3(x3,y3,z3,'b','linewidth',2);

zvm = zeros(9,length(z3));
zvm(n2,:)= ones(size(z3))/5;
zvm(n1,:)=zvm(n2,:)*4;
colselect=0;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%Now mid plates 2
x3 = [-0.1 0.1 0 0 -0.05 0.05]*2;
y3 = [-0.15 -0.15 -0.15 -0.05 -0.05 -0.05]+.025;
z3 = zeros(size(x3));

[x3,y3] = move(x,y,x3,y3);
   
h = plot3(x3,y3,z3,'b','linewidth',2);

zvm = zeros(9,length(z3));
zvm(n1,:)= ones(size(z3))*3/5;
zvm(n2,:)=ones(size(z3))*2/5;
colselect=0;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%Now mid plates 3
x3 = [-0.1 0.1 0 0 -0.05 0.05]*2;
y3 = [0.05 0.05 0.05 0.15 0.15 0.15]-.025;
z3 = zeros(size(x3));

[x3,y3] = move(x,y,x3,y3);
   
h = plot3(x3,y3,z3,'b','linewidth',2);

zvm = zeros(9,length(z3));
zvm(n1,:)= ones(size(z3))*2/5;
zvm(n2,:)=ones(size(z3))*3/5;
colselect=0;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});

%%%%%%%%%%Now mid plates 4
x3 = [-0.1 0.1 0 0 -0.05 0.05]*2;
y3 = [0.2 0.2 0.2 0.3 0.3 0.3]-.025;
z3 = zeros(size(x3));

[x3,y3] = move(x,y,x3,y3);
   
h = plot3(x3,y3,z3,'b','linewidth',2);

zvm = zeros(9,length(z3));
zvm(n1,:)= ones(size(z3))/5;
zvm(n2,:)=zvm(n1,:)*4;
colselect=0;

ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});



%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%

function ca = addwire(ca,n1,n2)

xes = [ 0 1 2 0 1 2 0 1 2];
ys =  [ 0 0 0 1 1 1 2 2 2];

x=[xes(n1),xes(n2)];
y=[ys(n1),ys(n2)];
z=zeros(size(x));

h = plot3(x,y,z,'k','linewidth',2);
hold on

zvm = zeros(9,length(z));
zvm(n1,1)=1;
zvm(n2,2)=1;
colselect=n1;

aux = zeros(0); %placeholder for future features, e.g. switches and exploding wires.
if length(ca)>0,
   ca = cat(2,ca,{h; [n1 n2]; zvm; colselect; aux});
else
   ca = {h; [n1 n2]; zvm; colselect; aux};
end

%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%

function ia = currentarrow(ia,n1,n2,cn,x)
offset=0.2;

xes = [ 0 1 2 0 1 2 0 1 2];
ys =  [ 0 0 0 1 1 1 2 2 2];

%Arrow

x0 = [ 0 0 -0.1 0 0.1];
y0 = [0.4 -0.4 -0.3 -0.4 -0.3];
z0 = zeros(size(x0));

x=[xes(n2),xes(n1)];
y=[ys(n2),ys(n1)];

[xplot,yplot] = move(x,y,x0+offset,y0);
   
h = plot3(xplot,yplot,z0,'color',[0 0.8 0],'linewidth',2);
hold on

zvm = zeros(9,length(z0));
zvm(n1,:)= ([x0]+0.5)/0.7;
zvm(n2,:)=-([x0]-0.35)/0.7;


if length(ia)>0,
   ia = cat(2,ia,{h; [n1 n2]; zvm;cn ;x0;x;y0;y});
else
   ia = {h; [n1 n2]; zvm;cn ;x0;x;y0;y};
end

% iac contains
% {graphics handle; [n1, n2];zv matrix; current number ;Arrx;xpos;Arry;ypos}
% ivec is a row of each numbered current.

%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%

function [xnew,ynew] = move(x,y,x1,y1)
% move component at 0,0, axis y, 
% feature locations x1,y1
% to new coordinates centered halfway between xes in
% x and ys in y, with axis whereever axis of 

v=[x1;y1];
rotate = [0 1;-1 0];

%keyboard

if y(1) > y(2) & x(1) == x(2)
   %need 180
   v = rotate*rotate*v;
elseif y(1) < y(2) & x(1) == x(2)
   %already oriented right
elseif y(1) == y(2) & x(1) > x(2)
   v = rotate^3*v;  
elseif y(1) == y(2) & x(1) < x(2)
   v = rotate*v;     
else
   error('Non Manhattan Geometry Specified')
end
  
  
xr = v(1,:);
yr = v(2,:);
  
xnew = xr + mean(x);
ynew = yr +mean(y);
 
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%

function updatec(cac,vvec,cvvec,switchm)
%keyboard

% cac contains
%{graphics handle; [n1, n2]; zv matrix; colselect; aux data}
% vvector is a row of voltage at each node.

h = cac{1};
%keyboard
zvm = cac{3};
zvalues = vvec*zvm;
aux = cac{5};
colselect = cac{4};

if colselect ~= 0
   colorvalue=cvvec(colselect);
   color = sqrt(abs(colorvalue))*([0 1 0]*(colorvalue < 0) + [1 0 0]);
   set(h,'zdata',zvalues,'color',color)
else
   % Color not to change, e.g. resistor
   set(h,'zdata',zvalues)
end

if any(aux)
   xvalues = aux(2,:)*(1-switchm(aux(1,1)))+aux(3,:)*switchm(aux(1));
   yvalues = aux(4,:)*(1-switchm(aux(1,1)))+aux(5,:)*switchm(aux(1));
   set(h,'xdata',xvalues,'ydata',yvalues)
   %keyboard
end
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function updatei(iac,vvec,ivec)

% iac contains
% {graphics handle; [n1, n2];zv matrix; current number ;Arrx;xpos;Arry;ypos}
% ivec is a row of each numbered current.
offset=0.2; %distance from branch

h = iac{1};
nodes = iac{2};
n1=nodes(1);
n2= nodes(2);

zvm = iac{3};

iselect = iac{4};
ivalue=ivec(iselect);
scale = sign(ivalue)*sqrt(abs(ivalue));
hypfact = 1/sqrt(1+(vvec(n1)-vvec(n2))^2);  %factor so length doesn't stretch vertically;

Arrx = iac{5};
Arrxs = Arrx*scale;
xpos = iac{6};
Arry = iac{7};
Arrys=Arry*scale*hypfact;
ypos = iac{8};
[xplot,yplot] = move(xpos,ypos,Arrxs+0.2,Arrys);

zvm(n1,:)= ([Arrys]+0.35)/0.7;
zvm(n2,:)=-([Arrys]-0.35)/0.7;

zvalues = vvec*zvm;

set(h,'zdata',zvalues,'xdata',xplot,'ydata',yplot,'linewidth',5*abs(scale)+eps);
%keyboard