In this unit we will learn to do some physics with this system.  This unit is about Newtonian mechanics.  We will get up to simple  problems of astronomy.  

Before we start we need to spend a little time on physical constants and units.  These can be obtained from within Mathematica. Let’s look at this quickly, but in most cases we will not use it.  

The command Quantity[magnitude, “item”] will give you some physical constant “item” with magnitude.

Quantity[2,"SpeedOfLight"]

myspeed= UnitConvert[Quantity[1/10,"SpeedOfLight"],"Meters/Sec"]

UnitConvert[Quantity[1,"ElectronMass"],"kg"]

You can extract both the magnitude and the units. This is useful when you want to get the numbers to setup or when you want to write a careful program. It gets cumbersome for day to day calculations.  And so I tend to set up the units in the beginning of a calculation into simple variables.

QuantityMagnitude[myspeed]

QuantityUnit[myspeed]

Lets calculate the escape velocity from the Earth

First we setup some constants.

myG = QuantityMagnitude[UnitConvert[Quantity["GravitationalConstant"],"Newtons*meter^2/kg^2"]]

radiusofearth = QuantityMagnitude[UnitConvert[Quantity["EarthMeanRadius"],"Meters"]]

massofearth = QuantityMagnitude[UnitConvert[Quantity["EarthMass"],"Kg"]]

masssun =  QuantityMagnitude[UnitConvert[Quantity["SunMass"],"Kg"]]

Now we setup some equations.  
1/2*  m v^2 = G* M *m / r  

v = Sqrt[ 2 G *M/r]

KineticEnergy = m * v^2/2;  
PotentialEnergy = -myG * massofearth * m/r; 
totalenergy = KineticEnergy+PotentialEnergy;

If a particle  barely escapes from any radius r from the center of the earth, then it must have total energy zero.  What is the velocity  ? Notice it is a quadratic equation and so it will have negative and positive

Solve[totalenergy == 0, v]

Let’s take the positive solution

answer  = %[[2]]

What is the answer at the Earth’s surface in units of meters/sec.  (Do you remember the answer from your courses  ?)

answer/.{r-> radiusofearth}

what about plotting it as a function of the height above the earth, up to 
10000 km.

Plot[v/.answer,{r,radiusofearth,radiusofearth+10000*10^3},PlotStyle->Bold,AxesLabel->{"radius (m)","escape velocity meters/sec"},LabelStyle->{Medium,Bold}]

moondistance = QuantityMagnitude[UnitConvert[PlanetaryMoonData["Moon","DistanceFromEarth"],"meters"]]

answer/.{r->moondistance}

Lets calculate the mass in the Galaxy that holds the solar system bound from the speed of the Sun around the Galaxy.  We need to get a bunch of numbers first.

massofsun =  QuantityMagnitude[UnitConvert[Quantity["SolarMass"],"kg"]]

suntogalaxylightyears =  QuantityMagnitude[UnitConvert[AstronomicalData["MilkyWay", "Distance"],"lightyear"]]

suntogalaxy =  QuantityMagnitude[UnitConvert[AstronomicalData["MilkyWay", "Distance"],"meters"]]

sunspeed= 220000;  (* meters/second *)

orbitperiod = 2*Pi*suntogalaxy/sunspeed   (* in seconds *)

orbitperiodyears = orbitperiod/(365*24*3600)

now use v^2/r = G*M/r^2  =>   M  = v^2*r/G

galaxymass = sunspeed^2*suntogalaxy/myG/massofsun

There is a mass of 90  billion sun’s of mass inside the sun’s radius from the center of the Galaxy !  Most of this is actually invisible in the form of Dark matter !  We are in territory that we have no instinct for !  
Galaxy halo extends out to at least 100000 lightyears.  This means there could be 90 billion X (100000/24000)^3 ~ 1.5 trillion  stars worth of mass in the galaxy !  Most of this invisible.

90*10^9*(100/24)^2

Now something about particles, dark matter, and neutrinos. 

The world is made of two types of particles: ones that stick together and form ordinary items such as people, children, planets, or food, and ones that do not stick together and float in space.  The non-stick particles also penetrate ordinary matter and occupy all available space as they diffuse from their origin.  Neutrinos are such particles.  They were found to exist because ordinary matter has the ability to  both emit and absorb such particles on rare occasions. We learn about this process as radioactivity. It is lucky for us that radioactivity is a rare process which keeps the Sun burning slowly for billions of years and allows life to evolve. It also causes the stars to burst when they reach a certain stage when all the matter is so compressed that it decays in a huge radioactive decay.  The universe is full of these neutrinos from both the big bang and the stars. They are nearly as common as light itself. 

Radioactive decay of a neutron particle

Let’s calculate what happens if we let the escape velocity become the speed of light.  We can turn the problem around and calculate the mass of the Earth if the escape velocity becomes the speed of light.  

 r =  2 G *M/c^2    
 
 Now this is actually the correct formula even when using Einstein’s theory of gravity (in 1916 from Karl Schwarzchild) 
 Even though this can be formulated from Newtonian mechanics.

lightspeed = QuantityMagnitude[UnitConvert[Quantity[1,"SpeedOfLight"],"Meters/Sec"]]; 
massofobject = lightspeed^2 *radiusofearth/myG/2

And so such an object has the mass of 2000 Suns !  Actually such an object was proposed in the 18th century by a fellow named John Michell and also by Laplace.  The idea was are there objects that are so dense that light bends backwards and cannot escape ?  It had to wait until the 1930’s that people became very serious that such objects might exists --- and they are called black holes

massofobject/masssun

What is the earth were to become a black hole; how big would it be ? This is in meters !   ~ 9 mm.

radiusofblackjholeearth  = 2*massofearth*myG/lightspeed^2

How about a projectile in a uniform gravitational field ? 
we are now going to solve Newton’s equations in a uniform field.  The units will be  meters, seconds, and kg. 
We have to learn a little bit about differential equations.

g=9.8;   (* gravity acceleration  in meters/sec^2 *)

?D

Differentiation is indicated by either   D[F[x],x]   or by F’’[t]. We use it for Newton’s equations of motion.    D[F[x],{x,n}] is the n’th derivative of F[x]. 

Let’s imagine a ball tossed.  z[t]  is the vertical motion and x[t] is the horizontal motion.  Remember that to define an equation you must use  ==.  An equation has a value of True or False that can be assigned to a variable !

equationsofmotion  =  {m * D[z[t],{t,2}] == -m g , m* D[x[t],{t,2}]==0}

vx0 = 20; vz0 = 20; 
initialconditions = {z[0]==0,  x[0]==0,z'[0]==vz0,x'[0]==vx0 }

alltogether = Join[equationsofmotion, initialconditions]

We use DSolve to solve the equations.  But we have to join all the equations together.

answer =  DSolve[alltogether, {z[t],x[t]},t]

when does the ball come back down to z=0 ? Solutions are in the form of a substitution rule.

times = Solve[z[t]==0/.answer,t]

totaltime = t/.times[[2]]

Lets make a plot of the trajectory !

ParametricPlot[{x[t],z[t]}/.answer, {t,0,totaltime},PlotStyle->{Red,Thick},Frame -> True,FrameLabel->{"X meters", "Z meters"},LabelStyle->{Medium,Bold},FrameStyle -> Bold]

Let’s imagine that there is a wall at x = xwall < 80 meters.  Then the time the ball hits the wall is 
twall = xwall/vx0.  The ball will reflect, but its upwards/downwards  motion is not disturbed.  
We can therefore replace x[t] with a If statement.  

This is getting pretty sophisticated !  But it is not difficult to do.  I even drew a line where the wall is.  
Try changing the parameter xwall to see what happens to the plot.

xwall = 30;  
twall = xwall/vx0;  
ParametricPlot[{If[t<xwall/vx0,x[t],2*xwall-x[t]],z[t]}/.answer, {t,0,totaltime},PlotStyle->{Red,Thick},Frame -> True,FrameLabel->{"X meters", "Z meters"},LabelStyle->{Medium,Bold},Epilog->{Directive[Thick],Line[{{xwall,0},{xwall,20}}]} ]

Lets make it even fancier.

Animate[ParametricPlot[{If[t<xwall/vx0,x[t],2*xwall-x[t]],z[t]}/.answer, {t,0,totaltime},PlotRange->{{-80,80},{0,25}},PlotStyle->{Red,Thick},Frame -> True,FrameLabel->{"X meters", "Z meters"},LabelStyle->{Medium,Bold},Epilog->{Directive[Thick],Line[{{xwall,0},{xwall,20}}]} ],{xwall,2,80}]

Lets review some of the many plotting functions.  

Plot  - we already used this  
LogPlot - makes a plot with vertical axes Log 
LogLinearPlot - makes a plot with horizontal axes log.  
LogLogPlot
Plot3D - plots a surface in 3D  
ParametricPlot -  plot a curve in 2D with a parameter (time)
ContourPlot - make a 3D plot with one dimension as contours
ListPlot - makes plots of lists of data points.

LogPlot[10^x,{x,0,100}]

LogLinearPlot[Sin[x^2],{x,0,3}]

LogLogPlot[ProductLog[x],{x,0,100}]  (* this is a weird function it is also called LambertW *)

?ProductLog

Plot3D[Sin[x+y],{x,-2Pi, 2Pi},{y,-2Pi,2Pi}]

ParametricPlot3D[{Sin[3 t],t Cos[5 t],t},{t,-Pi,Pi},PlotStyle->{Red,Thick},AxesStyle->Thick]

ContourPlot[Exp[-x^2]Exp[-(y-3)^2]+Exp[-(x-1)^2]+Exp[-(y)^2], {x,-5,5},{y,-5,5}]

ListPlot[{{-2,3},{4,5},{-3,1},{5,10},{6,4}},PlotMarkers->8,Joined->True]

Let’s do a  Rocket problem.  This is something that can be done in the class.  This one can be done with just algebra.  But it is not simple.  

A small rocket with mass m  has a thruster that operates for a very brief time (dt) and the rocket flies to a height h.  

How can we approximately figure out the average force (fave) in Newtons that was applied to the rocket by the thruster ? 

We can do the experiment with a bottle rocket filled with vinegar and soda.   We have to figure out both the time dt and the force fave and therefore we need two equations.  

What if we double the mass of the rocket by attaching some small weight  and measure the height  again ?

rocketmass  = 0.010 ;  (* in Kg *)  
height1 = 20;  (* meters with the first mass *)  
height2 = 4.9 ; (* meters with the mass doubled by attaching a coin for example *)  
Clear[g];

Now we do a little bit of algebra.  Remember that the rocket gets a short burst of thrust.  and so 
dt << time of flight.    Let’s leave everything in terms of variables,   F, mass, and g.  We will plug in numbers later.

v'[t] = F/mass - g;    (* this is just Newton's equation when F is on *)  
(* recall that F is on only for a short time dt threfore integrating both sides *)  
v0 = F*dt/mass - g dt;    (* we are trying to determine both F and dt. Here we use dt << t; otherwise this does not work *)  

vv  = v0 -g t;  (* velocity as a function of time after getting the initial velocity v0 *) ; 
sl = Solve[vv==0, t]//Flatten; (* Time to reach the top of the trjaectory *)  
ttop = t/.sl

Z[t_]:=  v0 t - g((t^2)/2);   (* formula for the height of the rocket as a function of time *)  
height = Z[ttop]//FullSimplify

Lets setup the equations to be solved in terms of variables. 
We will plug in the actual values at the end.  
remember that the mass is doubled for height 2.

eq1 = h1 ==height

eq2 = h2 == height/.{mass->2*mass}  (* mass is doubled *)

Now solve the two equations for two unknowns.  There should be 4 solutions !

solutions = Solve[eq1&&eq2,{dt,F}]//FullSimplify

Now we plug in numbers and see which solution has the correct form.

finalanswer  =  solutions/.{g->9.8, h1->height1,h2->height2,mass ->rocketmass}

Notice that two of the solutions have positive answers.  This is because these are quadratic equations.  Take the positive answer only.  Only one of the positive solutions satisfies the condition that dt is very small.  Therefore the correct answer is the second one.   F is in Newtons and dt is in seconds.  

Notice that it is possible to have no reasonable solution if height1,  height2 are not in the correct ratio.     Have fun thinking about this.  It is amazing that we can infer short time bursts from a simple measurement.  

Now what is the initial velocity obtained after the thrust.  This will be in meters per second. and the time to the top in seconds.

InitialVelocity = v0/.{g->9.8,mass->rocketmass}/.finalanswer[[2]]

TimeToTop =  ttop/.{g->9.8,mass->rocketmass}/.finalanswer[[2]]

Let’s now play with planet data that is in Mathematica.  Much of this  can be obtained by PlanetData[“planet”, “property”]. Get the details from Help.

?PlanetData

PlanetData["Jupiter","Mass"]

Let’s make a table of “name”, mass,  diameter,  orbit period,  distancefromsun, and eccentricity, Pluto is no longer a planet !

planets = {"Mercury", "Venus", "Earth",  "Mars", "Jupiter", "Saturn",  "Uranus", "Neptune"}

neededprop = {"Name",  "Mass",   "EquatorialDiameter", "OrbitPeriod", "DistanceFromSun", "Eccentricity" };

data =  Table[PlanetData[planets[[i]],neededprop[[j]] ] ,{i,1,8},{j,1,6}]

Prepend[data, neededprop]//TableForm

PlanetData[]

I will leave it to you to convert units into units that you like.

Graphics3D[PlanetData[PlanetData[],"OrbitPath"]]

Before we finish let’s learn how to do various operations with lists.  
Let’s first define a list.

somenumbers = {2.0, 3.4, 4.5,  1.1, 1.2, 2.3, 10.3}

morenumbers = {4.5,2.1, 90.0, 10^3,4.3,1.2,0.34}

Length of the list

listlength = Length[somenumbers]

Join two lists.

Join[somenumbers, morenumbers]

make a list of lists

nestedlist  = {{2,3},{4,5},{5,6},{2.1,1.1}}

Length[nestedlist]

Flatten out the list

flatlist = Flatten[nestedlist]

Length[flatlist]

Some arithmatic operations on lists

addedlists = somenumbers + morenumbers

substracted = somenumbers -morenumbers

squaredlist = somenumbers^2

Map a function on a list. How about exponential  ? Map[Exp,

Map[Exp, somenumbers]

Exp/@ somenumbers

definedF[x_,y_]:= x*y*Exp[x];

How do I apply this function to a list of two elements ?

xylist = {3,1}

definedF @@ xylist

This mapping and applying can get very complicated.  You have to try various ways. It is very helpful to make compact statement.

nestedlist

This will Map and apply to a nest list.

Apply[definedF,  nestedlist, 1]

definedF @@@ nestedlist

This will Thread the function over two lists: taking a number from each and joining them into a function.  x will come from the first list , and y from the second list.

MapThread[definedF,{somenumbers,morenumbers}]

list1 = {1, 2, 3} ; list2 = {4, 5, 6};  MapThread[definedF,{list1,list2}]

?MapThread