diff --git a/eos/options.py b/eos/options.py
index 1e43731..13ed1e3 100644
--- a/eos/options.py
+++ b/eos/options.py
@@ -517,6 +517,7 @@ class OutputFomatOption(StrEnum):
     orb = "orb"
     Rqz = "Rqz"
     Rlt = "Rlt"
+    OS_ort = "OS.ort"
 
 
 class PlotColormaps(StrEnum):
@@ -530,7 +531,7 @@ class PlotColormaps(StrEnum):
 @dataclass
 class ReflectivityOutputConfig(ArgParsable):
     outputFormats: List[OutputFomatOption] = field(
-            default_factory=lambda: ['Rqz.ort'],
+            default_factory=lambda: [],
             metadata={
                 'short': 'of',
                 'group': 'output',
@@ -595,6 +596,9 @@ class ReflectivityOutputConfig(ArgParsable):
                 or OutputFomatOption.Rlt_ort in outputFormat\
                 or OutputFomatOption.Rlt in outputFormat:
             format_list.append(OutputFomatOption.Rlt_ort)
+        if OutputFomatOption.ort in outputFormat\
+                or OutputFomatOption.OS_ort in outputFormat:
+            format_list.append(OutputFomatOption.OS_ort)
         if OutputFomatOption.orb in outputFormat\
                 or OutputFomatOption.Rqz_orb in outputFormat\
                 or OutputFomatOption.Rqz in outputFormat:
diff --git a/eos/projection.py b/eos/projection.py
index c7c3bac..e783bfa 100644
--- a/eos/projection.py
+++ b/eos/projection.py
@@ -13,6 +13,8 @@ from dataclasses import dataclass
 from .event_data_types import EventDatasetProtocol
 from .instrument import Detector, LZGrid
 from .normalization import LZNormalisation
+from .options import IncidentAngle
+
 
 class ProjectionInterface(ABC):
     @abstractmethod
@@ -92,6 +94,39 @@ class ProjectedReflectivity:
         plt.xlabel('Q / Å$^{-1}$')
         plt.ylabel('R')
 
+@dataclass
+class ProjectedOffSpecular:
+    mki_z: np.ndarray
+    kf_z: np.ndarray
+    I: np.ndarray
+    dI: np.ndarray
+
+    @property
+    def data(self):
+        """
+        Return combined data compatible with storing as columns in orso file.
+            Q, R, dR, dQ
+        """
+        ki_sort = np.argsort(self.mki_z[:,0])
+        kf_sort = np.argsort(self.kf_z[0])
+
+        return np.array([self.mki_z[ki_sort, :][:, kf_sort].flatten(),
+                         self.kf_z[ki_sort, :][:, kf_sort].flatten(),
+                         self.I[ki_sort, :][:, kf_sort].flatten(),
+                         self.dI[ki_sort, :][:, kf_sort].flatten()]).T
+
+    @property
+    def n_ki(self):
+        return self.I.shape[0]
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        from matplotlib.colors import LogNorm
+        plt.pcolormesh(self.kf_z-self.mki_z, self.mki_z+self.kf_z, self.I, norm=LogNorm(), **kwargs)
+        plt.xlabel('($k_{f,z}+k_{i,z}$) / Å$^{-1}$')
+        plt.ylabel('($k_{f,z}-k_{i,z}$) / Å$^{-1}$')
+        plt.colorbar()
+
 class LZProjection(ProjectionInterface):
     grid: LZGrid
     lamda: np.ndarray
@@ -110,6 +145,7 @@ class LZProjection(ProjectionInterface):
             ('qx', np.float64),
             ('norm', np.float64),
             ])
+    alphaI = None
 
     def __init__(self, tthh: float, grid: LZGrid):
         self.grid = grid
@@ -129,13 +165,16 @@ class LZProjection(ProjectionInterface):
         self.monitor = 0.
 
     @classmethod
-    def from_dataset(cls, dataset: EventDatasetProtocol, grid: LZGrid, has_offspecular=False):
+    def from_dataset(cls, dataset: EventDatasetProtocol, grid: LZGrid, alpha_i:IncidentAngle=IncidentAngle.alphaF):
         tthh  = dataset.geometry.nu - dataset.geometry.mu
         output = cls(tthh, grid)
         output.correct_gravity(dataset.geometry.detectorDistance)
-        if has_offspecular:
+        if alpha_i==IncidentAngle.mu:
             alphaI_lz = grid.lz()*(dataset.geometry.mu+dataset.geometry.kap+dataset.geometry.kad)
             output.calculate_q(alphaI_lz)
+        elif alpha_i==IncidentAngle.nu:
+            alphaI_lz = grid.lz()*((dataset.geometry.nu+dataset.geometry.kap+dataset.geometry.kad)/2.)
+            output.calculate_q(alphaI_lz)
         else:
             output.calculate_q()
         return output
@@ -144,6 +183,7 @@ class LZProjection(ProjectionInterface):
         self.alphaF += np.rad2deg( np.arctan( 3.07e-10 * detector_distance * self.lamda**2 ) )
 
     def calculate_q(self, alphaI=None):
+        self.alphaI = alphaI
         if alphaI is None:
             self.data.qz = 4.0*np.pi*np.sin(np.deg2rad(self.alphaF))/self.lamda
             self.data.qx = 0.*self.data.qz
@@ -280,6 +320,26 @@ class LZProjection(ProjectionInterface):
         dq_q      = np.where(fltr, dq_q/W_q, np.nan)
         return ProjectedReflectivity(R_q, dR_q, (q_q[1:]+q_q[:-1])/2., dq_q)
 
+    def extract_offspec(self, dalphaF=0.01):
+        if self.alphaI is None:
+            raise ValueError("Off-specular projection requires alphaI to be set")
+        if not self.is_normalized:
+            raise ValueError("Dataset needs to be normalized, first")
+        l = self.grid.lamda()
+        L = (l[1:]+l[:-1])/2.
+        alphaF_center = self.alphaF[0,:].mean()
+        norm_fltr = abs(self.alphaF[0,:]-alphaF_center) < dalphaF
+        lnorm = self.data.norm[:,norm_fltr].sum(axis=1)
+        lnorm = np.where(lnorm>0, lnorm, np.nan)
+        I_OS = self.data.I/lnorm[:,np.newaxis]*self.norm_monitor/self.monitor
+        dI_OS = np.sqrt(self.data.I)/lnorm[:,np.newaxis]*self.norm_monitor/self.monitor
+
+        k = 2.*np.pi/L
+
+        mkiz = k[:,np.newaxis]*np.sin(np.deg2rad(self.alphaI))
+        kfz = k[:,np.newaxis]*np.sin(np.deg2rad(self.alphaF))
+        return ProjectedOffSpecular(mkiz, kfz, I_OS, dI_OS)
+
     def plot(self, **kwargs):
         from matplotlib import pyplot as plt
         from matplotlib.colors import LogNorm
diff --git a/eos/reduction_reflectivity.py b/eos/reduction_reflectivity.py
index a4264ee..05603c8 100644
--- a/eos/reduction_reflectivity.py
+++ b/eos/reduction_reflectivity.py
@@ -84,6 +84,9 @@ class ReflectivityReduction:
 
         # transfer option to data reader
         AmorEventData.correctTimeFallback = self.config.reduction.correctTimeFallback
+        # default export format is specular
+        if self.config.output.is_default('outputFormats'):
+            self.config.output.outputFormats = ["Rqz.ort"]
 
     def reduce(self):
         if not os.path.exists(f'{self.config.output.outputPath}'):
@@ -112,6 +115,7 @@ class ReflectivityReduction:
         # load measurement data and do the reduction
         self.datasetsRqz = []
         self.datasetsRlt = []
+        self.datasetsOS = []
         for i, short_notation in enumerate(self.config.reduction.fileIdentifier):
             self.read_file_block(i, short_notation)
 
@@ -133,6 +137,9 @@ class ReflectivityReduction:
         if 'Rqz.ort' in self.config.output.outputFormats:
             self.save_Rqz()
 
+        if 'OS.ort' in self.config.output.outputFormats:
+            self.save_OS()
+
         if 'Rlt.ort' in self.config.output.outputFormats:
             self.save_Rtl()
 
@@ -257,6 +264,29 @@ class ReflectivityReduction:
                 # plot all reflectivity results in same graph
                 plt.figure(num=99)
                 result.plot(label=f'{self.config.reduction.fileIdentifier[i]}')
+        if 'OS.ort' in self.config.output.outputFormats:
+            result = proj.extract_offspec()
+            columns = [
+                fileio.Column('mki_z', '1/angstrom', 'negative z-component of incident wave vector'),
+                fileio.Column('kf_z', '1/angstrom', 'z-component of outgoing wave vector'),
+                fileio.Column('I', '', 'Scattering intensity'),
+                fileio.ErrorColumn(error_of='I', error_type='uncertainty', value_is='sigma'),
+                ]
+
+            headerOS = self.header.orso_header(columns=columns)
+            headerOS.data_set = f'dataset_{i}'
+            headerOS.data_source.measurement.instrument_settings.incident_angle = fileio.Value(
+                    proj.alphaI[0,0], 'degree')
+            headerOS.data_source.measurement.scheme = {'type': 'off-specular 2D',
+                'n_ki': result.n_ki}
+            orso_data = fileio.OrsoDataset(headerOS, result.data)
+            self.datasetsOS.append(orso_data)
+
+            if self.config.output.plot:
+                import matplotlib.pyplot as plt
+                plt.figure(num=i+100)
+                result.plot(label=f'{self.config.reduction.fileIdentifier[i]}')
+
         if 'Rlt.ort' in self.config.output.outputFormats:
             columns = [
                 fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'),
@@ -392,6 +422,39 @@ class ReflectivityReduction:
             theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | R(q_z)'
             fileio.save_orso(self.datasetsRqz, fname, data_separator='\n', comment=theSecondLine)
 
+    def save_OS(self):
+        fname = os.path.join(self.config.output.outputPath, f'{self.config.output.outputName}.OS.ort')
+        logging.warning(f'    {fname}')
+        theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | OS'
+        #fileio.save_orso(self.datasetsOS, fname, data_separator='\n', comment=theSecondLine)
+        # while there is not official off-specular 2D scheme support in orsopy, export manually
+        from orsopy.fileio.orso import ORSO_DESIGNATE
+
+        with open(fname, "w") as f:
+            header = f"{ORSO_DESIGNATE}\n"
+            if theSecondLine is not None:
+                header += f"# {theSecondLine}\n"
+
+            ds1:fileio.OrsoDataset = self.datasetsOS[0]
+            header += ds1.header()
+            n_kis = ds1.info.data_source.measurement.scheme['n_ki']
+            n_kfs = ds1.data.shape[0]//n_kis
+            np.savetxt(f, ds1.data[:n_kfs], header=header, fmt="%-22.16e")
+            for j in range(1, n_kis):
+                f.write(f"\n### j_ki = {j}\n")
+                np.savetxt(f, ds1.data[(j)*n_kfs:(j+1)*n_kfs],fmt="%-22.16e")
+
+            for dsi in self.datasetsOS[1:]:
+                # write an optional spacer string between dataset e.g. \n
+                f.write("\n\n")
+                hi = ds1.diff_header(dsi)
+                n_kis = dsi.info.data_source.measurement.scheme['n_ki']
+                n_kfs = dsi.data.shape[0]//n_kis
+                np.savetxt(f, dsi.data[:n_kfs], header=hi, fmt="%-22.16e")
+                for j in range(1, n_kis):
+                    f.write(f"\n### j_ki = {j}\n")
+                    np.savetxt(f, dsi.data[(j)*n_kfs:(j+1)*n_kfs], fmt="%-22.16e")
+
     def save_Rtl(self):
         fname = os.path.join(self.config.output.outputPath, f'{self.config.output.outputName}.Rlt.ort')
         logging.warning(f'    {fname}')
@@ -446,8 +509,7 @@ class ReflectivityReduction:
     def project_on_lz(self, dataset=None):
         if dataset is None:
             dataset=self.dataset
-        proj = LZProjection.from_dataset(dataset, self.grid,
-                                         has_offspecular=(self.config.experiment.incidentAngle!=IncidentAngle.alphaF))
+        proj = LZProjection.from_dataset(dataset, self.grid, alpha_i=self.config.experiment.incidentAngle)
 
         t0 = dataset.geometry.nu-dataset.geometry.mu
         if not self.config.reduction.is_default('thetaRangeR'):