function string9
% wave equation for driven damped string, trapezoidal method

close all;
% global vect
% vect = img2dat;

strlength = 0.32;       % length of string
dx = 0.01;              % size of x step
n = round(strlength/dx);    % number of x iterations

tfin = 50;              % length of time
dt = 0.01;              % size of t step, must be <= 0.02
m = round(tfin/dt);     % number of t iterations

% coefficients and other constant variables
z = 0.00063;            % linear mass density
k = 52;                 % tension
xf = 1/3*strlength;     % position of finger on string
xb = 1/2*strlength;     % position of bow on string

% damping due to finger
fingvec = zeros(n-1,1); % initialize empty vector for finger damping
fingloc = round(xf/dx); % index for location of damping on string
fingvec(fingloc) = 1;   % zero at all places except fingloc

% bow action
bowvec = zeros(n-1,1);  % initialize empty vector for bow action
bowloc = round(xb/dx);  % index for location of bow on string
bowvec(bowloc) = 1;     % zero at all places except bowloc

% initial conditions, t=0
x = (dx:dx:strlength-dx)'; % cut off ends of string for matching dimensions
V = zeros(2*(n-1),m);
V(1:n-1,1) = f(x);

Q = ones(n-1,3);
Q(:,2) = -2*Q(:,2);
A = sparse(zeros(2*(n-1),2*(n-1)));
A(1:n-1,n:2*(n-1)) = speye(n-1);
A(n:2*(n-1),1:n-1) = spdiags(Q,[-1 0 1],n-1,n-1)*(k/z*dt^2/dx^2);
A0 = A;
A0(n:2*(n-1),n:2*(n-1)) = -spdiags(fingvec,0,n-1,n-1)*c(0)*dt/(dx*z);
B = sparse(zeros(2*(n-1),1));
B0 = B;
B0(n:2*(n-1),1) = b(0,dt)*bowvec*dt^2/(dx*z);

for i=2:m
    t = (i-1)*dt;
    A(n-1+fingloc,n-1+fingloc) = -c(t)*dt/(dx*z);
    B(n-1+bowloc,1) = b(t,dt)*dt^2/(dx*z);
    D = speye(2*(n-1)) - dt/2*A;
    E = speye(2*(n-1)) + dt/2*A0;
    V(:,i) = D\(E*V(:,i-1) + dt/2*(B+B0));
    
    % reset "previous" values for next iteration
    A0 = A;
    B0 = B;
end
V = [zeros(1,m); V(1:n-1,:); zeros(2,m); V(n:2*(n-1),:); zeros(1,m)];
                % add rows of zeros for string's ends

figure(1);
rotate3d;                   % allow 3D plot to rotate
pinc = 100;
position = (0:dx:strlength)'*ones(1,m/pinc);
time = ones(n+1,1)*(0:pinc:m-1)*dt;    % length m/20 for a less dense plot
vplot = V(1:n+1,1:pinc:end);
hu = plot3(position,time,vplot,'Color',1/256*[139 131 120]);
hold on;

fingtime = dt:dt:tfin;
fingpos = xf*ones(1,length(fingtime));
fingplot = V(fingloc,1:length(fingtime));
hf = plot3(fingpos,fingtime,fingplot,'Color',1/256*[238 203 173], ...
                                                       'LineWidth',3);
bowtime = dt:dt:tfin;
bowpos = xb*ones(1,length(bowtime));
bowplot = V(bowloc,1:length(bowtime));
hb = plot3(bowpos,bowtime,bowplot,'Color',1/256*[139 69 19],'LineWidth',2);
grid on;
legend([hu(1) hf(1) hb(1)],'String''s Path','Finger Damping','Bow Action');
title('String''s Path Over Time');
xlabel('position on string (m)');
ylabel('time (sec)');
zlabel('vertical displacement of string (m)');
axis([0 strlength 0 tfin]);

% % sound
% figure(2);
% vsample = V(round((n+1)/2),1:end);
% %vol=1;         % increment to increase volume
% %soundvec=reshape(vol*U',q*w,1);  % reshape U matrix into a vector
% plot(vsample);
% soundsc(vsample,44000);     % play string's sound

% % movie
% figure(3);
% % find minimum and maximum vertical height of string to set axes
% U = V(1:n+1,:)';
% [a1,b1]=size(U);
% maxi=U(1,1);
% mini=maxi;
% for i=1:a1
%     for j=1:b1
%         if U(i,j)>maxi
%             maxi=U(i,j);
%         end
%         if U(i,j)<mini
%             mini=U(i,j);
%         end
%     end
% end
% ymax=maxi;
% ymin=mini;
% 
% numframes=m; % frame for each time increment
% A1=moviein(numframes); % create the movie matrix 
% set(gca,'NextPlot','replacechildren');  % set to replace the previous
% for i=1:m
%     axis([0 strlength ymin ymax]);      % set axis to minimum and maximum
%     plot((0:dx:strlength),U(i,:));   % plot string
%     title(['t = ',num2str(i*dt),' seconds']);
%     xlabel('position on string (m)');
%     ylabel('vertical displacement (m)');
%     A1(i)=getframe;              % capture plot of string
% end
% movie2avi(A1,'string7.avi');  % save movie as .avi file

end

function out = f(x)

%out = sin(pi*x);
%out = 0.01 - 0.04*abs(x-0.25);
out = 0;

end

function out=c(t)

%out = 1;
%out = 0;
out = 2*(exp(-3*t) - exp(-5*t));

end

function out=b(t,dt)

%out = 0;
%out = 1;
out = 4*(exp(-0.5*t) - exp(-2*t));
%out = imrend(t,dt);

end

% function out=imrend(t,dt)
% global vect
%     
% tlen=20;
% numsteps=tlen/dt;
% veclen=length(vect);
% if round(t./tlen*veclen)<1||t>tlen
%     out=0;
% else 
%     out=vect(round(t./tlen*veclen));
% end
% end
% 
% function vect=img2dat
%    
% B=imread('wavdat.bmp');
% B=sparse(1-B);
% A=imrotate(B,180);
% [rows,cols]=size(A);
% vect=zeros(cols,1);
% for i=1:cols
%     tot=sum(A(:,i));
%     if tot~=0
%         vect(i)=mean(find(A(:,i)));
%     end
% 
% end
% vect=fliplr(vect);
% vect=vect/max(vect);
% 
% figure(4);
% plot(vect);
%     
% end