In the lectures you should have learned:
Unix (Linux or Mac) terminals and how to work with directories and files
Python for scientific programming
Using open source software, using PDE athena as example
Things we didn't cover in Python:
In [ ]:
%matplotlib inline
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
from astropy.io import fits
from scipy.optimize import curve_fit
Create three variables a,b,c all representing the number 1 in the three basic python types, and print out their type using the type(a) function call:
a = ...
b = ...
c = ...
print...note the word "class" you should see, which we never really covered. But just remember: everything in python is an object (ie. class)
In [ ]:
#Q1
a = 1
b = 1.0
c = "1"
print(type(a),type(b),type(c))
In [ ]:
!echo "# x y yerr" > 123.tab
!echo 0.0 0.0 0.10 >> 123.tab
!echo 1.0 1.0 0.20 >> 123.tab
!echo 2.0 0.5 0.15 >> 123.tab
!echo 3.0 1.5 0.20 >> 123.tab
!echo 4.0 3.5 0.15 >> 123.tab
!pwd
!cat 123.tab
Use the "loadtxt" function from the numpy module to read this data, and create three arrays representing these three columns. Print out the shape of the data array and notice where the rows and colums are.
In [ ]:
#Q2
data = np.loadtxt('123.tab')
print(data.shape)
x = data[:,0]
y = data[:,1]
z = data[:,2]
One the many ipython magic %-commands we did not cover is the "%whos" command, which is like a who is who in memory. Try this out in the next cell. You can find more about the magic commands in e.g. https://ipython.org/ipython-doc/3/interactive/magics.html
In [ ]:
%whos
In [ ]:
#Q3
plt.errorbar(x,y,fmt='o',yerr=z)
plt.xlim(-0.5,4.5)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("My 123 Title")
In [ ]:
#Q4
def line(x, a, b):
return a*x + b
a,b = curve_fit(line,x,y)[0]
print(a,b)
yfit = a*x + b
plt.errorbar(x,y,fmt='o',yerr=z)
plt.plot(x,yfit)
plt.xlim(-0.5,4.5)
In [ ]:
#Q5
data = fits.open('../data/ngc6503.cube.fits')[0].data.squeeze()
for z in [1,40,50,87]:
print("z=%d rms=%g" % (z,data[z,:,:].flatten().std()))
Make a copy of the orbits-01 notebook to orbits-01b and change the potential to work in an ellipsoidal coordinate system emulating something like a (non-rotating) elliptical galaxy.
1) You would replace $$ r^2 = x^2 + y^2 $$ with $$ r^2 = x^2 + {y^2 \over b^2} $$ where $b$ is the effective axis ratio in the potential of the mass distribution. Then note that $$ {{\partial \Phi} \over { \partial x}} = {{\partial \Phi} \over { \partial r}}. {x \over r}, {{\partial \Phi} \over { \partial y}} = {{\partial \Phi} \over { \partial r}}. {y \over r b^2}, $$
2) Instead of launching the orbit from the X axis, now launch it from the Y axis. First confirm with b=1.0 that you get the same orbit as you got in the orbits-01 notebook, but use b=0.8 for the experiments. Take the same launch velocity of 0.3, and launch at different positions along the Y axis, stepping down from 1.0 to 0.9, 0.8 etc. Notice a transition? At what launch position does the transition occur?
3) Integrate the orbit at least 10 times longer, and use plt.plot() instead of plt.scatter(), since there would be too many points on the screen.
4) You can either use our handcrafted stepper, or feel free to use the scipy.integrate.odeint function, since that conserves the energy better. But be sure to check that energy is conserved. What about angular momentum?
In [ ]:
#Q6 (in a copy notebook)