MATLAB Answers

0

Index position exceeds array bounds(must not exceed 15)

Asked by Suraj Parasuram on 14 Jun 2019
Latest activity Answered by Guillaume
on 15 Jun 2019
clear all;
%% CONSTANTS AND PARAMETERS TO BE USED
M=[0 0 1]; % Magnetic moment of permanent magnet i.e M=Bp*Vp*b3
u=3.986004418e14; %Gravitational parameter of the Earth
mu=4e-7*pi; %Permeability of vacuum in free space
Hc=1.59; %Constant coercivity value of hyst rods
Bm=0.73; %Magnetization at saturation of hyst rods
Hr=1.696; %Magnetic remanence of hysteresis rods
h=0.1; % Step size for Rk4 solver in seconds
hours=0.05; % Simulation time in hours
time= hours*3600; % Simulation time in seconds
period=5746;
I=[0.026066,0,0;... %Inertia tensor assuming
0,0.026066,0;... %uniform mass distribution
0,0,0.00523];
rodsPerAxis=1; % No of hysteresis rods per axis (x&y) of the satellite
rodLength=0.07; %m %Hysteresis rod dimensions
rodRadius=0.0005; %m
rodVol=pi*rodRadius^2*rodLength; %Hysteresis rod volume
%% Orbital parameters for hot case - 31 JAN and calculation
a=6928.14e3; %m %Semimajor axis of orbit
e=0; %eccentricity
i=97.5897; %deg %inclination of orbit
Omega=98.9563; %deg %Longitude of Ascending Node
w_omega=0; %Argument of Periapsis
M_Epoch=0; % Mean Anomaly at t0
Kepler=[a;e;i;Omega;w_omega;M_Epoch]; %Initial orbital elements
[pos0,vel0]=KeplertoCartesian(Kepler); % Initial position and velocity
E0=[0;0;0];% Quaternion vector initially
n0=[1;0;0];
w0=[5;5;5]; % Initial Angular Velocity
%% Calculations
%For initial magnetization, we have that ECF and ECI are aligned so no need
%to rotate from ECI to ECF.
Rbi0=ECItoBody(n0,E0); % Rotation matrix from ECI to body frame
[~,H_E0,~]= Lor_ECI_Mag(pos0,vel0); % Magnetic field strength in ECI at t0
Hb0=Rbi0*H_E0; % Magnetic field strength in body frame
H1o=Hb0(1); % For hysteresis rod 1 in x direction
H2o=Hb0(2); % For hysteresis rod 2 in y direction
if H1o>0
B1o=2*(Bm/pi)*atan((H1o-Hc)/Hr);
else
B1o=2*(Bm/pi)*atan((H1o+Hc)/Hr);
end
if H2o>0
B2o=2*(Bm/pi)*atan((H2o-Hc)/Hr);
else
B2o=2*(Bm/pi)*atan((H2o+Hc)/Hr);
end
%% Integration using Rk4
g=[Hc,Hr,Bm,rodsPerAxis,rodVol,u,M]; %Matrix containing all constants
A=[w0;n0;E0;pos0;vel0;B1o;B2o]; % Matrix containing initial conditions
F_tx= @(t,x)propagate(t,x,I,g); % Main function to evaluate
[t,Z] = rk4(h,time,F_tx,A); % Orbit propogation using Rk4 solver
%% INITIALIZE RESULTS
H_E = zeros(3,time+1); %%Magnetic field strength in ECI
Hb = zeros(3,time+1); %%Magnetic field strength in Body Frame
B_E = zeros(3,time+1); %%Magnetic field flux density in ECI
Bb = zeros(3,time+1); %%Magnetic field flux density in Body Frame
H1 = zeros(1,time+1); %%Magnetic field strength along an x-axis rod
H2 = zeros(1,time+1); %%Magnetic field strength along a y-axis rod
B1 = zeros(1,time+1); %%Magnetic field generated along x-axis rod
B2 = zeros(1,time+1); %%Magnetic field generated along y-axis rod
Mp = zeros(3,time+1); %%Torque due to permanent magnet
M1 = zeros(3,time+1); %%Torque due to all x-axis rods
M2 = zeros(3,time+1); %%Torque due to all y-axis rods
Tgg = zeros(3,time+1); %%Torque due to the gravity gradient
Error = zeros(1,time+1); %%Angle between minor (z) axis and external magnetic field
timeDays = zeros(time+1,1); %%Time of simulation in days
%% Outputs
w = Z(:,1:3);
n = Z(:,4);
E = Z(:,5:7); % Vector component of quaterion..epsilon
r = Z(:,8:10);
v = Z(:,11:13);
B1o = Z(:,14);
B2o = Z(:,15);
%% Outputs over entire orbit
for k=1:time+1 %%Each index of k corresponds to the current time+1 (k=1->t=0)
%%Declare rotation cosine matrices
Rei = ECItoECF(k-1);
ReiDot = dECItoECFdt(k-1);
Rbi = ECItoBody(n(k,:),E(k,:));
%%Fill magnetic field elements
[B_E(:,k),H_E(:,k),~] = Lor_ECI_Mag(r(k,:),v(k,:));
Hb(:,k) = Rbi*Rei'*H_E(:,k);
Bb(:,k) = Rbi*Rei'*B_E(:,k);
H1(k) = Hb(1,k);
H2(k) = Hb(2,k);
B1(k) = Bb(1,k);
B2(k) = Bb(2,k);
%%Fill torque elements
Mp(:,k) = skew(M)*(Bb(:,k));
M1(:,k) = rodsPerAxis * skew([B1(k)*rodVol;0;0])*Hb(:,k);
M2(:,k) = rodsPerAxis * skew([0;B2(k)*rodVol;0])*Hb(:,k);
%Tgg(:,k) = gravity_gradient(u,r(k,:)',Rbi,I);
%Fill Error matrix
%Error(k) = pointerr(Rbi,Rei,H_E(:,k));
%Fill TimeDays matrix
timeDays(k) = t(k)/3600/24;
end
%% Graphs
%%Easily designate which results to graph
angVel = true;
errorAng = true;
torque = true;
hystCurves = true;
inertialOrbit = false;
qNorm = false;
if angVel==true
figure;
plot(time,w(1,:),'r'); hold on;
plot(time,w(2,:),'b')
plot(time,w(3,:),'k')
title('Body-referenced angular velocities','FontSize',14);
xlabel('time (days)','FontSize',12);
ylabel('(rad/s)','FontSize',12);
legend('w(x)','w(y)','w(z)');
end
if errorAng==true
figure;
plot(timeDays,Error)
title('Error angle','FontSize',14);
xlabel('time (days)','FontSize',12);
ylabel('degrees','FontSize',12);
end
if torque==true
figure;
plot(timeDays,sqrt(Mp(1,:).^2+Mp(2,:).^2+Mp(3,:).^2),'r'); hold on;
plot(timeDays,sqrt(M1(1,:).^2+M1(2,:).^2+M1(3,:).^2),'b');
plot(timeDays,sqrt(M2(1,:).^2+M2(2,:).^2+M2(3,:).^2),'k');
plot(timeDays,sqrt(Tgg(1,:).^2+Tgg(2,:).^2+Tgg(3,:).^2),'g');
title('Magnitudes of Torques on Satellite','FontSize',14);
xlabel('time (days)','FontSize',12);
ylabel('N*m','FontSize',12);
legend('Permanent Magnet','x-Axis Rods','y-Axis Rods','Gravity Gradient');
end
if hystCurves==true
figure;
plot(H1,B1);
title('Hysteresis loop of x-axis rod','FontSize',14);
xlabel('H','FontSize',12);
ylabel('B','FontSize',12);
figure;
plot(H2,B2);
title('Hysteresis loop of y-axis rod','FontSize',14);
xlabel('H','FontSize',12);
ylabel('B','FontSize',12);
end
if inertialOrbit==true
figure;
plot3(r(:,1),r(:,2),r(:,3),'k');
title('Orbit in Inertial Frame (m)','FontSize',12);
end
if qNorm==true
figure;
plot(timeDays,sqrt(q(:,1).^2+q(:,2).^2+q(:,3).^2+q(:,4).^2)) %normalization over time
end
And the error that pops is:
Index in position 1 exceeds array bounds (must not exceed 15).
Error in Main (line 100)
Rbi = ECItoBody(n(k,:),E(k,:));
I am unable to proceed as the value of n and E are NaN. They are called quaternions which are used to represent attitude.
Here is the source: http://arrow.utias.utoronto.ca/~damaren/papers/jgcdbehrad2016.pdf

  3 Comments

I cannot even begin to help you find the source of the error simply by looking at the objects in your workspace. It would be much more helpful to see the full code and full error list.
"...It would be much more helpful to see the full code and full error list."
See the attachment with the question.
I have edited the question. Please have a look.

Sign in to comment.

Products


Release

R2019a

1 Answer

Answer by Guillaume
on 15 Jun 2019

I'm not sure what you find unclear about the error and why you can't find the problem yourself.
Index in position 1 exceeds array bounds (must not exceed 15).
Error in Main (line 100)
Rbi = ECItoBody(n(k,:),E(k,:));
So, the indexing in this line is n(k, :) and E(k, :). Position 1 indexing is k in each case. As the error message tells us, n or E has only 15 rows, but k is greater than 15. So what is k? Let's look up
for k=1:time+1
k goes up to time+1. What is time?
hours=0.05; % Simulation time in hours
time= hours*3600; % Simulation time in seconds
So time is 180 and k will count up to 181. Unless both n and E have 181 or more rows, you're going to see the error.
Now you can go back, through the construction of n and E. They're column of Z that is created in rk4.m and set to the height of A,with A:
A=[w0;n0;E0;pos0;vel0;B1o;B2o];
I've not looked much further than that. You can do it yourself. But most of the inputs to A are small column vectors, so it's unlikely that A will have 181 or more elements.

  0 Comments

Sign in to comment.