diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md
index fc3962080..03eabf590 100644
--- a/CHANGELOG-unreleased.md
+++ b/CHANGELOG-unreleased.md
@@ -11,6 +11,7 @@ the released changes.
 ### Changed
 - Updated the `plot_chains` function in `event_optimize` so that the subplots are a fixed size to prevent the subplots from being condensed in the case of many fit parameters.
 ### Added
+- Added an option `nbin` to `photonphase` to decide how many phase bins to use for the phaseogram
 - Added an option `linearize_model` to speed up the photon phases calculation within `event_optimize` through the designmatrix.
 - Added AIC and BIC calculation to be written in the post fit parfile from `event_optimize`
 - When TCB->TDB conversion info is missing, will print parameter name
diff --git a/src/pint/scripts/photonphase.py b/src/pint/scripts/photonphase.py
index 7d72eec18..d21679138 100755
--- a/src/pint/scripts/photonphase.py
+++ b/src/pint/scripts/photonphase.py
@@ -106,6 +106,9 @@ def main(argv=None):
         help="Logging level",
         dest="loglevel",
     )
+    parser.add_argument(
+        "--nbin", help="Number of phase bins in the phaseogram", default=100, type=int
+    )
     parser.add_argument(
         "-v", "--verbosity", default=0, action="count", help="Increase output verbosity"
     )
@@ -254,7 +257,7 @@ def main(argv=None):
         print("Htest : {0:.2f} ({1:.2f} sigma)".format(h, h2sig(h)))
 
     if args.plot:
-        phaseogram_binned(mjds, phases, bins=100, plotfile=args.plotfile)
+        phaseogram_binned(mjds, phases, bins=args.nbin, plotfile=args.plotfile)
 
     # Compute orbital phases for each photon TOA
     if args.addorbphase: