%Doppler Effect Project Part 3:
%Source from left to right in front of observer
%srclen is the length of the source trajectory in meters
%v is the velocity of the source
%obsx is the x-position along that axis
%obsy is the distance of the source trajectory away from the observer
%f0-f2 are the starting frequencies
%Mike Battaglia

function doppler3(srclen, v, obsx, obsy, f0, f1, f2)
fs = 8000;
bitdepth = 16;
c = 341.17;
totaltime = srclen/v;
t = 1/fs:1/fs:totaltime;
srcx = v*t;
relx = srcx - obsx;
rad = sqrt(relx.^2 +obsy^2);
pressenv = 1./rad;
vrad = (relx*v)./sqrt(relx.^2 + obsy^2);
ff0 = f0*c./(c+vrad);
ff1 = f1*c./(c+vrad);
ff2 = f2*c./(c+vrad);
timeoffset = (obsx/srclen)*totaltime;
o0 = pressenv.*sin(2.*pi.*ff0.*(t-timeoffset));
o1 = pressenv.*sin(2.*pi.*ff1.*(t-timeoffset));
o2 = pressenv.*sin(2.*pi.*ff2.*(t-timeoffset));

output = (o0 + o1 + o2)./3;
output = output/max(output);

subplot(411); plot(t, vrad, 'LineWidth', 1); xlabel('Time'); ylabel('Radial Velocity'); title('Radial Velocity'); axis([min(t) max(t) min(vrad) max(vrad)]);
subplot(412); plot(t, pressenv, 'LineWidth', 1); xlabel('Time'); ylabel('Pressure'); title('Sound Pressure Envelope'); axis([min(t) max(t) min(pressenv) max(pressenv)]);
subplot(413); hold all;
ftot = cat(2, ff0, ff1, ff2);
fmin = min(ftot);
fmax = max(ftot);
plot(t, ff0, 'LineWidth', 2); xlabel('Time'); ylabel('Frequency'); title('Frequencies vs. Time');
plot(t, ff1, 'LineWidth', 2);
plot(t, ff2, 'LineWidth', 2);
axis([min(t) max(t) fmin fmax]);
hold off;
subplot(414); plot(t, output, 'LineWidth', 1); xlabel('Time'); ylabel('Amplitude'); title('Output Waveform'); axis([min(t) max(t) min(output) max(output)]);

sound(output, fs, bitdepth);
wavwrite(0.99*output, fs, bitdepth, 'doppler3.wav');