% SSi_6
clear all
clc

f = @(x) sqrt(4-x.^2);
j = @(x) sqrt(0.64-x.^2);
k = @(x) -j(x);
p = @(x) (-3*(x-1).^5)-8;

nd = 2000;

z = linspace(0,1,nd);
x = 2*cos(linspace(-pi,0,nd));
%x = linspace(-2,2,nd);
xh = x(1:ceil(length(x)/2));
x2 = [x(end:-1:1),xh];
f2 = [f(x(end:-1:1)),-f(xh)]+2;

% cartesian coordinates for SSi
xS1 = [x2,-fliplr(x2)];
fS1 = [f2,-fliplr(f2)];
xS2 = xS1-4;
fS2 = fS1-4;
xd = [x(end:-1:1),x]/2.5;
fd = real([j(x(end:-1:1)/2.5),k(x/2.5)])-2;
zi = z;
fi = p(zi);

% polar coordinates for SSi (from [-2 2])
rS1 = sqrt((fS1-2).^2 + (xS1+2).^2);
tS1 = atan2((fS1-2),(xS1+2));
rS2 = sqrt((fS2-2).^2 + (xS2+2).^2);
tS2 = atan2((fS2-2),(xS2+2));
rd = sqrt((fd-2).^2 + (xd+2).^2);
td = atan2((fd-2),(xd+2));
ri = sqrt((fi-2).^2 + (zi+2).^2);
ti = atan2((fi-2),(zi+2));

% random vectors for SSi explosion
xS1r = randn(1,length(xS1));
fS1r = randn(1,length(xS1));
xS2r = randn(1,length(xS2));
fS2r = randn(1,length(xS1));
xdr = randn(1,length(xd));
fdr = randn(1,length(fd));
zir = randn(1,length(zi));
fir = randn(1,length(fi));

i = 6*ones(1,10);
i = [i,i,9:3:45,45:-3:9,i,9:3:69,69:-3:9,i,9:6:200];

% balloon SSi
for l = 1:length(i)
figure(1)
plot(xS1,fS1,'.b',xS2,fS2,'.r',xd,fd,'.y',zi,fi,'.g', 'MarkerSize', i(l))
axis([-13,7,-11,7]);
axis off
SiSFrame(l)=getframe;
end

m = length(i);

t = linspace(0,pi/2,75);
i = 10*sin(t); j = i(4:end);

% explode SSi
for l = 1:length(j)
xS1p = xS1+j(l)*xS1r;
fS1p = fS1+j(l)*fS1r;
xS2p = xS2+j(l)*xS2r;
fS2p = fS2+j(l)*fS2r;
xdp = xd+j(l)*xdr;
fdp = fd+j(l)*fdr;
zip = zi+j(l)*zir;
fip = fi+j(l)*fir;

figure(1)
plot(xS1p,fS1p,'.b',xS2p,fS2p,'.r',xdp,fdp,'.y',zip,fip,'.g')
axis([-13,7,-11,7]);
axis off
SiSFrame(l+m) = getframe;
end

n = length(j);

t = linspace(pi/2,3/2*pi,75);
i = 1/2+sin(t)/2;

% reform SSi at -90 deg
theta_0 = -pi/2;

xS1t = rS1.*cos(tS1+theta_0)-2;
fS1t = rS1.*sin(tS1+theta_0)+2;
xS2t = rS2.*cos(tS2+theta_0)-2;
fS2t = rS2.*sin(tS2+theta_0)+2;
xdt = rd.*cos(td+theta_0)-2;
fdt = rd.*sin(td+theta_0)+2;
zit = ri.*cos(ti+theta_0)-2;
fit = ri.*sin(ti+theta_0)+2;

xS1r = xS1p - xS1t;
fS1r = fS1p - fS1t;
xS2r = xS2p - xS2t;
fS2r = fS2p - fS2t;
xdr = xdp - xdt;
fdr = fdp - fdt;
zir = zip - zit;
fir = fip - fit;

for l = 1:length(i)
xS1p = xS1t+i(l)*xS1r;
fS1p = fS1t+i(l)*fS1r;
xS2p = xS2t+i(l)*xS2r;
fS2p = fS2t+i(l)*fS2r;
xdp = xdt+i(l)*xdr;
fdp = fdt+i(l)*fdr;
zip = zit+i(l)*zir;
fip = fit+i(l)*fir;

figure(1)
plot(xS1p,fS1p,'.b',xS2p,fS2p,'.r',xdp,fdp,'.y',zip,fip,'.g')
axis([-13,7,-11,7]);
axis off
SiSFrame(l+m+n) = getframe;
end

o = length(i);

t_max = 12;
zeta = .2;
omeg = 2;

t = linspace(0,t_max);
theta = theta_0 * exp(-zeta*omeg.*t).*cos(omeg*sqrt(1-zeta^2).*t);
theta = [theta_0*ones(1,10), theta];

% make SSi pivot back to upright
for l = 1:length(theta)
xS1p = rS1.*cos(tS1+theta(l))-2;
fS1p = rS1.*sin(tS1+theta(l))+2;
xS2p = rS2.*cos(tS2+theta(l))-2;
fS2p = rS2.*sin(tS2+theta(l))+2;
xdp = rd.*cos(td+theta(l))-2;
fdp = rd.*sin(td+theta(l))+2;
zip = ri.*cos(ti+theta(l))-2;
fip = ri.*sin(ti+theta(l))+2;

figure(1)
plot(xS1p,fS1p,'.b',xS2p,fS2p,'.r',xdp,fdp,'.y',zip,fip,'.g')
axis([-13,7,-11,7]);
axis off
SiSFrame(l+m+n+o) = getframe;
end
%movie(SiSFrame,1,30)
%movie2avi(SiSFrame,'SiSFrame.avi', 'compression', 'none', 'fps', 20)