How to align my plots correctly with the x-axis?

First, I apologize for the title of the question. Essentially, when I plot two series individually (see P7.png and P25.png) and compare them with the dashed lines at 4, 7, 9.75 and 11 degrees (where the peaks should occur) they do pretty well. As soon as I try to plot them on the same graph they seem to be shifted and no longer align (see P7and25.png). This gets worse the more dataset that I have. I have included the data as test.mat. Thank you in advance for your help.

P7.png:

P25.png:

P7and25.png:

load(test.mat)

h = cell(1,length(d));
cmap = jet(m);
figure
%xlim([-5 20])
hold all
 for k=1:m 
    %Apply 50 degree elevation cutoff
    idx2 = find(All_oelev(:,k) > 50 & All_oelev(:,k) < 140); % indices to plot of All_oelev
   
    plot(xi2(k,idx2),VTEC2(idx2,k), 'Color', cmap(k, :))%, '--g','linewidth',2)% Fig5
    sz=size(idx2);
    idx2 = NaN(sz);
    %h{k} =sprintf('aircraft at %3.4f degrees latitude', parnam{1,k}.data(1,1)) ;
    end
     plot([4 4], [0 max(max(AllTEC))],'--k','LineWidth',2);
    plot([7 7], [0 max(max(AllTEC))],'--k','LineWidth',2);
    plot([9.75 9.75], [0 max(max(AllTEC))],'--k','LineWidth',2);
    plot([11 11], [0 max(max(AllTEC))],'--k','LineWidth',2);
%ylim([-0.4 0.8])
xlabel('Latitude (Degree)','fontsize',24, 'fontweight','bold'); 
ylabel('Total Electron Content (TECu = 10^{16} electrons/m^2)','fontsize',24, 'fontweight','bold')
 %legend(h)
 hold off

Best regards,

Alex

[Jakob B. Nielsen]

How do you generate your first two (individual) plots? Looking at your data, the peak that you want to have at 11 degrees is at 11,36 so I don't understand how you manage to make them align in your first plot.

[Alex]

Here is the data from the first two datasets (P7_test.mat and P25_test.mat).

Here is the complete script I am using. The first part is the engine used to open the raw datasets. It should be fine, It just parses the different datasets text files into parameters such as the quantitudy 'AllTEC', elevation angles 'Alloelev' and others that are not used.

clear all
close all
clc
help VTEC_mapping_AIS

Re =6371; %Earth radius (km)
h = 350;% (250 km in Sharma example) height (km) of airglow emission (or thin shell)

% the data
     proot='Plane*';% <- your FOLDER template
     froot='par_o.par'; % <- your FILENAME template
     hroot='Data.txt'; % <- your FILENAME template

% the engine 
     d=dir(proot);
     d={d([d.isdir]).name}.';
     fnam=cellfun(@(x) sprintf('%s%s%s',x,filesep,froot),d,'uni',false);%all data looks good (text and col headers not so much)
     gnam=cellfun(@(x) sprintf('%s%s%s',x,filesep,froot),d,'uni',false);
     hnam=cellfun(@(x) sprintf('%s%s%s',x,filesep,hroot),d,'uni',false);

%%%%variable for number of traces 
rays=importdata(hnam{1});%rays=4; use fixed value for FFTData or rewrite similar script
rays=rays(1,1);

% - processing each file
m=length(d);
for i=1:m
     cfile=fnam{i};
     gfile=gnam{i};
     disp(sprintf('working on file %s',cfile));
     disp(sprintf('working on file %s',gfile));%should be redundant
% code to read contents of CFILE
fnam{i}=importdata(cfile, ' ', 32); %delimeter and headerlines using import wizard 
%                                   to format unknown file(s) *.par
gnam{i}=importdata(cfile, ' ', 57);%171
p_line=rays+53;
parnam{i} =importdata(cfile, ' ', p_line);%p_line should give 113 for 60 rays, 171 for 118 rays

%extract starting latitude, longitude, and altitude to corresponding
%vectors
Latrayorigin(1:rays,i)=parnam{1,i}.data(1,1);%WRONG
Lonrayorigin(1:rays,i)=parnam{1,i}.data(1,2);
Altrayorigin(1:rays,i)=parnam{1,i}.data(1,3);

Latrayend(:,i)=fnam{i,1}.data(:,3);
Lonrayend(:,i)=fnam{i,1}.data(:,2);
Altrayend(:,i)=fnam{i,1}.data(:,4);
%rawTEC{i}=fnam{i}(1:12,:);
end

Alto=0.001*Altrayorigin;%convert origin to km to match endpoint
O(1,:)=reshape(Latrayorigin.',1,numel(Latrayorigin));%x=lat
O(2,:)=reshape(Alto.',1,numel(Alto));%y=alt
O(3,:)=reshape(Lonrayorigin.',1,numel(Lonrayorigin));%z=lon
%get direction unit vectors
numx=Latrayend-Latrayorigin;
numy=Altrayend-Alto;
numz=Lonrayend-Lonrayorigin;
numex=reshape(numx.',1,numel(numx));
numey=reshape(numy.',1,numel(numy));
numez=reshape(numz.',1,numel(numz));
denominator=sqrt(numex.^2 + numey.^2 + numez.^2);%magnitude of direction vector
D(1,:)=numex./denominator;%x=lat
D(2,:)=numey./denominator;%y=alt
D(3,:)=numez./denominator;%z=lon
%direction=   %direction unit vector
for i=1:m
     hfile=hnam{i};
     disp(sprintf('working on file %s',hfile));
% your code to read contents of CFILE
hnam{i}=importdata(hfile);

%satellite latitudes (convergence latittude)
rawsatlat{i}=hnam{i}(2:(rays+1));
ophase{i}=hnam{i}((rays+2):((2*rays)+1));
xphase{i}=hnam{i}(((2*rays)+2):((3*rays)+1));
oaspect{i}=hnam{i}(((3*rays)+2):((4*rays)+1));
xaspect{i}=hnam{i}(((4*rays)+2):((5*rays)+1));
altitude{i}=hnam{i}(((5*rays)+2):((6*rays)+1));
oelev{i}=hnam{i}(((6*rays)+2):((7*rays)+1));
xelev{i}=hnam{i}(((7*rays)+2):((8*rays)+1));
opath{i}=hnam{i}(((8*rays)+2):((9*rays)+1));
xpath{i}=hnam{i}(((9*rays)+2):((10*rays)+1));
otime{i}=hnam{i}(((10*rays)+2):((11*rays)+1));
xtime{i}=hnam{i}(((11*rays)+2):((12*rays)+1));
timediff{i}=hnam{i}(((12*rays)+2):((13*rays)+1));
FRone{i}=hnam{i}(((13*rays)+2):((14*rays)+1));
FRtwo{i}=hnam{i}(((14*rays)+2):((15*rays)+1));
rotrateone{i}=hnam{i}(((15*rays)+2):((16*rays)+1));
rotratetwo{i}=hnam{i}(((16*rays)+2):((17*rays)+1));
maxedens{i}=hnam{i}(((17*rays)+2):((18*rays)+1));
minedens{i}=hnam{i}(((18*rays)+2):((19*rays)+1));
TEC{i}=hnam{i}(((19*rays)+2):((20*rays)+1));
opower{i}=hnam{i}(((20*rays)+2):((21*rays)+1));
xpower{i}=hnam{i}(((21*rays)+2):((22*rays)+1));
orientation{i}=hnam{i}(((22*rays)+2):((23*rays)+1));
ellipticity{i}=hnam{i}(((23*rays)+2):((24*rays)+1));

%Combine TEC data from all iterations into a single array
rowrawTEC{i}=reshape(TEC{i}.',rays,1);
AllTEC(:,i)=rowrawTEC{i};

%Combine FR data from all iterations into a single array
rowrawFR{i}=reshape(FRone{i}.',rays,1);
AllFR(:,i)=rowrawFR{i};


% refine calculation by considering the elvation angle
row_oelev{i}=reshape(oelev{i}.',rays,1);
All_oelev(:,i)=abs(row_oelev{i});

% rawAlt{i}=hnam{i}(7,:);
% %Combine all Altitude data, from all iterations into a single Altitude Data
% %matrix
% rowrawAlt{i}=reshape(rawAlt{i}.',rays,1);
% AllAlt(:,i)=rowrawAlt{i};
end

TECfinal(1,:)=reshape(AllTEC.',1,numel(AllTEC));
FRfinal(1,:)=reshape(AllFR.',1,numel(AllFR));
latrayendfinal(1,:)=reshape(Latrayend.',1,numel(Latrayend));

% Find and remove any TEC ==0
for kk =1:1:m
Tidx = find(AllTEC(:,kk)==0); % indices
AllTEC(Tidx,kk)=NaN;
AllFR(Tidx,kk)=NaN;
All_oelev(Tidx,kk)=NaN;
Latrayend(Tidx,kk)=NaN;
Latrayorigin(Tidx,kk)=NaN;
end

%interpolate for NaN values
xi = linspace(min(Latrayend(:,1)), max(Latrayend(:,1)), numel(Latrayend(:,1))); % Interpolating Vector


 %2nd method of calculating VTEC
 zp =sind(90-All_oelev).*(Re/(Re+h));
 VTEC2 =AllTEC.*sqrt(1-sind(zp).^2);
 
% FIGURE 3
% plot VTEC as a function of satlat with lines of different color and style
figure
hold all
 for k=1:m
         VTEC2(:,k) = interp1(Latrayend(~isnan(VTEC2(:,k)),k), VTEC2(~isnan(VTEC2(:,k)),k), xi, 'linear');
     line(xi,VTEC2(:,k), 'Color', cmap(k, :), 'LineStyle', '--')% USE FOR LARGE NUMBER OF
    
end
     ylabel('Total Electron Content (TECu = 10^{16} electrons/m^2)','fontsize',24, 'fontweight','bold')
    xlabel('Satellite Latitude (Degree)','fontsize',24, 'fontweight','bold'); 
    title('FIG 4- VTEC (first and second method) as a function of satellite latitude') 
 legend('VTEC method 1', 'VTEC method 2')
 hold off


 
 
 
 
 
%Mapping the IPP subionopheric latitude (from VTEC_mapping_example.m)
%ground receiver located at:
gnd_lon = Lonrayorigin;%74.2;% 74.2 E longitude
L = Latrayorigin;%16.8; % gnd_lat 16.8 N latitude
% satellite transmitter located at:
sat_lon = Lonrayend;%71.2;%71.2 E longitude

% %TEST using latitude instead of longitude
for k=1:m
for ijk=1:length(VTEC)
    G(k,ijk) = gnd_lon(ijk,k) - sat_lon(ijk,k); %(3 degrees) difference between satellite longitude and ground receiver longitude
    Az(k,ijk) = 180 + atand(tand(G(k,ijk)./tand(L(ijk,k))));%(189.8 degrees) azimuth angle
    El(k,ijk) = All_oelev(ijk,k);%= atand((cosd(G(ijk)).*cosd(L(ijk)) - 0.1512)./(sqrt(1 - cosd(G(ijk)).^2.*cosd(L(ijk)).^2))) ;% (69.97 degrees) elevation angle
    Chi_pp(k,ijk) = 180/2 - El(k,ijk) - asind((Re/(Re + h)).*cosd(El(k,ijk))); % (0.79 degrees) Earth centred angle between ground and IPP 
    lat_pp(k,ijk) = asind(sind(L(ijk,k)).*cosd(Chi_pp(k,ijk)) + cosd(L(ijk,k)).*sind(Chi_pp(ijk)).*cosd(Az(ijk)));%(16.02 degrees) latitude of IPP
    lon_pp(k,ijk) = gnd_lon(ijk,k) + asind((sind(Chi_pp(k,ijk)).*sind(Az(k,ijk)))./cosd(lat_pp(k,ijk)));% (74.06 degrees) longitude of IPP
end
 %interpolate for NaN values in lat_pp
xi2(k,:) = linspace(min(lat_pp(k,:)), max(lat_pp(k,:)), numel(lat_pp(k,:))); % Interpolating Vector
 end

Previous code to plot the data goes here.