"""SolarSys1.py Simulates solar system formation """ # Defines the output file directory. OutputFile = "/home/tom/Desktop/Test1.txt" OutputFile2 = "/home/tom/Desktop/Test2.txt" BigOutput = "/home/tom/Desktop/PHYSICS/BIGOUTPUT.txt" # import necessary modules import math import random import numpy import scipy import Gnuplot # Set the Gravitational constant. Gravity = 100. # Set the density constant. Density = 10 # Set the number of initial particles ParticleMax = 1000 ParticleMass = 1. ParticleNumber = ParticleMax # Set the size of the field FieldSize = 1000 # Set the time limit and increment TimeMax = 1510 TimeSkip = 1. dT = 1. # Calculates the Escape Velocity. # Based on equation ESCVEL = sqrt(2*G*M/r) # M is the total mass of the particle cloud. # r is assumed to be 1/2 Fieldsize TotalMass = ParticleMass*ParticleMax #CloudRadius = FieldSize/2. ESCVEL = (math.sqrt(2.*Gravity*TotalMass/FieldSize)) # Define the cube size for splitting the field CubeSize = 1 CubeSide = FieldSize/CubeSize CubeNumber = math.pow(CubeSide,3) DeadParticle = 0 # Defines the Particle lists ParticleID = range(0,ParticleMax+1) XPOS = range(0,ParticleMax+1) YPOS = range(0,ParticleMax+1) RPOS = range(0,ParticleMax+1) ZPOS = range(0,ParticleMax+1) MASS = range(0,ParticleMax+1) XACCEL = range(0,ParticleMax+1) YACCEL = range(0,ParticleMax+1) ZACCEL = range(0,ParticleMax+1) XVELO = range(0,ParticleMax+1) YVELO = range(0,ParticleMax+1) ZVELO = range(0,ParticleMax+1) RADIUS = range(0,ParticleMax+1) # A range for randomly determining positive or negative. List L consists of 1 and -1. L = range(-1,2,2) # Loads the Particle lists with the initial particle data. for Particle in range(0,ParticleNumber): # Assigns an ID number to each particle. ParticleID[Particle] = 0 # Randomly determines X, Y, and radius R. XPOS[Particle] = (2*FieldSize +1)*random.random() - FieldSize YPOS[Particle] = (2*FieldSize +1)*random.random() - FieldSize ZPOS[Particle] = (2*FieldSize +1)*random.random() - FieldSize # If X^2 + Y^2 > R^2, X, Y, and R are redetermined until X^2 + Y^2 <= R^2. while XPOS[Particle]*XPOS[Particle] + YPOS[Particle]*YPOS[Particle] + ZPOS[Particle]*ZPOS[Particle] > 1000000.: XPOS[Particle] = (2*FieldSize +1)*random.random() - FieldSize YPOS[Particle] = (2*FieldSize +1)*random.random() - FieldSize ZPOS[Particle] = (2*FieldSize +1)*random.random() - FieldSize # Assigns a mass value to each particle according the predetermined particle mass. MASS[Particle] = (ParticleMass) # Sets initial force and acceleration for each particle to zero. XACCEL[Particle] = 0. YACCEL[Particle] = 0. ZACCEL[Particle] = 0. # Finds the escape velocity of each particle # Sets the velocity of each variable to a random value. (1.4) XVELO[Particle] = random.choice(L)*5*random.random() YVELO[Particle] = random.choice(L)*5*random.random() ZVELO[Particle] = random.choice(L)*5*random.random() # Sets the radius of each particle according to the Density term. RADIUS[Particle] = math.pow(MASS[Particle],.333)*Density # Checks for initial collisions. for Particle in range(0,ParticleNumber - 1): if ParticleID[Particle] == 0: for Other in range(Particle + 1, ParticleNumber,): if ParticleID[Other] == 0: # Calculates the distance between the particles. XDist = XPOS[Particle] - XPOS[Other] YDist = YPOS[Particle] - YPOS[Other] ZDist = ZPOS[Particle] - ZPOS[Other] Distance = math.sqrt(math.pow(XDist,2) + math.pow(YDist,2) + math.pow(ZDist,2)) # Checks if the distance is less than the collision distance, or the sum of the particles' radii. if Distance <= RADIUS[Particle] + RADIUS[Other]: # Averages the postions and adds the velocities of each particle. MASS2 = MASS[Particle] + MASS[Other] XPOS[Particle] = (XPOS[Particle]*MASS[Particle] + XPOS[Other]*MASS[Other])/MASS2 YPOS[Particle] = (YPOS[Particle]*MASS[Particle] + YPOS[Other]*MASS[Other])/MASS2 ZPOS[Particle] = (ZPOS[Particle]*MASS[Particle] + ZPOS[Other]*MASS[Other])/MASS2 XVELO[Particle] = (XVELO[Particle]*MASS[Particle] + XVELO[Other]*MASS[Other])/MASS2 YVELO[Particle] = (YVELO[Particle]*MASS[Particle] + YVELO[Other]*MASS[Other])/MASS2 ZVELO[Particle] = (ZVELO[Particle]*MASS[Particle] + ZVELO[Other]*MASS[Other])/MASS2 # Adds the masses of the particles. MASS[Particle] = MASS[Particle] + MASS[Other] # Calculates the new radius of the particle. RADIUS[Particle] = math.pow(MASS[Particle],.333)*Density # Removes the other particle from future calculations, keeping the particle with lower ID number. # The removed particle is given a particle ID that is above the range involved in calculations. MASS[Other] = 0 RADIUS[Other] = 0 XPOS[Other] = FieldSize YPOS[Other] = FieldSize ZPOS[Other] = FieldSize ParticleID[Other] = 1 #ParticleNumber = ParticleNumber - 1 DeadParticle = DeadParticle + 1 # Finds the center of Mass. XSUM1 = 0. YSUM1 = 0. ZSUM1 = 0. MASSSUM = 0. for Particle in range(0,ParticleNumber): XSUM1 = (XSUM1 + XPOS[Particle]*MASS[Particle]) YSUM1 = (YSUM1 + YPOS[Particle]*MASS[Particle]) ZSUM1 = (ZSUM1 + ZPOS[Particle]*MASS[Particle]) MASSSUM = (MASSSUM + MASS[Particle]) XCENTER = (XSUM1/MASSSUM) YCENTER = (YSUM1/MASSSUM) ZCENTER = (ZSUM1/MASSSUM) #print "The center of mass is ", XCENTER, YCENTER, ZCENTER,"." # Find the Average Momentum XSUM2 = 0. YSUM2 = 0. ZSUM2 = 0. for Particle in range(0,ParticleNumber): XSUM2 = (XSUM2 + XVELO[Particle]*MASS[Particle]) YSUM2 = (YSUM2 + YVELO[Particle]*MASS[Particle]) ZSUM2 = (ZSUM2 + ZVELO[Particle]*MASS[Particle]) XSUM2 = XSUM2/ParticleNumber YSUM2 = YSUM2/ParticleNumber ZSUM2 = ZSUM2/ParticleNumber #print "The average momentum is", XSUM2, YSUM2, ZSUM2,"." # Find the Average Angular Momentum XSUM3 = 0. YSUM3 = 0. ZSUM3 = 0. for Particle in range(0,ParticleNumber): XDIST = XPOS[Particle] YDIST = YPOS[Particle] ZDIST = ZPOS[Particle] XSUM3 = (XSUM3 + (YDIST*ZVELO[Particle] - ZDIST*YVELO[Particle])*MASS[Particle]) YSUM3 = (YSUM3 + (ZDIST*XVELO[Particle] - XDIST*ZVELO[Particle])*MASS[Particle]) ZSUM3 = (ZSUM3 + (XDIST*YVELO[Particle] - YDIST*XVELO[Particle])*MASS[Particle]) XSUM3 = XSUM3/ParticleNumber YSUM3 = YSUM3/ParticleNumber ZSUM3 = ZSUM3/ParticleNumber #print "The average angular momentum is", XSUM3, YSUM3, ZSUM3,"." # Finds the RMS of the Z varbiable ZSUM4 = 0. MASSSUM2 = 0. for Particle in range(0,ParticleNumber): ZSUM4 = ZSUM4 + MASS[Particle]*math.pow(ZPOS[Particle],2) MASSSUM2 = MASSSUM2 + MASS[Particle] ZRMS = math.sqrt(ZSUM4/MASSSUM2) #print "The RMS value of Z is ", ZRMS,"." # Finds the RMS of the velocity VELSUM = 0. for Particle in range(0,ParticleNumber): VELSUM = VELSUM + (XVELO[Particle]*XVELO[Particle]+YVELO[Particle]*YVELO[Particle]+ZVELO[Particle]*ZVELO[Particle]) VELRMS = math.sqrt(VELSUM/ParticleNumber) #print "The RMS value of velocity is ", VELRMS,"." print " " print " " # Corrects the position, momentum, and angular momentum of each particle. for Particle in range(0,ParticleNumber): if ParticleID[Particle] == 0: XPOS[Particle] = XPOS[Particle] - XCENTER YPOS[Particle] = YPOS[Particle] - YCENTER ZPOS[Particle] = ZPOS[Particle] - ZCENTER XVELO[Particle] = XVELO[Particle] - XSUM2/MASS[Particle] YVELO[Particle] = YVELO[Particle] - YSUM2/MASS[Particle] ZVELO[Particle] = ZVELO[Particle] - ZSUM2/MASS[Particle] # Calculates the force on each particle at each step open(OutputFile2,"w").write("Time\tPx\tPy\tPz\tLx\tLy\tLz\tParticles\n") open(BigOutput,"w").write("Time\tParticle\tLx/L\tLy/L\tLz/L\tX\tY\tZ\tMass\n") TimeCounter = 10. for TimeStep in range(0,TimeMax,TimeSkip): if TimeCounter != 10: TimeCounter = TimeCounter + 1 for Particle in range(0,ParticleNumber): XACCEL[Particle] = 0. YACCEL[Particle] = 0. ZACCEL[Particle] = 0. for Particle in range(0,ParticleNumber - 1): if ParticleID[Particle] == 0: for Other in range(Particle + 1,ParticleNumber): if ParticleID[Other] == 0: XDist = (XPOS[Particle] - XPOS[Other]) YDist = (YPOS[Particle] - YPOS[Other]) ZDist = (ZPOS[Particle] - ZPOS[Other]) R3 = math.pow((XDist*XDist+YDist*YDist+ZDist*ZDist),1.5) XACCEL[Particle] = XACCEL[Particle] + (-Gravity*MASS[Other]*XDist/R3) YACCEL[Particle] = YACCEL[Particle] + (-Gravity*MASS[Other]*YDist/R3) ZACCEL[Particle] = ZACCEL[Particle] + (-Gravity*MASS[Other]*ZDist/R3) XACCEL[Other] = XACCEL[Other] + (Gravity*MASS[Particle]*XDist/R3) YACCEL[Other] = YACCEL[Other] + (Gravity*MASS[Particle]*YDist/R3) ZACCEL[Other] = ZACCEL[Other] + (Gravity*MASS[Particle]*ZDist/R3) #print Particle, Other, XDist, YDist, ZDist for Particle in range(0,ParticleNumber): if ParticleID[Particle] == 0: #print "Forces", XACCEL[Particle], YACCEL[Particle], ZACCEL[Particle] XVELO[Particle] = XVELO[Particle] + XACCEL[Particle]*TimeSkip*dT YVELO[Particle] = YVELO[Particle] + YACCEL[Particle]*TimeSkip*dT ZVELO[Particle] = ZVELO[Particle] + ZACCEL[Particle]*TimeSkip*dT XPOS[Particle] = XPOS[Particle] + XVELO[Particle]*TimeSkip*dT YPOS[Particle] = YPOS[Particle] + YVELO[Particle]*TimeSkip*dT ZPOS[Particle] = ZPOS[Particle] + ZVELO[Particle]*TimeSkip*dT # Check for collisions. for Particle in range(0,ParticleNumber - 1): if ParticleID[Particle] == 0: for Other in range(Particle + 1, ParticleNumber,): if ParticleID[Other] == 0: # Calculates the distance between the particles. XDist = XPOS[Particle] - XPOS[Other] YDist = YPOS[Particle] - YPOS[Other] ZDist = ZPOS[Particle] - ZPOS[Other] Distance = math.sqrt(math.pow(XDist,2) + math.pow(YDist,2) + math.pow(ZDist,2)) # Checks if the distance is less than the collision distance, or the sum of the particles' radii. if Distance <= RADIUS[Particle] + RADIUS[Other]: # Averages the postions and adds the velocities of each particle. MASS2 = MASS[Particle] + MASS[Other] XPOS[Particle] = (XPOS[Particle]*MASS[Particle] + XPOS[Other]*MASS[Other])/MASS2 YPOS[Particle] = (YPOS[Particle]*MASS[Particle] + YPOS[Other]*MASS[Other])/MASS2 ZPOS[Particle] = (ZPOS[Particle]*MASS[Particle] + ZPOS[Other]*MASS[Other])/MASS2 XVELO[Particle] = (XVELO[Particle]*MASS[Particle] + XVELO[Other]*MASS[Other])/MASS2 YVELO[Particle] = (YVELO[Particle]*MASS[Particle] + YVELO[Other]*MASS[Other])/MASS2 ZVELO[Particle] = (ZVELO[Particle]*MASS[Particle] + ZVELO[Other]*MASS[Other])/MASS2 # Adds the masses of the particles. MASS[Particle] = MASS[Particle] + MASS[Other] # Calculates the new radius of the particle. RADIUS[Particle] = math.pow(MASS[Particle],.333)*Density # Removes the other particle from future calculations, keeping the particle with lower ID number. # The removed particle is given a particle ID that is above the range involved in calculations. #if Other != ParticleNumber: #print "collision", ParticleID[Particle],ParticleID[Other] #XPOS[Other] = XPOS[ParticleNumber] #YPOS[Other] = YPOS[ParticleNumber] #ZPOS[Other] = ZPOS[ParticleNumber] #XVELO[Other] = XVELO[ParticleNumber] #YVELO[Other] = YVELO[ParticleNumber] #ZVELO[Other] = ZVELO[ParticleNumber] MASS[Other] = 0 RADIUS[Other] = 0 XPOS[Other] = FieldSize YPOS[Other] = FieldSize ZPOS[Other] = FieldSize ParticleID[Other] = 1 #ParticleNumber = ParticleNumber - 1 DeadParticle = DeadParticle + 1 # Begins printing the output. # print "At time", TimeStep,"," # print "There are", ParticleNumber,"particles at distance", Distance,"," # Finds the center of Mass. XSUM1 = 0. YSUM1 = 0. ZSUM1 = 0. MASSSUM = 0. for Particle in range(0,ParticleNumber): XSUM1 = (XSUM1 + XPOS[Particle]*MASS[Particle]) YSUM1 = (YSUM1 + YPOS[Particle]*MASS[Particle]) ZSUM1 = (ZSUM1 + ZPOS[Particle]*MASS[Particle]) MASSSUM = (MASSSUM + MASS[Particle]) XCENTER = (XSUM1/MASSSUM) YCENTER = (YSUM1/MASSSUM) ZCENTER = (ZSUM1/MASSSUM) # print "The center of mass is ", XCENTER, YCENTER, ZCENTER,"," # Find the Average Momentum XSUM2 = 0. YSUM2 = 0. ZSUM2 = 0. for Particle in range(0,ParticleNumber): XSUM2 = (XSUM2 + XVELO[Particle]*MASS[Particle]) YSUM2 = (YSUM2 + YVELO[Particle]*MASS[Particle]) ZSUM2 = (ZSUM2 + ZVELO[Particle]*MASS[Particle]) XSUM2 = XSUM2/ParticleNumber YSUM2 = YSUM2/ParticleNumber ZSUM2 = ZSUM2/ParticleNumber # print "The average momentum is", XSUM2, YSUM2, ZSUM2,"," # Find the Average Angular Momentum XSUM3 = 0. YSUM3 = 0. ZSUM3 = 0. for Particle in range(0,ParticleNumber): XDIST = XPOS[Particle] YDIST = YPOS[Particle] ZDIST = ZPOS[Particle] XSUM3 = (XSUM3 + (YDIST*ZVELO[Particle] - ZDIST*YVELO[Particle])*MASS[Particle]) YSUM3 = (YSUM3 + (ZDIST*XVELO[Particle] - XDIST*ZVELO[Particle])*MASS[Particle]) ZSUM3 = (ZSUM3 + (XDIST*YVELO[Particle] - YDIST*XVELO[Particle])*MASS[Particle]) XSUM3 = XSUM3/ParticleNumber YSUM3 = YSUM3/ParticleNumber ZSUM3 = ZSUM3/ParticleNumber # print "The average angular momentum is", XSUM3, YSUM3, ZSUM3,"," # Finds the RMS of the Z varbiable ZSUM4 = 0. MASSSUM2 = 0. for Particle in range(0,ParticleNumber): ZSUM4 = ZSUM4 + MASS[Particle]*math.pow(ZPOS[Particle],2) MASSSUM2 = MASSSUM2 + MASS[Particle] ZRMS = math.sqrt(ZSUM4/MASSSUM2) # print "The RMS value of Z is ", ZRMS,"." # Finds the RMS of the velocity VELSUM = 0. for Particle in range(0,ParticleNumber): VELSUM = VELSUM + (XVELO[Particle]*XVELO[Particle]+YVELO[Particle]*YVELO[Particle]+ZVELO[Particle]*ZVELO[Particle]) VELRMS = math.sqrt(VELSUM/ParticleNumber) # print "The RMS value of velocity is ", VELRMS,"." #print "velocity of particle 10 is ", XVELO[1], XVELO[2], XVELO[3], TimeStep #print "X positions", XPOS[1], XPOS[2], XPOS[3] # print " " # print " " #Calculate the RMS radius for each time step. SUM = 0 MASSSUM4 = 0 for Particle in range(0,ParticleNumber): RADIUS2 = (XPOS[Particle]*XPOS[Particle]+YPOS[Particle]*YPOS[Particle]+ZPOS[Particle]*ZPOS[Particle]) SUM = SUM + MASS[Particle]*RADIUS2 MASSSUM4 = MASSSUM4 + MASS[Particle] RMSRADIUS = math.sqrt(SUM / MASSSUM4) WRITE = str(TimeStep) + "\t" + str(XSUM2) + "\t" + str(YSUM2) + "\t" + str(ZSUM2) +"\t" + str(XSUM3) + "\t" + str(YSUM3) + "\t" + str(ZSUM3) + "\t" + str(ParticleNumber - DeadParticle) + "\n" open(OutputFile2,"a").write(WRITE) print "RMS Radius is", RMSRADIUS print "There are", (ParticleNumber - DeadParticle), "particles." if TimeCounter == 10: Rval = range(0,101) Nval = range(0,101) for i in range(0,100): Rval[i] = 0. Nval[i] = 0. for Particle in range(0,ParticleNumber): RADIUS3 = int(math.sqrt((XPOS[Particle]*XPOS[Particle]+YPOS[Particle]*YPOS[Particle]+ZPOS[Particle]*ZPOS[Particle]))/10) if RADIUS3 > 100: RADIUS3 = 100 Rval[RADIUS3] = Rval[RADIUS3] + MASS[Particle] Nval[RADIUS3] = Rval[RADIUS3] + 1 OutputFile3 = "/home/tom/Desktop/PHYSICS/4Time" + str(TimeStep) open(OutputFile3,"w").write("Radius\tMass\n") for i in range(0,100): WRITE = str(i) + "\t" + str(Rval[i]) + "\t" + str(Nval[i]) + "\n" open(OutputFile3,"a").write(WRITE) TimeCounter = 0 #Plots the points #g = Gnuplot.Gnuplot(persist = 1) #g.clear() #d = Gnuplot.Data(XPOS,YPOS,ZPOS,title = str(TimeStep)) #g.splot(d) if (ParticleNumber - DeadParticle) < 50: for Particle in range(0,ParticleMax): if MASS[Particle] > 0: # Find the Average Angular Momentum XSUM3 = 0. YSUM3 = 0. ZSUM3 = 0. XDIST = XPOS[Particle] YDIST = YPOS[Particle] ZDIST = ZPOS[Particle] XSUM3 = (XSUM3 + (YDIST*ZVELO[Particle] - ZDIST*YVELO[Particle])*MASS[Particle]) YSUM3 = (YSUM3 + (ZDIST*XVELO[Particle] - XDIST*ZVELO[Particle])*MASS[Particle]) ZSUM3 = (ZSUM3 + (XDIST*YVELO[Particle] - YDIST*XVELO[Particle])*MASS[Particle]) LVector = math.sqrt(XSUM3*XSUM3 + YSUM3*YSUM3 + ZSUM3*ZSUM3) Lx = XSUM3/LVector Ly = YSUM3/LVector Lz = ZSUM3/LVector STRING = str(TimeStep) +"\t"+ str(Particle) +"\t"+ str(Lx) +"\t"+ str(Ly) +"\t"+ str(Lz) +"\t"+ str(XDIST) +"\t"+ str(YDIST) +"\t"+ str(ZDIST) +"\t"+ str(MASS[Particle]) + "\n" open(BigOutput,"a").write(STRING) # Creates an output file at directory specified above. open(OutputFile,"w").write("Particle\tRadius\tVelocity\tMass\n") for Particle in range(0,ParticleNumber): if ParticleID[Particle] == 0: # Calculates the velocity and radius for the output. ParticleX = ParticleID[Particle] Velocity = math.sqrt(XVELO[Particle]*XVELO[Particle]+YVELO[Particle]*YVELO[Particle]+ZVELO[Particle]*ZVELO[Particle]) Radius = math.sqrt(XPOS[Particle]*XPOS[Particle]+YPOS[Particle]*YPOS[Particle]+ZPOS[Particle]*ZPOS[Particle]) Mass = MASS[Particle] WRITE = str(ParticleX) + "\t" + str(Radius) + "\t" + str(Velocity) + "\t" + str(Mass) + "\n" open(OutputFile,"a").write(WRITE)