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

164 views (last 30 days)
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
  10 Comments

Sign in to comment.

Answers (1)

Meg Noah
Meg Noah on 1 Jan 2022
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
  1 Comment
Meg Noah
Meg Noah on 7 Jan 2022
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.

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!