Code covered by the BSD License  

Highlights from
Optimal Impulsive Orbital Transfer

Optimal Impulsive Orbital Transfer

by

 

27 Nov 2012 (Updated )

MATLAB script for the solution of the one and two impulse orbit transfer between two Earth orbits.

make_textfile(soln)
function make_textfile(soln)

% create disk text files of solutions

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global mu rtd top2eci filename fname fid1 fid2

global oev1 oevt1 oev2 output

global r1 v1 r2 v2 vt1 vt2 dv1 dv2

% create filename and open file

if (soln == 1)
       
    txt_filename1 = horzcat(fname, '_solutions', '.txt');
    
    fid1 = fopen(txt_filename1, 'w');
   
    txt_filename2 = horzcat(fname, '_summary', '.txt');
    
    fid2 = fopen(txt_filename2, 'w');

    % write name of input data file
    
    fprintf(fid1, '\ninput data file ==> ');
    
    fprintf(fid1, filename);
    
    fprintf(fid1, '\n');
    
    fprintf(fid2, '\ninput data file ==> ');
    
    fprintf(fid2, filename);
    
    fprintf(fid2, '\n\n');
    
    % write file header

    fprintf(fid2, '   oota      dv1 true      dv2 true      dv1 arg      dv2 arg      delta-v1      delta-v2        total\n');
    fprintf(fid2, ' solution     anomaly       anomaly      latitude     latitude     magnitude     magnitude      delta-v\n');
    fprintf(fid2, '  number     (degrees)     (degrees)     (degrees)    (degrees)      (m/s)         (m/s)         (m/s)\n');
    fprintf(fid2, ' --------    ---------     ---------     ---------    ---------    ---------     ---------     ---------\n\n');
end

fprintf(fid1, '\nsolution number  %3i \n', soln);

fprintf(fid1, '\nnumber of function evaluations  %2i \n', output.funcCount);

fprintf(fid1, '\nnumber of iterations            %2i \n\n', output.iterations);

fprintf(fid1, output.message);

fprintf(fid1, '\n');

% compute eci delta-v vectors

dv1eci = top2eci * dv1';

dv2eci = top2eci * dv2';

% compute eci delta-v magnitudes

dv1ecim = norm(dv1eci);

dv2ecim = norm(dv2eci);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute eci state vectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial orbit - first impulse

r1eci = top2eci * r1';

v1eci = top2eci * v1';

oev1 = eci2orb1 (mu, r1eci, v1eci);

upeci = dv1eci / dv1ecim;

[uplvlh, pitch1, yaw1] = eci2lvlh (r1eci, v1eci, upeci);

fprintf(fid1, '\ninitial orbit - prior to the first impulse');
fprintf(fid1, '\n------------------------------------------\n');

oeprint2(fid1, mu, oev1);

[r, v] = orb2eci(mu, oev1);

svprint1(fid1, r, v);

% transfer orbit - first impulse

rt1eci = top2eci * r1';

vt1eci = top2eci * vt1';

oevt1 = eci2orb1 (mu, rt1eci, vt1eci);

fprintf(fid1, '\ntransfer orbit - after the first impulse');
fprintf(fid1, '\n----------------------------------------\n');

oeprint2(fid1, mu, oevt1);

[r, v] = orb2eci(mu, oevt1);

svprint1(fid1, r, v);

% transfer orbit - second impulse

rt2eci = top2eci * r2';

vt2eci = top2eci * vt2';

oevt2 = eci2orb1 (mu, rt2eci, vt2eci);

upeci = dv2eci / dv2ecim;

[uplvlh, pitch2, yaw2] = eci2lvlh (rt2eci, vt2eci, upeci);

fprintf(fid1, '\ntransfer orbit - prior to the second impulse');
fprintf(fid1, '\n--------------------------------------------\n');

oeprint2(fid1, mu, oevt2);

[r, v] = orb2eci(mu, oevt2);

svprint1(fid1, r, v);

% final orbit - second impulse

r2eci = top2eci * r2';

v2eci = top2eci * v2';

oev2 = eci2orb1 (mu, r2eci, v2eci);

fprintf(fid1, '\nfinal orbit - after the second impulse');
fprintf(fid1, '\n--------------------------------------\n');

oeprint2(fid1, mu, oev2);

[r, v] = orb2eci(mu, oev2);

svprint1(fid1, r, v);

fprintf(fid1, '\nECI delta-v vectors, magnitudes and LVLH angles');
fprintf(fid1, '\n-----------------------------------------------\n');

fprintf(fid1, '\ndelta-v1x         %12.4f meters/second', 1000.0 * dv1eci(1));
fprintf(fid1, '\ndelta-v1y         %12.4f meters/second', 1000.0 * dv1eci(2));
fprintf(fid1, '\ndelta-v1z         %12.4f meters/second', 1000.0 * dv1eci(3));
fprintf(fid1, '\n\ndelta-v1          %12.4f meters/second\n', 1000.0 * dv1ecim);

fprintf(fid1, '\nLVLH pitch angle  %12.4f degrees', rtd * pitch1);
fprintf(fid1, '\nLVLH yaw angle    %12.4f degrees\n', rtd * yaw1);

fprintf(fid1, '\ndelta-v2x         %12.4f meters/second', 1000.0 * dv2eci(1));
fprintf(fid1, '\ndelta-v2y         %12.4f meters/second', 1000.0 * dv2eci(2));
fprintf(fid1, '\ndelta-v2z         %12.4f meters/second', 1000.0 * dv2eci(3));
fprintf(fid1, '\n\ndelta-v2          %12.4f meters/second\n', 1000.0 * dv2ecim);

fprintf(fid1, '\nLVLH pitch angle  %12.4f degrees', rtd * pitch2);

fprintf(fid1, '\nLVLH yaw angle    %12.4f degrees\n', rtd * yaw2);

fprintf(fid1, '\ntotal delta-v     %12.4f meters/second\n', ...
    1000.0 * (dv1ecim + dv2ecim));

% compute flight time (seconds)

dtof = tof1(mu, oevt1(1), oevt1(2), oevt1(6), oevt2(6));

fprintf(fid1, '\ntransfer time     %12.4f seconds', dtof);

fprintf(fid1, '\n                  %12.4f minutes', dtof / 60.0);

fprintf(fid1, '\n                  %12.4f hours\n\n', dtof / 3600.0);

% write information to summary text file

arglat1 = mod(oev1(6) + oev1(4), 2.0 * pi);

arglat2 = mod(oevt2(6) + oevt2(4), 2.0 * pi);

fprintf(fid2, '  %3i        %8.4f      %8.4f      %8.4f     %8.4f    %10.4f    %10.4f    %10.4f\n', ...
    soln, rtd * oev1(6), rtd * oevt2(6), rtd * arglat1, rtd * arglat2, ...
    1000.0 * dv1ecim, 1000.0 * dv2ecim, 1000.0 * (dv1ecim + dv2ecim));

Contact us