Barry-centric initial conditions KSP

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