% getting the climatology of non SSW years and climatology of the anomaly
% of the non SSW years.

load('YearlyConfidenceIntervalRothera82v.mat')
load('YearlyConfidenceIntervalRothera85v.mat')
load('YearlyConfidenceIntervalRothera88v.mat')
load('YearlyConfidenceIntervalRothera91v.mat')
load('YearlyConfidenceIntervalRothera94v.mat')
load('YearlyConfidenceIntervalRothera98v.mat')
load('YearlyDataRothera82v.mat')
load('YearlyDataRothera85v.mat')
load('YearlyDataRothera88v.mat')
load('YearlyDataRothera91v.mat')
load('YearlyDataRothera94v.mat')
load('YearlyDataRothera98v.mat')

nonSSWyearlyCI82 = YearlyConfidenceIntervalRothera82v;
nonSSWyearlyCI85 = YearlyConfidenceIntervalRothera85v;
nonSSWyearlyCI88 = YearlyConfidenceIntervalRothera88v;
nonSSWyearlyCI91 = YearlyConfidenceIntervalRothera91v;
nonSSWyearlyCI94 = YearlyConfidenceIntervalRothera94v;
nonSSWyearlyCI98 = YearlyConfidenceIntervalRothera98v;
nonSSWdata82 = YearlyDataRothera82v;
nonSSWdata85 = YearlyDataRothera85v;
nonSSWdata88 = YearlyDataRothera88v;
nonSSWdata91 = YearlyDataRothera91v;
nonSSWdata94 = YearlyDataRothera94v;
nonSSWdata98 = YearlyDataRothera98v;
nonSSWyearlyCI82([2,5,6,8,9,12],:,:) = [];
nonSSWyearlyCI85([2,5,6,8,9,12],:,:) = [];
nonSSWyearlyCI88([2,5,6,8,9,12],:,:) = [];
nonSSWyearlyCI91([2,5,6,8,9,12],:,:) = [];
nonSSWyearlyCI94([2,5,6,8,9,12],:,:) = [];
nonSSWyearlyCI98([2,5,6,8,9,12],:,:) = [];
nonSSWdata82([2,5,6,8,9,12],:,:) = [];
nonSSWdata85([2,5,6,8,9,12],:,:) = [];
nonSSWdata88([2,5,6,8,9,12],:,:) = [];
nonSSWdata91([2,5,6,8,9,12],:,:) = [];
nonSSWdata94([2,5,6,8,9,12],:,:) = [];
nonSSWdata98([2,5,6,8,9,12],:,:) = [];

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

HeightDatanonSSW = cat(4,nonSSWdata82,nonSSWdata85,nonSSWdata88,nonSSWdata91,nonSSWdata94,nonSSWdata98);
HeightConfidenceIntervalnonSSW = cat(4,nonSSWyearlyCI82,nonSSWyearlyCI85,nonSSWyearlyCI88,nonSSWyearlyCI91,nonSSWyearlyCI94,nonSSWyearlyCI98);
HeightData = cat(4,YearlyDataRothera82v,YearlyDataRothera85v,YearlyDataRothera88v,YearlyDataRothera91v,YearlyDataRothera94v,YearlyDataRothera98v);
HeightConfidenceInterval = cat(4,YearlyConfidenceIntervalRothera82v,YearlyConfidenceIntervalRothera85v,YearlyConfidenceIntervalRothera88v,YearlyConfidenceIntervalRothera91v,YearlyConfidenceIntervalRothera94v,YearlyConfidenceIntervalRothera98v);

for i=1:365
    for j=1:9
        for k=1:6
            NormalAvgnonSSW(k,i,j) = nanmean(HeightDatanonSSW(:,i,j,k));
            WeightedAvgnonSSW(k,i,j) = nansum(HeightDatanonSSW(:,i,j,k).*(HeightConfidenceIntervalnonSSW(:,i,j,k).^-2))/nansum(HeightConfidenceIntervalnonSSW(:,i,j,k).^-2);
        end
    end
end
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];

TripleAvgnonSSW = zeros(6,1095,5);
TripleAvg = zeros(6,1095,5);

for i = 1:6
    for j = 1:5
        TripleAvgnonSSW(i,:,j) = [WeightedAvgnonSSW(i,:,j) WeightedAvgnonSSW(i,:,j) WeightedAvgnonSSW(i,:,j)]';
    end
end
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
for i = 1:6
    AvgMeannonSSW(:,i) = feval(fit(TripleTime,TripleAvgnonSSW(i,:,1)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg48hnonSSW(:,i) = feval(fit(TripleTime,TripleAvgnonSSW(i,:,2)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg24hnonSSW(:,i) = feval(fit(TripleTime,TripleAvgnonSSW(i,:,3)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg12hnonSSW(:,i) = feval(fit(TripleTime,TripleAvgnonSSW(i,:,4)','smoothingspline','SmoothingParam',0.07),[366:730]);
    Avg8hnonSSW(:,i) = feval(fit(TripleTime,TripleAvgnonSSW(i,:,5)','smoothingspline','SmoothingParam',0.07),[366:730]); 
end
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

%make the average data to Rothera non SSW data years (9 years)
for i = 1:6
    LongAvgMeannonSSW(:,i) = repmat(AvgMeannonSSW(:,i),9,1);
    LongAvg48hnonSSW(:,i) = repmat(Avg48hnonSSW(:,i),9,1);
    LongAvg24hnonSSW(:,i) = repmat(Avg24hnonSSW(:,i),9,1);
    LongAvg12hnonSSW(:,i) = repmat(Avg12hnonSSW(:,i),9,1);
    LongAvg8hnonSSW(:,i) = repmat(Avg8hnonSSW(:,i),9,1);
end
for i = 1:6
    LongAvgMean(:,i) = repmat(AvgMean(:,i),9,1);
    LongAvg48h(:,i) = repmat(Avg48h(:,i),9,1);
    LongAvg24h(:,i) = repmat(Avg24h(:,i),9,1);
    LongAvg12h(:,i) = repmat(Avg12h(:,i),9,1);
    LongAvg8h(:,i) = repmat(Avg8h(:,i),9,1);
end

%Make the original data into a long vector 
for i = 1:6
    OrgMean(:,:,i) = HeightDatanonSSW(:,:,1,i);
    Org48h(:,:,i) = HeightDatanonSSW(:,:,2,i);
    Org24h(:,:,i) = HeightDatanonSSW(:,:,3,i);
    Org12h(:,:,i) = HeightDatanonSSW(:,:,4,i);
    Org8h(:,:,i) = HeightDatanonSSW(:,:,5,i);
    LongOrgMean(:,i) = reshape(OrgMean(:,:,i)',[3285,1]);
    LongOrg48h(:,i) = reshape(Org48h(:,:,i)',[3285,1]);
    LongOrg24h(:,i) = reshape(Org24h(:,:,i)',[3285,1]);
    LongOrg12h(:,i) = reshape(Org12h(:,:,i)',[3285,1]);
    LongOrg8h(:,i) = reshape(Org8h(:,:,i)',[3285,1]);
end

%substract the smoothed average from the original data and get the anomaly.
for i = 1:6
    AnomalyMean(:,i) = LongOrgMean(:,i) - LongAvgMean(:,i);
    Anomaly48h(:,i) = LongOrg48h(:,i) - LongAvg48h(:,i);
    Anomaly24h(:,i) = LongOrg24h(:,i) - LongAvg24h(:,i);
    Anomaly12h(:,i) = LongOrg12h(:,i) - LongAvg12h(:,i);
    Anomaly8h(:,i) = LongOrg8h(:,i) - LongAvg8h(:,i);
end

AnomalyMean2 = reshape(AnomalyMean, 365, 9, 6);
Anomaly48h2 = reshape(Anomaly48h, 365, 9, 6);
Anomaly24h2 = reshape(Anomaly24h, 365, 9, 6);
Anomaly12h2 = reshape(Anomaly12h, 365, 9, 6);
Anomaly8h2 = reshape(Anomaly8h, 365, 9, 6);

AnomalyMean3 = permute(AnomalyMean2,[2 1 3]);
Anomaly48h3 = permute(Anomaly48h2,[2 1 3]);
Anomaly24h3 = permute(Anomaly24h2,[2 1 3]);
Anomaly12h3 = permute(Anomaly12h2,[2 1 3]);
Anomaly8h3 = permute(Anomaly8h2,[2 1 3]);

AnomalyAvgMean = permute(nanmean(AnomalyMean3),[2 3 1]);
AnomalyAvg48h = permute(nanmean(Anomaly48h3),[2 3 1]);
AnomalyAvg24h = permute(nanmean(Anomaly24h3),[2 3 1]);
AnomalyAvg12h = permute(nanmean(Anomaly12h3),[2 3 1]);
AnomalyAvg8h = permute(nanmean(Anomaly8h3),[2 3 1]);

%plot the Climatology of non-event years

figure(1);
contourf(time,height,AvgMeannonSSW(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) Mean');
title('Climatology of non-SSW years');

figure(2);
contourf(time,height,Avg48hnonSSW(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 48h');
title('Climatology of non-SSW years');

figure(3);
contourf(time,height,Avg24hnonSSW(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 24h');
title('Climatology of non-SSW years');

figure(4);
contourf(time,height,Avg12hnonSSW(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 12h');
title('Climatology of non-SSW years');

figure(5);
contourf(time,height,Avg8hnonSSW(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 8h');
title('Climatology of non-SSW years');

%plot the Climatology of Anomaly of non-event years

figure(6);
contourf(time,height,AnomalyAvgMean(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) Mean');
title('Climatology of Anomaly (non-SSW)');

figure(7);
contourf(time,height,AnomalyAvg48h(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 48h');
title('Climatology of Anomaly (non-SSW)');

figure(8);
contourf(time,height,AnomalyAvg24h(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 24h');
title('Climatology of Anomaly (non-SSW)');

figure(9);
contourf(time,height,AnomalyAvg12h(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 12h');
title('Climatology of Anomaly (non-SSW)');

figure(10);
contourf(time,height,AnomalyAvg8h(:,:)');
xlabel('Day'),ylabel('Altitude(km)');
title(colorbar,'Meridional wind (m/s) 8h');
title('Climatology of Anomaly (non-SSW)');