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

[86]: # Example of calculations with variables a = 5 b = 7 # Make a new variable that multiplies the first two together: c = a * b # This is another way of including results in a print statement. Within the ␣ ↪ curly brackets {}, # the number gives the index of the list in 'format' to print in that location. print ( ' {0} x {1} = {2} ' . format(a,b,c)) 5 x 7 = 35 To practice doing calculations, set up your constants and calculate your value for 𝜌 𝑐 in the block below. You can also use the astropy package, which includes libraries of constants and units. If you’re careful to include the correct units for all your variables and constants, you can also use astropy to convert between units in your final answer. See more documentation here: astropy units 4

astropy_constants So our equation for 𝜌(𝑟) = 𝜌 𝑐 ∗ (1 − 𝑟/𝑅) . Since our radius array is normalized with 𝑥 = 𝑟/𝑅 , we can rewrite this as 𝜌(𝑥) = 𝜌 𝑐 ∗ (1 − 𝑥) . Run the block below to see how we can create an array for 𝜌(𝑟) in python. [232]: # Next, make arrays for P and T as a function of r/R # Follow the method for density, and use the normalized radius for your r/R ␣ ↪ terms # Watch units! mypress = 3.14 * G * rho_c **2 * rsun **2 * 5/36 * ((rrsun **2 ) * -0.67 ) + ( 0.78 * ␣ ↪ (rrsun **3 )) - ( 0.25 * (rrsun **4 )) mytemp = ( 15 * mmsun * G * mu * mp) / (rrsun * 36 * k) # It might help to print the first few values for your pressure and temperature ␣ ↪ arrays to check values & units: print ( 'First ten values of pressure: ' ) print (mypress[ 0 : 9 ]) print ( 'First ten values of temperature: ' ) print (mytemp[ 0 : 9 ]) First ten values of pressure: [-6.26481116e+11 -7.15011664e+11 -8.00569442e+11 -8.81704317e+11 -9.66754549e+11 -1.04564172e+12 -1.11720573e+12 -1.19113828e+12 -1.25639421e+12] First ten values of temperature: [5.11445578e-19 5.98421062e-19 6.78648940e-19 7.54449171e-19 9.26355804e-19 9.89696371e-19 1.05322107e-18 1.11273970e-18 1.17374517e-18] 1.5 Part 3: Plot your results and compare with the full stellar model In the next three cells, create a plot for each of your density, pressure, and temperature models as a function of the normalized radius. Plot your derived values with black lines. Be sure to label your axes. In the same figure, plot the full solar model results with red lines. Make both x and y axes in log scale. Add axis labels for each that show your units. The above plot for density is filled in the first cell to get you started. 1.6 Include these three plots in your assigment submission as Q2 part d) [190]: # Plot density plt . figure(figsize = ( 5 , 5 )) plt . plot(rnorm, my_rho_r, color = 'blue' ) # Always add axis labels! Include units when needed. 5

[245]: # Plot temperature plt . figure(figsize = ( 5 , 5 )) plt . plot(rrsun, mytemp) plt . xlabel( 'normalized radius (r/R)' ) plt . ylabel( 'tempurature (k)' ) plt . yscale( 'log' ) plt . xscale( 'log' ) 7

While our results aren’t perfect, you will hopefully find that we are generally within an order of magnitude of the solar model, which is impressive considering our very simple density assumption, and without including energy generation, opacity, etc.! [ ]: [102]: # Read in the full stellar model data by modifying the code below # Import the needed packages import numpy as np import matplotlib.pyplot as mpl # Change these array names to something more useful # Modify the filename # Modify the skiprows and usecols to read the correct column to the ␣ ↪ corresponding array mmsun, rrsun, t, rho, p = np . loadtxt( '222' , usecols = ( 0 , 1 , 2 , 3 , 4 ), ␣ ↪ skiprows =25 ,unpack = True ) 8