%% Simulation Project % Bradley Peckinpaugh % Section T4 % 4/15/2014 %% Main Function %the objective of this function is to import data from %a given file with information on a set of planets/moons %then create a simulation of how the gravity of each planet/moon %would effect each other with their given postions,velocitys, mass, etc function project3(FileName) if nargin ==0 %if no initial inputs are made it prompts for a txt file to be selected FileName =(uigetfile('*.txt', 'Select a file')); end out = importdata(FileName); %imports the txt file selected into matlab

FileTitle = out.textdata(1,1); %pulls out the title of the file from the text data for later use

mass =out.data(:,2); %pulls out the mass of each planet from the data file mass_scale = textscan(out.textdata{5,1},'%f'); %pulls out the mass scale factor from the text data mass = mass_scale{1}.*mass; %multiplies the mass scale factor by the mass

pos = out.data(:,3:5); %pulls out the initial position of the planets from the data file pos_scale = textscan(out.textdata{6,1},'%f'); %pulls out the position scale factor from the text data pos = pos_scale{1}.*pos; %multiplies the position scale factor by the positions

G = textscan(out.textdata{2,1},'%f'); %pulls out the gravitational constant from the text data G = G{1}; %changes the gravitational cell to a variable

vel = out.data(:,6:8); %pulls out the initial velocity data vel_scale = textscan(out.textdata{7,1},'%f'); %pulls out the velocity scale factor from the text data vel = vel_scale{1}.*vel; %multiplies the velocity scale factor by the initial velocities

colors = out.textdata(11:end,2); %pulls out the colors of the planets from the text data

max = textscan(out.textdata{8,1},'%f'); %pulls out the "max" axis limitations from the text data max = max{1}; max = pos_scale{1}.*max; %multiplies the positions scale factor by the axis limitations

a = accel(G,mass,pos); %finds the intial accelerations based on the given Gravitational constant,mass,and initial positions

dt = textscan(out.textdata{3,1},'%f'); %pulls out the change in time from the text data dt = dt{1};

path = zeros(length(mass),3,0); % creates an array of zeros for the same amount of rows as planets % and the 3 coloumns long for the path following the planets

step = textscan(out.textdata{4,1},'%f'); %pulls out the step in calculations from the text data step = step{1};

radius = out.data(:,1); %pulls out the radius size for the planets from the data

n=0; count=0; fignewton = figure; %opens a figure while ishandle(fignewton) %runs a while loop the entire time the figure is open. Ends when closed new_a = accel(G,mass,pos); %calls on the acceleration function to calc new accelerations vel = velocity(vel,new_a,a,dt); %calls on velocity function to calc new velocity based on previous calculations new_r = position(pos,vel,a,dt); %calls on position function to calc new position based on new calculations

pos = new_r; %assigns new postions to old to keep loop calculating correctly pos = centerofmass(mass,pos); %calls on the center of mass function to keep position of everything centered a = new_a; %assigns new acceleration to old to keep loop calculating correctly

n=n+dt; %adds the change in time up starting at 0 to appear in title later count =count +1; %counts the number of time the loop has run by adding 1 each time

if count==step %if the amount of iterations loop has ran is equal to step size it plots the position of the planets path(:,:,end+1) = pos ; %makes the path array equal to the position of the planets to keep line tracing path showplot(pos,path,colors,max,radius,n,FileTitle); %calls on the showplot function to plot all the planets and paths count=0; %starts the iteration count back at zero end end

end %% Acceleration Function %The objective of this function is calculate the accleration of %of each of the planets/moon given their mass and position (pos) %and the gravitational constant (G) function [a] = accel(G,mass,pos) a = zeros(length(mass),3); %creates an array of zeros for %for the starting acceleration %of the calculation for i = 1:length(mass)%runs a for loop the number of planets given for j = 1:length(mass) %runs a for loop for the same as previous if i==j %if the loops are in the same location they continue continue end r = pos(j,:)-pos(i,:); %calculates the differences %in vector position for all the %planets/moons posmag = norm(r); %finds the magnitude of distances %with respect to each other a(i,:) = a(i,:)+G.mass(j).r./posmag.^3; %uses the found calculations to plug into Newton's law of %Gravitation end

end end %% New Position Function %objective of function is to calculate the new position %of the planets based on the calculated acceleration, calculated/given %positions, the calculated/given velocity, and the given time step function [new_r] = position(pos,vel,a,dt) new_r = pos+vel.dt+.5.a.*dt.^2; %plugs in previous position (pos), velocity (vel), acceleration (a), and %time step (dt)

end

%% New Velocity Function %objective is to calculate the new velocity of the planets/moons based on %their previous veloctiy (vel), previous acceleration (a), %new acceleration (new_a), and the time step (dt) function [new_v] = velocity(vel,new_a,a,dt) new_v = vel+.5.(new_a+a).dt; %the given equation to solve for new velocity end

%% Plot Function %The objective of this function is plot all of the planets/moons %along with their paths they have traveled. function showplot(pos,path,colors,max,radius,n,FileTitle) [nr ns] = size(pos); %creates an array the same size as pos clf %clears the figure before reploting hold on %keeps both plots (planets&paths) on the same figure

for i = 1:nr %loop runs same size of pos and that plots objects

%plots the plantets with the given radius for markersize,
%the given color for markerfacecolor and markeredgecolor

h=plot(pos(i,1),pos(i,2),'o','MarkerSize',radius(i),'MarkerFaceColor',colors{i},'MarkerEdgeColor',colors{i});

%creates a string that uses the file title on one line %and the current time ran on the second line String =sprintf(' %s \n Time %.1f',FileTitle{1},n);

%sets the background color of the figure to black set(1,'color','k')

%uses the string created as the title written in white title(String,'color','w' ) ;

%sets the axis mins and max to given "max" from the file axis([-max max -max max]);

%sets x & y as the first and second % numbers of the array for paths created above x = path(i,1,:); y = path(i,2,:);

%plots the paths of the planets with the matching color

z= plot(squeeze(x),squeeze(y),colors{i});

%turns the grid on the plotting figure off grid off

%turns the axis on the plotting figure off axis off

end %updates the figure after the above for loop has ran drawnow end

%% Center of mass Function % Objective of this function is to create a center of % of mass for the planets to orbit around and be used to center %everything function pos = centerofmass(mass,pos) %calculates the x,y,and z postions for the center of mass by %taking the dot product of the mass and position and dividing %by the sum of the masses cx = dot(mass,pos(:,1))/sum(mass); cy = dot(mass,pos(:,2))/sum(mass); cz = dot(mass,pos(:,3))/sum(mass);

%calculates the differences of the x,y, and z of the postions %and the center of mass pos(:,1) = pos(:,1)-cx;
pos(:,2) = pos(:,2)-cy;
pos(:,3) = pos(:,3)-cz;
end