load('YearlyConfidenceIntervalRothera82u.mat')
load('YearlyConfidenceIntervalRothera85u.mat')
load('YearlyConfidenceIntervalRothera88u.mat')
load('YearlyConfidenceIntervalRothera91u.mat')
load('YearlyConfidenceIntervalRothera94u.mat')
load('YearlyConfidenceIntervalRothera98u.mat')
load('YearlyDataRothera82u.mat')
load('YearlyDataRothera85u.mat')
load('YearlyDataRothera88u.mat')
load('YearlyDataRothera91u.mat')
load('YearlyDataRothera94u.mat')
load('YearlyDataRothera98u.mat')

NormalAvg=zeros(6,365,9);
WeightedAvg=zeros(6,365,9);

HeightData = cat(4,YearlyDataRothera82u,YearlyDataRothera85u,YearlyDataRothera88u,YearlyDataRothera91u,YearlyDataRothera94u,YearlyDataRothera98u);
HeightConfidenceInterval = cat(4,YearlyConfidenceIntervalRothera82u,YearlyConfidenceIntervalRothera85u,YearlyConfidenceIntervalRothera88u,YearlyConfidenceIntervalRothera91u,YearlyConfidenceIntervalRothera94u,YearlyConfidenceIntervalRothera98u);

for i=1:365
    for j=1:9
        for k=1:6
            NormalAvg(k,i,j) = nanmean(HeightData(:,i,j,k));
            WeightedAvg(k,i,j) = nansum(HeightData(:,i,j,k).*(HeightConfidenceInterval(:,i,j,k).^-2))/nansum(HeightConfidenceInterval(:,i,j,k).^-2);
        end
    end
end

time = (1:365);
height = [82 85 88 91 94 98];
divtime = (-180:184);

TripleAvg = zeros(6,1095,5);

for i = 1:6
    for j = 1:5
        TripleAvg(i,:,j) = [WeightedAvg(i,:,j) WeightedAvg(i,:,j) WeightedAvg(i,:,j)]';
    end
end

TripleTime = (1:1095)';
 
%Smooth it and pick up the middle

for i = 1:6
    AvgMean(:,i) = feval(fit(TripleTime,TripleAvg(i,:,1)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg48h(:,i) = feval(fit(TripleTime,TripleAvg(i,:,2)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg24h(:,i) = feval(fit(TripleTime,TripleAvg(i,:,3)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg12h(:,i) = feval(fit(TripleTime,TripleAvg(i,:,4)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg8h(:,i) = feval(fit(TripleTime,TripleAvg(i,:,5)','smoothingspline','SmoothingParam',0.07),[366:730]); 
end

figure(1)
contourf(divtime,height,circshift(Avg48h,180,1)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Zonal wind (m/s) 48h');
hold on
contour(divtime,height,circshift(AvgMean,180,1)',[0 0],'LineColor','k','LineWidth',2);
hold off

figure(2)
contourf(time,height,Avg48h(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Zonal wind (m/s) 48h');
hold on
contour(time,height,AvgMean(:,:)',[0 0],'LineColor','k','LineWidth',2);
hold off

%smooth the wind in z-u plane
for i = 1:365
    Shear(i,:) = differentiate(fit(height',AvgMean(i,:)','smoothingspline','SmoothingParam',0.07),height);
end

figure(3)
contourf(time,height,Shear(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Zonal wind shear(1/s)');

figure(4)
contourf(divtime,height,circshift(Shear(:,:),180,1)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Zonal wind shear(1/s)');

