Merge pull request #195 from discsim/debris_interface

Debris interfacing

frank/default_parameters.json

@@ -26,7 +26,8 @@
     "pa"               : 0.0,
     "dra"              : 0.0,
     "ddec"             : 0.0,
-    "rescale_flux"     : true
+    "rescale_flux"     : true,
+    "scale_height"     : null
   },
 
   "hyperparameters" : {

frank/fit.py

@@ -152,6 +152,12 @@ def parse_parameters(*args):
         err = ValueError("method should be 'Normal' or 'LogNormal'")
         raise err
 
+    if model['geometry']['scale_height'] is not None and \
+       model['hyperparameters']['method'] == 'LogNormal':
+        logging.warning("WARNING: 'scale_height' is not 'None' AND 'method' is "
+                        "'LogNormal'. It is suggested to use the 'Normal' method " 
+                        "for vertical inference.")
+
     if model['hyperparameters']['nonnegative'] and \
        model['hyperparameters']['method'] == 'LogNormal':
         logging.warning("WARNING: 'nonnegative' is 'true' AND 'method' is "
@@ -373,6 +379,42 @@ def determine_geometry(u, v, vis, weights, model):
     return geom
 
 
+def get_scale_height(model):
+    """
+    Parse the functional form for disc scale-height in the parameter file
+
+    Parameters
+    ----------
+    model : dict
+        Dictionary containing model parameters the fit uses
+
+    Returns
+    -------
+    scale_height : function R --> H
+        Returns None if scale_height is 'null' in the input parameter file.
+        Else, returns a function for the vertical thickness of disc provided in 
+        the parameter file.
+
+    """
+
+    if model['geometry']['scale_height'] is None:
+        return
+
+    else:
+        if model['geometry']['rescale_flux']:
+            err = ValueError("scale_height should be 'null' if rescale_flux is 'true'")
+            raise err
+
+        def scale_height(R):
+            """Power-law with cutoff, unit=[arcsec]"""
+            vars = model['geometry']['scale_height']
+            h0, a, r0, b = vars['h0'], vars['a'], vars['r0'], vars['b']
+
+            return h0 * R ** a * np.exp(-(R / r0) ** b)
+
+        return scale_height
+
+
 def perform_fit(u, v, vis, weights, geom, model):
     r"""
     Deproject the observed visibilities and fit them for the brightness profile
@@ -404,6 +446,8 @@ def perform_fit(u, v, vis, weights, geom, model):
     need_iterations = model['input_output']['iteration_diag'] or \
         model['plotting']['diag_plot']
 
+    scale_height = get_scale_height(model)
+
     t1 = time.time()
     FF = radial_fitters.FrankFitter(Rmax=model['hyperparameters']['rout'],
                                     N=model['hyperparameters']['n'],
@@ -416,8 +460,9 @@ def perform_fit(u, v, vis, weights, geom, model):
                                     I_scale=model['hyperparameters']['I_scale'],
                                     max_iter=model['hyperparameters']['max_iter'],
                                     store_iteration_diagnostics=need_iterations,
-                                    assume_optically_thick=model['geometry']['rescale_flux'],
                                     convergence_failure=model['hyperparameters']['converge_failure'],
+                                    scale_height=scale_height,
+                                    assume_optically_thick=model['geometry']['rescale_flux']
                                     )
 
     sol = FF.fit(u, v, vis, weights)

frank/parameter_descriptions.json

@@ -26,7 +26,8 @@
     "pa"               : "Position angle, defined east of north, unit=[deg]",
     "dra"              : "Delta (offset from 0) right ascension, defined positive for offsets toward east, unit=[arcsec]",
     "ddec"             : "Delta declination, unit=[arcsec]",
-    "rescale_flux"     : "Whether to rescale the total flux to account for the source's inclination (see frank.geometry.rescale_total_flux)"
+    "rescale_flux"     : "Whether to rescale the total flux to account for the source's inclination (see frank.geometry.rescale_total_flux)",
+    "scale_height"     : "Parameter values for calcuating the vertical thickness of the disk in terms of its (assumed Gaussian) scale-height: 'h0', 'a', 'r0', 'b' in H(r) = h0 * r**a * np.exp(-(r/r0)**b). 'r0' should be in [arcsec]. Example (replace single- with double-quotes): {'h0': -1.0, 'a': 1.0, 'r0': 1.0, 'b': 1.0}",
   },
 
   "hyperparameters" : {
@@ -62,7 +63,7 @@
   "analysis" : {
     "compare_profile"  : "Path of file with comparison profile to be overplotted with frank fit (.txt or .dat). Columns: r [arcsec], I [Jy / sr], optional 3rd column with a single error [Jy / sr] or optional 3rd and 4th columns with lower and upper errors",
     "clean_beam"       : "Dictionary of B_major [arcsec], B_minor [arcsec], position angle [deg] of beam to convolve frank profile",
-    "bootstrap_ntrials": "Number of trials (dataset realizations) to perform in a bootstrap analysis",
+    "bootstrap_ntrials": "Number of trials (dataset realizations) to perform in a bootstrap analysis"
 
   }
 }