Simulation in Matlab running too slow

My simulation in matlab is running too slow because it is asking me to preallocate. The problem is that I do not know the dimension of the matrix because I do not know how many iterations will require (and requires a lot). What can I do? Also, if someone knows how to keep instead of whole matrix only the last two values and updating in each iteration, due to I am using double integrator in my loop. Thanks!

[per isakson]

See Preallocation performance section "Non-deterministic pre-allocation"

[Walter Roberson]

How does this question differ from https://www.mathworks.com/matlabcentral/answers/413431-matlab-is-running-too-slow#comment_596397 ?

[Alba Cuevas]

It is basically the same, because is my question

[KALYAN ACHARJYA]

Do not post the multiple questions, if your question having sufficient information and relevant, definitely you will get the answer.

Now the question part, remove the loops by using other alternatives method.

[Alba Cuevas]

I know, has been a mistake because I thought that have deleted the previous question because was too confused. I have been thinking in another methods but I need a for and if, I believe that the problem is the huge memory that Matlab has to use because the number of iterations is too large, so after while, Matlab begins to slow down.

[Walter Roberson]

I see the other question is gone now. But if I recall correctly, the original version had code for us to look at. When there is no code for us to look at, we can only reply in generalities, such as per's reference to Non-deterministic pre-allocation.

[Alba Cuevas]

I post the whole code, if you increase the "h" you will see that Matlab slows down because at each iteration there are more data to keep

%%CIRCULAR ORBIT%% clear all clc

%Definition of the trajectory rE=6371e3; %(m) h=500e3; %(m) r0=rE+h; %(m) % H=8.5e3; %(m) lim=60e3; %(m) %At this altitude can be considered that the debris has re-entry

G=6.67259e-11; %(N*m^2/kg^2) M=5.972e24; %(kg)

%% Space Debris m=[0.5 0.1]; %(kg) rob=[0.2 0.1]; %(m)

%Gravitational force mu=G*M; %(N*m^2/kg) vorb=sqrt(mu/r0); %(m/s) % w=vorb/r; %(rad/s) % T=2*pi/w; %(s) % f=1/T; %(1/s) % rho0=1.07e-16; %(kg/m^3)

%Aerodynamic drag A=pi*rob.^2; % sphere(m^2) Cd=1.2; v=vorb; %(m/s) %ballistic coefficient BC=(A./m)*Cd;

a=mu/r0^2; % a=Fg/m (m/s^2)

%% Atmosphere % alt(1)=r(1)-rE; %(m)

% %Name of the method % rho(1)=rho0*exp(-(r-rE)/H); %(kg/m^3)

% %Name of the method % Temp(1)=-131.21+0.00299*alt(1); %(C) % p(1)=2.488*((T(1)+273.1)/216.6)^(-11.388); %(K Pa) % rho(1)=p(1)/(0.2869*(T(1)+273.1)); %(kg/m^3)

%Solar Cycle xalt=0:10:1000; %(km) xx=0:1:1000; %(km)

%F10.7=72.8 y01=[1.168 4.200E-01 9.533E-02 1.796E-02 3.938E-03 1.053E-03 3.209E-04 8.546E-05 1.721E-05 3.505E-06 6.070E-7 7.369E-8 1.451E-8 5.540E-9 2.743E-9 1.544E-9 9.418E-10 6.067E-10 4.065E-10 2.805E-10 1.981E-10 1.425E-10 1.042E-10 7.716E-11 5.781E-11 4.377E-11 3.344E-11 2.576E-11 2.000E-11 1.563E-11 1.229E-11 9.714E-12 7.719E-12 6.162E-12 4.939E-12 3.973E-12 3.208E-12 2.597E-12 2.109E-12 1.717E-12 1.401E-12 1.145E-12 9.384E-13 1.251E-13 6.336E-13 5.220E-13 4.308E-13 3.561E-13 2.949E-13 2.446E-13 2.033E-13 1.827E-13 1.527E-13 1.279E-13 1.073E-13 9.019E-14 7.598E-14 6.416E-14 5.432E-14 4.610E-14 3.924E-14 3.350E-14 2.869E-14 2.465E-14 2.126E-14 1.841E-14 1.600E-14 1.397E-14 1.224E-14 1.078E-14 9.537E-15 8.476E-15 7.569E-15 6.790E-15 6.120E-15 5.542E-15 5.041E-15 4.605E-15 4.225E-15 3.892E-15 3.599E-15 3.340E-15 3.110E-15 2.905E-15 2.721E-15 2.557E-15 2.408E-15 2.273E-15 2.150E-15 2.037E-15 1.934E-15 1.839E-15 1.752E-15 1.670E-15 1.595E-15 1.524E-15 1.458E-15 1.396E-15 1.338E-15 1.284E-15 1.232E-15]; yy01=interp1(xalt,y01,xx,'spline'); %F10.7=93.6 y02=[1.167 4.196E-01 9.524E-02 1.794E-02 3.934E-03 1.052E-03 3.206E-04 8.528E-05 1.722E-05 3.425E-06 5.841E-7 7.436E-8 1.476E-8 5.631E-9 2.810E-9 1.601E-9 9.917E-10 6.502E-10 4.441E-10 3.128E-10 2.257E-10 1.660E-10 1.240E-10 9.385E-11 7.186E-11 5.558E-11 4.336E-11 3.410E-11 2.700E-11 2.152E-11 1.725E-11 1.390E-11 1.125E-11 9.150E-12 7.469E-12 6.119E-12 5.028E-12 4.145E-12 3.426E-12 2.839E-12 2.357E-12 1.962E-12 1.636E-12 1.366E-12 1.143E-12 9.582E-13 8.043E-13 6.761E-13 5.691E-13 4.797E-13 4.049E-13 3.820E-13 3.247E-13 2.765E-13 2.356E-13 2.011E-13 1.719E-13 1.471E-13 1.260E-13 1.082E-13 9.298E-14 8.005E-14 6.903E-14 5.963E-14 5.161E-14 4.476E-14 3.889E-14 3.387E-14 2.957E-14 2.587E-14 2.270E-14 1.997E-14 1.761E-14 1.559E-14 1.383E-14 1.231E-14 1.100E-14 9.857E-15 8.864E-15 7.998E-15 7.242E-15 6.580E-15 6.000E-15 5.490E-15 5.040E-15 4.642E-15 4.290E-15 3.977E-15 3.698E-15 3.449E-15 3.225E-15 3.024E-15 2.842E-15 2.678E-15 2.529E-15 2.393E-15 2.268E-15 2.154E-15 2.049E-15 1.953E-15 1.863E-15]; yy02=interp1(xalt,y02,xx,'spline'); %F10.7=139.0 y03=[1.165 4.190E-01 9.510E-02 1.792E-02 3.928E-03 1.050E-03 3.201E-04 8.493E-05 1.684E-05 3.554E-06 6.224E-7 6.940E-8 1.474E-8 5.641E-9 2.823E-9 1.631E-9 1.029E-9 6.900E-10 4.826E-10 3.484E-10 2.577E-10 1.944E-10 1.489E-10 1.156E-10 9.073E-11 7.189E-11 5.743E-11 4.622E-11 3.744E-11 3.051E-11 2.499E-11 2.057E-11 1.700E-11 1.411E-11 1.176E-11 9.826E-12 8.237E-12 6.925E-12 5.837E-12 4.931E-12 4.175E-12 3.543E-12 3.011E-12 2.564E-12 2.187E-12 1.868E-12 1.598E-12 1.369E-12 1.174E-12 1.008E-12 8.668E-13 1.261E-12 1.103E-12 9.652E-13 8.455E-13 7.414E-13 6.506E-13 5.715E-13 5.024E-13 4.420E-13 3.892E-13 3.430E-13 3.025E-13 2.670E-13 2.358E-13 2.085E-13 1.845E-13 1.634E-13 1.448E-13 1.284E-13 1.140E-13 1.013E-13 9.014E-14 8.026E-14 7.153E-14 6.382E-14 5.700E-14 5.096E-14 4.562E-14 4.089E-14 3.670E-14 3.297E-14 2.967E-14 2.673E-14 2.412E-14 2.180E-14 1.973E-14 1.789E-14 1.625E-14 1.478E-14 1.347E-14 1.230E-14 1.125E-14 1.031E-14 9.466E-15 8.708E-15 8.028E-15 7.415E-15 6.863E-15 6.365E-15 5.915E-15]; yy03=interp1(xalt,y03,xx,'spline'); %F10.7=170.0 y04=[1.163 4.181E-01 9.490E-02 1.788E-02 3.920E-03 1.048E-03 3.194E-04 8.460E-05 1.693E-05 3.369E-06 5.699E-7 7.169E-8 1.527E-8 5.843E-9 2.971E-9 1.748E-9 1.126E-9 7.697E-10 5.492E-10 4.045E-10 3.052E-10 2.347E-10 1.834E-10 1.451E-10 1.160E-11 9.361E-11 7.614E-11 6.237E-11 5.139E-11 4.258E-11 3.546E-11 2.965E-11 2.490E-11 2.099E-11 1.775E-11 1.505E-11 1.281E-11 1.092E-11 9.337E-12 8.000E-12 6.868E-12 5.908E-12 5.092E-12 4.395E-12 3.800E-12 3.290E-12 2.852E-12 2.476E-12 2.151E-12 1.872E-12 1.631E-12 1.454E-12 1.268E-12 1.108E-12 9.685E-13 8.474E-13 7.421E-13 6.504E-13 5.705E-13 5.008E-13 4.400E-13 3.869E-13 3.405E-13 2.999E-13 2.643E-13 2.332E-13 2.059E-13 1.819E-13 1.609E-13 1.425E-13 1.262E-13 1.120E-13 9.940E-14 8.834E-14 7.860E-14 7.000E-14 6.242E-14 5.573E-14 4.981E-14 4.459E-14 3.996E-14 3.587E-14 3.224E-14 2.902E-14 2.617E-14 2.363E-14 2.138E-14 1.937E-14 1.759E-14 1.600E-14 1.458E-14 1.331E-14 1.218E-14 1.117E-14 1.026E-14 9.443E-15 8.712E-15 8.054E-15 7.461E-15 6.927E-15 6.445E-15]; yy04=interp1(xalt,y04,xx,'spline'); %F10.7=156.3 y05=[1.167 4.196E-01 9.525E-02 1.795E-02 3.935E-03 1.052E-03 3.206E-04 8.498E-05 1.703E-05 3.390E-06 5.743E-7 7.233E-8 1.518E-8 5.809E-9 2.943E-9 1.722E-9 1.103E-9 7.501E-10 5.327E-10 3.908E-10 2.938E-10 2.253E-10 1.756E-10 1.386E-10 1.106E-11 8.917E-11 7.245E-11 5.929E-11 4.882E-11 4.044E-11 3.366E-11 2.815E-11 2.364E-11 1.992E-11 1.685E-11 1.430E-11 1.217E-11 1.038E-11 8.882E-12 7.615E-12 6.543E-12 5.632E-12 4.857E-12 4.196E-12 3.631E-12 3.146E-12 2.730E-12 2.372E-12 2.063E-12 1.797E-12 1.567E-12 1.532E-12 1.343E-12 1.180E-12 1.037E-12 9.118E-13 8.027E-13 7.072E-13 6.236E-13 5.503E-13 4.860E-13 4.296E-13 3.800E-13 3.363E-13 2.979E-13 2.641E-13 2.343E-13 2.080E-13 1.848E-13 1.643E-13 1.462E-13 1.302E-13 1.160E-13 1.035E-13 9.241E-14 8.258E-14 7.385E-14 6.612E-14 5.925E-14 5.314E-14 4.772E-14 4.290E-14 3.860E-14 3.478E-14 3.138E-14 2.834E-14 2.563E-14 2.321E-14 2.105E-14 1.912E-14 1.739E-14 1.584E-14 1.445E-14 1.321E-14 1.209E-14 1.108E-14 1.018E-14 9.369E-15 8.637E-15 7.976E-15 7.380E-15]; yy05=interp1(xalt,y05,xx,'spline'); %F10.7=209.3 y06=[1.172 4.213E-01 9.562E-02 1.802E-02 3.950E-03 1.056E-03 3.218E-04 8.506E-05 1.704E-05 3.304E-06 5.482E-7 7.241E-8 1.569E-8 6.039E-9 3.114E-9 1.866E-9 1.226E-9 8.560E-10 6.245E-10 4.705E-10 3.634E-10 2.863E-10 2.290E-10 1.856E-10 1.521E-11 1.257E-10 1.048E-10 8.787E-11 7.415E-11 6.290E-11 5.361E-11 4.588E-11 3.942E-11 3.398E-11 2.939E-11 2.549E-11 2.217E-11 1.933E-11 1.689E-11 1.479E-11 1.297E-11 1.140E-11 1.004E-11 8.850E-12 7.815E-12 6.911E-12 6.119E-12 5.425E-12 4.815E-12 4.278E-12 3.805E-12 3.375E-12 3.009E-12 2.685E-12 2.399E-12 2.144E-12 1.919E-12 1.718E-12 1.539E-12 1.380E-12 1.238E-12 1.112E-12 9.989E-13 8.979E-13 8.077E-13 7.269E-13 6.546E-13 5.899E-13 5.318E-13 4.797E-13 4.330E-13 3.910E-13 3.533E-13 3.193E-13 2.888E-13 2.614E-13 2.367E-13 2.144E-13 1.944E-13 1.763E-13 1.600E-13 1.452E-13 1.319E-13 1.199E-13 1.091E-13 9.928E-14 9.041E-14 8.238E-14 7.511E-14 6.854E-14 6.258E-14 5.717E-14 5.227E-14 4.783E-14 4.380E-14 4.014E-14 3.681E-14 3.379E-14 3.104E-14 2.854E-14 2.626E-14]; yy06=interp1(xalt,y06,xx,'spline'); %F10.7=134.6 y07=[1.165 4.188E-01 9.506E-02 1.791E-02 3.927E-03 1.050E-03 3.200E-04 8.492E-05 1.716E-05 3.325E-06 5.561E-7 7.450E-8 1.517E-8 5.790E-9 2.928E-9 1.703E-9 1.081E-9 7.275E-10 5.110E-10 3.705E-10 2.752E-10 2.085E-10 1.604E-10 1.251E-10 9.859E-11 7.847E-11 6.297E-11 5.090E-11 4.141E-11 3.389E-11 2.788E-11 2.304E-11 1.912E-11 1.594E-11 1.333E-11 1.118E-11 9.411E-12 7.942E-12 6.719E-12 5.698E-12 4.842E-12 4.123E-12 3.517E-12 3.006E-12 2.573E-12 2.205E-12 1.893E-12 1.627E-12 1.400E-12 1.206E-12 1.041E-12 8.980E-13 7.765E-13 6.721E-13 5.823E-13 5.050E-13 4.384E-13 3.810E-13 3.314E-13 2.885E-13 2.514E-13 2.193E-13 1.914E-13 1.673E-13 1.464E-13 1.282E-13 1.124E-13 9.866E-14 8.670E-14 7.628E-14 6.720E-14 5.928E-14 5.236E-14 4.632E-14 4.104E-14 3.642E-14 3.237E-14 2.883E-14 2.572E-14 2.298E-14 2.058E-14 1.847E-14 1.661E-14 1.497E-14 1.352E-14 1.224E-14 1.111E-14 1.010E-14 9.214E-15 8.424E-15 7.721E-15 7.094E-15 6.535E-15 6.035E-15 5.587E-15 5.186E-15 4.825E-15 4.499E-15 4.206E-15 3.941E-15 3.700E-15]; yy07=interp1(xalt,y07,xx,'spline'); %F10.7=109.1 y08=[1.165 4.189E-01 9.510E-02 1.792E-02 3.928E-03 1.050E-03 3.201E-04 8.507E-05 1.717E-05 3.393E-06 5.754E-7 7.415E-8 1.489E-8 5.679E-9 2.846E-9 1.634E-9 1.021E-9 6.760E-10 4.665E-10 3.320E-10 2.420E-10 1.799E-10 1.358E-10 1.038E-10 8.031E-11 6.273E-11 4.941E-11 3.922E-11 3.134E-11 2.520E-11 2.038E-11 1.656E-11 1.352E-11 1.108E-11 9.117E-12 7.528E-12 6.235E-12 5.180E-12 4.314E-12 3.602E-12 3.014E-12 2.527E-12 2.123E-12 1.787E-12 1.506E-12 1.272E-12 1.075E-12 9.105E-13 7.719E-13 6.553E-13 5.571E-13 6.882E-13 5.914E-13 5.088E-13 4.381E-13 3.777E-13 3.260E-13 2.817E-13 2.436E-13 2.109E-13 1.829E-13 1.587E-13 1.379E-13 1.200E-13 1.045E-13 9.115E-14 7.962E-14 6.965E-14 6.103E-14 5.356E-14 4.708E-14 4.146E-14 3.658E-14 3.234E-14 2.865E-14 2.544E-14 2.263E-14 2.019E-14 1.805E-14 1.618E-14 1.454E-14 1.310E-14 1.184E-14 1.072E-14 9.744E-15 8.879E-15 8.114E-15 7.436E-15 6.835E-15 6.300E-15 5.823E-15 5.397E-15 5.015E-15 4.673E-15 4.366E-15 4.088E-15 3.838E-15 3.610E-15 3.404E-15 3.216E-15 3.044E-15]; yy08=interp1(xalt,y08,xx,'spline'); %F10.7=95.8 y09=[1.167 4.194E-01 9.521E-02 1.794E-02 3.933E-03 1.052E-03 3.204E-04 8.524E-05 1.714E-05 3.474E-06 5.984E-7 7.313E-8 1.467E-8 5.599E-9 2.788E-9 1.587E-9 9.819E-10 6.433E-10 4.390E-10 3.089E-10 2.226E-10 1.635E-10 1.220E-10 9.224E-11 7.053E-11 5.448E-11 4.244E-11 3.333E-11 2.636E-11 2.098E-11 1.679E-11 1.351E-11 1.093E-11 8.873E-12 7.233E-12 5.918E-12 4.857E-12 3.999E-12 3.301E-12 2.732E-12 2.266E-12 1.883E-12 1.568E-12 1.308E-12 1.093E-12 9.154E-13 7.674E-13 6.443E-13 5.417E-13 4.561E-13 3.846E-13 4.534E-13 3.863E-13 3.296E-13 2.815E-13 2.408E-13 2.062E-13 1.768E-13 1.518E-13 1.305E-13 1.124E-13 9.690E-14 8.369E-14 7.240E-14 6.275E-14 5.449E-14 4.740E-14 4.132E-14 3.610E-14 3.162E-14 2.775E-14 2.442E-14 2.155E-14 1.907E-14 1.692E-14 1.506E-14 1.345E-14 1.204E-14 1.082E-14 9.755E-15 8.824E-15 8.008E-15 7.292E-15 6.663E-15 6.108E-15 5.618E-15 5.183E-15 4.797E-15 4.453E-15 4.146E-15 3.871E-15 3.624E-15 3.401E-15 3.200E-15 3.017E-15 2.851E-15 2.699E-15 2.560E-15 2.433E-15 2.316E-15 2.207E-15]; yy09=interp1(xalt,y09,xx,'spline'); %F10.7=78.3 y10=[1.168 4.200E-01 9.534E-02 1.796E-02 3.938E-03 1.053E-03 3.209E-04 8.544E-05 1.716E-05 3.529E-06 6.139E-7 7.296E-8 1.449E-8 5.536E-9 2.741E-9 1.545E-9 9.451E-10 6.110E-10 4.112E-10 2.851E-10 2.024E-10 1.464E-10 1.075E-10 8.008E-11 6.032E-11 4.591E-11 3.526E-11 2.730E-11 2.130E-11 1.672E-11 1.321E-11 1.050E-11 8.380E-12 6.721E-12 5.412E-12 4.374E-12 3.547E-12 2.885E-12 2.353E-12 1.925E-12 1.578E-12 1.296E-12 1.067E-12 8.796E-13 7.266E-13 6.013E-13 4.984E-13 4.138E-13 3.441E-13 2.866E-13 2.391E-13 2.752E-13 2.321E-13 1.960E-13 1.658E-13 1.405E-13 1.192E-13 1.014E-13 8.631E-14 7.364E-14 6.296E-14 5.394E-14 4.632E-14 3.987E-14 3.441E-14 2.978E-14 2.584E-14 2.249E-14 1.964E-14 1.721E-14 1.513E-14 1.335E-14 1.183E-14 1.052E-14 9.387E-15 8.413E-15 7.570E-15 6.840E-15 6.204E-15 5.650E-15 5.165E-15 4.740E-15 4.366E-15 4.036E-15 3.743E-15 3.483E-15 3.250E-15 3.042E-15 2.855E-15 2.687E-15 2.534E-15 2.395E-15 2.268E-15 2.152E-15 2.045E-15 1.947E-15 1.856E-15 1.772E-15 1.694E-15 1.620E-15 1.552E-15]; yy10=interp1(xalt,y10,xx,'spline'); %F10.7=76.1 y11=[1.170 4.205E-01 9.546E-02 1.798E-02 3.943E-03 1.054E-03 3.213E-04 8.556E-05 1.731E-05 3.446E-06 5.894E-7 7.507E-8 1.467E-8 5.599E-9 2.783E-9 1.574E-9 9.658E-10 6.267E-10 4.235E-10 2.950E-10 2.104E-10 1.530E-10 1.131E-10 8.467E-11 6.415E-11 4.910E-11 3.793E-11 2.954E-11 2.318E-11 1.830E-11 1.454E-11 1.162E-11 9.325E-12 7.519E-12 6.088E-12 4.947E-12 4.033E-12 3.298E-12 2.704E-12 2.223E-12 1.832E-12 1.512E-12 1.251E-12 1.037E-12 8.612E-13 7.163E-13 5.967E-13 4.979E-13 4.160E-13 3.482E-13 2.918E-13 2.580E-13 2.172E-13 1.831E-13 1.547E-13 1.308E-13 1.109E-13 9.412E-14 8.005E-14 6.823E-14 5.828E-14 4.989E-14 4.282E-14 3.684E-14 3.179E-14 2.751E-14 2.388E-14 2.080E-14 1.817E-14 1.594E-14 1.403E-14 1.240E-14 1.100E-14 9.800E-15 8.766E-15 7.873E-15 7.101E-15 6.432E-15 5.849E-15 5.340E-15 4.895E-15 4.504E-15 4.159E-15 3.854E-15 3.583E-15 3.342E-15 3.126E-15 2.932E-15 2.758E-15 2.599E-15 2.456E-15 2.325E-15 2.205E-15 2.095E-15 1.993E-15 1.900E-15 1.813E-15 1.732E-15 1.657E-15 1.586E-15 1.520E-15]; yy11=interp1(xalt,y11,xx,'spline');

years=[yy01;yy02;yy03;yy04;yy05;yy06;yy07;yy08;yy09;yy10;yy11]; %% Initial conditions % ellipse orbit % vt=vorb; % beta(1)=0; % r(1)=r; % x(1)=r*cos(beta); % vx=0; % ax(1)=-(mu/(r)^2)*cos(beta); % xp=x; % y(1)=r*sin(beta); % vy=sqrt(vt^2-vx^2)+1000; % ay(1)=-(mu/(r)^2)*sin(beta); % yp=y;

% cicular orbit beta=[0 0]; r=[r0 r0]; alt=round((r0-rE*ones(1,length(r)))/1e3); %(km)(Solar Cycle) x=r.*cos(beta); vx=sqrt(mu./r).*sin(beta); ax=-(mu./(r).^2).*cos(beta); xp=x; y=r.*sin(beta); vy=sqrt(mu./r).*cos(beta); ay=-(mu./(r).^2).*sin(beta); yp=y;

dt=10; %(s) % t=0:dt:1*T; % t(1)=dt*1;

i=2; rho=years(1,alt+1); %(vary the initial year of the Solar Cycle modifying the number of rows in the years' matrix) j=0; %(number of solar cycle+1) % for i=2:size(t,2) added(1:5)=false;

%% PROCESS while max(r(i-1,:))>rE+lim && i<1e8 tic %%Name of the method % rho(i)= rho0*exp(-(r(i)-rE)/H);

    %%Name of the method
    % Temp(i)=-131.21+0.00299*alt(i);
    % p(i)=2.488*((T(i)+273.1)/216.6)^(-11.388);
    % rho(i)=p(i)/(0.2869*(T(i)+273.1));

    time(i)=(dt*i)/86400; %(days)
    t=(dt*i)/31536000; %(years)

    p=t-11*j; %(years 1-11 years)

% if (p>2)&&(p<3) && ~added(1) % m_add=[1 1]; % rob_add=[0.2 0.2]; % r_add=[r0 r0]; % beta_add=[0 0]; % num=2; % number of objects that have been added % [beta,r,alt,rho,x,vx,ax,xp,y,vy,ay,yp,BC,m,A,rob]=spacedebris(m,rob,beta_add,beta,num,r,alt,r_add,x,vx,ax,y,vy,ay,Cd,rE,years,mu,m_add,rob_add); % added(1)=true; % elseif (p>4)&&(p<5) && ~added(2) % m_add=[1 1]; % rob_add=[0.2 0.2]; % r_add=[r0 r0]; % beta_add=[0 0]; % num=2; % number of objects that have been added % [beta,r,alt,rho,x,vx,ax,xp,y,vy,ay,yp,BC,m,A,rob]=spacedebris(m,rob,beta_add,beta,num,r,alt,r_add,x,vx,ax,y,vy,ay,Cd,rE,years,mu,m_add,rob_add); % added(2)=true; % elseif (p>5)&&(p>6) && ~added(3) % m_add=[1 1]; % rob_add=[0.2 0.2]; % r_add=[r0 r0]; % beta_add=[0 0]; % num=2; % number of objects that have been added % [beta,r,alt,rho,x,vx,ax,xp,y,vy,ay,yp,BC,m,A,rob]=spacedebris(m,rob,beta_add,beta,num,r,alt,r_add,x,vx,ax,y,vy,ay,Cd,rE,years,mu,m_add,rob_add); % added(3)=true; % elseif (p>7)&&(p<8) && ~added(4) % m_add=[1 1]; % rob_add=[0.2 0.2]; % r_add=[r0 r0]; % beta_add=[0 0]; % num=2; % number of objects that have been added % [beta,r,alt,rho,x,vx,ax,xp,y,vy,ay,yp,BC,m,A,rob]=spacedebris(m,rob,beta_add,beta,num,r,alt,r_add,x,vx,ax,y,vy,ay,Cd,rE,years,mu,m_add,rob_add); % added(4)=true; % elseif (p>9)&&(p<11) && ~added(5) % m_add=[1 1]; % rob_add=[0.2 0.2]; % r_add=[r0 r0]; % beta_add=[0 0]; % num=2; % number of objects that have been added % [beta,r,alt,rho,x,vx,ax,xp,y,vy,ay,yp,BC,m,A,rob]=spacedebris(m,rob,beta_add,beta,num,r,alt,r_add,x,vx,ax,y,vy,ay,Cd,rE,years,mu,m_add,rob_add); % added(5)=true; % end

    for k=1:size(r,2)
        if r(i-1,k)<=rE+lim
            r(i,k)=rE; %stop the calculation because has arrive to the Earth 
        else
            axd(i,k)=((1/2)*rho(k)*BC(k)*vx(i-1,k)^2)*sin(beta(k));
            ayd(i,k)=((1/2)*rho(k)*BC(k)*vy(i-1,k)^2)*cos(beta(k));
            ax(i,k)=-(mu./(r(i-1,k).^2))*cos(beta(k))+axd(i,k);
            ay(i,k)=-(mu./(r(i-1,k).^2))*sin(beta(k))-ayd(i,k);
            vx(i,k)=vx(i-1,k)+ax(i,k)*dt;
            vy(i,k)=vy(i-1,k)+ay(i,k)*dt;
            x(i,k)=x(i-1,k)+vx(i,k)*dt;
            y(i,k)=y(i-1,k)+vy(i,k)*dt;
            r(i,k)=sqrt(x(i,k).^2+y(i,k).^2);
            beta(k)=atan((y(i,k)/x(i,k)));             
    % re-ajusting the angles to obtain the correct position
    % MATLAB tends to calculate the smaller angle with the same result
            if x(i,k)<0
                beta(k)=beta(k)+pi;
            end
            alt(i,k)=round((r(i,k)-rE*ones(1,length(r(i,k))))/1e3); %(km) (Solar Cycle)   

            %Solar Cycle
            if (p<=1)
                rho(k)=years(1,alt(i,k)+1);
            elseif (p>1)&&(p<=2)
                rho(k)=years(2,alt(i,k)+1);
            elseif (p>2)&&(p<=3)
                rho(k)=years(3,alt(i,k)+1);
            elseif (p>3)&&(p<=4)
                rho(k)=years(4,alt(i,k)+1);
            elseif (p>4)&&(p<=5)
                rho(k)=years(5,alt(i,k)+1);
            elseif (p>5)&&(p<=6)
                rho(k)=years(6,alt(i,k)+1);
            elseif (p>6)&&(p<=7)
                rho(k)=years(7,alt(i,k)+1);
            elseif (p>7)&&(p<=8)
                rho(k)=years(8,alt(i,k)+1);
            elseif (p>8)&&(p<=9)
                rho(k)=years(9,alt(i,k)+1);
            elseif (p>9)&&(p<=10)
                rho(k)=years(10,alt(i,k)+1);
            elseif (p>10)&&(p<=11)
                rho(k)=years(11,alt(i,k)+1);
            else
                j=j+1;
                rho(k)=years(1,alt(i,k)+1);
            end
        end
    end

    % plot(x,y,'b')
    % % axis ([-r r -r r])
    % viscircles([0,0],rE)
    % axis equal
    %
    % hold on
    % plot(x(i),y(i), 'or','MarkerSize',3,'MarkerFaceColor','r')
    % plot(0,0,'*') %centre
    %
    % hold off
    % pause(0.0000000001)

    i=i+1;
    toc
end

% plot(x,y,'b'); % hold on % scatter (0,0,300,'blue','filled') % % viscircles([0,0],rE); % axis equal

plot(time(1:20:end),alt(1:20:end,:)) title('Variation altitude vs. time in the Solar Cycle'),xlabel('days'), ylabel('km')