<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Transit-Network
Making plots for my SETI idea, dreamed up on the airplane home from AAS 227
SHELVED
I put this idea on the backburner and removed it from the paper because it has an extra layer of confusion
Step1:  The idea is simple
Step2:  a simple 3D model
Step3:  OK, so our toy model works... but how do we actually detect these beacons among the noise of naturally occuring exoplanets we've detected? 
if we had the Kepler exoplanet database
and "injected" this signal (assuming the Kepler cadence was 100% effective in detecting transiting systems), would we see an "over density" or concentration of orbital periods as a funcrtion of location?
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.cm as cm
import matplotlib
matplotlib.rcParams.update({'font.size':18})
matplotlib.rcParams.update({'font.family':'serif'})
Explanation: Transit-Network
Making plots for my SETI idea, dreamed up on the airplane home from AAS 227
SHELVED
I put this idea on the backburner and removed it from the paper because it has an extra layer of confusion:
if the ETI goal is to simply send a 1-bit message ("we're here!") then we don't need an extra unknown relation (i.e. the relation of orbital period and separation). Instead it is easier to find (and probably build) a beacon where all the periods are the same in a region (or if not periods, transit depths, or time of transit, etc)
End of explanation
Alpha = 3./2.
ab_dist = np.array([2., 3.,3.5, 6., 6.7, 7., 8.2, 14., 18., 20.])
ab_per = ab_dist ** (Alpha)
# a figure of the 1-d projection, the key for SETI.
plt.figure(figsize=(5,4))
plt.plot(ab_dist, ab_per, 'k')
plt.scatter(ab_dist, ab_per,color='k')
plt.xlabel('Home $-$ Beacon Separation (pc)')
plt.ylabel('P$_{orb}$ of Beacon (days)')
plt.ylim(0,100)
plt.savefig('dist_per.pdf', dpi=150, bbox_inches='tight', pad_inches=0.25)
# figure angle on sky to plot
# theta = np.random.random(len(ab_dist)) * 2. * np.pi
 # freeze this random config, i like it
theta = np.array([ 4.52448995,  3.46489278,  0.33872438,  1.6891746 ,  2.37611205,
        2.72516744,  5.41764719,  4.01860732,  1.72938583,  0.60279578])
x = ab_dist * np.cos(theta)
y = ab_dist * np.sin(theta)
# a figure of the 2-d observed (sky) plane 
plt.figure(figsize=(6,5))
plt.axes()
# the central red cicrle
circ = plt.Circle((0,0), radius=1.4, fc='r', zorder=0)
plt.gca().add_patch(circ)
# make the concentric circles
for k in range(5,29,3):
    circ = plt.Circle((0,0), radius=k, fc='none', alpha=0.35, color='k')
    plt.gca().add_patch(circ)
plt.scatter(x,y, c=ab_per, cmap=cm.viridis_r, s=90, alpha=0.7, edgecolors='k', zorder=2)
plt.xlim(-20,20)
plt.ylim(-20,20)
plt.xlabel('RA (deg)')
plt.ylabel('Dec (deg)')
cbar = plt.colorbar()
cbar.set_label('P$_{orb}$ of Beacon (days)')
plt.savefig('sky_per.pdf', dpi=150, bbox_inches='tight', pad_inches=0.25)
plt.show()
Explanation: The idea is simple: ET's would place artificial satellites in orbit around stars surrounding their home star system. The orbit period of the remote artificial satellites (or beacons) would be proportional to the distance from the remote star to the home star system. Each remote beacon system would be part of a network map leading to the home system. If the orbital period vs beacon distance relationship was known, the exact location of the home system could be triangulated using only a subset of the beacons.
First toy model (2D)
End of explanation
n_b = 1000 # number of beacon systems
# n_b = 1e6
i_max = 90.
# i_max = 10. # what if we forced beacons to be roughly aligned w/ Galactic plane?
d_min = 2 # min distance from ET home world to place beacons (pc)
d_max = 50 # max distance from ET home world to place beacons (pc)
d_home = 1000 # distance from Earth to ET Home world (in pc)
alpha = 3. / 2. # the coefficient to encode the period-distance relationship - in this case Kepler 3rd law
R_star = 6.955e10 # cm (R_Sun)
R_planet = 7149235 # cm (R_Jupiter)
AU = 1.49598e13 # cm (1 AU)
#__ the part to repeat __
rad = np.random.random(n_b) * (d_max - d_min) + d_min # dist from ET Home to Beacons (pc)
per = rad**Alpha # the period, in days by arbitrary construction
a_AU = (per / 365.)**(2./3.) # the orbital semimajor axis (in AU), assuming solar scaling
incl = np.random.random(n_b) * i_max # orbit plane inclination (deg)
#__
# plt.scatter(a_AU, per, s=90, alpha=0.6)
# plt.xlabel('a (AU)')
# plt.ylabel('Period (days)')
# determine if beacon is "visible", i.e. does it Transit?
b = a_AU * AU * np.sin(incl / 180. * np.pi)
Transit = b < (R_star + R_planet)
no_Transit = b >= (R_star + R_planet)
print(sum(Transit), n_b, float(sum(Transit)) / n_b)
# plt.scatter(a_AU[no_Transit], per[no_Transit], s=20, alpha=0.6, c='blue', lw=0)
# plt.scatter(a_AU[Transit], per[Transit], s=100, alpha=0.6, c='red', lw=0)
# plt.xlabel('a (AU)')
# plt.ylabel('Period (days)')
# plt.xlim(0,2)
# make a plot of fraction of systems that transit as a function of orbital semimajor axis (a)
yy, aa = np.histogram(a_AU[Transit], bins=25, range=(min(a_AU),1))
nn, aa = np.histogram(a_AU[no_Transit], bins=25, range=(min(a_AU),1))
plt.plot((aa[1:] + aa[0:-1])/2., np.array(yy, dtype='float') / nn)
plt.xlabel('a (AU)')
plt.ylabel('Fraction that transit')
# now put beacons in random places in space to illustrate on the sky
theta = np.random.random(n_b) * 2 * np.pi
phi = np.random.random(n_b) * np.pi
x = rad * np.cos(theta)
y = rad * np.sin(theta)
plt.figure(figsize=(5,5))
plt.scatter(x[no_Transit], y[no_Transit], s=10, alpha=0.1)
plt.scatter(x[Transit], y[Transit], s=100, alpha=0.5, c='red')
plt.xlim(-60,60)
plt.ylim(-60,60)
plt.xlabel('RA (deg)')
plt.ylabel('Dec (deg)')
plt.savefig('3d_model.pdf', dpi=150, bbox_inches='tight', pad_inches=0.25)
'''
# repeat this 3d toy model 1000 times to get smooth recovery fraction
hist_smooth = np.zeros_like(yy)
num_smooth  = np.zeros_like(nn)
num_transit = np.zeros(1000)
for k in range(1000):
    rad = np.random.random(n_b) * (d_max - d_min) + d_min # dist from ET Home to Beacons (pc)
    per = rad**alpha # the period, in days by arbitrary construction
    a_AU = (per / 365.)**(2./3.) # the orbital semimajor axis (in AU), assuming solar scaling
    incl = np.random.random(n_b) * i_max # orbit plane inclination (deg)
    b = a_AU * AU * np.sin(incl / 180. * np.pi)
    Transit = b < (R_star + R_planet)
    no_Transit = b >= (R_star + R_planet)
    yy, aa = np.histogram(a_AU[Transit], bins=25, range=(0,2))
    nn, aa = np.histogram(a_AU[no_Transit], bins=25, range=(0,2))
    hist_smooth = hist_smooth + np.array(yy, dtype='float')
    num_smooth = num_smooth +  np.array(nn, dtype='float')
    
#     plt.plot((aa[1:] + aa[0:-1])/2., np.array(yy, dtype='float') / nn, alpha=0.1, c='k')
    num_transit[k] = (float(sum(Transit)) / n_b)
plt.plot((aa[1:] + aa[0:-1])/2., hist_smooth / num_smooth, lw=2);
plt.xlabel('a (AU)');
plt.ylabel('Fraction that transit');
# plt.savefig('recov_fraction.pdf', dpi=150, bbox_inches='tight', pad_inches=0.25)
print(np.mean(num_transit), np.std(num_transit))
''';
Explanation: a simple 3D model
End of explanation
plt.figure(figsize=(5,5))
plt.scatter(x[no_Transit], y[no_Transit], s=10, alpha=0.1)
plt.scatter(x[Transit], y[Transit], s=100, alpha=0.5, c=per[Transit], edgecolors='k', cmap=cm.viridis_r)
plt.xlim(-60,60)
plt.ylim(-60,60)
plt.xlabel('RA (deg)')
plt.ylabel('Dec (deg)')
plt.savefig('3dcolor.pdf', dpi=150, bbox_inches='tight', pad_inches=0.25)
Nother = 200
plt.figure(figsize=(5,5))
# plt.scatter(x[no_Transit], y[no_Transit], s=10, alpha=0.1)
plt.scatter(x[Transit], y[Transit], alpha=0.5, c=per[Transit], cmap=cm.viridis_r)
plt.scatter(np.random.random(Nother)*100-50,np.random.random(Nother)*100-50, 
            c=np.random.random(Nother)*250+5, alpha=.5, s=10)
plt.xlim(-60,60)
plt.ylim(-60,60)
plt.xlabel('RA (deg)')
plt.ylabel('Dec (deg)')
plt.savefig('3dcolor_bkgd.pdf', dpi=150, bbox_inches='tight', pad_inches=0.25)
Explanation: OK, so our toy model works... but how do we actually detect these beacons among the noise of naturally occuring exoplanets we've detected? 
if we had the Kepler exoplanet database
and "injected" this signal (assuming the Kepler cadence was 100% effective in detecting transiting systems), would we see an "over density" or concentration of orbital periods as a funcrtion of location?
End of explanation
<END_TASK>