Index position exceeds array bounds(must not exceed 15)
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
madhan ravi
2019-6-15
"...It would be much more helpful to see the full code and full error list."
See the attachment with the question.
回答(1 个)
Guillaume
2019-6-15
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Satellite Mission Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!