calc_sbp_excess.py: new; to calculate the central brightness excess value and ratio

1 files changed, 117 insertions, 0 deletions

diff --git a/calc_sbp_excess.py b/calc_sbp_excess.py
new file mode 100755
index 0000000..2b472ca
--- /dev/null
+++ b/calc_sbp_excess.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Aaron LI
+# Created: 2016-04-26
+# Updated: 2016-04-26
+#
+# TODO:
+#   * to estimate errors for the excess value and ratio (Monte Carlo:
+#     repeatedly fit the randomized SBP data and calculate excess)
+#
+
+"""
+Calculate the central brightness excess value and ratio with respect to the
+fitted SBP model (i.e., single-beta model).
+
+NOTE:
+* excess value: brightness_observed - brightness_model_predicted
+* excess ratio: excess_value / brightness_observed
+"""
+
+__version__ = "0.1.0"
+__date__    = "2016-04-26"
+
+
+import sys
+import argparse
+import json
+from collections import OrderedDict
+
+import numpy as np
+from configobj import ConfigObj
+
+from sbp_fit import make_model
+
+
+def calc_excess(data, model):
+    """
+    Calculate the central brightness excess value and ratio with respect
+    to the fitted SBP model (single-beta).
+
+    Arguments:
+      * data: raw 4-column SBP data
+      * model: fitted SBP model
+    """
+    radii            = data[:, 0]
+    radii_err        = data[:, 1]
+    brightness       = data[:, 2]
+    brightness_model = model.f(radii)
+    rin   = radii - radii_err
+    rout  = radii + radii_err
+    areas = np.pi * (rout**2 - rin**2)
+    bsum_obs   = np.sum(brightness * areas)
+    bsum_model = np.sum(brightness_model * areas)
+    excess_value = bsum_obs - bsum_model
+    excess_ratio = excess_value / bsum_obs
+    excess = {
+            "brightness_obs":   bsum_obs,
+            "brightness_model": bsum_model,
+            "excess_value":     excess_value,
+            "excess_ratio":     excess_ratio,
+    }
+    return excess
+
+
+def main():
+    # default arguments
+    default_outfile = "excess.json"
+
+    parser = argparse.ArgumentParser(
+            description="Calculate the central brightness excess value",
+            epilog="Version: %s (%s)" % (__version__, __date__))
+    parser.add_argument("-V", "--version", action="version",
+            version="%(prog)s " + "%s (%s)" % (__version__, __date__))
+    parser.add_argument("config", help="Config file for SBP fitting")
+    parser.add_argument("outfile", nargs="?", default=default_outfile,
+            help="Output json file to save the results " + \
+                 "(default: %s)" % default_outfile)
+    args = parser.parse_args()
+
+    config = ConfigObj(args.config)
+    modelname = config["model"]
+    config_model = config[modelname]
+
+    model = make_model(config, modelname=modelname)
+    print("SBP model: %s" % model.name, file=sys.stderr)
+
+    sbpfit_outfile = config.get("outfile")
+    sbpfit_outfile = config_model.get("outfile", sbpfit_outfile)
+    sbpfit_results = json.load(open(sbpfit_outfile),
+                               object_pairs_hook=OrderedDict)
+
+    # Load fitted parameters for model
+    for pname, pvalue in sbpfit_results["params"].items():
+        model.set_param(name=pname, value=pvalue[0])
+
+    sbpdata = np.loadtxt(config["sbpfile"])
+    excess = calc_excess(data=sbpdata, model=model)
+
+    excess_data = OrderedDict([
+            ("name",             config["name"]),
+            ("obsid",            int(config["obsid"])),
+            ("model",            modelname),
+            ("brightness_obs",   excess["brightness_obs"]),
+            ("brightness_model", excess["brightness_model"]),
+            ("excess_value",     excess["excess_value"]),
+            ("excess_ratio",     excess["excess_ratio"]),
+    ])
+    excess_data_json = json.dumps(excess_data, indent=2)
+    print(excess_data_json)
+    open(args.outfile, "w").write(excess_data_json+"\n")
+
+
+if __name__ == "__main__":
+    main()
+
+#  vim: set ts=4 sw=4 tw=0 fenc=utf-8 ft=python: #