# Introduction to Computational Physics Lab¶

*name*: ____________

This is a *Jupyter Notebook*. It is organized in "cells", some containing
Python code (those labelled with `[]:`

) and others containing text (such as this one).

When you start this file with Jupyter, you can edit and execute the code
blocks. Click in a code block, type or edit some code, then type `Shift-Enter`

to execute the code. You can also edit text blocks like this one; type
`Shift-Enter`

to render the text after editing. If the address of this page in
your browser ends in `.html`

then you are viewing a static rendering of the
notebook and you will not be able to edit it.

To complete this lab, you must

**Get Jupyter running**. See How to Get Started with the Jupyter Notebook.**Download the code files**and save them in a folder where you plan to do your computational work this semester. The files may be found on the lab web page.**Type your name**in the text block at the top of this document. (Double-click the block to open it for editing.)**Read through**this entire notebook to see how it's organized.**Execute the import lines**without error.**Answer every "Question" in B. Exercises**by writing and executing the appropriate code. That section starts with "Example 1" and "Question 1". You can answer the "Question 1" by following the model in the "Example 1". Repeat for all Questions.*Section C. Importing Packages*is optional.**Study the methods illustrated in D. Vector Sum Example**.**Answer all the E. Vector Problems**as directed by writing and executing the code required.**Save this notebook as a .ipynb file and submit it**on the course web page when the is lab due.

## A. Set Up¶

To be able to execute much of the code we use, you will need to first execute
these two `import`

statements.

```
from numpy import *
import scipy.constants
```

The NumPy package (short for *numerical python*)
will allow you to use vectors, trig functions and much more.

The `scipy.constants`

package contains a lot of constants and unit conversion values.
The SciPy package (short for *scientific python*)
contains many advanced computational physics libraries.

## B. Exercises¶

Below are a series of examples, each followed by a similar question. Copy the example code and paste it into the question's answer cell. Then edit as needed to answer the question.

### Example 1¶

Mt. Everest is 29 000 ft high. How many meters is that?

```
h_ft = 29000
h_m = h_ft * scipy.constants.foot
print(f'height is {h_m} meters')
```

height is 8839.199999999999 meters

Note: The

`scipy.constants`

conversion factors all convert to SI units.

### Question 1¶

The Capitol is 22 800 ft from Truax. How many meters is that?

```
```

### Example 2¶

The SI unit of atmospheric pressure is the Pascal (Pa). One Pascal is one Newton of force per square meter of area:

$$1\;\rm{Pa} = 1 \frac{\rm{N}}{\rm{m}^2}$$

Standard Atmospheric Pressure is $1.013 \times 10^5\;\rm{Pa}$. Express this in pounds per square inch.

Notes:

- to convert
from(notto) SI units, you have todivideby the conversion factor.- to square a quantity, python uses a double-asterisk, as in
`x**2`

```
p = 101300 # p is atmospheric pressure in N per m²
p = p / scipy.constants.lbf # p is now in lb per m²
p = p * scipy.constants.foot**2 # lb per ft²
p = p / 12**2 # lb per in²
print(f'Standard Atmospheric Pressure is {p} lb/in²')
```

Standard Atmospheric Pressure is 14.692322832070191 lb/in²

### Question 2¶

A tire needs to be inflated to a pressure of $2.07 \times 10^5$ Pa. What is this pressure in lb/in²?

Hint: enter numbers in scientific notation with an

`e`

: $2.07 \times 10^5$ =`2.07e5`

.

```
```

### Example 3¶

You measure the diameter of a sphere to be $4.02 \pm 0.02$ cm.

Find the sphere's volume in cm$^3$ along with its absolute uncertainty.

```
d = 4.02 # diameter, cm
unc_d = 0.02/d # relative uncertainty in d
r = 0.5*d # find radius
unc_r = unc_d
V = 4/3*pi * r**3 # use ** to raise to a power
unc_V = 3*unc_r
unc_V = unc_V * V # convert to absolute uncertainty
print(f'volume is {V} ± {unc_V} cm^3')
```

volume is 34.01549392577856 ± 0.5076939391907248 cm^3

You can control the number of digits after the decimal that print out
with *string formatting* as follows:

```
print(f'volume is {V:.2f} ± {unc_V:.2f} cm^3')
```

volume is 34.02 ± 0.51 cm^3

The `:.2f`

directs python to show 2 decimal places of a "floating point" (ie,
decimal) value.
If you would like to use scientific notation, use an `e`

instead of `f`

as follows:

```
print(f'volume is {V:.2e} ± {unc_V:.2e} cm^3')
```

volume is 3.40e+01 ± 5.08e-01 cm^3

You can learn more about fancier printing in Python by a web search.

### Question 3¶

Astronomers report the diameter of a planet to be $(6.37 \pm 0.02) \times 10^6$ m.

Find the planet's volume in m$^3$ along with its absolute uncertainty.

```
```

### Example 4¶

The coordinates of point $P$ are measured to be $x = 2.55 \pm 0.02$ m, $y = 1.24 \pm 0.04$ m. What is the angle $\theta$ (in degrees) and its uncertainty?

Finding $\theta$ requires the inverse tangent function
(also known as the *arctangent*) as follows.

```
x = 2.55 # m
dx = 0.02 # m, absolute uncertainty
y = 1.24 # m
dy = 0.04 # m, absolute uncertainty
Q = arctan(y/x)
print(f'theta = {Q*180/pi} degrees')
```

theta = 25.93247363696144 degrees

Computers always work in radians,
so the last line converts the radians to degrees (`*180/pi`

).

We have no rule about propagating the uncertainty through an `arctan`

function.
But we can find the range of values we would get when $x$ and $y$ are at the
edge of their uncertainty.

```
print(f'the angle is theta = {Q*180/pi} degrees')
print(f'other theta values within our uncertainty are...')
print( arctan((y+dy)/(x+dx)) *180/pi )
print( arctan((y+dy)/(x-dx)) *180/pi )
print( arctan((y-dy)/(x+dx)) *180/pi )
print( arctan((y-dy)/(x-dx)) *180/pi )
```

the angle is theta = 25.93247363696144 degrees other theta values within our uncertainty are... 26.47580548677303 26.836164829561802 25.02909705182668 25.375373720723683

So a good measure of the uncertainty in $\theta$ is 1/2 of the difference $\max (\theta) - \min (\theta)$. Looking at the four values printed above, you need to find the maximum and minimum.

```
Q1 = arctan((y+dy)/(x-dx))*180/pi
Q2 = arctan((y-dy)/(x+dx))*180/pi
dQ = 0.5*( Q1 - Q2 )
print(f'theta = {Q*180/pi:.2f} ± {dQ:.2f} degrees')
```

theta = 25.93 ± 0.90 degrees

We've answered the question, but in case you're curious,
a better way to find the `min`

and `max`

would be to put the results in a
Python list
(comma separated between square brackets) and then use the `min`

and `max`

functions:

```
Q_list = [ arctan((y+dy)/(x+dx)) *180/pi,
arctan((y+dy)/(x-dx)) *180/pi,
arctan((y-dy)/(x+dx)) *180/pi,
arctan((y-dy)/(x-dx)) *180/pi ]
dQ = 0.5*( max(Q_list) - min(Q_list) )
Q_ave = 0.5*( max(Q_list) + min(Q_list) )
print(f'theta = {Q*180/pi:.2f} ± {dQ:.2f} degrees')
```

theta = 25.93 ± 0.90 degrees

### Question 4¶

A distance along a slope is measured to be $d=2.00 \pm 0.02$ m. The slope lies at an angle of $\theta = 30.0 \pm 0.1$ degrees above the horizontal. The vertical height of this slope is given by

$$y = d \sin \theta$$

Find the value and uncertainty of $y$ (ie, find $y \pm dy$).

```
```

## C. (Optional) Importing packages¶

If you have a computational problem to solve, there's a good chance someone has already solved it. Python's FLOSS (free-libre open source software) philosophy lends itself to community-shared solutions.

For example, you might wonder if anyone has written code to automatically
handle uncertainty calculations. Answering this kind of question is best done
with an internet search for, in this case,
"python uncertainty".
You would quickly learn that the `uncertainties`

package is popular.

However, you cannot import it since it does not come with the standard
Python/Jupyter installation. You will receive an `import error`

if you execute
the `import`

line below before installing the package.
In this case (and for many Python packages),
installation is easy. Execute the following at a command line:

`pip install uncertainties`

Note: You can obtain a command line in Windows by running app

`cmd`

, in a Mac by launching a`terminal`

. This is not available in`jupyter.org/try`

.

Now you may import and use the package. Check its documentation for how to use it.

```
import uncertainties
```

Now uncertainty calculations can be done automatically (and more accurately)!

Consider Example 3, above, using this package.

```
d = uncertainties.ufloat(3.14, 0.02) # from Example 3, above
r = 0.5*d
V = 4/3*pi * r**3
print(f'volume is {V} cm^3')
```

volume is 16.21+/-0.31 cm^3

## D. Vectors in Python¶

Suppose we know the vector sum $\vec{F_1} +\vec{F_2} + \vec{F_3} = \vec{0}$, and we are given the following information.

- vector $\vec{F_1}$ lies along the positve $x$ axis, with magnitude 200 N
- vector $\vec{F_2}$ makes an angle of $107^{\circ}$ with the $+x$ axis, and has magnitude 300 N

Answer the following questions:

- What is the direction and magnitude of vector $\vec{F_{4}} = \vec{F_1} + \vec{F_2}$ ?
- What must be the direction and magnitude of vector $\vec{F_3}$ ?

**Solution**

First we need to define vectors $\vec{F_1}$ and $\vec{F_2}$.
In Python, this is done with the `array`

statement.
Here are two equivalent ways we could define $\vec{F_1}$ :

```
F1 = array(( 200, 0 )) # specify the components of F1
F1 = 200 * array(( 1, 0 )) # write F1 as a scalar times a unit vector
```

For $\vec{F_2}$ it's easier to use the unit vector form. All angles must be in radians.

```
Q = 107*pi/180 # angle in radians
F2 = 300 * array(( cos(Q), sin(Q) ))
```

If you add the variables `F1`

and `F2`

,
it will (correctly) add the components of the vectors.
Let's add them and call the sum $\vec{F_{4}}$:

```
F4 = F1 + F2
print('F4 = ',F4)
```

F4 = [112.28848858 286.89142679]

This means that $\vec{F_{4}} = 112 \hat{\imath} + 287 \hat{\jmath}$ N.

Let's find the magnitude of this vector: $F = \sqrt{ F_x^2 + F_y^2}$.

- to square a Python array use
`F4**2`

; this will square each component individually `sum`

will add up all the (squared) array components`sqrt`

is the square root

```
mag = sqrt(sum(F4**2))
print(f"magnitude of F4 is {mag} N")
```

magnitude of F4 is 308.08342284724057 N

To find the direction of $\vec{F}$, we use $\theta = \arctan \left (\frac{F_y}{F_x} \right )$.

The components of a python array are accessed by indexing them, starting from zero, like this:

$F_x =$ `F4[0]`

$F_y =$ `F4[1]`

```
angle = arctan( F4[1] / F4[0] )
print(f'angle is {angle*180/pi} degrees') # convert rad to degrees
```

angle is 68.62476644579117 degrees

The vector $\vec{F_3}$ should be in the opposite direction to $\vec{F_{4}}$. So it should have the same magnitude with an angle $\theta+180^{\circ}$ or $\theta+\pi$ radians.

```
F3 = mag * array(( cos(angle+pi), sin(angle+pi) ))
print("F3 =",F3)
```

F3 = [-112.28848858 -286.89142679]

We know that $\vec{F}_3 = -\left (\vec{F_1} +\vec{F_2} \right ) = -\vec{F_{4}}$. Look above and note that the components of $\vec{F_3}$ are the negative of the components of $\vec{F_{4}}$.

That's it, we're done. The following code block summarizes the operations required to answer the question.

```
F1 = 200 * array(( 1, 0 ))
F2 = 300 * array(( cos(107*pi/180), sin(107*pi/180) ))
F4 = F1 + F2
mag = sqrt(sum(F4**2))
angle = arctan(F4[1]/F4[0])
print(f'magnitude of F4 is {mag:0.2f} N and its direction is {angle*180/pi:0.2f} degrees')
```

magnitude of F4 is 308.08 N and its direction is 68.62 degrees

### A Note about `arctan`

¶

The `arctan`

function above has a problem:
it can give only give angles between $-\pi/2$ and $+\pi/2$,
in the first and fourth quadrants.
This is because the computer cannot distinguish between, for example,
`arctan(-1/1)`

, which should be $- \pi/4$ and
`arctan(1/-1)`

, which should be $3 \pi/4$.

The solution is to use the `arctan2(y,x)`

function,
giving it the $y$ and $x$ values individually
(not their ratio) so it knows which angle to return.

```
print('arctan(1/-1)=',arctan(1/-1)*180/pi)
print(' wrong quatrant!')
print('arctan2(1,-1)=',arctan2(1,-1)*180/pi)
print(' correct!')
```

arctan(1/-1)= -45.0 wrong quatrant! arctan2(1,-1)= 135.0 correct!

Note the

`arctan2`

function requires`y,x`

, not`y/x`

.

## E. Vector Problems¶

Now it's your turn.

Find the magnitude and direction of the following vectors. Give your answers as follows:

- all values are decimals expressed with three significant figures
- directions are angles in degrees off the $+x$ axis (postive angles for counterclockwise and negative angles for clockwise)

(a) $\vec{A} = - \frac{1}{\sqrt{2}} \hat{\imath} - \frac{1}{\sqrt{2}} \hat{\jmath}$

(b) $\vec{B} = - 10 \hat{\imath} + 10 \hat{\jmath}$

(c) $\vec{C} = 1000 \hat{\imath} - \hat{\jmath}$

(d) $10\vec{A} - \frac{\vec{B}}{2}$

```
```

Hint: for question (d), you can write

`10*A - B/2`

(assuming you defined`A`

and`B`

correctly).