% Author: Ankur Kumar (ak364)
% References: Jon del Gaizo and Stephen Jhun and Prof Hutchins
%coupled spring mass in parallel coordinates using numerical method in 1D
%
%  |                              |
%  |/\/\/\/\|M|/\/\/\/\|M|/\/\/\/\|
%  |                              |
% 
%premise: each spring has a different k, each mass can be different, the
%system can be started in any setting to create displacements of the two masses (x1 and x2). 
%Analytically, the equation is cumbersome, but numerically (using Newton + Hooke) life can be easier
%applet: http://www.lon-capa.org/~mmp/applist/coupled/osc2.htm

%Governing equations (i think): 
%M_1*x" = -K_1*x_1 - K_2*x_1 + K_2*x_2
%M_2*x" = -K_3*x_2 - K_2*x_2 + K_2*x_2

clf
map = winter(100);
m1=input('give value for first mass ');
m2=input('give value for second mass ');
k1=input('give value for first k ');
k2=input('give value for second k ');
k3=input('give value for third k ');
x1=input('where is first mass stretched to? ');
x2=input('what is second mass stretched to? ');
xx1=[];
xx2=[];

v1=0;
v2=0;

T=.1; %deltaT
for n=1:100
    a1=(-k1*x1-k2*x1+k2*x2)/m1; %acceleration of x1
    a2=(-k3*x2-k2*x2+k2*x1)/m2; %acceleration of x2
    v1=v1+a1*T;
    v2=v2+a2*T;
    x1=x1+v1*T;
    x2=x2+v2*T;
    xx1(n)=x1;
    xx2(n)=x2;
end
% tempx1 = max(xx1)/100;
% tempx2 = max(xx2)/100;
tempx1 = 1/100; %scale output
tempx2 = 1/100; %scale output
figure(1)
subplot(2,2,1)
plot([0:99].*tempx1, xx1)
subplot(2,2,2)
plot([0:99].*tempx2, xx2, 'r')

%parallel coordinates
subplot(2,2,3);
hold on
ndim = 2;
%the axes for parallel coords
axisx = [1;ndim];
axisy = [0;0];
label = [];
for i = 1:ndim
    axisx = [axisx,[i;i]];
    axisy = [axisy,[0;1]];
end
label = [{'t'},{'x'}];
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:ndim]);
set(gca, 'xticklabel', label);
vx1=vertcat([0:99].*tempx1,xx1)'; %data points x-dim
for i=1:100
    plot(vx1(i,:),'color',map(i,:));
end

subplot(2,2,4);
hold on
ndim = 2;
%the axes for parallel coords
axisx = [1;ndim];
axisy = [0;0];
label = [];
for i = 1:ndim
    axisx = [axisx,[i;i]];
    axisy = [axisy,[0;1]];
end
label = [{'t'},{'x'}];
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:ndim]);
set(gca, 'xticklabel', label);
vx2=vertcat([0:99].*tempx2,xx2)'; %data points x-dim
for i=1:100
    plot(vx2(i,:),'color',map(i,:));
end


