Add files via upload

animateMission.m

@@ -0,0 +1,64 @@
+function animateMission(N1, N2, T1, T2, startTime, maxDuration, proceedHandle, axHandle)
+
+planets = unique([N1 N2]);
+numPlanets = length(planets);
+
+phase = 1;
+
+% animation step in days -- must evenly divide maxDuration
+step = 5;
+
+lag = 0*step;
+
+h = zeros((maxDuration/step)+1,numPlanets);
+
+%set(proceedHandle,'String','Animate')
+%disp(get(proceedHandle,'Value'))
+%waitfor(proceedHandle,'Value',1)
+input('Press Return key to watch mission animation.');
+
+hold on
+cla(axHandle)
+axes(axHandle)
+axis equal off
+
+% plot sun and orbits
+plot3(axHandle, 0,0,0,'yo')
+for body = 1:max(planets)
+    drawOrbit(axHandle, body,startTime);
+end
+drawnow
+
+% 1 day/frame
+for day = 0:step:maxDuration
+
+    % check phase
+    if day>T2(phase)
+        phase = phase + 1;
+    end
+
+    % plot planets
+    for body = 1:numPlanets
+        if planets(body) == N1(phase) && planets(body) == N2(phase)
+            drawPoint_Planet(axHandle, planets(body),startTime+day/36525,'ro');
+        else
+            h(day/step+1,body) = drawPoint_Planet(axHandle, planets(body),startTime+day/36525,'ko');
+        end
+        if day>lag
+            if h((day-lag)/step,body)
+                delete(h((day-lag)/step,body));
+            end
+            %h((day-lag)/step,body) = drawPoint_Planet(planets(body),startTime+(day-lag)/36525,'k.');
+        end 
+    end
+
+
+    % plot transit
+    if N1(phase) ~= N2(phase)
+        drawPoint_Transit(axHandle, N1(phase),N2(phase),startTime+T1(phase)/36525,T2(phase)-T1(phase),day-T1(phase),'b.');
+    end
+
+    drawnow
+end
+
+end

calculateDV.m

@@ -0,0 +1,46 @@
+function [deltaV] = calculateDV(dt, p1Index, p2Index, T, dtOption, c3Opt)
+
+p1Struct = planetState(p1Index,T);
+p2Struct = planetState(p2Index,T+dt/36525);
+
+% [V1, V2] = lambertSolver(p1Struct.r, p2Struct.r, dt*24*3600, 'pro');
+% deltaV = norm(V1-p1Struct.v)+norm(V2-p2Struct.v);
+
+[V1a, V2a, ~, errora] = lambert(p1Struct.r', p2Struct.r', dt, 0, 132712440018.9);
+V1a = V1a';
+V2a = V2a';
+
+[V1b, V2b, ~, errorb] = lambert(p1Struct.r', p2Struct.r', -dt, 0, 132712440018.9);
+V1b = V1b';
+V2b=V2b';
+
+if errora<1 || errorb<1
+    disp('Error calculating deltaV.')
+end
+
+%%% deltaV cost methods
+
+% sum of c3 and dv
+if c3Opt
+    deltaVa = norm(V1a-p1Struct.v)^2+norm(V2a-p2Struct.v);
+    deltaVb = norm(V1b-p1Struct.v)^2+norm(V2b-p2Struct.v);
+else
+    % sum of individual dv
+    deltaVa = norm(V1a-p1Struct.v)+norm(V2a-p2Struct.v);
+    deltaVb = norm(V1b-p1Struct.v)+norm(V2b-p2Struct.v);
+end
+
+% flyby hack
+%deltaVa = abs(norm(V1a)-norm(p1Struct.v))+abs(norm(V2a)-norm(p2Struct.v));
+%deltaVb = abs(norm(V1b)-norm(p1Struct.v))+abs(norm(V2b)-norm(p2Struct.v));
+
+%%% cost increases with time of flight
+switch dtOption
+    case 1
+        dtFactor = 1;
+    case 2
+        dtFactor = sqrt(dt);
+end
+deltaV = min([deltaVa deltaVb])*dtFactor;
+
+end

centPastJ2000.m

@@ -0,0 +1,7 @@
+function T = centPastJ2000(date)
+
+J0 = 367*date.y - fix(7*(date.y+fix((date.m+9)/12))/4) + fix(275*date.m/9) + date.d + 1721013.5;
+JD = J0 + date.UT/24;
+T = (JD-2451545)/36525;
+
+end

defineConstraints.m

@@ -0,0 +1,108 @@
+
+%%% Indequalities
+
+%%% 1 Minimum number of Transits equal to number of planets
+%%% 2  Depart from each Planet at least once
+%%% 3 Arrive at each Planet at least once
+%%% 4 Only One Event Can Begin at Each Time
+%%% 5 Only One Event Can End at Each Time
+%%% 6 Limit Orbit Time to Zero or One
+
+
+
+%%% Equalities
+
+%%% 1 Home Planet Must Be First
+%%% 2 Target Planet Must Be Last
+%%% 3 Time Pairs Must Match
+
+function [Ai_x, Ae_x] = defineConstraints(mi,me,t1,t2,n1_index,n2_index,t1_index,startPlanetIndex,endPlanetIndex,arriveTimeVector,departTimeVector,numberOfPlanets)
+
+Ai_x = zeros(sum(mi),1);
+Ae_x = zeros(sum(me),1);
+t2_c = find(arriveTimeVector==t2);
+
+nonReturn = n1_index~=n2_index;
+
+%%% Indequalities
+
+%%% 1. Minimum number of Transits equal to number of planets (or 1 less if
+%%% non-return)
+% Greater than or Equal to
+
+if n1_index ~= n2_index
+    Ai_x(1) = -1;
+end
+
+
+
+%%%2.  Depart from each Planet at least once (unless end planet of
+%%%nonreturn mission)
+% Greater than or Equal to
+if n1_index ~= n2_index
+    Ai_x(mi(1)+n1_index) = -1;
+end
+
+
+
+%%% 3. Arrive at each Planet at least once (unless start planet of
+%%% nonreturn mission)
+% Greater than or Equal to
+
+if n1_index ~= n2_index
+    Ai_x(sum(mi(1:2))+n2_index) = -1;
+end
+
+
+
+%%% 4. Only One Event Can Begin at Each Time
+% Less than or Equal to
+Ai_x(sum(mi(1:3))+t1_index) = 1;
+
+
+
+%%% 5. Only One Event Can End at Each Time
+% Less than or Equal to
+Ai_x(sum(mi(1:4))+t2_c) = 1;
+
+
+
+%%% 6. Limit Orbit Time to Zero or One
+% Less than or Equal to
+
+%if n1_index == n2_index
+    %if ~nonReturn && t2~=arriveTimeVector(end)
+    %    Ai_x(sum(mi(1:5))+numberOfPlanets) = 1;
+%elseif  ~nonReturn && t2~=arriveTimeVector(end)
+    %else
+%        Ai_x(sum(mi(1:5))+n1_index) = 1;
+    %end
+%end
+
+
+
+
+%%% Equalities
+
+
+%%% 1 Home Planet Must Be First
+if n1_index == startPlanetIndex && t1 == departTimeVector(1)
+    Ae_x(1) = 1;
+end
+
+%%% 2 Target Planet Must Be Last
+if n2_index == endPlanetIndex && t2 == arriveTimeVector(end)
+    Ae_x(2) = 1;
+end
+
+%%% 3 Time Pairs Must Match
+if t1_index ~= 1
+    Ae_x(me(1)+(length(departTimeVector)-1)*(n1_index-1)+(t1_index-1)) = 1;
+end
+
+if t2~=arriveTimeVector(end)
+    Ae_x(me(1)+(length(departTimeVector)-1)*(n2_index-1)+t2_c) = -1;
+end
+
+
+end

displayProgress.m

@@ -0,0 +1,115 @@
+function displayProgress(progressTracker,iterationCounter,totalIterations)
+
+
+% Define strings used in DisplayProgress
+str = '............................................................';
+str2= '                                                            ';
+
+%{
+% Only Consider Nonzero Entries
+progressTracker = progressTracker(progressTracker~=0);
+
+clc;
+disp(' ')
+
+disp('Data Point:')
+disp(iterationCounter)
+
+if length(progressTracker) >= 25
+    disp('Current Speed (ms per calculation):')
+    disp(1000*mean(progressTracker(end-24:end)))
+else
+    disp('Current Speed (ms per calculation):')
+    disp('...')
+end
+
+disp('Average Speed (ms per calculation):')
+disp(1000*mean(progressTracker))
+
+disp('Total Time Elapsed:')
+s = sum(progressTracker);
+hours = floor(s/3600);
+s = s - hours*3600;
+mins = floor(s/60);
+secs = s-mins*60;
+hms = sprintf('%02d:%02d:%05.2f\n', hours, mins, secs);
+disp(hms)
+
+
+disp('Percent Complete:')
+disp(iterationCounter/nn_eff*100)
+
+disp('Time Remaining:')
+s = (nn_eff-iterationCounter)*mean(progressTracker);
+hours = floor(s/3600);
+s = s - hours*3600;
+mins = floor(s/60);
+secs = s-mins*60;
+hms = sprintf('%02d:%02d:%05.2f\n', hours, mins, secs);
+disp(hms)
+
+mark = fix(length(str)*iterationCounter/nn_eff);
+
+bar = str(1:mark);
+blank = str2(1:length(str)-mark);
+
+progress = sprintf('|%s%s|',bar,blank);
+
+disp(progress)
+%}
+
+% Only Consider Nonzero Entries
+progressTracker = progressTracker(progressTracker~=0);
+
+clc;
+disp(' ')
+
+disp('Data Point:')
+disp(iterationCounter)
+
+if length(progressTracker) >= 25
+    disp('Current Speed (ms per calculation):')
+    disp(1000*mean(progressTracker(end-24:end)))
+    nextUpdate = 1/mean(progressTracker(end-24:end));
+else
+    disp('Current Speed (ms per calculation):')
+    disp('...')
+    nextUpdate = 100;
+end
+
+disp('Average Speed (ms per calculation):')
+disp(1000*mean(progressTracker))
+
+disp('Total Time Elapsed:')
+%disp(sum(time))
+s = sum(progressTracker);
+hours = floor(s/3600);
+s = s - hours*3600;
+mins = floor(s/60);
+secs = s-mins*60;
+hms = sprintf('%02d:%02d:%05.2f\n', hours, mins, secs);
+disp(hms)
+
+
+disp('Percent Complete:')
+disp(iterationCounter/totalIterations*100)
+
+disp('Time Remaining:')
+s = (totalIterations-iterationCounter)*mean(progressTracker);
+hours = floor(s/3600);
+s = s - hours*3600;
+mins = floor(s/60);
+secs = s-mins*60;
+hms = sprintf('%02d:%02d:%05.2f\n', hours, mins, secs);
+disp(hms)
+
+mark = fix(length(str)*iterationCounter/totalIterations);
+
+bar = str(1:mark);
+blank = str2(1:length(str)-mark);
+
+progress = sprintf('|%s%s|',bar,blank);
+
+disp(progress)
+
+end

drawOrbit.m

@@ -0,0 +1,18 @@
+function drawOrbit(axHandle, pIndex,T0)
+
+tau = 2*pi;
+
+[pStruct] = planetState(pIndex,T0);
+
+nu = linspace(0,tau,361);
+r = zeros(361,3);
+
+for deg = 0:360
+    pStruct.nu = nu(deg+1);
+    pStruct = oe2sv(pStruct);
+    r(deg+1,:) = pStruct.r;
+end
+
+plot3(axHandle, r(:,1), r(:,2), r(:,3),'k')
+
+end

drawPoint_Planet.m

@@ -0,0 +1,7 @@
+function h = drawPoint_Planet(axHandle, pIndex,T,form)
+
+pStruct = planetState(pIndex,T);
+
+h = plot3(axHandle, pStruct.r(1),pStruct.r(2),pStruct.r(3),form);
+
+end