diff --git a/binning_ddxc_given_c_from_raw_domain.py b/binning_ddxc_given_c_from_raw_domain.py
new file mode 100644
index 0000000..ac9524a
--- /dev/null
+++ b/binning_ddxc_given_c_from_raw_domain.py
@@ -0,0 +1,133 @@
+import argparse
+import datetime
+import sys
+import os
+import numpy as np
+import dnstool
+from pycompact import CompactScheme
+
+
+program_description = '''\
+Read all DNS data files and collect c and ddx(c) at all x coordinates.
+Compute conditional mean of ddx(c) given condition c by data binning
+'''
+
+# Commandline argument parser
+parser = argparse.ArgumentParser(description=program_description)
+parser.add_argument("-c", "--case", help="target case name, required")
+parser.add_argument("-l", "--list", help="list all case names", action='store_true')
+parser.add_argument("-t", "--test", help="test run, use only 10 timesteps", action='store_true')
+args = parser.parse_args()
+params = vars(args)
+
+cases = dnstool.case_library()
+
+if params["list"] :
+    for c in cases:
+        print(c, cases[c].case_root)
+    sys.exit(0)
+else:
+    if params["case"] is None :
+        sys.exit("--case CASE not given")
+
+casename = params["case"]
+case = cases[casename]
+nx, ny, nz = case.shape
+
+n_timestep = len(case.data_files)
+if params["test"] :
+    n_timestep = 10
+
+y_idx = 4
+
+x_indices = np.arange(152, 408, dtype=np.int)
+
+print("Calculating Conditional Mean of ddx(c) given c")
+print(" *    Case:", casename)
+print(" * x index:", x_indices)
+
+
+class ConditionalMeanBinning:
+    def __init__ (self, name='', nbins=100, boundary=0.1):
+        self.eps = np.finfo(np.float).eps
+
+        self.name = name
+
+        # Point grid containing interval boundaries
+        self.grid = np.hstack(
+                (np.linspace(0, boundary, nbins),
+                 np.linspace(boundary, 1-boundary, nbins)[1:],
+                 np.linspace(1-boundary, 1, nbins)[1:],
+                )
+                )
+        self.grid[0] = -self.eps # to include c = 0 samples to the first bin
+
+        # interval center points
+        self.cstar = np.zeros(self.grid.size + 1)
+        self.cstar[1:-1] = (self.grid[1:]+self.grid[:-1])/2.
+        self.cstar[-1] = 1.0
+
+        self.v_sums = np.zeros(self.grid.size)
+        self.v_counts = np.zeros(self.grid.size, dtype=int)
+
+    def feed_samples (self, c, v):
+        ''' Feed arrays c, v of samples and conditions '''
+
+        bin_indices = np.digitize(c, self.grid, right=True)
+        # assign bin indices. bin intervals are right inclusive
+        # i-th bin: c_bin[i-1] < c <= c_bin[i]
+
+        for i in range(1, self.grid.size):
+            binned = v[bin_indices == i]
+            self.v_sums[i] += binned.sum()
+            self.v_counts[i] += binned.size
+
+    def average (self):
+        v_avg = np.zeros(self.v_sums.size+1)
+        v_avg[:-1] = self.v_sums / (self.v_counts + self.eps)
+        # add eps to prevent nan values from divide by zero
+
+        d = self.sum_count()
+        d["bin_cstar"] = self.cstar
+        d["bin_avg"] = v_avg
+        return d
+
+    def sum_count (self):
+        return {
+            "bin_grid" : self.grid, 
+            "bin_sum" : self.v_sums, 
+            "bin_count" : self.v_counts
+        }
+
+cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2)
+
+# c mean at sampling x coordinates (verification purpose)
+cmean_verify = np.zeros(len(x_indices))
+
+avg_objs = [ConditionalMeanBinning(name=str(xi)) for xi in x_indices]
+
+for i, fname in enumerate(case.data_files[:n_timestep]):
+    print(datetime.datetime.now(), '{:5d}/{}'.format(i+1, n_timestep), fname)
+
+    fpath = os.path.join(case.case_root, fname)
+    time, y = case.read_single_field(fpath, y_idx)
+    c = (1. - y)
+    ddxc = cs.ddx(c)
+
+    for j, xi in enumerate(x_indices):
+        avg_objs[j].feed_samples(c[:,:,xi], ddxc[:,:,xi])
+        cmean_verify[j] += np.sum(c[:,:,xi])
+print(datetime.datetime.now(), 'Finished')
+
+print()
+print("Verify mean(c) at sampling x coordinates")
+cmean_verify /= (n_timestep * ny * nz)
+for xi_cmean_pair in zip(x_indices, cmean_verify):
+    print (xi_cmean_pair)
+
+print()
+print("Save the result")
+for j, xi in enumerate(x_indices):
+    filename = "cavg_ddxc_c_{}_{}".format(casename, avg_objs[j].name)
+    np.savez(filename, **avg_objs[j].average())
+    print("saved ", filename)