Index position exceeds array bounds(must not exceed 15)

 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

[James Browne]

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.

[madhan ravi]

"...It would be much more helpful to see the full code and full error list."

See the attachment with the question.

[Suraj Parasuram]

I have edited the question. Please have a look.

Answer 1

[Guillaume]

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.