%{
Hello intrepid simulator. 
This is an n-body simulator. Running it is pretty easy, and it has a bunch
of features

Running it the way it is now will begin a simple simulation.
If you want to model a collision, you will need to 
1)generate a 'Set' by running a code once to make some object.
2) Storing that set with the first commented line.
3) For 2 unique objects, make 'Set' & 'Set2'
4) Turn 'crashing' to True.
5) Modify code inside the crashing-dependant conditional to suit your
needs. Right now it uses two identical objects, copied from 'Set' only.

Thansk for downloading, and be patient for the disks! Or turn maxC up.
%}


%Set = horzcat(xPos,yPos, zPos, xVel, yVel, zVel, mass , minCMtx);
clearvars -except Set Set2
close all;
crashing = false;


%----Phyical & Simulation constants
G = 5;                  % Gravitational constant
colDis = 2;             % Maximum distance of inelastic interaction
maxC = 0.0006;          % Maximum strength of inelastic interaction (distance=0)
dT = 0.0001;            % Timestep
macC = maxC*(dT/1e-3);  % Adjustment ofinteraction based on timestep


nPlanets = 1400;         % Number of objects; 500 = test; 2000 = HQ video.
temp = 200;             % Velocities of those objects- relative & modified by creation code.

%----To set up look. 'field'  wil help you zoom in and out. Bigger = out.
camtarget([0,0,0]);
set(gcf,'color','black');
set(gcf,'InvertHardcopy','off')
colormap autumn;
spotsize = 1;
viewpoint = [10 10 0];
field = 8;

% % z = 0;              %For "weaving" up and down to see the disk better
% % dz = 0.1;

%----Initiate the information vectors. The reason they are all stored as
%separate vectos rather than one big matrix is speed - you don't need to
%call the big huge matrix every time you want to do anything. In n-body
%sims like these, the number of calculations is very high. Inelegancies
%arise.

mass=zeros(nPlanets,1);

xPos=zeros(nPlanets,1);
yPos=zeros(nPlanets,1);
zPos=zeros(nPlanets,1);

xVel=zeros(nPlanets,1);
yVel=zeros(nPlanets,1);
zVel=zeros(nPlanets,1);

xFor=zeros(nPlanets,1);
yFor=zeros(nPlanets,1);
zFor=zeros(nPlanets,1);

minCMtx = zeros(nPlanets,1);

minCMtx(1:0.5*nPlanets) = maxC; % Half the particles have no inelastic interaction. This is 'dark matter".

%-----Creation loop. Currently sets up a rotating cylinder.
for createIter = 1:nPlanets
    mass(createIter) = 800;
    
    randfacTh = 360*rand();
    randfacR = .5+7*rand();
    
    xPos(createIter) = randfacR*(sind(randfacTh));
    yPos(createIter) = randfacR*(cosd(randfacTh));
    zPos(createIter) = 6*(rand()-.5);
    
    xVel(createIter) = temp*(yPos(createIter))*(randfacR/randfacR^2+0.05);
    yVel(createIter) = temp*(-xPos(createIter))*(randfacR/randfacR^2+0.05);
    zVel(createIter) = .1*(rand()-.5);
end

%Specialized code for collision simulation. Sets are generated on earlier
%runs and copied to disposable matrices & modified to set up a collision.
%Then they are copied back into the data vectors. Order doesn't matter, of course.
if crashing
    Set_cp = Set; %#ok<UNRCH>
    Set_cp(:,1:3) =  Set(:,1:3) * [1 0 0 ; 0 cosd(40) -sind(40) ; 0 sind(40) cosd(40)];
    Set_cp(:,4:6) =  Set(:,4:6) * [1 0 0 ; 0 cosd(40) -sind(40) ; 0 sind(40) cosd(40)];
    Set_cp(:,1) = Set(:,1)+12;
    Set_cp(:,2) = Set(:,2)-2;
    Set_cp(:,4) = Set(:,4)-2e2; 
    
    Set2_cp = Set;
    Set2_cp(:,1) = Set(:,1)-12;
    Set2_cp(:,2) = Set(:,2)+2;
    Set2_cp(:,4) = Set(:,4)+2e2;
    
    xPos = [Set2_cp(:,1); Set_cp(:,1)];
    yPos = [Set2_cp(:,2); Set_cp(:,2)];
    zPos = [Set2_cp(:,3); Set_cp(:,3)];
    
    xVel = [Set2_cp(:,4); Set_cp(:,4)];
    yVel = [Set2_cp(:,5); Set_cp(:,5)];
    zVel = [Set2_cp(:,6); Set_cp(:,6)];
    
    mass=[Set2_cp(:,7); Set_cp(:,7)];
    minCMtx = [Set2_cp(:,8); Set_cp(:,8)];
end


%---- Useful code for testing 2-body systems & types of interaction.
% mass(1) = 100;
% mass(2) = 100;
% 
% xPos(1) = 0;
% xPos(2) = .25;
% 
% yVel(1) = -15;
% yVel(2) = 15; 
% 
% VelMtx = zeros(1,1000);
% DisMtx = VelMtx;

%----Simulation loop.
for timeIter = 1:1e6;       
    xFor=zeros(nPlanets,1); % Forces have no memory
    yFor=zeros(nPlanets,1);
    zFor=zeros(nPlanets,1);
    
    %-----Force Calc Loop
    for Cur = 1:nPlanets-1; %B/c of Newton's 3rd, only half of the "square" of gravities needs to be calculated.
        for Tar = (Cur+1):nPlanets
            
            massCur = mass(Cur); %Strip into individual variables for speed
            massTar = mass(Tar);
            massTot = massCur+massTar;
            
            xVelCur = xVel(Cur);
            yVelCur = yVel(Cur);
            zVelCur = zVel(Cur);
            
            
            xVelTar = xVel(Tar);
            yVelTar = yVel(Tar);
            zVelTar = zVel(Tar);
            
            xDis = xPos(Tar)-xPos(Cur); % Distance calc.
            yDis = yPos(Tar)-yPos(Cur);
            zDis = zPos(Tar)-zPos(Cur);
            tDis = sqrt(xDis^2 + yDis^2 +zDis^2);
            if tDis < colDis && minCMtx(Cur)*minCMtx(Tar) ~= 0 
                % Inelastic interaction only within maximum distance, and 
                % only if neither particle is Dark.
                
                C = -maxC/colDis*tDis+maxC; 
                %Actual interaction strength, decreases linearly with distance.
                
                xRelVel = (xVelTar-xVelCur); % Relative velocities
                yRelVel = (yVelTar-yVelCur);
                zRelVel = (zVelTar-zVelCur);
                
                projmag = (xDis*xRelVel+yDis*yRelVel+zDis*zRelVel)/(tDis^2); 
                % component of relative velocities that projects onto 
                % relative position vector, i.e. proportion of velocity radial
                
                parX = projmag*xDis;
                parY = projmag*yDis;
                parZ = projmag*zDis;
                % Actual radial velocity.
                
                inelX = (massTar*parX)/(massTot);
                inelY = (massTar*parY)/(massTot);
                inelZ = (massTar*parZ)/(massTot);
                %Inelastic velocity shift. Assumes equal masses, bleh.
                
                xVel(Cur) =  xVelCur + C *inelX;
                yVel(Cur) =  yVelCur + C *inelY;
                zVel(Cur) =  zVelCur + C *inelZ;
                
                xVel(Tar) =  xVelTar - C *inelX;
                yVel(Tar) =  yVelTar - C *inelY;
                zVel(Tar) =  zVelTar - C *inelZ;
                %Partial velocity traded.
                
%----Useful code for analyzing changes in tangential velocity, looking at various dissipation methodds.              
%                 perpVel = sqrt((xRelVel-parX)^2+(yRelVel-parY)^2+(zRelVel-parZ)^2);
%                 VelMtx(timeIter) = perpVel;
%                 DisMtx(timeIter) = tDis;                
                
%----Pure inelastic collisions. Tends to collapse to singularity because of
%    angular momentum loss.
% % %                 inelX = (massCur*xVelCur + massTar*xVelTar)/(massTot);
% % %                 inelY = (massCur*yVelCur + massTar*yVelTar)/(massTot);
% % %                 inelZ = (massCur*zVelCur + massTar*zVelTar)/(massTot);
% % %                 
% % %                 xVel(Cur) =  C*xVelCur +(1-C)*inelX;
% % %                 yVel(Cur) =  C*yVelCur +(1-C)*inelY;
% % %                 zVel(Cur) =  C*zVelCur +(1-C)*inelZ;
% % %                 
% % %                 xVel(Tar) =  C*xVelTar +(1-C)*inelX;
% % %                 yVel(Tar) =  C*yVelTar +(1-C)*inelY;
% % %                 zVel(Tar) =  C*zVelTar +(1-C)*inelZ;
                
            end
                
                Force = (G*massCur*massTar)*  (1/ (tDis^2)); %Universal Gravitation
                if Force > 1e6;                 %Flatten force up close to remove artifacts of simulation.
                    Force = 1e6;
                end
                xForCalc = xDis/tDis * Force;
                yForCalc = yDis/tDis * Force;
                zForCalc = zDis/tDis * Force;
                
                xFor(Cur) = xFor(Cur)+ xForCalc ;
                yFor(Cur) = yFor(Cur)+ yForCalc ;
                zFor(Cur) = zFor(Cur)+ zForCalc ;
                
                xFor(Tar) = xFor(Tar)- xForCalc ;
                yFor(Tar) = yFor(Tar)- yForCalc ;
                zFor(Tar) = zFor(Tar)- zForCalc ;
                
        end
    end
    
    %---- Integration loop
    for sumIter = 1:nPlanets        
        interMass = mass(sumIter);
        
        xVel(sumIter) = xVel(sumIter) + dT*(xFor(sumIter)/interMass);
        xPos(sumIter) = xPos(sumIter) + dT*xVel(sumIter);
        
        yVel(sumIter) = yVel(sumIter) + dT*(yFor(sumIter)/interMass);
        yPos(sumIter) = yPos(sumIter) + dT*yVel(sumIter);
        
        zVel(sumIter) = zVel(sumIter) + dT*(zFor(sumIter)/interMass);
        zPos(sumIter) = zPos(sumIter) + dT*zVel(sumIter);
        
    end
    
    % Visualize as scatterplot. Color determined by interaction strength.
    subplot(3,2,1);

    scatter3(xPos,yPos,zPos, spotsize, minCMtx, 'filled');
    axis([-field,field,-field,field, -field,field]);
    axis off;
    campos([10 0 0]);
    
    subplot(3,2,2);
    scatter3(xPos,yPos,zPos, spotsize, minCMtx, 'filled');
    axis([-field,field,-field,field, -field,field]);
    axis off;
    campos([0 0 10]);
    
    
    subplot(3,2,[3,6]);
    scatter3(xPos,yPos,zPos, spotsize, minCMtx, 'filled');
    axis([-field,field,-field,field, -field,field]);
    axis off;
    campos([10 10 20]);
    
    
    % Camera that goes up & Down.
% %     campos([10 10 z]);
% %     z = z + dz;
% %     if mod(timeIter,25) == 0
% %        dz = -dz; 
% %     end
    
%---- Make images!
%     print(gcf,'-dpng','-r120',sprintf('disk%03d.png',timeIter))
    

    pause(.00000001)   %Don't know why, but you need this.
end

%-----Helpful charting tools if you're nutty about those interactions.

% subplot(1,3,1);
% plot(VelMtx);
% axis([0,4000,0,40])
% title('Relative Velocities');
% subplot(1,3,2);
% plot(DisMtx,'r');
% axis([0,4000,0,.5])
% title('Relative Positions');
% subplot(1,3,3);
% plot(VelMtx.*DisMtx);
% axis([0,4000,0,10])
% title('Angular momentum');