AST221_2023F_Assignment2_Q2

AST221_2023F_Assignment2_Q2 October 22, 2023 1 Assignment 2 Q2 d): A simple stellar model 1.1 Before you begin Assuming you have loaded this file into your Jupyter Notebooks workspace, make sure to press the “play” button at the top of the page in each box. This will render the Markdown into nicely- formatted text and execute all of the sections of Python code. You can also hit “shift”+“enter” for the same result. You will need to upload the text file “bs05_agsop.txt’ to your Jupyter hub folder or wherever you are running your jupyter notebook. 1.2 Introduction This Jupyter Notebook is Q2d) of Assignment 2. You will be working through the content it contains and then filling in your own work in the spaces provided. In this notebook, you will use the theoretical equations you derived for density, pressure, and temperature as a function of radius given our simple density model. You will read in values for 𝑟 , 𝜌 , 𝑃 , and 𝑇 as a function of 𝑟/𝑅 that were derived from a full stellar model, and see how well we are able to reproduce the full stellar model with a simple approximation for the density. 1.2.1 For Assignment 2, include your final three plots in your solution for Q2 d)! 1.3 PART 1: Read in the full stellar model data Follow your previous tutorial and assignment solutions to read in the needed columns from the stellar model data file. You will need to set the following parameters in np.loadtxt: usecols, skiprows. Be sure to set unpack = True. Most error messages stem from incorrect mapping of parameters to columns, or an incorrect number of rows at the top of the file to skip. (Other issues include unclosed brackets and the like, so check your lines carefully!) Open the file in your favourite text editor (or on the Jupyter hub) to figure out which columns you will need and understand their units. The units for the pressure in this file are dyne/cm 2 . The units for temperature are K. The units for density are g/cm 3 . You can either perform your calculations below in cgs to obtain the correct units, or you can convert the values in the stellar model to SI. 1 dyne = 1 g cm/s 2 1 dyne/cm 2 = 1 N/m 2 = 0.1 Pa 1

1 g/cm 3 = 1000 kg/m 3 If you prefer to convert the stellar model data into SI units, do so in the next step. [126]: # Set up constants, e.g. M_Sun, R_Sun, etc. # numpy already has pi defined as np.pi # If you want to try astropy units and constants, you need to import the ␣ ↪ libraries: import astropy.constants as c import astropy.units as u # Here setting mu 0.61 for the Sun, where mu = the mean particle mass mu = 0.61 # Set up your own constants here. Add a comment for each that shows units. msun = 199 rsun = 696000000 G = 6.67e-11 x = 'r/R' mp = 1.7e-27 k = 1.4e-23 # One parameter you'll want to calculate, for example, is the central density ␣ ↪ (change this!): rho_c = 16e+04 # If you are converting the solar model data to different units, do so here: # (You can make new arrays, i.e., array1_new = array1 * factor, # or simply say array1 = array1 * factor. Probably one of these is better in ␣ ↪ terms of # proper coding practice, but do what works for you.) print ( ' {0} ' . format (rsun)) 696000000.0 1.4 Part 2: Calculate the stellar parameters as a function of radius from the simple stellar model From your simple stellar model, we have equations for 𝜌(𝑟) , 𝑚(𝑟) , 𝑃(𝑟) , and 𝑇(𝑟) . Note that your equations all have terms in 𝑟/𝑅 , rather than just 𝑟 . In the complete stellar model data, the radius is normalized and therefore also given in 𝑟/𝑅 . The means that we can use the normalized radius array as an input to calculate our values of the stellar parameters as a function of radius. Python is very useful for doing calculations like this, where we can use an array in 𝑟 to calculate 𝜌(𝑟) at every 𝑟 in the array. Some of the mathematical operators you will need: x to the power of y: x**y square root of a number x: np.sqrt{number} or x**(1/2) x times y: x * y 2

x divided by y: x/y Recall the order of operations - python will perform brackets, exponents, division and multiplication, addition and subtraction in that order! So use brackets to ensure your calculations are performed the way you want them to be done. When doing calculations, it is much easier to set variables for the constants you will use frequently. This allows you to see more clearly what parameters you are using in your calculations. Then when doing mathematical operations, you can use those variables in place of the numbers. As an example, run the cell below: [109]: # Use the equation for density in your simple stellar model to calculate rho(r) # Here I'm just reloading the normalized radius (r/R) column so you don't have ␣ ↪ to edit the cell. rnorm = np . loadtxt( '222' , skiprows =26 , unpack = True , usecols = ( 1 )) # By multiplying rho_c by the rnorm array, we create an array in my_rho_r my_rho_r = rho_c * ( 1 - rnorm) # Print the first ten values to show you have an array! print (my_rho_r[ 0 : 9 ]) # Plot the density as a function of radius, see that it's linear! # First, we need to import the matplotlib.pyplot library import matplotlib.pyplot as plt # Create a figure to plot in. Inside the brackets, you can set things like ␣ ↪ figsize=(6,4), # or the dpi (dots per inch). plt . figure(figsize = ( 5 , 5 )) plt . plot(rnorm, my_rho_r, color = 'black' ) # Always add axis labels! Include units when needed. plt . xlabel( 'Normalized radius (r/R)' ) plt . ylabel( 'Density (g/cm^3)' ) # Units here will depend on what you used in ␣ ↪ your rho_c calculation! [159724.8 159708.8 159694.4 159680. 159667.2 159656. 159644.8 159635.2 159624. ] [109]: Text(0, 0.5, 'Density (g/cm^3)') 3