Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- function [X, V, body] = initialConditions(barry)
- % When barry is True then position and velocity will be done relative to
- % the barry center. Otherwise the barry center is assumed to be equal to
- % the center of mass of the primary body, for instance the sun for the planets.
- A = zeros(17,8);
- % 1 2 3 4 5 6 7 8
- % MU SMA ECC INC LPE LAN MNA REF Name
- A( 1,:) = [1.17233279483249e18 0 0 0 0 0 0.00 0]; % Sun
- A( 2,:) = [1.68609378654509e11 5.2631383040e09 0.20000 7.000 15 70 3.14 1]; % Moho
- A( 3,:) = [8.17173022921085e12 9.8326845440e09 0.01000 2.100 0 15 3.14 1]; % Eve
- A( 4,:) = [8.28944981471635e06 3.1500000000e07 0.55000 12.00 10 80 0.90 3]; % Gilly
- A( 5,:) = [3.53160000000000e12 1.3599840256e10 0.00000 0.000 0 0 3.14 1]; % Kerbin
- A( 6,:) = [6.51383975207806e10 1.2000000000e07 0.00000 0.000 0 0 0.90 5]; % Mun
- A( 7,:) = [1.76580002631247e09 4.7000000000e07 0.00000 6.000 38 78 1.70 5]; % Minmus
- A( 8,:) = [3.01363211975098e11 2.0726155264e10 0.05100 0.060 0 135 3.14 1]; % Duna
- A( 9,:) = [1.85683685731441e10 3.2000000000e06 0.03000 0.200 0 0 1.70 8]; % Ike
- A(10,:) = [2.14844886000000e10 4.0839348203e10 0.14500 5.000 90 280 3.14 1]; % Dres
- A(11,:) = [2.82528004209995e14 6.8773560320e10 0.05000 1.304 0 52 0.10 1]; % Jool
- A(12,:) = [1.96200002923608e12 2.7184000000e07 0.00000 0.000 0 0 3.14 11]; % Laythe
- A(13,:) = [2.07481499473751e11 4.3152000000e07 0.00000 0.000 0 0 0.90 11]; % Vall
- A(14,:) = [2.82528004209995e12 6.8500000000e07 0.00000 0.025 0 0 3.14 11]; % Tylo
- A(15,:) = [2.48683494441491e09 1.2850000000e08 0.23500 15.00 25 10 0.90 11]; % Bop
- A(16,:) = [7.21702080000000e08 1.7989000000e08 0.17085 4.250 15 2 0.90 11]; % Pol
- A(17,:) = [7.44108145270496e10 9.0118820000e10 0.26000 6.150 260 50 3.14 1]; % Eeloo
- M = A(:,1);
- A(13,7) = 3.14 - pi; % Change relative initial phase of Vall to something maybe more stable
- A(13,2) = A(12,2) * (2 * (M(11) + M(12)) / (M(11) + M(13)))^(2/3); % Orbital period twice Laythes
- A(14,2) = A(12,2) * (4 * (M(11) + M(12)) / (M(11) + M(14)))^(2/3); % Orbital period four times Laythes
- name = {'Sun' 'Moho' 'Eve' 'Gilly' 'Kerbin' 'Mun' 'Minmus' 'Duna' 'Ike' ...
- 'Dres' 'Jool' 'Laythe' 'Vall' 'Tylo' 'Bop' 'Pol' 'Eeloo' 'Asteroid 1' 'Asteroid 2'};
- color = [1.0000 0.9216 0.3882;
- 0.7294 0.5765 0.4353;
- 0.4784 0.1686 0.8235;
- 0.3843 0.2549 0.1961;
- 0.5059 0.8196 0.7765;
- 0.3255 0.3412 0.3765;
- 0.1961 0.3882 0.4431;
- 0.8157 0.3882 0.2588;
- 0.3176 0.3451 0.4078;
- 0.2824 0.2235 0.1608;
- 0.4510 0.6471 0.1137;
- 0.1647 0.2431 0.5255;
- 0.3294 0.4588 0.5294;
- 0.7176 0.6000 0.6000;
- 0.4039 0.3412 0.2627;
- 0.4627 0.4824 0.3294;
- 0.3255 0.3255 0.3255];
- body.m = A(:,1);
- body.name = name;
- N = numel(body.m);
- for i = 1 : N
- body.col{i} = color(i,:);
- end
- % reference direction: x-axis
- % reference plane: x-y plane
- pos = zeros(3, N);
- vel = zeros(3, N);
- for i = 2 : N
- n = [sind(A(i,4)) * sind(A(i,6)); -sind(A(i,4)) * cosd(A(i,6)); cosd(A(i,4))]; % Normal of the orbital plane
- eAN = [cosd(A(i,6)); sind(A(i,6)); 0]; % Unit vector pointing to the ascending node
- temp = cross(n, eAN); % Unit vector perpendicular to the two previous vectors
- angle = A(i,5)*pi/180 + A(i,7); % Total angle away from ascending node
- er = cos(angle) * eAN + sin(angle) * temp; % Unit vector of the position
- et = cos(angle) * temp - sin(angle) * eAN; % Unit vector perpendicular to the position and in the orbital plane
- pos(:,i) = A(i,2) * (1 - A(i,3)^2) / (1 + A(i,3) * cos(A(i,7))) * er; % Position relative to barry center
- mu = barry * M(A(i,8))^3 / (M(A(i,8)) + M(i))^2 + ~barry * M(A(i,8)); % Effective gravitaional parameter
- vel(:,i) = sqrt(mu / (A(i,2) * (1 - A(i,3)^2))) * ...
- (A(i,3) * sin(A(i,7)) * er + (1 + A(i,3) * cos(A(i,7))) * et); % Velocity relative to barry center
- end
- moons = [4 6 7 9 12:16];
- X = zeros(3, N);
- V = zeros(3, N);
- for j = moons
- X(:,j) = pos(:,j) + pos(:,A(j,8));
- V(:,j) = vel(:,j) + vel(:,A(j,8));
- end
- if barry
- X(:,2) = pos(:,2);
- V(:,2) = vel(:,2);
- X(:,3) = pos(:,3) - A(4,1) * pos(:,4) / A(3,1);
- V(:,3) = vel(:,3) - A(4,1) * vel(:,4) / A(3,1);
- X(:,5) = pos(:,5) - (A([6,7],1)' * pos(:,[6,7])' / A(5,1))';
- V(:,5) = vel(:,5) - (A([6,7],1)' * vel(:,[6,7])' / A(5,1))';
- X(:,8) = pos(:,8) - A(9,1) * pos(:,9) / A(8,1);
- V(:,8) = vel(:,8) - A(9,1) * vel(:,9) / A(8,1);
- X(:,10) = pos(:,10);
- V(:,10) = vel(:,10);
- X(:,11) = pos(:,11) - (A(12:16,1)' * pos(:,12:16)' / A(11,1))';
- V(:,11) = vel(:,11) - (A(12:16,1)' * vel(:,12:16)' / A(11,1))';
- X(:,17) = pos(:,17);
- V(:,17) = vel(:,17);
- X(:,1) = pos(:,1) - (A(2:17,1)' * X(:,2:17)' / A(1,1))';
- V(:,1) = vel(:,1) - (A(2:17,1)' * V(:,2:17)' / A(1,1))';
- else
- X(:,[1 2 3 5 8 10 11 17]) = pos(:,[1 2 3 5 8 10 11 17]);
- V(:,[1 2 3 5 8 10 11 17]) = vel(:,[1 2 3 5 8 10 11 17]);
- end
- end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement