Repaso de programación científica en Python¶

Variables¶

Tipos fundamentales

In [1]:
i = 4 # integer
x = -0.4 # float
z = 3 + 4j # complex
t = "text" # string
b = True # boolean

Operaciones fundamentales con los tipos fundamentales

In [2]:
x = -0.5
y = 2

print(f"x+y = {x+y}, x-y={x-y}, x*y = {x*y}, x/y = {x/y}, x//y = {x//y}, y**x = {y**x}")

x+y = 1.5, x-y=-2.5, x*y = -1.0, x/y = -0.25, x//y = -1.0, y**x = 0.7071067811865476

Organización de la información¶

Listas

In [3]:
l = [1,2,'a',-1.5] # list
In [4]:
# access an element of the list
l[2]
Out[4]:
'a'
In [5]:
# access the last element of the list
l[-1]
Out[5]:
-1.5
In [6]:
# access some elements of the list
l[0:2]
Out[6]:
[1, 2]
In [7]:
# add an element to a list
l.append('new')
In [8]:
l
Out[8]:
[1, 2, 'a', -1.5, 'new']
In [9]:
# check whether an element is in a list or not
print(3 in l)
print(2 in l)

False
True

Diccionarios

In [10]:
d = {'physics': 'science that studies the principles of the universe',
    'astronomy': 'science that studies the '}
In [11]:
# access an element
d['physics']
Out[11]:
'science that studies the principles of the universe'
In [12]:
# define a new element
d[130] = 4
d
Out[12]:
{'physics': 'science that studies the principles of the universe',
 'astronomy': 'science that studies the ',
 130: 4}
In [13]:
# get only the keys
d.keys()
Out[13]:
dict_keys(['physics', 'astronomy', 130])

Clases y objetos¶

In [14]:
class Particle:
    mass = 0
    position = 0
In [15]:
ball = Particle()
ball.mass = 4
ball.position = -0.4

Condiciones¶

In [16]:
energy = -4.0

if energy < 0:
    print("energy is leaving the system")
elif energy == 0:
    print("energy is conserved")
else:
    print("energy is entering the system")

energy is leaving the system

Bucles¶

In [17]:
l
Out[17]:
[1, 2, 'a', -1.5, 'new']
In [18]:
for element in l:
    print("the element",element,"is a",type(element))

the element 1 is a <class 'int'>
the element 2 is a <class 'int'>
the element a is a <class 'str'>
the element -1.5 is a <class 'float'>
the element new is a <class 'str'>
In [19]:
for i,element in enumerate(l):
    print(i,element)

0 1
1 2
2 a
3 -1.5
4 new
In [20]:
# Imagine you want to plot several lines, each one has different colors, and different start and end points
line_colors = ['blue','red','orange','black']
intervals = [[1,2],[-3,4.4],[-1.2,0.4],[5,6]]

for line_color, interval in zip(line_colors,intervals):
    print(f"this line extends from {interval[0]} to {interval[1]} and has color {line_color}")

this line extends from 1 to 2 and has color blue
this line extends from -3 to 4.4 and has color red
this line extends from -1.2 to 0.4 and has color orange
this line extends from 5 to 6 and has color black

One-liners¶

In [21]:
# List comprehensions
new_list = [i for i in range(45,50)]
new_list
Out[21]:
[45, 46, 47, 48, 49]
In [22]:
positions = [3,-4,-1,2]

# get only the positive positions (filter a list)
positions = [pos for pos in positions if pos > 0]
positions
Out[22]:
[3, 2]
In [23]:
frequency = 50e5 # Hz

is_a_radio_wave = True if frequency > 1e4 or frequency < 1e7 else False
print(is_a_radio_wave)

True

Numpy¶

In [24]:
import numpy as np
import matplotlib.pyplot as plt
In [25]:
# array of integers
a = np.array([1,2,3,4,6],dtype=int)
a
Out[25]:
array([1, 2, 3, 4, 6])
In [26]:
# array of floats
a = np.array([3,4,3,5],dtype=float)
a
Out[26]:
array([3., 4., 3., 5.])
In [27]:
# bidimensional array
x_points = np.linspace(-10,10,1000)
y_points = np.linspace(-10,10,1000)

x, y = np.meshgrid(x_points, y_points, indexing='ij')
x_points.shape, y_points.shape, x.shape, y.shape
Out[27]:
((1000,), (1000,), (1000, 1000), (1000, 1000))
In [28]:
# calculation of the distance to the origin
r = np.sqrt(x**2 + y**2)

# plot of the distance as a pcolormesh
plt.pcolormesh(x, y, r, cmap='gray_r')
plt.colorbar()
Out[28]:
<matplotlib.colorbar.Colorbar at 0x111002510>
No description has been provided for this image
In [29]:
# delete all distances larger than a threshold
threshold = 8

# first we make a copy of the original array
modif_r = np.array(r)
modif_r[ r > 8] = 0.


# plot of the distance as a pcolormesh
plt.pcolormesh(x, y, modif_r, cmap='gray_r')
plt.colorbar()
Out[29]:
<matplotlib.colorbar.Colorbar at 0x111188550>
No description has been provided for this image
In [30]:
# calculation of a force
force = 0.1/r**2

# plot of the distance as a pcolormesh
plt.pcolormesh(x, y, force, vmin=0,vmax=10,cmap='gray_r')
plt.colorbar()
Out[30]:
<matplotlib.colorbar.Colorbar at 0x11126cb90>
No description has been provided for this image

Example with simulation data¶

This simulation data has been produced by PLUTO and then interpolated into a Cartesian grid. All quantities are expressed in Gaussian-CGS units. The Cartesian grid extends from x = 0 to x = 200 au, and from z = 0 to z = 200 au and has 1000 cells in each direction.

At the instant of this simulation data, the star has a mass of 2.90 solar masses.

In [31]:
x_points = np.linspace(0,200,1000)
z_points = np.linspace(0,200,1000)

x, z = np.meshgrid(x_points, z_points, indexing='ij')
In [32]:
rho =  np.load('sim1data/rho.npy')
v_r =  np.load('sim1data/vx1.npy')
v_th = np.load('sim1data/vx2.npy')
v_ph = np.load('sim1data/vx3.npy')
In [33]:
# density
plt.pcolormesh(x,z, np.log10(rho))
plt.colorbar()
Out[33]:
<matplotlib.colorbar.Colorbar at 0x111327750>
No description has been provided for this image
In [34]:
# radial velocity in km/s
plt.pcolormesh(x,z, v_r/1e5, vmin=-10,vmax=10, cmap='seismic_r')
plt.colorbar()
Out[34]:
<matplotlib.colorbar.Colorbar at 0x11140fc50>
No description has been provided for this image
In [35]:
# Separating the outflow from the inflow
outflow = np.array(v_r)/1e5 # km/s
outflow[outflow < 0] = 0
plt.pcolormesh(x,z, outflow, cmap='plasma')
plt.colorbar()
Out[35]:
<matplotlib.colorbar.Colorbar at 0x1114ffb10>
No description has been provided for this image

Let's compute the magnitude of the gravitational force that the star exerts on the circumstellar material.

In [36]:
au_in_cm = 14950000000000.0
r = np.sqrt(x**2 + z**2)*au_in_cm
G = 6.67e-08 # gravitational constant G in cgs units
Msun = 2e33 # solar mass in g

Mstar = 2.9*Msun # (given from the simulation)

F_grav = rho*G*Mstar/r**2 # in dynes/cm^3

plt.pcolormesh(x,z,F_grav, cmap='plasma')
plt.colorbar()

/var/folders/1r/kzh10jnx48jcn2l34mt99ctw0000gn/T/ipykernel_2776/1293598483.py:8: RuntimeWarning: divide by zero encountered in divide
  F_grav = rho*G*Mstar/r**2 # in dynes/cm^3
Out[36]:
<matplotlib.colorbar.Colorbar at 0x1115ea0d0>
No description has been provided for this image

The result is somewhat meaningless because we don't have anyting to compare it against. Let's compute the centrifugal force per volume (in the co-rotating frame), that is, $\rho v_\phi^2/r \, \mathbf e_x$, and then let's compare it against the horizontal component of gravity, that is, $\rho GM/r^2 \sin\theta \, \mathbf e_x$.

In [37]:
th = np.atan2(x,z)

F_centrif_x = rho*v_ph**2/r
F_grav_x = F_grav*np.sin(th)

/var/folders/1r/kzh10jnx48jcn2l34mt99ctw0000gn/T/ipykernel_2776/3700062375.py:3: RuntimeWarning: divide by zero encountered in divide
  F_centrif_x = rho*v_ph**2/r
/var/folders/1r/kzh10jnx48jcn2l34mt99ctw0000gn/T/ipykernel_2776/3700062375.py:4: RuntimeWarning: invalid value encountered in multiply
  F_grav_x = F_grav*np.sin(th)

Now, let's divide the magnitude of both components and plot the result

In [38]:
plt.pcolormesh(x,z, np.abs(F_centrif_x)/np.abs(F_grav_x), vmin=0, vmax=2, cmap='seismic_r')
plt.colorbar()

/var/folders/1r/kzh10jnx48jcn2l34mt99ctw0000gn/T/ipykernel_2776/46531648.py:1: RuntimeWarning: divide by zero encountered in divide
  plt.pcolormesh(x,z, np.abs(F_centrif_x)/np.abs(F_grav_x), vmin=0, vmax=2, cmap='seismic_r')
Out[38]:
<matplotlib.colorbar.Colorbar at 0x1116ccf50>
No description has been provided for this image

The red regions (the ratio < 1) mean that the gravity wins. The blue regions (ratio > 1) mean that the centrifugal force wins. The centrifugal force, together with the magnetic field, is the cause the jet is being launched. We'll see this later in the course.