Astronomy.py

# -*- coding: utf-8 -*-

from math import pi, exp, log


AU = 1.49597870691e11 #meters (m)
# light year (ly) ........................... 9.460536207 × 1012 km = 63,240 A.U.
ly = 63240*AU
# parsec (pc) .............................. 3.08567802 × 1013 km = 206,265 A.U.
pc = 206265*AU
# sidereal year (yr) ..................... 365.2564 days
# elementary charge ................... 1.6021765 × 10−19 coulombs (C)
# atomic mass unit (mu) ............. 1.6605389 × 10−27 kg
atom_mass = 1.6605389e-27
year_coeff = 365.* 24 * 60 * 60
# Physical Constants
G = 6.673e-11 #m3 kg−1 s−2
g = 9.82 #m/s2
#speed of light in vacuum (c) ....... 2.99792458 × 105 km/s
c = 2.9979e8
kB = 1.380648e-23 #joules/kelvin (J/K)
# electron mass (me) .................. 9.109382 × 10−31 kg
# proton mass (mp) .................... 1.6726217 × 10−27 kg
mp = 1.6726217e-27
# proton-to-electron mass ratio ..... 1.836152672 × 103
# Planck constant (h) ................. 6.626069 × 10−34 joule-seconds (J⋅sec)
# Rydberg constant (R∞) ............ 1.09737315685 × 107 m−1
# Stefan-Boltzmann constant (σ) .. 5.6703 × 10−8 W m−2 K−4
sigma = 5.6703e-8
# Wien's Constant (b) ................. 0.0028977685 m K
# Solar constant ........................ 1.37 × 103 W m−2

# Astronomical Data
# Earth mass ............................. 5.97219 × 1024 kilograms (kg)
Me = 5.97219e24 #kg
# Earth radius (mean) .................. 6.371 × 103 kilometers
Re = 6.371e6 #m
Msun = 1.9891e30 #kg
# Sun radius (mean) ................... 6.95508 × 105 km
Rsun = 6.95508e8 # m
# Sun luminosity ........................ 3.83 × 1026 watts (W)
Lsun = 3.83e26
# Sun surface temperature .......... 5,777 K
Tsun = 5777.

Bsun = 1400


def hw3_1():
  # 1. We estimated the temperature at which a cloud becomes unstable to gravitational collapse to be given by kBT∼GMmR0 where M is the total mass of the cloud, R0 its size, and m the mass of a typical particle. For a cloud of molecular hydrogen (formed from two hydrogen atoms), m= 3.34e-27 kg. If the cloud that formed our solar nebula had a mass of M=1.5M⊙ and a radius of R0= 1e4 AU, calculate its temperature in kelvin from the equation above. Round to two significant figures.
	# kBT∼GMmR0
	# M is the total mass of the cloud
	# R0 its size
	# m the mass of a typical particle

	m= 3.34e-27 #kg.
	M=1.5 * Msun #and a radius of 
	R0= 1e4 * AU #AU, calculate its temperature in kelvin from the equation above. Round to two significant figures.

	T = G*M*m/R0/kB
	print T, " HW3.1"


	#2

	#G*M^2*(1/R-1/R_0)
	# 3. For M=M⊙, R0= 1e4 AU and R= 6.95508e8 meters (the Sun’s radius) find the amount of gravitational energy released in forming the Sun, in joules. Round to one significant figure.
	M = Msun
	R_0 = 1e4 * AU
	R= 6.95508e8 
	Energy = G*M**2*(1./R-1./R_0)

	print Energy, " HW3.3"

	#4 The Sun’s luminosity is L⊙= 3.83e26 watts (1 watt = 1 joule/second). How long does it take for the Sun to lose the amount of energy generated by collapse, by radiating it away at this rate? Give your answer in years and round to one significant figure. For your convenience, a year is approximately 3.16e7 seconds.
	flow = 3.83e26 # watts (1 watt = 1 joule/second)
	year = 3.16e7 #seconds.
	collapse_time = Energy /flow / year
	print collapse_time, " HW3.4"


	#5 
	# 5. In this problem we will explore how gravitational forces, while microscopically far weaker than chemical forces, dominate them at large distances. 

	# The key idea is that chemical bonds are, at a microscopic level, much stronger than gravitational attraction. A rock is held together by chemical bonds and not its own gravity. But chemical bonds are short-range forces. A given atom in a rock is bound chemically only to its nearest neighbors. Cutting a rock in half only requires breaking the bonds near the cut. Gravitational attraction is a long-range force, so for a large enough object it dominates. 

	# To get a sense of this, we will ask a related question. To lift a rock out of the Earth, we need to break the chemical bonds at its edge as well as overcome Earth’s gravity to lift it. As the rock gets bigger, the latter force (to overcome gravity) will eventually need to be bigger than the former (to overcome chemical bonds). Our question: How big does a mountain need to be for its weight to be larger than its chemical adhesion to bedrock? 

	# To keep calculations simple, we will consider a mountain that is a pyramid (with four smooth equilateral triangle faces) of side length a, made of rock which has a density of ρ= 3e3 kg/m3, and a tensile strength (the force it takes to break the bonds for a cut of unit area) of 1e7 N/m2. Find the side length a of a mountain for which its weight is equal to the force with which it adheres to bedrock. Write your answer in meters and round to two significant digits. 

	# Useful formulae: The area of a pyramid’s side is A=3√4a2, its volume is 2√12a3, and mass = density * volume.

	p= 3e3 # density , kg/m3

	chem_force =  1.e7 # N/m2
	# A = 3√4a2
	# v = 2√12a3

	# mass = p * v
	# A*F=m*G
	a = 3.**0.5/4.*chem_force / (2.**0.5 / 12. * p * g)
	print a, " HW3.5"

	#7. In fact, the observed temperature of Saturn’s clouds is 93K. This means Saturn radiates more power than it absorbs from the Sun, the difference being made up mostly by Kelvin-Helmholtz processes and latent heat of helium condensation. What fraction of the total power radiated by Saturn is due to these processes? Note that an understanding of the details of Kelvin-Helmholtz processes or helium condensation is not required to solve this problem. Give your answer as a decimal fraction of the total power and round to two significant digits.
	T1 = 93.
	T2 = 77.
	#(P1−P2)/P1
	answer = (T1**4-T2**4)/T1**4
	print answer, " HW3.7"

def hw3_2():
	"1. The planet Kepler-11f orbits its star (Kepler-11) at a radius of 0.25 AU with a period of 46.7 (Earth) days. Use this to estimate the mass of the star it orbits, in terms of the solar mass M⊙. Your answer should be a number - the ratio of the Kepler-11’s mass to that of the Sun - rounded to two significant figures."
	R = 0.25
	P = 46.7
	M_star_k = (R*AU)**3*4*3.14**2/G/(P*60*60*24)**2
	delta = M_star_k / Msun
	print delta, " HW3_1"
	

	"10. The planet Kepler-11e orbits at a radius of 0.194 AU. Predict the period of this planet's orbit, in Earth days. Round your answer to two significant digits."
	R_11e = 0.194 *AU
	"P?"
	P_11e = (4.*3.14**2/G/M_star_k*R_11e**3 ) ** 0.5 /(60*60*24.)
	print P_11e," HW3_2"

	"11. The star Kepler-11 has a surface temperature of 5680 K and a radius of 1.1 R⊙. Find the luminosity of Kepler-11 in terms of the solar luminosity L⊙. Round your answer to two significant figures."
	"L"
	T_11 = 5680.
	R_11 = 1.1
	"L = 4 piR^2b and b = sigma T^4"
	L_11 = R_11**2* T_11 **4 / Tsun**4
	print L_11, " HW3_3"



	"8.Suppose we want to directly image Kepler-11f in visible light. Find the ratio of the planet's luminosity in reflected starlight to that of the star Kepler-11. The planet has an albedo of 0.5."
	"Lp/Lstar = .5Rp^2/4D^2 = Rp^2/8D^2"

	"Rp = 2.61Re andD = .25 AU???"
	R_11f = 2.61 * Re
	D_11f = 0.25 * AU
	answer8 = R_11f**2/8./D_11f**2
	print answer8, " HW3_8"

	"16-B. Also find the maximum angular separation, in arcseconds, between the planet and star as seen from Earth. The star Kepler-11 lies about 1.264×108 AU away from Earth. Round to two significant digits. This problem and the preceding problem should give you a sense of the difficulty faced by direct imaging methods."
	D_11 = 1.264e8 * AU
	AB = D_11f
	alpha = AB/D_11*206265
	print alpha, " HW3_9"


	"10. Vk=(G*M/R)**0.5(Mb*M)"
	R =.091 * AU
	Mb = 4.3 * Me
	
	v_11= (G*M_star_k/R)**0.5*(Mb/M_star_k)
	print v_11, " HW3_10"

	"11. 18. If we were to use a hydrogen absorption line with a wavelength of λ0=656.3 nm in the spectrum of Kepler-11 to detect the motion of its orbit due to the gravitational attraction of Kepler-11b as it circles the center of gravity, what is the longest wavelength, in nanometers at which we would detect this line? Keep as many significant digits as you need to represent the first deviation from 656.3 nm."
	l0=656.3 #nm
	l= l0*(1+v_11/c)
	print l, " HW3_11"
#hw3_1()
#hw3_2()

def hw4_1():
	"2"
	T = 2000000.
	R = 2 * Rsun
	m = mp # proton
	v_termal = (3 * kB * T / m)**0.5
	v_escape = ((2*G*Msun)/R)**0.5
	print v_escape, v_termal, "HW4_2"



	"4"
	D = 3.9
	u = 8.67
	v_t=4.74 * u * D
	print v_t, "HW4_4"

	"5"
	l0 = 592 
	l2 =  592.48
	
	v = (l2-l0/l0)*c
	print v, "HW4_5"

# hw4_1()

def hw4_2():
	"""1 The α (Alpha) Centauri binary star system mentioned in Clip 11 has a period of 79.9 years and an apparent semimajor axis of 17.57 arcsecond. The system has a parallax angle of .747 arcsecond. Find the total mass of the system in terms of the <b>Solar mass</b>. Round to two significant figures."""
	p=79.9 # years
	sm_axis = 17.57 # arcsec
	parallax = .747 # arcsec

	"""Mtotal - ?"""

	Mtotal = ( sm_axis / parallax )**3 * p**-2
	print Mtotal, "HW4_B1"

	"""2. mAvA = mBvB """
	vv = 1.213
	mB = Mtotal / (vv + 1)
	mA = mB * vv
	print mA, mB, " HW4_B2"

	"4. "
	T1 = 9700. 
	T2 = 5800.
	R1 = 3.04 
	R2 = 0.90

	"R/Rsun =( L/Lsun) ^1/2 * (T/Tsun)^-2"

	ratio = pi* (R2)**2/ ( pi*(R1)**2 )
	# print ratio

	L1 = (R1*(T1/Tsun) **2)**2
	# print L1
	L2 = (R2*(T2/Tsun) **2)**2
	# print L2
	Ltotal = L1 + L2
	print (L1-L2)/Ltotal, (L2+(1-ratio)*L1)/Ltotal, "HW4_B4"


	"5. "
	v1 = 35000.
	v2 = 99000.
	vsum = v1 + v2

	t_decrease = (R2*2)*Rsun / vsum
	t_stable = (R1*2 - R2*2)*Rsun / vsum
	print t_stable, t_decrease, "HW4_B5"


# hw4_2()


def hw5_1():
	"1. The Sun, in its T-Tauri phase, may have been losing mass at a rate of 10−8M⊙/yr for up to ten million years, ending with a mass M⊙. As a main sequence star, it loses mass at a rate of about 2×10−14M⊙/yr for ten billion years. As a red giant, it may lose up to 28% of the remaining mass. Estimate, in terms of M, the mass at the start of the T-Tauri phase, the mass of the remaining star at the end of the red giant phase. Round to two significant figures."
	rate_ti_tauri = 1e-8*Msun
	time_ti_tauri = 10000000
	rate_main_sequence = 2e-14*Msun
	time_main_sequence = 10000000000
	red_coeff = 1-0.28
	m_ti_tauri = (Msun + rate_ti_tauri * time_ti_tauri)/Msun
	m_main_sequence = (Msun - rate_main_sequence * time_main_sequence)/Msun
	m_red  = m_main_sequence * red_coeff

	print m_red/m_ti_tauri, "HW5_A1"

	"2."
	Pe = 1.8e15
	Pt = 1.95e16
	print Pe/Pt, "HW5_A2"

	"3."
	Pe = 1.5e18
	Pt = 3.1e18
	print Pe/Pt, "HW5_A3"
	"5. The Cat’s Eye nebula is about 213 pc from Earth and has an angular radius of 8 arc minutes. It observed (from Doppler shifts) to be expanding at 20 km/s. Assuming the nebula is expanding at a constant rate away from a central star, find its age in years. Round to two significant figures."
	d = 213. #pc
	angular_r = 8.*60 #arcsec
	v = 20.*1000 # m/sec
	"in years"
	radius = d*360*60*60/2/pi*AU*angular_r/206265

	age = radius / v /60/60/24/365
	print age, "HW5_A5"

	"6. Find the escape velocity in km/s at:"
# a) The surface of the Sun as an AGB giant with mass 0.7M⊙ and radius 213R⊙.
# b) The surface of an O-type star with a mass 30 M⊙ and a radius 5R⊙.
# c) The Sun’s remnant white dwarf with mass 0.6M⊙ and a radius 0.014R⊙.
# d) A neutron star with a mass of 1M⊙ and radius 20km.
# Your answer should be four numbers, separated by one or more spaces. Round to two significant figures."
	"v = sqrt(2*G*M/R)"
	def v(M,R):
		return (2*G*M*Msun/(R*Rsun))**0.5 /1000.
	a = v(0.7,213.)
	b = v(30.,5.)
	c = v(0.6,0.014)
	d = v(1.,20000./Rsun)
	print a,b,c,d, "HW5_A6"

def hw5_2():
	"""1
		The Crab nebula is observed via Doppler shifts to be expanding at 1500km/s. Comparisons of the angular size of the nebula over years of observation show that the nebula’s angular radius (more precisely, the semi-minor axis of an elliptical approximation to its shape, which is expected to correspond to the expansion in the direction of our line of sight) is growing at a rate of 0.15 arcsec/yr. The apparent semimajor axis was 113.1 arcseconds in 1987. Using these data, find the following. Your answer should be three numbers, separated by one or more spaces.
		a) The distance to the Crab nebula in pc. Round to two significant figures.
		b) The year (AD) in which we would have observed the supernova, assuming the expansion has been uniform since then. Round to three significant figures.
		c) The year (BC) in which the explosion actually occurred, assuming your result from b) were the actual date on which it was observed. You may find the conversion 1pc = 3.26ly to be useful here (a light-year is the distance light travels in a year and is a common measure of astronomical distances). Enter your answer as a positive number for the date BC of the explostion. Round to two significant figures.
	"""

	v = 1500000. # m / sec
	g_rate = 0.15 # arcsec / year
	r0 = 113.1 #arcsec in 1987
	parsec =  30856775814671900. # m
	a = v*206265/(g_rate/(365*24*60*60))/parsec
	b = 2013-r0/g_rate
	year_shift = 3.26*a
	c = year_shift-b
	print a, b, c, "HW5_B1"


	"2"
	p = 5.4 #days
	rate=5.4e-13 # brigtness of sun

	l = 1.e3
	d = l * AU / parsec
	d = (l/rate)**0.5 *AU/ parsec
	print d, "HW5_B2"

	"3"
	ext =  0.37
	l = 100.
	b = 3.5e-15
	b_correct = b/ext
	d = (l / b_correct)**0.5 *AU/parsec
	print d, "HW5_B2"


	"4"
	R = 10*Rsun
	t=3*60*60.
	v = R/t
	print v, "HW5_B4"

	"5"
	d = 427. # l years
	L =  3.e8
	# b = L*Lsun/4/pi/(d*ly)**2/Bsun
	b = L/(d*ly/AU)**2
	print b, "HW5_B5"



# hw5_1()


# hw5_2()

def hw6_1():
	"2"
	"E = mc**2"
	m_muon = 106e6 #eV/c2
	m_electron = 0.51e6 #ev/c2
	j_coeffs = 1.7e-11/106e6
	m_loss = m_muon - m_electron
	e = m_loss * j_coeffs
	print e, "hw6_A2"

	"3"
	e = 10e9 #ev/c2
	e_kin = e-m_muon
	v = e_kin /m_muon
	v = (1-(m_muon/e)**2)**0.5
	print v, "hw6_A3"

	"4"
	h = 20000. # m
	v = v*c
	t = h/v
	print t, "hw6_A4"

	"5"
	t_half = 1.525e-6 # s
	n = 2**(-t/t_half)
	print n, "hw6_A5"

	"6"
	v=0.995*c
	t_rel=(t-v*h/(c**2))/ ((1-v**2/c**2)**0.5)
	print t_rel, "hw6_A6"

	"7"

	n_rel = 2**(-t_rel/t_half)
	print n_rel, "hw6_A7"

	"8"
	x_rel=((h-v*t)/ ((1-v**2/c**2)**0.5))/1000.#km
	print x_rel,"hw6_A8"


def hw6_2():
	"1"
	j = 4.2e6
	m = j/c**2
	print m, "hw6_B1"

	"2"
	m_h = 1.00783 # amu
	m_coeff = 1.66053886e-27

	m_he = 4.002602 # amu

	x = 1/ ( m_h*4/m_he)
	remain = 1-x
	print remain, "hw6_B2"

	"3"
	m_h_in_sun = Lsun /c**2 / remain
	print m_h_in_sun, "hw6_B3"

	"4"
	t = 1e10*365*24*60*60. # years
	m_h_loss = m_h_in_sun * t
	m_h_loss_rate = m_h_loss / Msun
	print m_h_loss_rate, "hw6_B4"

	"5"

	m_loss = m_h_loss_rate * remain
	print m_loss, "hw6_B5"

	"6"
	M = 5 * Msun
	density1 = M / ( 4*pi/3 * (2*G*M/c**2)**3)
	M = 4e6 * Msun
	density2 = M / ( 4*pi/3 * (2*G*M/c**2)**3)
	print round(density1,2), round(density2,2), "hw6_B6"

	"9"
	m = 10 * Msun
	rs = 2.*G*m/c**2
	r = 3. * rs
	l=2.068e-10
	"l0 - ?"
	l0 = l*(1-1/3.)**0.5
	print l0,"hw6_B9"

	"10"
	v=(G*m/r)**0.5
	print v/c, "hw6_B10"

	"11"

	l0=l0*(1-(v/c)**2)**0.5
	print l0, "HW6_b11"





# hw6_2()

def hw7_1():
	"1"
	density=0.054*Msun/pc**3#Msun/pc3
	print density, "hw7_A1"

	"2"
	m_bulge = 1e10*Msun
	r = 4000.*pc
	den_bulge = m_bulge /  (4/3.*pi*r**3)
	print den_bulge, "hw7_a2"

	"3"

	m_dm = 2e12*Msun
	r = 8000.*pc
	m_inside = 1e10*Msun
	r_dm = (m_dm / m_inside * r**3)**0.3333333
	den_dm = m_dm / ( 4./3.*pi*r_dm**3)
	print r_dm/pc/1000, den_dm, "hw7_a3"

	"4"
	m_elem_dm = 100. * mp
	den_elem_dm = den_dm/m_elem_dm
	print den_elem_dm, "hw7_a4"

	"5"

	# v = (G*m_inside/r)**0.5
	v = 220000.
	num_elem_dm = v*den_elem_dm
	print num_elem_dm, "hw7_a5"

	"6"
	e = m_elem_dm*v**2/2.
	print e, "hw7_a6"

	"7"
	m = 6.8e10*Msun


def hw7_2():
    "1"
    l = 403.2  # nm
    l0 = 393.3  # nm
    z = (l - l0) / l0
    print z, "hw7_b1"

    "2"
    v = (c * (1 + z) ** 2 - c) / (1 + (1 + z) ** 2)
    print v / 1000, "hw7_b2"

    "3"
    h0 = 71  # km/sec / Mpc
    d = v / 1000 / h0
    print d, "hw7_b3"

    "4"
    f_shift = 1 + z
    print f_shift, "hw7_b4"

    "5"
    h0y_minus_one = 1.37e10 
    t = - z * h0y_minus_one
    print t, "hw7_b5"

    "6"
    d = c * -t * 365 * 24 * 60 * 60
    print d / pc / 1.e6, "hw7_b6"

    "8"

    m_gal = 5.e11 * Msun
    r_gal = 50000. * pc
    t = (2 * r_gal / (G * m_gal / r_gal ** 2)) ** 0.5
    print t / year_coeff, "hw7_b8"

    "9"
    m_lmc = 2e10 * Msun
    r_lmc = 2.2 * 1000 * pc
    r_orbital = 50 * 1000 * pc

    m_milky = 2e12 * Msun
    a_g = G * m_lmc / r_lmc ** 2
    a_t = (2 * r_lmc) * G * m_milky / r_orbital ** 3
    print a_g / a_t, "hw7_b9"







   
# hw7_2()
def hw8_1():
	z = 5.34
	"1"
	a = 1 / (1+z)
	print a, "hw8_a1"

	"2"
	h0 = 7.26e-11
	t0 = 2/3./h0
	t = a**(3/2.) /3.*2./h0
	print t0-t, "hw8_a2"

	"3"
	km_to_ly = 9.460536207e12
	d =  3*c*t0*year_coeff*(1-a**0.5)
	print d/ly, "hw8_a3"

	"4"
	d0 = d * a
	print d0/ly, "hw8_a4"

	"5"
	r = 50*1000.*pc
	ang_r = 206265*r/(d0)
	print ang_r, "hw8_a5" 

	"6"
	L = 1.e10
	d_l = d * (1+z)
	b = L/(d_l**2/AU**2)
	print b, "hw8_a6"

	"7"
	f = 6
	days = 40
	print days*(1+z)

	"8"
	t0 = 2.7260
	t_em = t0 * (1+z)
	print t_em, "hw8_a8"


def hw8_2():
    "1"
    e_ion = 2.17896E-18
    t = e_ion / kB
    print t, "1"
    "2"
    z = 1100
    t0 = 2.7260
    t_em = t0 * (1 + z)
    print t_em, "2"
    
    "3"
    z1 = 1400
    d = 4.17e-28  # kg/m3
    coeff = 0.75
    h_mass = 1.007
    d_z1 = d*(1+z1)**3

    num_h_density = d_z1*coeff/h_mass/atom_mass
    
    t_z1 = t0*(1+z1)
    e_density = 2/4./c*t_z1**4*sigma
    e_photon = kB*t_z1
    num_p_density = e_density/e_photon

    x = e_ion / kB/t_z1

    fraction = x**5*exp(-x)
    num_p_ion_den  = fraction * num_p_density

    print num_h_density, num_p_density, num_p_ion_den, "3"

    d = 201  # kpc
    d_now = d*(1+z)
    print d_now/1000, "4"

    d = 4.17e-28
    t = (d*c**3/4/sigma)**0.25
    print t, "5"


    z = 5.34
    H_0 = 72.58/1.e12

    d = 1./H_0*z
    print d, "7"
	# c/H_0*z    6
	# c/H_0*exp(-H_0*t)   8
    d = 2.5e6
    t = log(d/1.*H_0)/-H_0-1/H_0
    print t, '9'







hw8_2()