Skip to content

Commit

Permalink
Final set of code updates for VVmax paper
Browse files Browse the repository at this point in the history
  • Loading branch information
Clancy James committed Jul 18, 2024
1 parent 9616401 commit 8d83f02
Show file tree
Hide file tree
Showing 16 changed files with 3,175 additions and 0 deletions.
36 changes: 36 additions & 0 deletions papers/Vvmax/171020_dist.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
"""
Calculates primary beam distance for this FRB. Use this for beam values B!
Although likely this is incorrect... (because the primary beam value is
highly uncertain)
"""

import numpy as np

def main():
#detection beam was beam 00
bra = 333.216312589
bdec = -19.989420334
bx,by,bz = calc_xyz_deg(bra,bdec)

hra = 22*15. +15*15/60. + 24.61*15/3600
hdec = -19 - 35/60. - 4/3600.
hx,hy,hz = calc_xyz_deg(hra,hdec)

ctheta = bx*hx + by*hy + bz*hz
dtheta = np.arccos(ctheta)
dtheta_deg = dtheta * 180./np.pi
print(dtheta_deg)

def calc_xyz_deg(ra,dec):
"""
"""
ra *= np.pi/180.
dec *= np.pi/180.

x = np.cos(ra)*np.cos(dec)
y = np.sin(ra)*np.cos(dec)
z = np.sin(dec)
return x,y,z

main()
30 changes: 30 additions & 0 deletions papers/Vvmax/Data/all_loc_frb_data_w171020.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
# FRB S/N DM DMgal Freq dt w theta_d Fnu Zloc Nants RF
20171020 19.5 114.1 38.4 1297.5 1.260 1.70 0.722 200 0.00867 1 200
20180924 21.1 362.4 40.5 1297.5 0.864 1.76 0.23 16 0.32140 24 18.4
20181112 19.3 589.0 40.2 1297.5 0.864 2.10 0.31 26 0.47550 12 28
20190102 14.0 364.5 57.3 1271.5 0.864 1.70 0.23 14 0.29000 23 16
20190608 16.1 339.5 37.2 1271.5 1.728 6.00 0.36 26 0.11780 25 28
20190611 9.5 322.2 57.6 1271.5 1.728 2.00 0.48 10 0.37800 25 10
20190711 23.8 594.6 56.6 1271.5 1.728 6.50 0.08 34 0.52200 29 36
20190714 10.7 504.7 38.5 1271.5 1.728 2.90 0.46 12 0.20900 28 12
20191001 37.1 506.9 44.2 919.5 1.728 4.20 0.58 143 0.23000 30 120
20191228 22.9 297.5 32.9 1271.5 1.728 2.30 0.49 40 0.24300 28 67
20200430 13.9 380.1 27.0 864.5 1.728 6.50 0.48 35 0.16100 26 35
20200906 10.5 577.8 35.9 864.5 1.728 6.00 0.65 59 0.36879 7 53
20210117 27.1 730.0 34.4 1271.5 1.182 3.2 0.464 36 0.21400 25 36
20210320 15.3 384 42 864.5 1.182 5.4 1.12 -1 0.28 26 59
20210807 47.1 251.9 121.2 920.5 1.182 10.00 0.453 113 0.12969 23 100
20211127 37.9 234.8 42.5 1271.5 1.182 1.41 0.16 31 0.04695 24 35
20211203 14.2 636.2 63.4 920.5 1.182 9.60 0.212 28 0.34386 24 30
20211212 12.8 206.0 27.1 1632.5 1.182 2.70 0.77 36 0.07150 24 131
20220105 9.8 583.0 22.0 1632.5 1.182 2.00 0.443 16 0.27850 22 19
20220501 16.1 449.5 30.6 863.5 1.182 6.50 0.516 30 0.38100 23 32
20220610 29.8 1458.1 31.0 1271.5 1.182 5.60 0.073 43 1.01600 22 47
20220725 12.7 290.4 30.7 920.5 1.182 4.10 1.272 100 0.19260 25 72
20220918 26.4 657 41 1271.5 1.182 7.1 0.457 -1 0.491 25 55
20230526 22.1 361.4 50 1271.5 1.182 4.7 0.383 -1 0.157 22 34
20230708 31.5 411.5 50.2 920.5 1.182 23.6 0.657 -1 0.105 23 111
20230718 10.9 477 395.6 1271.5 1.182 3.5 0.317 -1 0.035 22 14
20230731 16.6 701 547.1 1271.5 1.182 3.5 0.51 -1 -1 25 25
20230902 11.8 440 34.3 831.5 1.182 5.9 0.63 -1 0.3619 22 23
20231226 36.7 329.9 38 863.5 1.182 11.8 0.739 -1 0.1569 22 78
64 changes: 64 additions & 0 deletions papers/Vvmax/Data/all_unloc_frb_data.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
# FRB DM Freq S/N S/Nth Dmg dt w theta_d Fnu Zmac Nants RF
20170107 610 1320.5 16 9.5 37 1.26 2.4 0.163 58 0.449 1 58
20170416 523 1320.5 13.1 9.5 40 1.26 5 0.332 96 0.369 1 96
20170428 992 1320.5 10.5 9.5 40 1.26 4.4 0.041 34 0.786 1 34
20170712 313 1296.5 12.7 9.5 39 1.26 3.5 0.281 52 0.174 1 52
20170707 235 1296.5 9.5 9.5 36 1.26 1.4 0.133 54 0.101 1 54
20170906 390 1296.5 17 9.5 39 1.26 2.5 0.237 74 0.247 1 74
20171003 463 1297.5 13.8 9.5 41 1.26 2 0.342 82 0.313 1 82
20171004 304 1297.5 10.9 9.5 39 1.26 2 0.203 44 0.166 1 44
20171019 461 1297.5 23.4 9.5 37 1.26 5.4 0.379 219 0.314 1 219
20171020 114 1297.5 19.5 9.5 38 1.26 1.7 0.630 200 -0.024 1 200
20171116 619 1297.5 11.8 9.5 36 1.26 3.2 0.346 64 0.458 1 64
20171213 159 1297.5 25.1 9.5 37 1.26 1.5 0.513 133 0.023 1 133
20171216 203 1297.5 8 8 37 1.26 1.9 0.491 40 0.068 1 40
20180110 716 1297.5 35.6 9.5 39 1.26 7.88 0.430 422 0.542 1 422
20180119 403 1297.5 15.9 9.5 36 1.26 2.7 0.298 110 0.262 1 110
20180128 441 1297.5 12.4 9.5 32 1.26 2.9 0.158 51 0.301 1 51
20180128 496 1297.5 9.6 9.5 41 1.26 2.3 0.396 66 0.343 1 66
20180130 344 1297.5 10.3 9.5 39 1.26 4.1 0.380 95 0.203 1 95
20180131 658 1297.5 13.8 9.5 40 1.26 4.5 0.391 100 0.490 1 100
20180212 168 1297.5 18.3 9.5 31 1.26 1.81 0.391 96 0.038 1 96
20180315 479 1297.5 10.5 9.5 101 1.26 2.4 0.395 56 0.272 1 56
20180324 431 1297.5 9.8 9.5 64 1.26 4.3 0.494 71 0.262 1 71
20180417 475 1297.5 17.5 9.5 26 1.26 2.3 0.458 49 0.337 1 49
20180430 264 1297.5 28.2 9.5 169 1.26 1.2 0.392 177 -0.005 1 177
20180515 355 1297.5 12.1 9.5 33 1.26 1.9 0.191 46 0.220 1 46
20180525 388 1297.5 27.4 9.5 31 1.26 3.8 0.510 300 0.253 1 300
20180924 362 1297.5 21.1 9.5 41 0.864 1.76 0.234 16 0.219 24 18.4
20181112 589 1297.5 19.3 9.5 40 0.864 2.1 0.312 26 0.428 12 28
20190102 365 1271.5 14 9.5 57 0.864 1.7 0.229 14 0.205 23 16
20190608 340 1271.5 16.1 9.5 37 1.728 6 0.362 26 0.201 25 28
20190611 322 1271.5 9.5 9.5 58 1.728 2 0.481 10 0.165 25 10
20190711 595 1271.5 23.8 9.5 57 1.728 6.5 0.079 34 0.418 29 36
20190714 505 1271.5 10.7 9.5 39 1.728 2.9 0.458 12 0.353 28 13
20191001 507 919.5 37.1 9.5 44 1.728 4.2 0.579 143 0.350 30 120
20191228 298 1271.5 22.9 9.5 33 1.728 2.3 0.486 40 0.165 28 67
20200430 380 864.5 13.85 9.5 27 1.728 6.5 0.475 35 0.249 26 35
20200627 294 920.5 11.0 9.5 40 1.728 11 0.34 -1 -1 31 28
20200906 578 864.5 10.5 9.5 36 1.728 6 0.57 59 0.422 7 53
20210117 730 1271.5 27.1 9.5 34 1.182 3.2 0.292 36 0.559 25 36
20210214 398 1271.5 11.6 9.5 32 1.182 3.5 0.442 -1 -1 26 13
20210320 384 864.5 15.3 9.5 42 1.182 5.4 1.146 -1 -1 24 59
20210407 1785 1271.5 19.1 9.5 154 1.182 6.6 0.37 -1 -1 24 36
20210807 252 920.5 47.1 9.5 121 1.182 10 0.602 113 0.032 23 100
20210809 652 920.5 16.8 9.5 190 1.182 14 0.190 45 0.349 23 45
20210912 1235 1271.5 31.7 9.5 31 1.182 5.5 0.475 90 1.010 23 70
20211127 235 1271.5 37.9 9.5 43 1.182 1.4 0.161 31 0.094 24 35
20211203 636 920.5 14.2 9.5 63 1.182 9.6 0.212 28 0.449 24 30
20211212 206 1632.5 12.8 9.5 27 1.182 2.7 0.772 36 0.081 24 131
20220105 583 1632.5 9.8 9.5 22 1.182 2.0 0.443 16 0.439 22 19
20220501 450 863.5 16.1 9.5 31 1.182 6.5 0.516 30 0.310 23 32
20220531 727 1271.5 9.7 9.5 70 1.182 11.0 0.790 30 0.525 23 30
20220610 1458 1271.5 29.8 9.5 31 1.182 5.6 0.073 43 1.212 22 47
20220725 290 920.5 12.7 9.5 31 1.182 4.1 1.272 100 0.160 25 72
20220918 657 1271.5 26.4 9.5 41 1.182 7.1 0.457 55 0.488 25 55
20221106 344 1631.5 35.1 9.5 35 1.182 5.7 0.361 80 0.207 21 80
20230521 640.2 831.5 15.2 9.5 41.8 1.182 16.5 0.401 -1 -1 23 34
20230526 361.4 1271.5 22.1 9.5 50 1.182 4.7 0.383 -1 -1 22 34
20230708 411.5 920.5 31.5 9.5 50.2 1.182 23.6 0.657 -1 -1 23 111
20230718 477 1271.5 10.9 9.5 395.6 1.182 3.5 0.317 -1 -1 22 14
20230731 701 1271.5 16.6 9.5 547.1 1.182 3.5 0.51 -1 -1 25 25
20230902 440 831.5 11.8 9.5 34.3 1.182 5.9 0.63 -1 -1 22 23
20231006 509.7 863.5 15.2 9.5 67.5 1.182 8.3 0.534 -1 -1 24 25
20231226 329.9 863.5 36.7 9.5 38 1.182 11.8 0.739 -1 -1 22 78
29 changes: 29 additions & 0 deletions papers/Vvmax/Data/unbiased_loc_data.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
# FRB S/N DM DMgal Freq dt w theta_d Fnu Zloc Nants RF
20180924 21.1 362.4 40.5 1297.5 0.864 1.76 0.23 16 0.32140 24 18.4
20181112 19.3 589.0 40.2 1297.5 0.864 2.10 0.31 26 0.47550 12 28
20190102 14.0 364.5 57.3 1271.5 0.864 1.70 0.23 14 0.29000 23 16
20190608 16.1 339.5 37.2 1271.5 1.728 6.00 0.36 26 0.11780 25 28
20190611 9.5 322.2 57.6 1271.5 1.728 2.00 0.48 10 0.37800 25 10
20190711 23.8 594.6 56.6 1271.5 1.728 6.50 0.08 34 0.52200 29 36
20190714 10.7 504.7 38.5 1271.5 1.728 2.90 0.46 12 0.20900 28 12
20191001 37.1 506.9 44.2 919.5 1.728 4.20 0.58 143 0.23000 30 120
20191228 22.9 297.5 32.9 1271.5 1.728 2.30 0.49 40 0.24300 28 67
20200430 13.9 380.1 27.0 864.5 1.728 6.50 0.48 35 0.16100 26 35
20200906 10.5 577.8 35.9 864.5 1.728 6.00 0.65 59 0.36879 7 53
20210117 27.1 730.0 34.4 1271.5 1.182 3.2 0.464 36 0.21400 25 36
20210320 15.3 384 42 864.5 1.182 5.4 1.12 -1 0.28 26 59
20210807 47.1 251.9 121.2 920.5 1.182 10.00 0.453 113 0.12969 23 100
20211127 37.9 234.8 42.5 1271.5 1.182 1.41 0.16 31 0.04695 24 35
20211203 14.2 636.2 63.4 920.5 1.182 9.60 0.212 28 0.34386 24 30
20211212 12.8 206.0 27.1 1632.5 1.182 2.70 0.77 36 0.07150 24 131
20220105 9.8 583.0 22.0 1632.5 1.182 2.00 0.443 16 0.27850 22 19
20220501 16.1 449.5 30.6 863.5 1.182 6.50 0.516 30 0.38100 23 32
20220610 29.8 1458.1 31.0 1271.5 1.182 5.60 0.073 43 1.01600 22 47
20220725 12.7 290.4 30.7 920.5 1.182 4.10 1.272 100 0.19260 25 72
20220918 26.4 657 41 1271.5 1.182 7.1 0.457 -1 0.491 25 55
20230526 22.1 361.4 50 1271.5 1.182 4.7 0.383 -1 0.157 22 34
20230708 31.5 411.5 50.2 920.5 1.182 23.6 0.657 -1 0.105 23 111
20230718 10.9 477 395.6 1271.5 1.182 3.5 0.317 -1 0.035 22 14
20230731 16.6 701 547.1 1271.5 1.182 3.5 0.51 -1 -1 25 25
20230902 11.8 440 34.3 831.5 1.182 5.9 0.63 -1 0.3619 22 23
20231226 36.7 329.9 38 863.5 1.182 11.8 0.739 -1 0.1569 22 78
37 changes: 37 additions & 0 deletions papers/Vvmax/README.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
This code is the calculation of V/Vmax for [paper link]

It runs with four sets of data:
- All localised FRBs
- All FRBs used z(DM)
- Unbiased localised FRBs (max z = 0.7) ("UnbiasedLocalised")
- Unlocalised FRBs setting a minimum z ("Minz")


###### calc_vvmax.py #######

This will generate V and Vmax values for FRBs, and store them in files
in a defined output directory (e.g. v2Output)
the various combinations are:
- Localised or Macquart DM
- Min redshift z
- Max redshift z
- value of spectral index alpha
- value of population evolution scaling

Output in "[label]Output/"

##### calc_lf_errors.py #####

Runs a MC rountine to calculate errors in the luminosity function
Inputs from above, outputs in another directory
("[label]LumData")


##### plotting! ######

Plots for the paper are:


- make_paper_unbiased_lum_plot
Plots Minz and UnbiasedLocalised data
Shows effects of systematics (not strong! Interesting...)
125 changes: 125 additions & 0 deletions papers/Vvmax/RLFs.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
import numpy as np
from scipy.integrate import quad

global E0
# normalisation energy - no physical meaning, just makes units more sensible
E0=1e22


def make_first_guess(E,L, params,method):
"""
Makes a first guess by setting the amplitude
to that of the first bin
Params are always gamma,Emax
"""
if method==0:
mag = Schechter(E[0],1.,params[0],params[1])
guess = [L[0]/mag,params[0],params[1]]
elif method==1:
mag = integrateSchechter(E[0],1.,params[0],params[1])
guess = [L[0]/mag,params[0],params[1]]
elif method==2:
# returns the log10 magnitude expected with an amplitude of unity
mag = logIntegrateSchechter(np.log10(E[0]),1.,params[0],params[1])
guess = [L[0]/10.**mag,params[0],params[1]]
else:
mag = power_law(E[0],1.,params[0],params[1])
guess = [L[0]/mag,params[0],params[1]]


return guess




def logIntegrateSchechter(log10E,A,gamma,log10Emax):
"""
Returns the logarithms of the IntegrateSchechter function
Requires log10 E to be sent
"""
res = integrateSchechter(10.**log10E,A,gamma,log10Emax)
res = np.log10(res)
return res

def logIntegratePowerLaw(log10E,A,gamma):
"""
Returns the logarithms of the IntegrateSchechter function
Requires log10 E to be sent
"""
res = integratePowerLaw(10.**log10E,A,gamma)
res = np.log10(res)
return res

def integrateSchechter(E,A,gamma,log10Emax):
"""
Schechter function
Gamma is differential
Relative to E0, which is arbitrary
An extra E/E0 due to log binning of histogram
"""
Emax = 10**log10Emax
bw = 10**0.5
if isinstance(E,np.ndarray):
y=np.zeros(E.size)
for i,bc in enumerate(E):
res = quad(Schechter,bc/bw,bc*bw,args=(A,gamma,Emax))
y[i]=res[0]
else:
res = quad(Schechter,E/bw,E*bw,args=(A,gamma,Emax))
y = res[0]
return y

def integratePowerLaw(E,A,gamma):
"""
Schechter function
Gamma is differential
Relative to E0, which is arbitrary
An extra E/E0 due to log binning of histogram
"""
bw = 10**0.5
if isinstance(E,np.ndarray):
y=np.zeros(E.size)
for i,bc in enumerate(E):
res = A * ((bc*bw)**(gamma+1) - (bc/bw)**(gamma+1))/(gamma+1)
#res = quad(Schechter,bc/bw,bc*bw,args=(A,gamma,Emax))
y[i]=res
else:
#res = quad(Schechter,E/bw,E*bw,args=(A,gamma,Emax))
res = A * ((E*bw)**(gamma+1) - (E/bw)**(gamma+1))/(gamma+1)
y = res
return y

def Schechter(E,A,gamma,Emax):
"""
Schechter function
Gamma is differential
Relative to E0, which is arbitrary
An extra E/E0 due to log binning of histogram
"""
global E0
return A*(E/E0)**gamma * np.exp(-E/Emax)

def best_fits(JHz = True):
"""
Returns best-fit parameters of the Schechter function according to
different papers.
Only gives gamma and Emax
JHz means divide by 1e7 (erg) and 1e9 (GHz bandwidth)
"""
JamesFit = [-1.95,41.26]
RyderFit = [-1.95,41.7]
ShinFit = [-2.3,41.38]
LuoFit = [-1.79,41.46]
for item in [JamesFit,RyderFit,ShinFit,LuoFit]:
#item[1] = 10.**item[1]
if JHz:
# convert from erg to J Hz
item[1] = item[1] - 16
return JamesFit,RyderFit,ShinFit,LuoFit

Loading

0 comments on commit 8d83f02

Please sign in to comment.