# Low-Thrust Orbital Transfer with Solar-Electric Propulsion

### David Eagle (view profile)

29 Apr 2008 (Updated )

orbital transfer using solar-electric propulsion

sep_ltot.m
```% sep_ltot.m      April 29, 2008

% low-thrust orbit transfer between
% non-coplanar circular orbits using
% solar-electric propulsion (sep)

% Edelbaum equations

% Orbital Mechanics with Matlab

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

% conversion factors

dtr = pi/180;               % degrees to radians

rtd = 180/pi;               % radians to degrees

% astrodynamic constants

req = 6378.14;              % Earth equatorial radius (kilometers)

mu = 398600.5;              % Earth gravitational constant (km^3/sec^2)

% request inputs

clc; home;

fprintf('\n   SEP Low-thrust Orbit Transfer Analysis\n');

while (1)
fprintf('\n\nplease input the initial altitude (kilometers)\n');

alt1 = input('? ');

if (alt1 > 0)
break;
end
end

while (1)
fprintf('\nplease input the final altitude (kilometers)\n');

alt2 = input('? ');

if (alt2 > 0)
break;
end
end

while (1)
fprintf('\nplease input the initial orbital inclination (degrees)');
fprintf('\n(0 <= inclination <= 180)\n');

inc1 = input('? ');

if (inc1 >= 0 && inc1 <= 180)
break;
end
end

while (1)
fprintf('\nplease input the final orbital inclination (degrees)');
fprintf('\n(0 <= inclination <= 180)\n');

inc2 = input('? ');

if (inc2 >= 0 && inc2 <= 180)
break;
end
end

while (1)
fprintf('\nplease input the initial spacecraft mass (kilograms)\n');

xmass0 = input('? ');

if (xmass0 > 0)
break;
end
end

while (1)
fprintf('\nplease input the SEP propulsive efficiency (non-dimensional)\n');

teff = input('? ');

if (teff > 0 && teff <= 1)
break;
end
end

while (1)
fprintf('\nplease input the SEP input power (kilowatts)\n');

tpower = input('? ');

if (tpower > 0)
break;
end
end

while (1)
fprintf('\nplease input the SEP specific impulse (seconds)\n');

xisp = input('? ');

if (xisp > 0)
break;
end
end

% compute thrust (newtons)

thrust = 1000.0d0 * (2.0d0 * teff * tpower / (9.80665d0 * xisp));

% compute the thrust acceleration to km/sec^2

thracc = (thrust / xmass0) / 1000;

inc1 = inc1 * dtr;

inc2 = inc2 * dtr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the orbit transfer problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculate total inclination change

dinct = abs(inc2 - inc1);

% check for coplanar orbits

if (dinct == 0)
dinct = 1.0e-8;
end

% compute geocentric radii of initial and final orbits (kilometers)

r1 = req + alt1;

r2 = req + alt2;

% compute "local circular velocity" of initial and final orbits

v1 = sqrt(mu / r1);

v2 = sqrt(mu / r2);

% initial yaw angle

beta0 = atan(sin(0.5 * pi * dinct) / ((v1/v2) - cos(0.5 * pi * dinct)));

% delta-v

dvt = v1 * cos(beta0) - v1 * sin(beta0) / tan(0.5 * pi * dinct + beta0);

% thrust duration

tdur = dvt / thracc;

if (tdur < 3600)
% minutes
tdflag = 1;
tdur = tdur / 60;
elseif (tdur < 86400)
% hours
tdflag = 2;
tdur = tdur / 3600;
else
% days
tdflag = 3;
tdur = tdur / 86400;
end

dtstep = tdur / 100;

tsim = - dtstep;

for i = 1:1:101
tsim = tsim + dtstep;

if (tdflag == 1)
tsec = 60 * tsim;
elseif (tdflag == 2)
tsec = 3600 * tsim;
else
tsec = 86400 * tsim;
end

t(i) = tsim;

beta(i) = rtd * atan(v1 * sin(beta0) / (v1 * cos(beta0) ...
- thracc * tsec));

tmp1 = atan((thracc * tsec - v1 * cos(beta0))/(v1 * sin(beta0)));

dinc = rtd * (2/pi) * (tmp1 + 0.5 * pi - beta0);

inc(i) = rtd * inc1 - dinc;

v(i) = 1000 * sqrt(v1 * v1 - 2 * v1 * thracc * tsec * cos(beta0) ...
+ thracc * thracc * tsec * tsec);

sma(i) = 1.0e6 * mu / (v(i) * v(i));
end

xmprop = xmass0 * (1.0d0 - exp(-1000.0d0 * dvt / (9.80665d0 * xisp)));

xmassf = xmass0 - xmprop;

% print results

clc; home;

fprintf('\n   SEP Low-thrust Orbit Transfer Analysis \n\n\n');

fprintf('initial orbit altitude      %10.4f kilometers \n\n', alt1);

fprintf('initial orbit inclination   %10.4f degrees \n\n', inc1 * rtd);

fprintf('initial orbit velocity      %10.4f meters/second \n\n\n', 1000 * v1);

fprintf('final orbit altitude        %10.4f kilometers \n\n', alt2);

fprintf('final orbit inclination     %10.4f degrees \n\n', inc2 * rtd);

fprintf('final orbit velocity        %10.4f meters/second \n\n', 1000 * v2);

fprintf('\npropulsive efficiency      %10.4f \n\n', teff);

fprintf('input power                 %10.4f kilowatts\n\n', tpower);

fprintf('specific impulse            %10.4f seconds\n\n', xisp);

fprintf('thrust                      %10.4f newtons\n\n', thrust);

fprintf('initial spacecraft mass     %10.4f kilograms\n\n', xmass0);

fprintf('final spacecraft mass       %10.4f kilograms\n\n', xmassf);

fprintf('propellant mass             %10.4f kilograms\n\n', xmprop);

fprintf('\ntotal inclination change    %10.4f degrees\n\n', rtd * dinct);

fprintf('total delta-v               %10.4f meters/second \n\n', 1000 * dvt);

if (tdflag == 1)
fprintf('thrust duration             %10.4f minutes \n\n', tdur);
elseif (tdflag == 2)
fprintf('thrust duration             %10.4f hours \n\n', tdur);
else
fprintf('thrust duration             %10.4f days \n\n', tdur);
end

fprintf('initial yaw angle           %10.4f degrees \n\n', rtd * beta0);

fprintf('thrust acceleration         %10.6f meters/second^2 \n\n', 1000 * thracc);

keycheck;

% display graphics

subplot(2,1,1);

plot(t, beta);

title('SEP Low-thrust Orbit Transfer', 'FontSize', 16);

ylabel('Yaw Angle (deg)', 'FontSize', 12);

grid;

subplot(2,1,2);

plot(t, inc);

if (tdflag == 1)
xlabel('Simulation Time (minutes)', 'FontSize', 12);
elseif (tdflag == 2)
xlabel('Simulation Time (hours)', 'FontSize', 12);
else
xlabel('Simulation Time (days)', 'FontSize', 12);
end

ylabel('Inclination (deg)', 'FontSize', 12);

grid;

print -depsc -tiff -r300 sep_ltot1.eps

keycheck;

subplot(2,1,1);

plot(t, v);

title('SEP Low-thrust Orbit Transfer', 'FontSize', 16);

ylabel('Velocity (m/s)', 'FontSize', 12);

grid;

subplot(2,1,2);

plot(t, sma);

if (tdflag == 1)
xlabel('Simulation Time (minutes)', 'FontSize', 12);
elseif (tdflag == 2)
xlabel('Simulation Time (hours)', 'FontSize', 12);
else
xlabel('Simulation Time (days)', 'FontSize', 12);
end

ylabel('Semimajor Axis (kilometers)', 'FontSize', 12);

grid;

print -depsc -tiff -r300 sep_ltot2.eps
```