After running the code, the arrays that should have the calculated values for all the rows are still empty?

2 views (last 30 days)
Arpan Ojha
Arpan Ojha on 29 Nov 2021
Answered: Prachi Kulkarni on 3 Dec 2021
I have this code that I am trying to loop for each row. After running the code, the arrays that should have the calculated values for all the rows are still empty. When I open the arrays from the workplace they are still empty with a value of 0 in one cell.
RA_NGP = deg2rad(15*(12+(51/60)+(26.28/3600)));
DEC_NGP = deg2rad(27+(7/60)+(41.7/3600));
l_NCP = deg2rad(123+(55/60)+(55.2/3600));
data_path = 'Data\';
file_name = 'Arpan_BINOCULAR-20200119-fft-0.csv';
file_path = strcat(data_path,file_name);
data_i = readmatrix(file_path);
rows = size(data_i, 1);
all_rows_longitudes = zeros(length(rows),1);
all_rows_latitudes = zeros(length(rows),1);
all_rows_peak_velocities = zeros(length(rows),1);
all_rows_peak_velocities_errors = zeros(length(rows),1);
all_rows_intensities = zeros(length(rows),1);
for row = 1:rows
data = data_i(rows, :);
RA_H = data(4);
RA_M = data(5);
RA_S = data(6);
RA = deg2rad(15*(RA_H+(RA_M/60)+(RA_S/3600)));
DEC = deg2rad(data(9));
Div1 = cos(DEC)*sin(RA-RA_NGP);
Div2 = cos(DEC_NGP)*sin(DEC)-sin(DEC_NGP)*cos(DEC)*cos(RA-RA_NGP);
l = l_NCP - atan2(Div1, Div2);
longitude = rad2deg(l);
if longitude < 0
longitude = longitude + 360;
end
Add = sin(DEC)*sin(DEC_NGP)+cos(DEC)*cos(DEC_NGP)*cos(RA-RA_NGP);
b = asin(Add);
latitude = rad2deg(b);
if ((longitude >= 0) && (longitude <= 90) || (longitude >=270) && (longitude <= 360)) && ((latitude >= -20) && (latitude <= 20))
x_data = (1:8192);
y_data = data(:, 10:end);
linear_y = 10.^(T./10);
frequency_rest = data(:, 7);
x_frequency = linspace(frequency_rest-1.25, frequency_rest+1.25, 8192);
x_1 = x_frequency(1500);
x_2 = x_frequency(6500);
y_1 = linear_y(1500);
y_2 = linear_y(6500);
m = ((y_1-y_2)/(x_1-x_2));
c = y_1 - m*x_1;
y_fit = m.*x_frequency + c;
y_fitted = linear_y-y_fit;
radial_vel = -((x_frequency-frequency_rest)./frequency_rest).*(3*10^5);
condition = (-150 < radial_vel) & (radial_vel < 150);
radial_vel_focussed = radial_vel(condition);
y_focussed = y_fitted(condition);
y_smooth = smooth(radial_vel_focussed, y_focussed, 0.1, 'rloess');
[radial_vel_focussed, order] = sort(radial_vel_focussed);
y_smooth = y_smooth(order);
[pks, locs, w, p] = findpeaks(y_smooth, radial_vel_focussed, 'MinPeakProminence', 0.02e-10);
peak_velocity = max(locs);
peak_data = [pks, locs', w', p];
if isempty(peak_data)
peak_data = zeros(1,4)
peak_velocity = 0
end
[a b] = find(peak_data == peak_velocity);
rightmost_peak_data = peak_data(a, :);
rightmost_pks = rightmost_peak_data(1);
rightmost_locs = rightmost_peak_data(2);
rightmost_w = rightmost_peak_data(3);
rightmost_p = rightmost_peak_data(4);
all_rows_latitudes(row) = latitude;
all_rows_longitudes(row) = longitude;
all_rows_peak_velocities(row) = peak_velocity;
all_rows_peak_velocities_errors(row) = rightmost_w;
all_rows_intensities(row) = rightmost_w*rightmost_p;
end
end

Answers (1)

Prachi Kulkarni
Prachi Kulkarni on 3 Dec 2021
Hi,
Here is a slight modification of your code that can help.
Hi,
Here is a slight modification of your code that can help.
RA_NGP = deg2rad(15*(12+(51/60)+(26.28/3600)));
DEC_NGP = deg2rad(27+(7/60)+(41.7/3600));
l_NCP = deg2rad(123+(55/60)+(55.2/3600));
data_path = 'Data\';
file_name = 'Arpan_BINOCULAR-20200119-fft-0.csv';
file_path = strcat(data_path,file_name);
data_i = readmatrix(file_path);
rows = size(data_i, 1);
all_rows_longitudes = zeros(rows,1); % Modified
all_rows_latitudes = zeros(rows,1); % Modified
all_rows_peak_velocities = zeros(rows,1); % Modified
all_rows_peak_velocities_errors = zeros(rows,1); % Modified
all_rows_intensities = zeros(rows,1); % Modified
for row = 1:rows
data = data_i(row, :); % Modified
RA_H = data(4);
RA_M = data(5);
RA_S = data(6);
RA = deg2rad(15*(RA_H+(RA_M/60)+(RA_S/3600)));
DEC = deg2rad(data(9));
Div1 = cos(DEC)*sin(RA-RA_NGP);
Div2 = cos(DEC_NGP)*sin(DEC)-sin(DEC_NGP)*cos(DEC)*cos(RA-RA_NGP);
l = l_NCP - atan2(Div1, Div2);
longitude = rad2deg(l);
if longitude < 0
longitude = longitude + 360;
end
Add = sin(DEC)*sin(DEC_NGP)+cos(DEC)*cos(DEC_NGP)*cos(RA-RA_NGP);
b = asin(Add);
latitude = rad2deg(b);
if ((longitude >= 0) && (longitude <= 90) || (longitude >=270) && (longitude <= 360)) && ((latitude >= -20) && (latitude <= 20))
x_data = (1:8192);
y_data = data(:, 10:end);
linear_y = 10.^(T./10);
frequency_rest = data(:, 7);
x_frequency = linspace(frequency_rest-1.25, frequency_rest+1.25, 8192);
x_1 = x_frequency(1500);
x_2 = x_frequency(6500);
y_1 = linear_y(1500);
y_2 = linear_y(6500);
m = ((y_1-y_2)/(x_1-x_2));
c = y_1 - m*x_1;
y_fit = m.*x_frequency + c;
y_fitted = linear_y-y_fit;
radial_vel = -((x_frequency-frequency_rest)./frequency_rest).*(3*10^5);
condition = (-150 < radial_vel) & (radial_vel < 150);
radial_vel_focussed = radial_vel(condition);
y_focussed = y_fitted(condition);
y_smooth = smooth(radial_vel_focussed, y_focussed, 0.1, 'rloess');
[radial_vel_focussed, order] = sort(radial_vel_focussed);
y_smooth = y_smooth(order);
[pks, locs, w, p] = findpeaks(y_smooth, radial_vel_focussed, 'MinPeakProminence', 0.02e-10);
peak_velocity = max(locs);
peak_data = [pks, locs', w', p];
if isempty(peak_data)
peak_data = zeros(1,4)
peak_velocity = 0
end
[a b] = find(peak_data == peak_velocity);
rightmost_peak_data = peak_data(a, :);
rightmost_pks = rightmost_peak_data(1);
rightmost_locs = rightmost_peak_data(2);
rightmost_w = rightmost_peak_data(3);
rightmost_p = rightmost_peak_data(4);
all_rows_latitudes(row) = latitude;
all_rows_longitudes(row) = longitude;
all_rows_peak_velocities(row) = peak_velocity;
all_rows_peak_velocities_errors(row) = rightmost_w;
all_rows_intensities(row) = rightmost_w*rightmost_p;
end
end

Community Treasure Hunt

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

Start Hunting!

Translated by