%importData();
%load('helios.mat')
close all;
figure;
startTime = timing.clock(8861);
endTime = timing.clock(8929);
%plot(timing.clock, gps.altitude)
plot(clock, altitude)
hold on
plot(timing.clock(8861), gps.altitude(8861), 'rx')
plot(timing.clock(8929), gps.altitude(8929), 'rx')
title 'Hieght versus Time'
%%
figure; 
plot(clock, speed);
hold on;
plot(timing.clock(8861), gps.speed(8861), 'rx')
plot(timing.clock(8929), gps.speed(8929), 'rx')
title 'Total Speed versus Time'
%Note how it was supposed to take five seconds to open the valve, and it
%only took three!
%% speed vs time
figure; 
len =length(altitude);
%aspeed = 1:len-1;
%for i=1:len-1
%    aspeed(i) = (altitude(i+1)-altitude(i))/delay(i+1);
%end
aspeed = 1:20:len-20;
newClock = 1:20:len-20;
for i=1:20:len-20
    change = delay(i+20)+delay(i+19)+delay(i+18)+delay(i+17)+delay(i+16)+delay(i+15)+delay(i+14)+delay(i+13)+delay(i+12)+delay(i+11)+...
      delay(i+10)+delay(i+9)+delay(i+8)+delay(i+7)+delay(i+6)+delay(i+5)+delay(i+4)+delay(i+3)+delay(i+2)+delay(i+1);
    aspeed((i-1)/20+1) = (altitude(i+20)-altitude(i))/change;
    newClock((i-1)/20+1) = clock(i+20);
end
plot(newClock, aspeed);
hold on;
plot(timing.clock(8861), (gps.altitude(8861)-gps.altitude(8860))/delay(8861), 'rx')
plot(timing.clock(8929), (gps.altitude(8929)-gps.altitude(8928))/delay(8929), 'rx')
title 'Vertical Speed versus Time'
%%
%figure;
%newAltitude = unique(altitude, 'stable');
%aspeed = 1:length(newAltitude);
%for i=1:length(newAltitude)
%    aspeed(i) = (altitude(i+1)-altitude(i))/delay(i+1);
%end
%plot(1:length(newAltitude), aspeed);
%%
figure; hold on;
%plotyy(clock, presIn, clock, presOut)
plot(clock, presIn);
plot(clock, presOut, 'g');
legend Inside Outside
plot(timing.clock(8861), altimeter.presIn(8861), 'ro')
plot(timing.clock(8929), altimeter.presIn(8929), 'ro')
title 'Pressure versus Time'
%%
figure; 
plot(clock, speed);
hold on;
plot(timing.clock(8861), gps.speed(8861), 'rx')
plot(timing.clock(8929), gps.speed(8929), 'rx')
title 'Total Speed versus Time'
%Note how it was supposed to take five seconds to open the valve, and it
%only took three!
%% speed vs time
figure; hold on;
len=length(altitude);
balt = 1:len;
bclock = clock(1:len);
slope1 = (gps.altitude(10484)-gps.altitude(1))/(timing.clock(10484)-timing.clock(1));
slope2 = (gps.altitude(8929)-gps.altitude(8861))/(timing.clock(8929)-timing.clock(8861));
line1 = 1:len; line1(1) = altitude(1);
line2 = 1:len; line2(1) = altitude(1)+2000;
pos = 1;
balt(1)=altitude(1);
for i=2:len
    balt(i) = altitude(i);
    line1(i) = line1(i-1)+slope1*(bclock(i)-bclock(i-1));
    line2(i) = line2(i-1)+slope2*(bclock(i)-bclock(i-1));
    if (altitude(i) ~= altitude(pos))
        pos = i;
        aChange = altitude(i)-altitude(pos);
        tChange = double(0);
        for j=pos+1:i
            tChange = tChange+(bclock(j)-bclock(j-1));
        end
        perChange = aChange/tChange;
        for j=pos+1:i
            balt(j)=balt(j-1)+perChange;
        end
    end
end
plot(bclock, balt, 'b');
plot(bclock(1:7000), line1(1:7000), 'g');
plot(bclock(1:7000), line2(1:7000), 'c');
legend altitude 'average speed' 'speed while open'
hold on;
plot(timing.clock(8861), gps.altitude(8861), 'rx')
plot(timing.clock(8929), gps.altitude(8929), 'rx')
title 'Height Filled In versus Time'
%Now will use this data to make a speed graph
figure; 
len =length(altitude);
cspeed = 1:30:len-30;
newClock = 1:30:len-30;
newAltitude = altitude(1:30:len-30);
for i=1:30:len-30
    change = clock(i+30)-clock(i);
    cspeed((i-1)/30+1) = (altitude(i+30)-altitude(i))/change;
    newClock((i-1)/30+1) = clock(i+30);
end
plot(newClock, cspeed);
hold on;
plot(timing.clock(8861), (gps.altitude(8861)-gps.altitude(8860))/delay(8861), 'rx')
plot(timing.clock(8929), (gps.altitude(8929)-gps.altitude(8928))/delay(8929), 'rx')
title 'Vertical Speed Filled In versus Time'
%Now to plot the speed versus altitude
figure;
plot(newAltitude, cspeed)
hold on;
plot(gps.altitude(8861), (gps.altitude(8861)-gps.altitude(8860))/delay(8861), 'rx')
plot(gps.altitude(8929), (gps.altitude(8929)-gps.altitude(8928))/delay(8929), 'rx')
title 'Speed Filled In versus Altitude'
%%
figure; hold on;
%plot(clock, presIn);
f = fit(timing.clock(5573:10484), altimeter.presIn(5573:10484), 'exp2');
plot(f, clock, presIn, '-', '-');
%legend 'Best Fit' Inside
plot(timing.clock(8861), altimeter.presIn(8861), 'ro')
plot(timing.clock(8929), altimeter.presIn(8929), 'ro')
title 'Pressure Inside versus Time with Best Fit'
%%
figure; hold on;
%plot(clock, presIn);
f = fit(timing.clock(5573:10484), altimeter.presOut(5573:10484), 'exp1');
plot(f, clock, presOut, '-', '-');
f = fit(timing.clock(5573:10484), altimeter.presOut(5573:10484), 'exp2');
plot(f, clock, presOut, 'c-', 'c-');
%legend 'Best Fit' Inside
plot(timing.clock(8861), altimeter.presOut(8861), 'ro')
plot(timing.clock(8929), altimeter.presOut(8929), 'ro')
title 'Pressure Outside versus Time with Best Fit'
%%
figure; hold on;
%plot(clock, presIn);
f = fit(timing.clock(5573:10484), altimeter.presIn(5573:10484), 'exp1');
g = fit(timing.clock(5573:8861), altimeter.presIn(5573:8861), 'exp1');
%h = fit(timing.clock(5573:8861), altimeter.presIn(5573:8861), 'exp1');
h = fit(timing.clock(8929:10484), altimeter.presIn(8929:10484), 'exp1');
plot(g, timing.clock(5573:10484), altimeter.presIn(5573:10484), 'c-', 'c-');
%plot(h, timing.clock(5573:8861), altimeter.presIn(5573:8861), 'c-', 'c-');
plot(h, timing.clock(5573:10484), altimeter.presIn(5573:10484), 'g-', 'g-');
plot(f, clock, presIn, '-');
plot(timing.clock(8861), altimeter.presIn(8861), 'ro')
plot(timing.clock(8929), altimeter.presIn(8929), 'ro')
plot([startTime, startTime], [-50000, 150000], 'k--');
plot([endTime, endTime], [-50000, 150000], 'k--');
axis([0 6500000 -50000 150000]); 
title 'Pressure Inside versus Time with Best Fit Separated'
%The two best fit lines in this graph actually intersect at the beginning
%of the jet stream. The next intersection has little to do with the jet
%stream, but is also at the point where total speed levels off. This is for
%a best fit two term exponential, and may not be currently shown.
%Conclusion, as much as I would like to say that we altered the pressure in
%the balloon, that is not what the data says. The exponential for before
%opening follows the exponential for the whole trip very well. The
%exponential closing has too few data points to be accurate. Note that a
%two term exponential is actually more accurate than a one term, but is
%used only in figure nine, not in this grpah.
%%
figure; hold on;
plot(gps.altitude(5573:10484), altimeter.presIn(5573:10484)/1000, 'b-')
plot(gps.altitude(5573:10484), altimeter.presOut(5573:10484)/1000, 'g-')
hold on;
legend Inside Outside
plot([gps.altitude(8861), gps.altitude(8861)], [0, 120000], 'k--');
plot([gps.altitude(8929), gps.altitude(8929)], [0, 120000], 'k--');
xlabel 'Altitude (m)'; ylabel 'Pressure (kPa)'
axis([0 30000 0 120]); 
title 'Pressure versus Altitude'
%% bno vertical speed very simplified
figure; hold on;
plot([startTime, startTime], [0, 12000000], 'k--');
plot([endTime, endTime], [0, 12000000], 'k--');
vel=5573:10484;
vel(1)=0;
%The sensor rose while shaking but primarily with the x axis being vertical
for i=5574:10484
    vel(i-5572) = vel(i-5573)+abs((imu.accelX(i)-imu.gravityX(i))*(timing.clock(i)-timing.clock(i-1)));
end
plot(timing.clock(5573:10484), vel);
title 'Vertical Speed versus Time BNO Simplified'
%% now plotting movement on all three axes, again this assumes no rotation which is not true
%note that these very simplified models may not be too innaccurate as the
%data shows fairly minimal rotation on the upward flight. Down is much
%worse though.
figure;
subplot(1,3,1); hold on;
plot([startTime, startTime], [0, 1200000], 'k--');
plot([endTime, endTime], [0, 1200000], 'k--');
vel=5573:10484;
vel(1)=0;
%The sensor rose while shaking but primarily with the x axis being vertical
for i=5574:10484
    vel(i-5572) = vel(i-5573)+(imu.accelX(i)-imu.gravityX(i))*(timing.clock(i)-timing.clock(i-1));
end
plot(timing.clock(5573:10484), vel);
title 'velocityX'
subplot(1,3,2); hold on; title velocityY
plot([startTime, startTime], [0, 200000], 'k--');
plot([endTime, endTime], [0, 200000], 'k--');
vely=5573:10484;
vely(1)=0;
%The sensor rose while shaking but primarily with the x axis being vertical
for i=5574:10484
    vely(i-5572) = vely(i-5573)+(imu.accelY(i)-imu.gravityY(i))*(timing.clock(i)-timing.clock(i-1));
end
plot(timing.clock(5573:10484), vely);
subplot(1,3,3); hold on; title velocityZ
plot([gps.altitude(8861), gps.altitude(8861)], [0, 200000], 'k--');
plot([gps.altitude(8929), gps.altitude(8929)], [0, 200000], 'k--');
velz=5573:10484;
velz(1)=0;
%The sensor rose while shaking but primarily with the x axis being vertical
for i=5574:10484
    velz(i-5572) = velz(i-5573)+(imu.accelZ(i)-imu.gravityZ(i))*(timing.clock(i)-timing.clock(i-1));
end
plot(timing.clock(5573:10484), velz);
%%
figure; hold on;
newTempIn = tempIn;
newTempOut = tempOut;
diffIn = (tempIn(10)-850)*ones(length(tempIn), 1);
diffOut = (tempOut(10)-850)*ones(length(tempOut), 1);
len = length(tempIn);
for i=2:len
    if (abs(tempIn(i)-tempIn(i-1)) > 200)
        diff=0;
    else
        diff=tempIn(i)-tempIn(i-1);
    end
    newTempIn(i) = newTempIn(i-1)+diff;
    
    if (abs(tempOut(i)-tempOut(i-1)) > 200)
        diff=0;
    else
        diff=tempOut(i)-tempOut(i-1);
    end
    newTempOut(i) = newTempOut(i-1)+diff;
end
newTempIn = newTempIn - diffIn;
newTempOut = newTempOut - diffOut;
plot(clock, newTempIn./100, 'b');
plot(clock, newTempOut./100, 'g');
%Now comparing to raw data
dT = rawD2In-31846*256;
tempRawIn=2000+(dT.*27581)./(2^23);
plot(clock, tempRawIn./100, 'c');
dT = rawD2Out-32593*256;
tempRawOut=2000+(dT.*27590)./(2^23);
plot(clock, tempRawOut./100, 'k');
plot(clock, tempBNO, 'm');
legend Inside Outside BNO RawIn RawOut
title 'Temperature versus Time'
axis([0 9000000 -50 50])
%Conclusion - Temperature inside the balloon is much more stable, and the
%box is much warmer. On the descent, the temperature inside became much
%closer to that inside the box because of the latex cover. However, this
%likely added little protection once the parachute opened.
%The method of differences actually worked rather well given the comparison
%to raw data. I've left both there, but I suggest that actual diagrams use
%only the raw
%% temperature versus altitude
figure; hold on;
newTempBNO = tempBNO - (tempBNO(10)-8.50)*ones(length(tempBNO), 1);
plot(altitude(1:(10484-5573)), tempRawIn(1:(10484-5573))./100, 'b');
plot(altitude(1:(10484-5573)), tempRawOut(1:(10484-5573))./100, 'g');
plot(altitude(1:(10484-5573)), newTempBNO(1:(10484-5573)), 'm');
legend Inside Outside Payload;
xlabel 'Altitude (m)'; ylabel 'Temperature (C)' 
title 'Temperature versus Altitude on Ascent'
%% altitude vs time
figure; hold on;
plot(timing.clock(5573:10484), gps.altitude(5573:10484), 'b.');
g = fit(timing.clock(5573:8861), gps.altitude(5573:8861), 'lin');
h = fit(timing.clock(8929:10484), gps.altitude(8929:10484), 'poly1');
plot(g, timing.clock(5573:10484), gps.altitude(5573:10484), 'b-', 'b-');
plot(h, timing.clock(5573:10484), gps.altitude(5573:10484), 'k-', 'k-');
title 'Altitude versus Time'


%Discoveries - The total speed starts vastly decreasing right below the jet
%stream, as expected. The speed starts vastly increasing actually before we
%reach the jet stream, at 6000m instead of 8000m. As soon as the speed
%starts vastly increasing, the pressure readings in the balloon become much
%more anomalous.