Simulation in Matlab running too slow

12 views (last 30 days)
Alba Cuevas
Alba Cuevas on 4 Aug 2018
Commented: Alba Cuevas on 5 Aug 2018
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!
  7 Comments
Walter Roberson
Walter Roberson on 5 Aug 2018
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
Alba Cuevas on 5 Aug 2018
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')

Sign in to comment.

Answers (0)

Categories

Find more on Gravitation, Cosmology & Astrophysics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!