MovingUniverseLab · Kerianne28 · Jan 29, 2020 · Jan 31, 2020 · Feb 4, 2020 · Feb 5, 2020
diff --git a/README.rst b/README.rst
@@ -10,7 +10,10 @@ PopSyCLE
 
 Dependencies
 ------------
-`Galaxia <http://galaxia.sourceforge.net>`_
+`galaxia <http://galaxia.sourceforge.net>`_
+PopSyCLE requires a custom version of galaxia in order to support
+user selected galaxy models. Please follow the installation instructions
+found at our galaxia GitHub repo: https://github.com/jluastro/galaxia_.
 
 `PyPopStar <https://pypopstar.readthedocs.io/en/latest/>`_
 

diff --git a/popsycle/add_pbh_plots.py b/popsycle/add_pbh_plots.py
@@ -0,0 +1,189 @@
+import matplotlib
+matplotlib.use('PDF')
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.backends.backend_pdf import PdfPages
+
+def print_plots(output_root, galactic_lin_distance, galactic_lin_b, galactic_lin_l, galactocen_lin_spherical_distance,
+                galactocen_lin_spherical_b, galactocen_lin_spherical_l, rho_lin, r_max, n_lin, cdf_los,
+                x_cyl, y_cyl, r_cyl, r_proj_los_cyl, n_pbh, d_galac, b_galac, l_galac, area_proj_los_cyl,
+                mask_obs_cone, field_of_view_radius, l_radian, b_radian, f_cdf_d, pbh_mass):
+
+    ##################################################################################
+    #Line of Sight - Galactic Coordinates
+    ##################################################################################
+    fig_arr = []
+    fig1, axs = plt.subplots(3, sharex=True, gridspec_kw={'hspace': 0}, figsize=(8,8))
+    fig_arr.append(fig1)
+    axs[0].plot(galactic_lin_distance, '.', label='Distance', c='C0')
+    axs[1].plot(galactic_lin_b, '.', label='Latitude',c='C1')
+    axs[2].plot(galactic_lin_l, '.', label='Longitude', c='C2')
+
+
+    plt.legend(loc='best')
+    # Hide x labels and tick labels for all but bottom plot.
+    for ax in axs:
+        ax.label_outer()
+        ax.legend(loc='upper left')
+
+    axs[0].set_title('Line of Sight Linear Space Galactic Coordinates')
+    axs[0].set_ylabel('[kpc]')
+
+    axs[1].set_ylabel('[deg]')
+
+    axs[2].set_ylabel('[deg]')
+    axs[2].set_xlabel('Line-of-Sight Linear Space Coordinate Index')
+
+    axs[0].set_ylim(0, axs[0].set_ylim()[1])
+
+    ##################################################################################
+    #Line of Sight - Galactocentric Coordinates
+    ##################################################################################
+    fig2, axs = plt.subplots(3, sharex=True, gridspec_kw={'hspace': 0}, figsize=(8,8))
+    fig_arr.append(fig2)
+    axs[0].plot(galactic_lin_distance,
+                galactocen_lin_spherical_distance, '.', label='Distance', c='C0')
+    axs[1].plot(galactic_lin_distance,
+                galactocen_lin_spherical_b, '.', label='Latitude', c='C1')
+    axs[2].plot(galactic_lin_distance,
+                galactocen_lin_spherical_l, '.', label='Longitude',c='C2')
+
+    plt.legend(loc='best')
+    # Hide x labels and tick labels for all but bottom plot.
+    for ax in axs:
+        ax.label_outer()
+        ax.legend(loc='upper left')
+
+    axs[0].set_title('Line of Sight Linear Space\nGalactocentric Coordinates')
+    axs[0].set_ylabel('[kpc]')
+
+    axs[1].set_ylabel('[deg]')
+
+    axs[2].set_ylabel('[deg]')
+    axs[2].set_xlabel('Galactic Coordinate Distance [kpc]')
+
+    axs[0].set_ylim(0, axs[0].set_ylim()[1])
+
+    ##################################################################################
+    #Denisty along Line of Sight
+    ##################################################################################
+    fig3, (ax1, ax2) = plt.subplots(1, 2, figsize=(9,4),
+                                   sharey=True, gridspec_kw={'wspace': 0})
+    fig_arr.append(fig3)
+    # 
+    ax1.semilogy(galactic_lin_distance,
+                 rho_lin, '.',c='k')
+    ax1.set_ylabel('DENSITY Msun*pc^-3')
+    ax1.set_xlabel('Galactic Coordinate Distance [kpc]')
+    ax1.set_xlim(0, ax1.set_xlim()[1])
+    ax2.semilogy(galactocen_lin_spherical_distance,
+                 rho_lin, '.',c='k')
+    ax2.set_xlabel('Galactocentric Distance [kpc]')
+    ax2.set_xlim(0, ax2.set_xlim()[1])
+
+    plt.suptitle('Milky Way Density along Line-of-Sight')
+
+    ##################################################################################
+    #Cumulative Projected Desnity
+    ##################################################################################
+    fig4 = plt.figure(figsize=(10,7))
+    fig_arr.append(fig4)
+
+    plt.plot(galactic_lin_distance,
+             np.cumsum(rho_lin) * r_max / n_lin,
+             '.', label='Density',c='k')
+    plt.legend(loc='best')
+    plt.ylabel('Cumulative Projected Density} [M_sun*pc^{-2}]$')
+    plt.xlabel('Galactic Distance [kpc]')
+    plt.xlim(0, plt.xlim()[1])
+
+    ##################################################################################
+    #Discrete CDF
+    ##################################################################################
+    fig5 = plt.figure(figsize=(10,7))
+    fig_arr.append(fig5)
+    plt.plot(cdf_los, lw=3)
+    plt.ylabel('CDF')
+    plt.xlabel('Line-of-Sight Linear Space Coordinate Index')
+
+    ##################################################################################
+    #Interpolated CDF
+    ##################################################################################
+    fig6 = plt.figure(figsize=(10,7))
+    fig_arr.append(fig6)
+    x = np.linspace(0,1,100)
+    y = f_cdf_d(x)
+    plt.plot(y, x, lw=3)
+    plt.ylabel('CDF')
+    plt.xlabel('Galactic Coordinate Distance along LOS [kpc]')
+
+    ##################################################################################
+    #Cylindrical distribution of PBHs
+    #Infered Radii Distribution
+    #Galactic Distance Distribution
+    ##################################################################################
+    fig7, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14,4))
+    fig_arr.append(fig7)
+    # Create a plot to look at the projected cylindrical distribution of PBH
+    ax1.plot(x_cyl * 1000, y_cyl * 1000, ',')
+    ax1.set_aspect('equal', 'datalim')
+    ax1.set_ylabel('Cylindrical LOS Samples l-axis [pc]')
+    ax1.set_xlabel('Cylindrical LOS Samples b-axis [pc]')
+    # Create a histogram of the infered radii to make sure that it shows the expected linear trend
+    bins = 20
+    ax2.hist(r_cyl * 1000, bins=bins, label='Sampled Number')
+    ax2.plot((0, r_proj_los_cyl * 1000), (0, 2 * n_pbh / bins), c='C2', label='Expected Number')
+    ax2.legend(loc='upper left')
+    ax2.set_ylabel('Number of PBH per LOS Cylindrical Annuli')
+    ax2.set_xlabel('Cylindrical Radius [pc]')
+
+    # Create a histogram of the galactic distance distribution
+    bins=100
+    n_per_bin, _, _ = ax3.hist(d_galac, bins=bins, label='Sampled Number')
+    ax3.plot(galactic_lin_distance,
+             rho_lin * area_proj_los_cyl / pbh_mass * r_max / bins * 1000**3,c='C2', label='Approximate Expectation')
+    ax3.legend(loc='best')
+    ax3.set_ylabel('Number of PBH')
+    ax3.set_xlabel('Galactic Distance [kpc]')
+
+    fig7.tight_layout(pad=2.0)
+
+    plt.suptitle('Distribution of PBHs in Cylindrical Line-of-Sight Tube (i.e. not light cone)')
+
+    ##################################################################################
+    #Checking Light Cone Boundaries
+    ##################################################################################
+    fig8, (ax1, ax2) = plt.subplots(2, 1, figsize=(5,8),
+                                   sharex=True, gridspec_kw={'hspace': 0})
+    fig_arr.append(fig8)
+    ax1.plot(d_galac[~mask_obs_cone], r_cyl[~mask_obs_cone] * 1000,
+             '.', alpha=0.1, label='Outside Light Cone')
+    ax1.plot(d_galac[mask_obs_cone], r_cyl[mask_obs_cone] * 1000,
+             '.', alpha=0.1, label='Inside Light Cone')
+    ax1.plot((0,r_max),
+             (0, field_of_view_radius * np.pi / 180 * r_max *1000),
+             label='Light Cone Boundry')
+    ax1.legend(loc='upper left')
+    ax1.set_ylabel('Radius on LOS Cylinder [pc]')
+
+    ax2.plot(d_galac[~mask_obs_cone], np.sqrt((l_galac[~mask_obs_cone] - l_radian)**2 + 
+                                     (b_galac[~mask_obs_cone] - b_radian)**2) * 180/np.pi,
+             '.', alpha=0.1, label='Outside Light Cone' )
+    ax2.plot(d_galac[mask_obs_cone], np.sqrt((l_galac[mask_obs_cone] - l_radian)**2 +
+                                    (b_galac[mask_obs_cone] - b_radian)**2) * 180/np.pi,
+             '.', alpha=0.1, label='Inside Light Cone')
+    ax2.plot((0,r_max),
+             (field_of_view_radius, field_of_view_radius), label='Light Cone Boundry')
+    ax2.set_ylim(0,5*field_of_view_radius)
+    ax2.legend(loc='best')
+    ax2.set_ylabel('Field of View Radius [deg]')
+    ax2.set_xlabel('Galactic Coordinate Distance [kpc]')
+
+    ##################################################################################
+    #Saving all plots to PDF
+    ##################################################################################
+    pp = PdfPages(output_root+'_pbh_diagnostic_plots.pdf')
+    for fig in fig_arr:
+        pp.savefig(fig)
+        plt.close(fig)
+    pp.close()
diff --git a/popsycle/converter.py b/popsycle/converter.py
@@ -47,6 +47,7 @@ def convert_h5_array_dtype_to_compound_dtype(hdf5_file):
     hdf5_file_new = hdf5_file.replace('.h5', '.compound_dtype.h5')
     f_in = h5py.File(hdf5_file, 'r')
     f_out = h5py.File(hdf5_file_new, 'w')
+    f_out['add_pbh'] = False
 
     # Looping over all of the datasets
     for key in f_in:

diff --git a/popsycle/data/pbh_config.yaml b/popsycle/data/pbh_config.yaml
@@ -0,0 +1,9 @@
+# Example configuration parameters for PBHs
+
+fdm: 1 
+pbh_mass: 40
+r_max: 16.6 
+r_s: 18.6 
+gamma: 1 
+v_esc: 550 
+rho_0 : 0.0093
diff --git a/popsycle/data/popsycle_config.yaml b/popsycle/data/popsycle_config.yaml
@@ -6,6 +6,7 @@ n_obs: 101
 theta_frac: 2
 blend_rad: 0.75
 isochrones_dir: /Users/myself/popsycle_isochrones
+galaxia_galaxy_model_filename: /Users/myself/galaxia_galaxy_model_filename
 bin_edges_number: None
 BH_kick_speed_mean: 50
 NS_kick_speed_mean: 400

diff --git a/popsycle/data/radial_velocity_profile_middle.csv b/popsycle/data/radial_velocity_profile_middle.csv
@@ -0,0 +1,62 @@
+r,v
+0.010792642776929309, 111.71180517031452
+0.012590321940393834, 118.88165941281024
+0.014687431970009386, 126.54067684534101
+0.01713384764066691, 134.0005466632754
+0.01993184060733583, 142.11168012051704
+0.023448011680987736, 150.78286434254127
+0.027455988964315874, 159.3571491494456
+0.03213404568533922, 168.4661816124074
+0.03764430047370784, 177.73468233692108
+0.043859894720648336, 186.54417763189218
+0.05108586704256386, 196.01349797292718
+0.05934511962521002, 204.67635482054385
+0.07083192093472244, 215.1824749096273
+0.08253367166054373, 223.90379647052103
+0.09684442649610563, 232.571184061413
+0.11350401149180162, 242.54850927105662
+0.13427738741960332, 251.89354878059038
+0.1586056689369005, 261.0497880274381
+0.186179188876288, 269.5935014641335
+0.21617856025114549, 276.95879153868964
+0.2526914812832889, 283.38932791498735
+0.2962593087128887, 290.36094394136614
+0.35792137580603484, 297.77197612125576
+0.4309060654759373, 304.89151975943804
+0.5059594624841087, 309.244695083941
+0.5778483816986392, 312.37804867591933
+0.6816933312568133, 314.7048135621079
+0.927699047785235, 320.24095532420716
+1.0671697449159139, 319.8132280123578
+1.4727675574625494, 314.96530138192946
+1.7718422284532012, 311.9846498190986
+2.1197411537756987, 307.83085241102907
+2.4598627413134806, 301.78571980162417
+2.8847016095688685, 295.5628794127563
+3.369909155131431, 289.9952727173605
+4.037236963747592, 283.1343596020251
+4.820932000343341, 279.00220082094125
+5.650248971118326, 275.5083826876096
+6.703101332022676, 271.60158332331133
+7.770476112914143, 265.4265828672343
+9.096562802640236, 259.22858778705637
+10.611733993556404, 252.9216664462353
+12.3792800416127, 246.90957929498413
+14.44123782613885, 240.62655413150787
+16.846646109472406, 235.17599223163336
+19.652711807301024, 229.3228663966612
+22.766181747138315, 223.97434631744974
+26.744873675700347, 218.17494681231219
+31.068822348129252, 212.83026716581725
+36.396415667350134, 206.61996451113927
+43.05764653615527, 199.86831942509124
+50.22955388134484, 193.3372953786207
+58.59605171407531, 186.69858418697999
+68.54786206569099, 179.3094384837924
+79.24703527090092, 171.75301799473309
+93.02409457669958, 163.93495203048826
+108.51867546639197, 155.9656915244779
+126.5941149824405, 147.0301382654459
+147.68029446830633, 138.0518741736095
+172.27869855775518, 128.43932707311797
+192.7048730125345, 121.55334325539347
diff --git a/popsycle/data/radial_velocity_profile_shallow.csv b/popsycle/data/radial_velocity_profile_shallow.csv
@@ -0,0 +1,66 @@
+r,v
+0.010792642776929309, 204.26643023097807
+0.012590321940393834, 210.92487127516114
+0.014687431970009386, 218.17494681231219
+0.01713384764066691, 225.4158429807519
+0.019987751131241467, 232.89705353032554
+0.023317015749352955, 240.21341604790845
+0.02720082013657757, 247.6177427638289
+0.031731531344141896, 255.39654951292243
+0.03701690156357718, 262.5159561489671
+0.04318263075634799, 269.98842786464206
+0.05037535612850712, 277.3556810971669
+0.05876613954792163, 284.43477324360026
+0.0685545358439959, 291.3605743829068
+0.07997333874472501, 298.11329956878745
+0.0932941173219056, 304.6732953099007
+0.10883367461568658, 311.0211336366553
+0.1269615820414861, 316.5932048851325
+0.14810896876722457, 322.0805612681069
+0.17277877509531006, 326.913144325661
+0.20155771370168304, 331.43832299459757
+0.235130223201544, 335.44920768526305
+0.2742947458941446, 338.73163728292803
+0.3199827168140894, 342.046186032515
+0.3732807156985778, 343.8140575474295
+0.4354563087024429, 346.3837933152392
+0.5079881944447303, 346.3837933152392
+0.5926013713388398, 346.3837933152392
+0.6913081625775506, 346.3837933152392
+0.8064560744546262, 343.61717701711524
+0.9407836262193813, 341.654561829232
+1.0974854792445174, 338.15005999198166
+1.2802884144496132, 334.10688186636213
+1.4935399649225864, 329.73408677356304
+1.7423118116552716, 324.4878564913184
+2.03252040141478, 318.41185653421763
+2.3710676553598873, 312.4496290180272
+2.7660051148222395, 306.07263557365945
+3.2267254280694635, 299.1396191737312
+3.764185732107431, 292.02890586026
+4.391168242126657, 285.9048791942064
+5.122584246093038, 278.7891986399261
+5.975828734271417, 271.383867634092
+6.971194097701775, 265.692780180841
+8.132352734462563, 258.3391991378944
+9.48692004136333, 253.79240406567155
+11.067111180485346, 247.90157653062977
+12.970922127990187, 241.5318118152786
+15.060949011985775, 236.1208479753841
+17.618863818498795, 229.54628085055066
+20.496060936083154, 223.2314835842758
+23.909992300202614, 217.3019029730845
+27.892565970532186, 211.4086969647845
+32.53849799081271, 206.147098218244
+37.958280805597894, 199.9828369269898
+44.28081106028432, 193.7807781680185
+51.65645510129202, 186.37803645873726
+60.26062508202531, 179.36080012164197
+70.29795072379233, 172.80561961785477
+82.0071459470937, 165.2548659942107
+95.6666861145325, 158.0340430777453
+111.60143090506897, 149.06535054405288
+130.19035032893146, 139.90937104204863
+151.8755376280755, 130.51589873409813
+177.1727234125981, 121.49765534647908
+195.4228471568372, 116.70016866901982