Index exceeds the number of array elements. Index must not exceed 1.

I'm trying to write this program for radialvelocities, but when I try to get my peaks.

This is the code I am running. But it gives me Index exceeds the number of array elements. Index must not exceed 0.

Error in Radialvelocity (line 95)

if pks(1) > pks(2).

I am starting i at 3 since the first two files in my dir are '.' and '..', is this the reason why?

And if so is there an easier way to skip '.' and '..' in the dir?

%Data
files = dir("C:\Users\Jeffs\OneDrive\Documents\MATLAB\KIC4054905");
data_reference = load('5750_45_p00p00_template.dat');


%Test
Test=load(files(3).name);
TestWave=Test(:,1); TestSolarFlux=Test(:,2);%Trial waves
%Plotting
fig1=figure(1);
hold on;
plot(TestWave,TestSolarFlux,'b')
plot(data_reference(:,1),data_reference(:,2),'r')
title('Original data')
xlabel('Wavelength[Aangstrom]')
ylabel('Solar flux')
hold off;

%Start = 5000.35Å; end = 5497.52Å
LambdaInterpolate=lambdavector(5e03, 5.497e3,1)

%New tests
New_test_flux = interp1(Test(:,1),Test(:,2),LambdaInterpolate);
New_comparison = interp1(data_reference(:,1),data_reference(:,2),LambdaInterpolate);

%Plotting
fig2=figure(2);
hold on;
title('Interpoleret Data')
plot(LambdaInterpolate,New_test_flux,'b')
plot(LambdaInterpolate,New_comparison,'r')
xlabel('Wavelength[Aangstrom]')
ylabel('Solar flux')
hold off;

%Crossrelations
[crossrelationTest,lagTest]=xcorr(New_test_flux,New_comparison,200)

%Plotting
fig3 = figure(3);
hold on
title('Crosskorrelation as func of lag/radialvelocity)')
xlabel('Lag/radialvelocity[km/s]')
ylabel('Crosskorrelation')
grid on;
plot(lagTest,crossrelationTest,'r')

%Peaks
[pks locs] = findpeaks(crossrelationTest,lagTest,'MinPeakProminence',4,'Annotate','extents');%Det lader til at virke ved 4.
% plot(-108,1.3497*1.0e+04,'.b','markersize',20)
% plot(1,1.3485*1.0e+04,'.b','markersize',20)

%Constant:
LambdaInterpolate=lambdavector(5334,5.6124e+03,1);
New_comparison=interp1(data_reference(:,1),data_reference(:,2),LambdaInterpolate);

%We need arrays for data
v1=[];
v2=[];
Times=[];
for i=3:15
    
    DataAll=load(files(i).name);
    Time=files(i).datenum;
    Times=[Times Time];
    
    LambdaAll=DataAll(:,1);
    FluxerAll=DataAll(:,2);
    
    
    New_flux=interp1(LambdaAll,FluxerAll,LambdaInterpolate);
    
    %Crossrelations
    [crossrelationDataAll,lagDataAll]=xcorr(New_flux,New_comparison,200,'coef');
    
    %Determining radialvelocities
    
    [pks locs widths proms] = findpeaks(crossrelationDataAll,lagDataAll,'MinPeakProminence',0.004,'Annotate','extents');
    
    fig2=figure(2);
    hold on;
    plot(lagDataAll,crossrelationDataAll,'r')
    
    
    if length(pks) == 1 %Her overlapper hastighedderne, altså de er ens/krydser.
        v1=[v1 locs(1)];
        v2=[v2 locs(1)];
    else
        if pks(1) > pks(2)
            v1=[v1 locs(1)];
            v2=[v2 locs(2)];
        else
            v1=[v1 locs(2)];
            v2=[v2 locs(1)];
        end
    end
    
    
    
    
end

%Function
function vector = lambdavector(lstart,lslut,dv)
c=2.9979e5;%Enhed km/s.
n=log10(lslut/lstart)/(log10(1+dv/c))+1;
vector=[];
for i=1:n
    li=lstart*(1+dv/c).^(i-1);
    vector=[vector,li];
end
vector;
end

[Walter Roberson]

I am starting i at 3 since the first two files in my dir are '.' and '..', is this the reason why?

No, but you should not do that.

files = dir("C:\Users\Jeffs\OneDrive\Documents\MATLAB\KIC4054905");
files([files.isdir]) = [];   %remove folders including . and ..
data_reference = load('5750_45_p00p00_template.dat');
%Test
Test=load(files(1).name);

%[...]
for i=1:13

[Walter Roberson]

[pks locs] = findpeaks(crossrelationTest,lagTest,'MinPeakProminence',4,'Annotate','extents');%Det lader til at virke ved 4.

What if the prominence is not at least 4?

[Jeff Sejr]

It definetly changed my data a lot, so that definetly helped. But even if I do that and change my for loop to i = 1:12, I get the same error: Index exceeds the number of array elements. Index must not exceed 0.

Error in Radialvelocity (line 98)

if pks(1) > pks(2)

[Walter Roberson]

You correlate relative to the first file you read in, which was offset 3 in your case (but above I show how to make it offset 1). But you start from that same file for the correlation, so the data would be identical to your reference. I would be concerned about whether that would throw things off. Perhaps the first data file should be the reference and you loop over the other files ?

[Walter Roberson]

You should not assume you are getting more than one peak back from findpeaks() . You should be checking the length(pks) before you index the output.

[Jeff Sejr]

I'm actuallly not getting any peaks anymore, although, there are peaks in my graph. I tried moving the template file out of the zip file, and ran i = 1:12

[Jeff Sejr]

They seems to disappear when I start my for loop, if I comment out that section I get a list of 19 peaks and 19 locs

[Jeff Sejr]

Might it be due to:

DataAll=load(files(i).name);

Time=files(i).datenum;

[Jan]

This line is a waste of time only:

vector;

You can replace:

vector=[];
for i=1:n
    li=lstart*(1+dv/c).^(i-1);
    vector=[vector,li];
end

by

vector = lstart * (1 + dv / c) .^ (0:n-1);

[Meg Noah]

to remove . and .. see:

https://www.mathworks.com/matlabcentral/answers/750639-how-to-remove-empty-directory-names-such-as-returned-by-dir-function

Answer 1

[Meg Noah]

Hi,

I was trying to figure out what your data are to help, but I couldn't find it online. I did find a spectrum at simbad. Is this the data you are working with?

% general way to read the file information
% download from simbad at 130.79.128.5 then unzip it
% at /0/more/splib120/5750_45_p02p00.ms.fits.gz
info = fitsinfo('5750_45_p02p00.ms.fits');
myTable = cell2table(info.PrimaryData.Keywords);

% general way to read the file information
fitsdisp('5750_45_p02p00.ms.fits','Index',1,'Mode','full')

% general way to read the data
spc = fitsread('5750_45_p02p00.ms.fits');

% reading just the data of interest
import matlab.io.*
fptr = fits.openFile('5750_45_p02p00.ms.fits');
n = fits.getNumHDUs(fptr);
fits.movAbsHDU(fptr,1);
[CRVAL1,comment] = fits.readKeyDbl(fptr,'CRVAL1');
[CDELT1,comment] = fits.readKeyDbl(fptr,'CDELT1');
[NAXIS1,comment] = fits.readKeyDbl(fptr,'NAXIS1');
[TEFF,comment] = fits.readKeyDbl(fptr,'TEFF');
fits.closeFile(fptr);

expected_lambda_of_max_Angstroms = 1e10*(2.897771955e-3/TEFF); 
lambda_Angstroms = CRVAL1:CDELT1:(NAXIS1*CDELT1 + CRVAL1);
lambda_Angstroms = lambda_Angstroms(1:NAXIS1);
figure()
subplot(2,1,1)
plot(1e-4*lambda_Angstroms,spc(1,:))
xlabel('Wavelength [\mum]')
title('KIC4054905');
grid on
subplot(2,1,2)
plot(1e-4*lambda_Angstroms,spc(2,:));
hold on
plot([1e-4*expected_lambda_of_max_Angstroms 1e-4*expected_lambda_of_max_Angstroms], ...
    [0 max(spc(2,:))],'r');
xlabel('Wavelength [\mum]')
grid on

[Meg Noah]

My point is that maybe you're trying to find the peaks that coorespond to a catalog in the wrong spectral reason because of the amount of doppler to shift where you searched may have been incorrect.

Answer 2

[美盈 康]

我出现了一样的问题，然后我发现在执行for循环下的两个if语句时跳过了其中一个，导致数组出现间断。我发现你的代码中也存在for循环下嵌套多个if语句。或许你可以研究一下怎么绕过这种嵌套方法。据说2009以前的matlab版本不会出现这种问题，当然我也没试过。