clear 
close all

%% Define Colormap MPG
m = 300; % steps in RGB colors
mpg = zeros(m , 3);
T = [0,     0,   0        %// black
     0,    85,  85        %// MP dark green
     199, 212,  40        %// MP light green
     255, 255, 255]./255; %// white
x = [0
     100
     200
     300];
mpg = interp1(x/255,T,linspace(0,1,255));
clear m x T

%% General folder
folder = ['Z:\data_analysis\231212_substrate\OpenSourceData\']; %path to global data folder
freq = [600:2:1050]; %measured frequency range [lowest : stepsize : highest]
nfreq = length(freq);
pdate = '231004'; %measurement date of p-polarized data
sdate = '231212'; %measurement date of s-polarized data
cos = 1; %cosmic ray correction, script will be faster without
% Note: data is already averaged and background corrected

%% Import Raw Data (FEL in p)
for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISpIRp_cx']
    filefolder = [folder,name,'\'];
    end
    filename = [pdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.PPx(:,:,kfreq) = image_out;
end

for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISpIRp_cy']
    filefolder = [folder,name,'\'];
    end
    filename = [pdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.PPy(:,:,kfreq) = image_out;
end

for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISsIRp_cx']
    filefolder = [folder,name,'\'];
    end
    filename = [pdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.SPx(:,:,kfreq) = image_out;
end

for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISsIRp_cy']
    filefolder = [folder,name,'\'];
    end
    filename = [pdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.SPy(:,:,kfreq) = image_out;
end

%% Import Raw Data (FEL in s)
for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISpIRs_cx']
    filefolder = [folder,name,'\'];
    end
    filename = [sdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.PSx(:,:,kfreq) = image_out;
end

for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISpIRs_cy']
    filefolder = [folder,name,'\'];
    end
    filename = [sdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.PSy(:,:,kfreq) = image_out;
end

for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISsIRs_cx']
    filefolder = [folder,name,'\'];
    end
    filename = [sdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.SSx(:,:,kfreq) = image_out;
end

for kfreq = 1:nfreq
    if kfreq == 1
    name = ['VISsIRs_cy']
    filefolder = [folder,name,'\'];
    end
    filename = [sdate,'_SFGmic_ALN_substrate_',name,'_',num2str(kfreq-1,'%.3d'),'.mat'];
    load([filefolder,filename]);        
    raw.SSy(:,:,kfreq) = image_out;
end

clear image_out name filename filefolder

%% Cosmic Ray Correction and Integration
if cos
    
% PPx
for kfreq = 1:nfreq;
    datamax.PPx = max(max(raw.PPx(:,:,kfreq))); %Calculate maximum
    RGBdata.PPx = raw.PPx(:,:,kfreq)./datamax.PPx;  %Normalize
    rawdata.PPx = im2uint8(RGBdata.PPx); %convert to uint8
    [rows columns numberOfColorBands] = size(rawdata.PPx); % Get the dimensions of the image.  numberOfColorBands should be = 1. 
    medianFiltered.PPx = medfilt2(rawdata.PPx, [4 4]); % Median Filtering
    noise.PPx = (rawdata.PPx < 0 | rawdata.PPx >= (180/2+mean(rawdata.PPx,'all'))); %Find the noise
    NoiseFreeData.PPx = rawdata.PPx; %Inititalize
    NoiseFreeData.PPx(noise.PPx) = medianFiltered.PPx(noise.PPx);   % Replace with median
    doubledata.PPx = im2double(NoiseFreeData.PPx);  % Reconvert to double
    cosmic.PPx(:,:,kfreq)= doubledata.PPx(:,:)*datamax.PPx; % rescale to previous scaling
    display([kfreq, sum(sum(noise.PPx))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered
end
raw = rmfield(raw, 'PPx') %note: has to be removed because raw and cosmic struct together take too much memory
clear noise

% PPy
for kfreq = 1:nfreq;
    datamax.PPy = max(max(raw.PPy(:,:,kfreq)));
    RGBdata.PPy = raw.PPy(:,:,kfreq)./datamax.PPy;
    rawdata.PPy = im2uint8(RGBdata.PPy);
    medianFiltered.PPy = medfilt2(rawdata.PPy, [4 4]);
    noise.PPy = (rawdata.PPy < 0 | rawdata.PPy >= (180/2+mean(rawdata.PPy,'all')));
    NoiseFreeData.PPy = rawdata.PPy;
    NoiseFreeData.PPy(noise.PPy) = medianFiltered.PPy(noise.PPy);
    doubledata.PPy = im2double(NoiseFreeData.PPy);   
    cosmic.PPy(:,:,kfreq)= doubledata.PPy(:,:)*datamax.PPy;
    display([kfreq, sum(sum(noise.PPy))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered
end
raw = rmfield(raw, 'PPy')
clear noise
       
% SPx
for kfreq = 1:nfreq;
    datamax.SPx = max(max(raw.SPx(:,:,kfreq)));
    RGBdata.SPx = raw.SPx(:,:,kfreq)./datamax.SPx;
    rawdata.SPx = im2uint8(RGBdata.SPx);
    medianFiltered.SPx = medfilt2(rawdata.SPx, [4 4]); 
    noise.SPx = (rawdata.SPx < 0 | rawdata.SPx >= (180/2+mean(rawdata.SPx,'all')));    
    NoiseFreeData.SPx = rawdata.SPx;   
    NoiseFreeData.SPx(noise.SPx) = medianFiltered.SPx(noise.SPx);
    doubledata.SPx = im2double(NoiseFreeData.SPx);
    cosmic.SPx(:,:,kfreq)= doubledata.SPx(:,:)*datamax.SPx;
    display([kfreq, sum(sum(noise.SPx))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered
end
raw = rmfield(raw, 'SPx')
clear noise

% SPy
for kfreq = 1:nfreq;
    datamax.SPy = max(max(raw.SPy(:,:,kfreq)));
    RGBdata.SPy = raw.SPy(:,:,kfreq)./datamax.SPy;
    rawdata.SPy = im2uint8(RGBdata.SPy);
    medianFiltered.SPy = medfilt2(rawdata.SPy, [4 4]);
    noise.SPy = (rawdata.SPy < 0 | rawdata.SPy >= (180/2+mean(rawdata.SPy,'all')));
    NoiseFreeData.SPy = rawdata.SPy;
    NoiseFreeData.SPy(noise.SPy) = medianFiltered.SPy(noise.SPy);    
    doubledata.SPy = im2double(NoiseFreeData.SPy);     
    cosmic.SPy(:,:,kfreq)= doubledata.SPy(:,:)*datamax.SPy;
    display([kfreq, sum(sum(noise.SPy))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered
end
raw = rmfield(raw, 'SPy')
clear noise

% PSx
for kfreq = 1:nfreq;    
    datamax.PSx = max(max(raw.PSx(:,:,kfreq)));
    RGBdata.PSx = raw.PSx(:,:,kfreq)./datamax.PSx;
    rawdata.PSx = im2uint8(RGBdata.PSx);
    medianFiltered.PSx = medfilt2(rawdata.PSx, [4 4]); 
    noise.PSx = (rawdata.PSx < 0 | rawdata.PSx >= (180/2+mean(rawdata.PSx,'all')));
    NoiseFreeData.PSx = rawdata.PSx;
    NoiseFreeData.PSx(noise.PSx) = medianFiltered.PSx(noise.PSx);
    doubledata.PSx = im2double(NoiseFreeData.PSx); 
    cosmic.PSx(:,:,kfreq)= doubledata.PSx(:,:)*datamax.PSx;
    display([kfreq, sum(sum(noise.PSx))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered
end
raw = rmfield(raw, 'PSx')
clear noise

% PSy
for kfreq = 1:nfreq; 
    datamax.PSy = max(max(raw.PSy(:,:,kfreq)));
    RGBdata.PSy = raw.PSy(:,:,kfreq)./datamax.PSy;
    rawdata.PSy = im2uint8(RGBdata.PSy);
    medianFiltered.PSy = medfilt2(rawdata.PSy, [4 4]);
    noise.PSy = (rawdata.PSy < 0 | rawdata.PSy >= (180/2+mean(rawdata.PSy,'all')));
    NoiseFreeData.PSy = rawdata.PSy;
    NoiseFreeData.PSy(noise.PSy) = medianFiltered.PSy(noise.PSy);
    doubledata.PSy = im2double(NoiseFreeData.PSy);
    cosmic.PSy(:,:,kfreq)= doubledata.PSy(:,:)*datamax.PSy;
    display([kfreq, sum(sum(noise.PSy))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered
end
raw = rmfield(raw, 'PSy')
clear noise

% SSx
for kfreq = 1:nfreq; 
    datamax.SSx = max(max(raw.SSx(:,:,kfreq)));
    RGBdata.SSx = raw.SSx(:,:,kfreq)./datamax.SSx;
    rawdata.SSx = im2uint8(RGBdata.SSx);
    medianFiltered.SSx = medfilt2(rawdata.SSx, [4 4]);
    noise.SSx = (rawdata.SSx < 0 | rawdata.SSx >= (180/2+mean(rawdata.SSx,'all')));
    NoiseFreeData.SSx = rawdata.SSx;
    NoiseFreeData.SSx(noise.SSx) = medianFiltered.SSx(noise.SSx);
    doubledata.SSx = im2double(NoiseFreeData.SSx);
    cosmic.SSx(:,:,kfreq)= doubledata.SSx(:,:)*datamax.SSx;
    display([kfreq, sum(sum(noise.SSx))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered
end
raw = rmfield(raw, 'SSx')
clear noise

% SSy
for kfreq = 1:nfreq; 
    datamax.SSy = max(max(raw.SSy(:,:,kfreq)));
    RGBdata.SSy = raw.SSy(:,:,kfreq)./datamax.SSy;
    rawdata.SSy = im2uint8(RGBdata.SSy);
    medianFiltered.SSy = medfilt2(rawdata.SSy, [4 4]);
    noise.SSy = (rawdata.SSy < 0 | rawdata.SSy >= (180/2+mean(rawdata.SSy,'all')));
    NoiseFreeData.SSy = rawdata.SSy;
    NoiseFreeData.SSy(noise.SSy) = medianFiltered.SSy(noise.SSy);
    doubledata.SSy = im2double(NoiseFreeData.SSy);
    cosmic.SSy(:,:,kfreq)= doubledata.SSy(:,:)*datamax.SSy;
    display([kfreq, sum(sum(noise.SSy))])
    clear datamax doubledata NoiseFreeData rawdata RGBdata medianFiltered   
end
raw = rmfield(raw, 'SSy')
clear noise rows columns numberOfColorBands

pxc = length(1:1024)^2;
int.PPx =  reshape(sum(sum(cosmic.PPx(:,:,:))),1,226)/pxc;
int.PPy =  reshape(sum(sum(cosmic.PPy(:,:,:))),1,226)/pxc;
int.SPx =  reshape(sum(sum(cosmic.SPx(:,:,:))),1,226)/pxc;
int.SPy =  reshape(sum(sum(cosmic.SPy(:,:,:))),1,226)/pxc;
int.PSx =  reshape(sum(sum(cosmic.PSx(:,:,:))),1,226)/pxc;
int.PSy =  reshape(sum(sum(cosmic.PSy(:,:,:))),1,226)/pxc;
int.SSx =  reshape(sum(sum(cosmic.SSx(:,:,:))),1,226)/pxc;
int.SSy =  reshape(sum(sum(cosmic.SSy(:,:,:))),1,226)/pxc;  

else
    
% Calculate Integrals without Cosmic Ray
pxc = length(1:1024)^2;
int.PPx =  reshape(sum(sum(raw.PPx(:,:,:))),1,226)/pxc;
int.PPy =  reshape(sum(sum(raw.PPy(:,:,:))),1,226)/pxc;
int.SPx =  reshape(sum(sum(raw.SPx(:,:,:))),1,226)/pxc;
int.SPy =  reshape(sum(sum(raw.SPy(:,:,:))),1,226)/pxc;
int.PSx =  reshape(sum(sum(raw.PSx(:,:,:))),1,226)/pxc;
int.PSy =  reshape(sum(sum(raw.PSy(:,:,:))),1,226)/pxc;
int.SSx =  reshape(sum(sum(raw.SSx(:,:,:))),1,226)/pxc;
int.SSy =  reshape(sum(sum(raw.SSy(:,:,:))),1,226)/pxc;    
end

%% Import wavenumber
w1 = load([folder,'VISpIRp_cx\AlN_substrate_VISpIRp_cx_wn.txt']);
w.PPx = w1(1,:);
clear w1
  
w1 = load([folder,'VISpIRp_cy\AlN_substrate_VISpIRp_cy_wn.txt']);
w.PPy = w1(1,:);
clear w1

w1 = load([folder,'VISsIRp_cx\AlN_substrate_VISsIRp_cx_wn.txt']);
w1(1,127) = 0.5*w1(1,126)+0.5*w1(1,128);
w.SPx = filloutliers(w1(1,:),'center','movmedian', 3);
clear w1

w1 = load([folder,'VISsIRp_cy\AlN_substrate_VISsIRp_cy_wn.txt']);
w.SPy = w1(1,:);
clear w1

w1 = load([folder,'VISpIRs_cx\AlN_substrate_VISpIRs_cx_wn.txt']);
w.PSx = w1(1,:);
clear w1
  
w1 = load([folder,'VISpIRs_cy\AlN_VISpIRs_cy_wn.txt']);
w1(1,121) = 0.6*w1(1,120)+0.4*w1(1,123); %two outliers --> remove manually
w1(1,122) = 0.4*w1(1,120)+0.6*w1(1,123);
w.PSy = w1(1,:);
clear w1

w1 = load([folder,'VISsIRs_cx\AlN_substrate_VISsIRs_cx_wn.txt']);
w.SSx = w1(1,:);
clear w1
  
w1 = load([folder,'VISsIRs_cy\AlN_substrate_VISsIRs_cy_wn.txt']);
w.SSy = w1(1,:);
clear w1

%% Calculate Interpolation
intp.PPx = interp1(w.PPx,int.PPx,freq,'pchip','extrap');
intp.PPy = interp1(w.PPy,int.PPy,freq,'pchip','extrap');
intp.SPx = interp1(w.SPx,int.SPx,freq,'pchip','extrap');
intp.SPy = interp1(w.SPy,int.SPy,freq,'pchip','extrap');
intp.PSx = interp1(w.PSx,int.PSx,freq,'pchip','extrap');
intp.PSy = interp1(w.PSy,int.PSy,freq,'pchip','extrap');
intp.SSx = interp1(w.SSx,int.SSx,freq,'pchip','extrap');
intp.SSy = interp1(w.SSy,int.SSy,freq,'pchip','extrap');

%% Filter Correction
filter = importdata([folder,'FilterAllInterpolated.txt']);
f45 = permute(filter.data(201:426,7),[2 1])/max(permute(filter.data(201:426,7),[2 1])); %SPy, PPy
f44 = permute(filter.data(201:426,6),[2 1])/max(permute(filter.data(201:426,6),[2 1])); 
f445 = (f44+f45)/2; %all the others

fi.PPx = intp.PPx./(f445);
fi.PPy = intp.PPy./f45;
fi.SPx = intp.SPx./f445;
fi.SPy = intp.SPy./f45;
fi.PSx = intp.PSx./f445;
fi.PSy = intp.PSy./f445;
fi.SSx = intp.SSx./f445;
fi.SSy = intp.SSy./f445;

clear filter f45 f44 f445 intp_PPx

%% Air absorption Correction
Air = importdata('IR_transmission_in_air.dat');
L = permute(Air.data(11050:24050,1),[2 1]); %wavelength in um (equidisant)
TL = permute(Air.data(11000:24000,3),[2 1]); %Transmission for wavelength

%Convert to correct length
wind = 0.005*(10000/mean(freq))*2;
Abs = log(TL);
Air9 = exp(Abs*9.4); %change to 9.4 meter
Tconv9 = conv_alex(Air9,L,wind); %transmission convoluted

%Convert to wavenumber and interpolate
w = 10000./(L(end:-1:1)); %coonvert to wavenumber
T9 = Tconv9(end:-1:1); %coonvert to wavenumber
AirA9off = interp1(w,T9,freq); %interpoaltion to freq
AirA9 = AirA9off + (1-max(AirA9off)); %add offset

%Actual correction
fA.PPx = fi.PPx./AirA9;
fA.PPy = fi.PPy./AirA9;
fA.SPx = fi.SPx./AirA9;
fA.SPy = fi.SPy./AirA9;
fA.PSx = fi.PSx./AirA9;
fA.PSy = fi.PSy./AirA9;
fA.SSx = fi.SSx./AirA9;
fA.SSy = fi.SSy./AirA9;

clear Air L TL wind Abs Air9 Air4 Tconv4 Tconv9 w T9 T4 AirA4off AirA9off

%% Plot Final Spectra
close all
f = figure(1)

%Paper format
    set(gcf, 'color', 'white');
    set(gcf, 'PaperUnits', 'centimeters');
    set(gcf, 'PaperType', 'A4');
    size = get(gcf,'PaperSize');
    width = 800;     % Initialize a variable for width.
    height = 800;      % Initialize a varible for height.
    left = 500;
    bottom = 150;
    f.Position = [left, bottom, width, height];
    set(0, 'DefaultAxesLineWidth', 1)
    
%subplot labels
    A = uicontrol('style','text')
    set(A,'String','(a) VISpIRp')
    set(A,'Position',[220,750,80,18])%[xpos,ypos,xsize,ysize]
    set(A, 'BackgroundColor', 'white');
    set(A, 'FontSize', 11);
    set(A, 'fontweight', 'bold'); 
    B = uicontrol('style','text')
    set(B,'String','(b) VISsIRp')
    set(B,'Position',[550,750,80,18])%[xpos,ypos,xsize,ysize]
    set(B, 'BackgroundColor', 'white');
    set(B, 'FontSize', 11);
    set(B, 'fontweight', 'bold'); 
    C = uicontrol('style','text')
    set(C,'String','(c) VISpIRs')
    set(C,'Position',[220,380,80,18])%[xpos,ypos,xsize,ysize]
    set(C, 'BackgroundColor', 'white');
    set(C, 'FontSize', 11);
    set(C, 'fontweight', 'bold'); 
    D = uicontrol('style','text')
    set(D,'String','(d) VISsIRs')
    set(D,'Position',[550,380,80,18])%[xpos,ypos,xsize,ysize]
    set(D, 'BackgroundColor', 'white');
    set(D, 'FontSize', 11);
    set(D, 'fontweight', 'bold'); 

%subplot 1
    sp1=subplot(2,2,1) 
    hold on
    scatter(freq, fA.PPx, 'MarkerFaceColor',mpg(1,:),'MarkerEdgeColor','none','DisplayName','c || x')
    scatter(freq, fA.PPy, 'MarkerFaceColor',mpg(120,:),'MarkerEdgeColor','none','DisplayName','c || y')
    box on
    set(gca,'XTickLabel',[]);
    ylabel('Average Intensity per Pixel')
    ylim([0 330])
    xlim([600 1050])

%subplot 2    
    sp2=subplot(2,2,2) 
    hold on
    scatter(freq, fA.SPx, 'MarkerFaceColor',mpg(1,:),'MarkerEdgeColor','none','DisplayName','c || x')
    scatter(freq, fA.SPy, 'MarkerFaceColor',mpg(120,:),'MarkerEdgeColor','none','DisplayName','c || y')
    box on
    set(gca,'XTickLabel',[]);
    set(gca,'YTickLabel',[]);
    ylim([0 330])
    xlim([600 1050])
    
    
%subplot 3
    sp3=subplot(2,2,3) 
    hold on
    scatter(freq, fA.PSx, 'MarkerFaceColor',mpg(1,:),'MarkerEdgeColor','none','DisplayName','c || x')
    scatter(freq, fA.PSy, 'MarkerFaceColor',mpg(120,:),'MarkerEdgeColor','none','DisplayName','c || y')
    box on
    ylim([0 2200])
    legend('location','west')
    xlabel('FEL frequency [cm-1]')
    ylabel('Average Intensity per Pixel')
    xlim([600 1050])
     
%subplot 4
    sp4=subplot(2,2,4) 
    hold on
    scatter(freq, fA.SSx, 'MarkerFaceColor',mpg(1,:),'MarkerEdgeColor','none','DisplayName','c || x')
    scatter(freq, fA.SSy, 'MarkerFaceColor',mpg(120,:),'MarkerEdgeColor','none','DisplayName','c || y')
    box on
    set(gca,'YTickLabel',[]);
    ylim([0 2200])
    xlabel('FEL frequency [cm-1]')
    xlim([600 1050])
       
%Subplot positions
    pos1 = get(sp1, 'Position'); 
    new_pos1 = pos1 + [-0.055 -0.066 0.12 0.12];
    set(sp1, 'Position',new_pos1)
    pos2 = get(sp2, 'Position');
    new_pos2 = pos2 + [-0.043 -0.066 0.12 0.12];
    set(sp2, 'Position',new_pos2)
    pos3 = get(sp3, 'Position');
    new_pos3 = pos3 + [-0.055 -0.052 0.12 0.12];
    set(sp3, 'Position',new_pos3)
    pos4 = get(sp4, 'Position');
    new_pos4 = pos4 + [-0.043 -0.052 0.12 0.12];
    set(sp4, 'Position',new_pos4)
  
%Clear variables
clear new_pos1 new_pos2 new_pos3 new_pos4 pos1 pos2 pos3 pos4 sp1 sp2 sp3 sp4 A B C D f bottom height width left size

%% Plot Images
f = figure(2)

%Paper format
    set(gcf, 'color', 'white');
    set(gcf, 'PaperUnits', 'centimeters');
    set(gcf, 'PaperType', 'A4');
    size = get(gcf,'PaperSize');
    width = 800;     % Initialize a variable for width.
    height = 400;      % Initialize a varible for height.
    left = 500;
    bottom = 150;
    f.Position = [left, bottom, width, height];
    set(0, 'DefaultAxesLineWidth', 1)
    colormap(mpg)
    
%subplot labels
    E = uicontrol('style','text')
    set(E,'String','(e) 948 cm-1')
    set(E,'Position',[130,370,150,18])%[xpos,ypos,xsize,ysize]
    set(E, 'BackgroundColor', 'white');
    set(E, 'FontSize', 11);
    set(E, 'fontweight', 'bold'); 
    F = uicontrol('style','text')
    set(F,'String','(f) 810 cm-1')
    set(F,'Position',[490,370,150,18])%[xpos,ypos,xsize,ysize]
    set(F, 'BackgroundColor', 'white');
    set(F, 'FontSize', 11);
    set(F, 'fontweight', 'bold'); 

%subplot 1
    sp1=subplot(1,2,1) 
    imagesc(cosmic.PSx(:,:,175))
    axis image
    
%subplot 2    
    sp2=subplot(1,2,2) 
    imagesc(cosmic.PSx(:,:,106))
    axis image
    colorbar
    
%Subplot positions
    pos1 = get(sp1, 'Position'); 
    new_pos1 = pos1 + [-0.08 -0.05 0.06 0.06];
    set(sp1, 'Position',new_pos1)
    pos2 = get(sp2, 'Position');
    new_pos2 = pos2 + [-0.06 -0.05 0.06 0.06];
    set(sp2, 'Position',new_pos2)
  
%Clear variables
clear new_pos1 new_pos2 new_pos3 new_pos4 pos1 pos2 pos3 pos4 sp1 sp2 sp3 sp4 A B C D f bottom height width left size