SHARE
TWEET

Barry-centric initial conditions KSP

a guest Sep 21st, 2014 186 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. function [X, V, body] = initialConditions(barry)
  2. % When barry is True then position and velocity will be done relative to
  3. % the barry center. Otherwise the barry center is assumed to be equal to
  4. % the center of mass of the primary body, for instance the sun for the planets.
  5.     A = zeros(17,8);
  6.     %          1                   2                       3       4     5   6   7    8
  7.     %          MU                  SMA                     ECC     INC   LPE LAN MNA  REF    Name
  8.     A( 1,:) = [1.17233279483249e18 0                       0       0       0   0 0.00  0]; % Sun
  9.     A( 2,:) = [1.68609378654509e11 5.2631383040e09 0.20000 7.000  15  70 3.14  1]; % Moho
  10.     A( 3,:) = [8.17173022921085e12 9.8326845440e09 0.01000 2.100   0  15 3.14  1]; % Eve
  11.     A( 4,:) = [8.28944981471635e06 3.1500000000e07 0.55000 12.00  10  80 0.90  3]; % Gilly
  12.     A( 5,:) = [3.53160000000000e12 1.3599840256e10 0.00000 0.000   0   0 3.14  1]; % Kerbin
  13.     A( 6,:) = [6.51383975207806e10 1.2000000000e07 0.00000 0.000   0   0 0.90  5]; % Mun
  14.     A( 7,:) = [1.76580002631247e09 4.7000000000e07 0.00000 6.000  38  78 1.70  5]; % Minmus
  15.     A( 8,:) = [3.01363211975098e11 2.0726155264e10 0.05100 0.060   0 135 3.14  1]; % Duna
  16.     A( 9,:) = [1.85683685731441e10 3.2000000000e06 0.03000 0.200   0   0 1.70  8]; % Ike
  17.     A(10,:) = [2.14844886000000e10 4.0839348203e10 0.14500 5.000  90 280 3.14  1]; % Dres
  18.     A(11,:) = [2.82528004209995e14 6.8773560320e10 0.05000 1.304   0  52 0.10  1]; % Jool
  19.     A(12,:) = [1.96200002923608e12 2.7184000000e07 0.00000 0.000   0   0 3.14 11]; % Laythe
  20.     A(13,:) = [2.07481499473751e11 4.3152000000e07 0.00000 0.000   0   0 0.90 11]; % Vall
  21.     A(14,:) = [2.82528004209995e12 6.8500000000e07 0.00000 0.025   0   0 3.14 11]; % Tylo
  22.     A(15,:) = [2.48683494441491e09 1.2850000000e08 0.23500 15.00  25  10 0.90 11]; % Bop
  23.     A(16,:) = [7.21702080000000e08 1.7989000000e08 0.17085 4.250  15   2 0.90 11]; % Pol
  24.     A(17,:) = [7.44108145270496e10 9.0118820000e10 0.26000 6.150 260  50 3.14  1]; % Eeloo
  25.  
  26.     M = A(:,1);
  27.    
  28.     A(13,7) = 3.14 - pi; % Change relative initial phase of Vall to something maybe more stable
  29.    
  30.     A(13,2) = A(12,2) * (2 * (M(11) + M(12)) / (M(11) + M(13)))^(2/3); % Orbital period twice Laythes
  31.     A(14,2) = A(12,2) * (4 * (M(11) + M(12)) / (M(11) + M(14)))^(2/3); % Orbital period four times Laythes
  32.  
  33.     name = {'Sun' 'Moho' 'Eve' 'Gilly' 'Kerbin' 'Mun' 'Minmus' 'Duna' 'Ike' ...
  34.             'Dres' 'Jool' 'Laythe' 'Vall' 'Tylo' 'Bop' 'Pol' 'Eeloo' 'Asteroid 1' 'Asteroid 2'};
  35.     color = [1.0000    0.9216    0.3882;
  36.              0.7294    0.5765    0.4353;
  37.              0.4784    0.1686    0.8235;
  38.              0.3843    0.2549    0.1961;
  39.              0.5059    0.8196    0.7765;
  40.              0.3255    0.3412    0.3765;
  41.              0.1961    0.3882    0.4431;
  42.              0.8157    0.3882    0.2588;
  43.              0.3176    0.3451    0.4078;
  44.              0.2824    0.2235    0.1608;
  45.              0.4510    0.6471    0.1137;
  46.              0.1647    0.2431    0.5255;
  47.              0.3294    0.4588    0.5294;
  48.              0.7176    0.6000    0.6000;
  49.              0.4039    0.3412    0.2627;
  50.              0.4627    0.4824    0.3294;
  51.              0.3255    0.3255    0.3255];
  52.  
  53.     body.m = A(:,1);
  54.     body.name = name;
  55.     N = numel(body.m);
  56.     for i = 1 : N
  57.         body.col{i} = color(i,:);
  58.     end
  59.    
  60.     % reference direction: x-axis
  61.     % reference plane:     x-y plane
  62.     pos = zeros(3, N);
  63.     vel = zeros(3, N);
  64.     for i = 2 : N
  65.         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
  66.         eAN = [cosd(A(i,6)); sind(A(i,6)); 0]; % Unit vector pointing to the ascending node
  67.         temp = cross(n, eAN); % Unit vector perpendicular to the two previous vectors
  68.         angle = A(i,5)*pi/180 + A(i,7); % Total angle away from ascending node
  69.         er = cos(angle) * eAN + sin(angle) * temp; % Unit vector of the position
  70.         et = cos(angle) * temp - sin(angle) * eAN; % Unit vector perpendicular to the position and in the orbital plane
  71.         pos(:,i) = A(i,2) * (1 - A(i,3)^2) / (1 + A(i,3) * cos(A(i,7))) * er; % Position relative to barry center
  72.         mu = barry * M(A(i,8))^3 / (M(A(i,8)) + M(i))^2 + ~barry * M(A(i,8)); % Effective gravitaional parameter
  73.         vel(:,i) = sqrt(mu / (A(i,2) * (1 - A(i,3)^2))) * ...
  74.                (A(i,3) * sin(A(i,7)) * er + (1 + A(i,3) * cos(A(i,7))) * et); % Velocity relative to barry center
  75.     end
  76.    
  77.     moons = [4 6 7 9 12:16];
  78.     X = zeros(3, N);
  79.     V = zeros(3, N);
  80.      
  81.     for j = moons
  82.         X(:,j) = pos(:,j) + pos(:,A(j,8));
  83.         V(:,j) = vel(:,j) + vel(:,A(j,8));
  84.     end
  85.     if barry
  86.         X(:,2) = pos(:,2);
  87.         V(:,2) = vel(:,2);
  88.        
  89.         X(:,3) = pos(:,3) - A(4,1) * pos(:,4) / A(3,1);
  90.         V(:,3) = vel(:,3) - A(4,1) * vel(:,4) / A(3,1);
  91.        
  92.         X(:,5) = pos(:,5) - (A([6,7],1)' * pos(:,[6,7])' / A(5,1))';
  93.         V(:,5) = vel(:,5) - (A([6,7],1)' * vel(:,[6,7])' / A(5,1))';
  94.        
  95.         X(:,8) = pos(:,8) - A(9,1) * pos(:,9) / A(8,1);
  96.         V(:,8) = vel(:,8) - A(9,1) * vel(:,9) / A(8,1);
  97.        
  98.         X(:,10) = pos(:,10);
  99.         V(:,10) = vel(:,10);
  100.        
  101.         X(:,11) = pos(:,11) - (A(12:16,1)' * pos(:,12:16)' / A(11,1))';
  102.         V(:,11) = vel(:,11) - (A(12:16,1)' * vel(:,12:16)' / A(11,1))';
  103.        
  104.         X(:,17) = pos(:,17);
  105.         V(:,17) = vel(:,17);
  106.        
  107.         X(:,1) = pos(:,1) - (A(2:17,1)' * X(:,2:17)' / A(1,1))';
  108.         V(:,1) = vel(:,1) - (A(2:17,1)' * V(:,2:17)' / A(1,1))';
  109.     else
  110.         X(:,[1 2 3 5 8 10 11 17]) = pos(:,[1 2 3 5 8 10 11 17]);
  111.         V(:,[1 2 3 5 8 10 11 17]) = vel(:,[1 2 3 5 8 10 11 17]);
  112.     end
  113.  
  114. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top