diff --git a/README.md b/README.md
index efdc345fcb6f5050e05854fc24dabe67f76a0f10..22d5de5c8d7e85ef801afdaf33d70f66c2ba5771 100644
--- a/README.md
+++ b/README.md
@@ -9,4 +9,21 @@
 [![PyPI](https://img.shields.io/pypi/wheel/firanka.svg)]()
 [![license](https://img.shields.io/github/license/mashape/apistatus.svg)]()
 
-Calculations on real functions
+Calculations on continuous, domain-being-a-single-interval, real domain
+functions.
+
+Functions of index: float -> any;
+
+
+
+You can:
+
+* join two series with a single operation
+* use a function on each value
+
+```python
+from firanka.series import DiscreteSeries, FunctionBasedSeries
+
+ds = DiscreteSeries(list of tuple(index, value), '<-20;4)')
+
+```
\ No newline at end of file
diff --git a/firanka/series/__init__.py b/firanka/series/__init__.py
index 63374d300700521420ba7dca1348c165c7c06385..433511507577aa789e26da50156158c3ee5a429d 100644
--- a/firanka/series/__init__.py
+++ b/firanka/series/__init__.py
@@ -3,15 +3,15 @@ from __future__ import print_function, absolute_import, division
 import six
 import logging
 
-from .exceptions import OutOfRangeError, EmptyDomainError
+from .exceptions import NotInDomainError, FirankaError
 from .range import Range, REAL_SET
 from .series import DiscreteSeries, FunctionBasedSeries, ModuloSeries
 
 
 __all__ = [
     'REAL_SET',
-    'OutOfRangeError',
-    'EmptyDomainError',
+    'FirankaError',
+    'NotInDomainError',
     'Range',
     'FunctionBasedSeries',
     'DiscreteSeries',
diff --git a/firanka/series/exceptions.py b/firanka/series/exceptions.py
index fd588a8578fbf08cd4464a4755eb0545d6bc57a5..1144a1df43d7c6c1870b7e49e18250a3c0b8850e 100644
--- a/firanka/series/exceptions.py
+++ b/firanka/series/exceptions.py
@@ -3,16 +3,12 @@ from __future__ import print_function, absolute_import, division
 import six
 import logging
 
-logger = logging.getLogger(__name__)
 
-
-class FirankaException(Exception):
-    pass
-
-
-class OutOfRangeError(FirankaException):
+class FirankaError(Exception):
     pass
 
 
-class EmptyDomainError(FirankaException):
-    pass
\ No newline at end of file
+class NotInDomainError(FirankaError, ValueError):
+    """
+    Requested index is beyond this domain
+    """
diff --git a/firanka/series/series.py b/firanka/series/series.py
index cad10115e2e5fa261e8983f3347122c97e263da3..c819859a162b033b987e1a1583b4d2bcac5b7cee 100644
--- a/firanka/series/series.py
+++ b/firanka/series/series.py
@@ -5,7 +5,8 @@ import six
 import functools
 import itertools
 
-from .range import Range, REAL_SET
+from .range import Range, REAL_SET, EMPTY_RANGE
+from .exceptions import NotInDomainError
 
 
 class Series(object):
@@ -16,16 +17,22 @@ class Series(object):
         self.domain = domain
 
     def __getitem__(self, item):
+        """
+        Return a value for given index, or a subslice of this series
+        :param item: a float, or a slice
+        :return: Series instance or a value
+        :raises NotInDomainError: index not in domain
+        """
         if isinstance(item, slice):
             item = Range(item)
             if item not in self.domain:
-                raise ValueError('slicing beyond series domain')
+                raise NotInDomainError('slicing beyond series domain')
 
             newdom = self.domain.intersection(item)
             return SlicedSeries(self, newdom)
         else:
             if item not in self.domain:
-                raise ValueError('item not in domain')
+                raise NotInDomainError('item not in domain')
 
             return self._get_for(item)
 
@@ -33,15 +40,67 @@ class Series(object):
         raise NotImplementedError
 
     def eval_points(self, points):
+        """
+        Return values for given points. A mass [] one could say
+        :param points: iterable of indices
+        :return: a list of values
+        """
         return [self[p] for p in points]
 
-    def apply(self, series, fun):
-        return AppliedSeries(self, series, fun)
+    def apply(self, fun):
+        """
+        Return this series with a function applied to each value
+        :param fun: callable/1 => 1
+        :return: Series instance
+        """
+        return AppliedSeries(self, fun)
+
+    def discretize(self, points, domain=None):
+        """
+        Return this as a DiscreteSeries, sampled at points
+        :return: a DiscreteSeries instance
+        """
+        if len(points) == 0:
+            return DiscreteSeries([])
+
+        if domain is None:
+            domain = Range(points[0], points[-1], True, True)
+
+        if domain not in self.domain:
+            raise NotInDomainError('points not inside this series!')
+
+        data = [(i, self[i]) for i in points]
+        return DiscreteSeries(data, domain)
+
+    def join(self, series, fun):
+        """
+        Return a new series with values of fun(v1, v2)
+
+        :param series: series to join against
+        :param fun: callable/2 => value
+        :return: new Series instance
+        """
+        return JoinedSeries(self, series, fun)
 
     def translate(self, x):
+        """
+        Translate the series by some distance
+        :param x: a float
+        :return: new Series instance
+        """
         return TranslatedSeries(self, x)
 
 
+class AppliedSeries(Series):
+    def __init__(self, series, applyfun):
+        super(AppliedSeries, self).__init__(series.domain)
+        self.fun = applyfun
+        self.series = series
+
+    def _get_for(self, item):
+        return self.fun(self._get_for(item))
+
+
 class TranslatedSeries(Series):
     def __init__(self, series, x):
         super(TranslatedSeries, self).__init__(self.domain.translate(x))
@@ -64,12 +123,17 @@ class SlicedSeries(Series):
 class DiscreteSeries(Series):
 
     def __init__(self, data, domain=None):
-        if domain is None:
+        if len(data) == 0:
+            domain = EMPTY_RANGE
+        elif domain is None:
             domain = Range(data[0][0], data[-1][0], True, True)
 
         self.data = data
         super(DiscreteSeries, self).__init__(domain)
 
+    def apply(self, fun):
+        return DiscreteSeries([(k, fun(v)) for k, v in self.data], self.domain)
+
     def _get_for(self, item):
         for k, v in reversed(self.data):
             if k <= item:
@@ -80,7 +144,7 @@ class DiscreteSeries(Series):
     def translate(self, x):
         return DiscreteSeries([(k+x, v) for k, v in self.data], self.domain.translate(x))
 
-    def apply(self, series, fun):
+    def join_discrete(self, series, fun):
         new_domain = self.domain.intersection(series.domain)
 
         if isinstance(series, DiscreteSeries):
@@ -156,16 +220,24 @@ class DiscreteSeries(Series):
 
 
 class FunctionBasedSeries(Series):
+    """
+    Series with values defined by a function
+    """
     def __init__(self, fun, domain):
         super(FunctionBasedSeries, self).__init__(domain)
         self.fun = fun
-        self._get_for = fun
 
+    def _get_for(self, item):
+        return self.fun(item)
 
-class AppliedSeries(Series):
+
+class JoinedSeries(Series):
+    """
+    Series stemming from performing an operation on two series
+    """
     def __init__(self, ser1, ser2, op):
-        super(AppliedSeries, self).__init__(
-            ser1.domain.intersection(ser2.domain))
+        domain = ser1.domain.intersection(ser2.domain)
+        super(JoinedSeries, self).__init__(domain)
         self.ser1 = ser1
         self.ser2 = ser2
         self.op = op
@@ -173,13 +245,18 @@ class AppliedSeries(Series):
     def _get_for(self, item):
         return self.op(self.ser1._get_for(item), self.ser2._get_for(item))
 
+
 class ModuloSeries(Series):
+
     def __init__(self, series):
         super(ModuloSeries, self).__init__(REAL_SET)
 
         self.series = series
         self.period = self.series.domain.length()
 
+        if self.period == 0:
+            raise ValueError('Modulo series cannot have a period of 0')
+
     def _get_for(self, item):
         if item < 0:
             item = -(item // self.period) * self.period + item
@@ -189,4 +266,3 @@ class ModuloSeries(Series):
             item = 0
 
         return self.series._get_for(self.series.domain.start + item)
-
diff --git a/tests/test_series/test_series.py b/tests/test_series/test_series.py
index 15a1517703857c7e7ad7ea206a26210210519c12..7036336a8e0f04e340e92efee39a978e55569479 100644
--- a/tests/test_series/test_series.py
+++ b/tests/test_series/test_series.py
@@ -48,26 +48,30 @@ class TestDiscreteSeries(unittest.TestCase):
         sa = DiscreteSeries([[0, 0], [1, 1], [2, 2]])
         sb = DiscreteSeries([[0, 1], [1, 2], [2, 3]])
 
-        sc = sa.apply(sb, lambda a, b: a+b)
+        sc = sa.join_discrete(sb, lambda a, b: a+b)
         self.assertIsInstance(sc, DiscreteSeries)
         self.assertEqual(sc.eval_points([0,1,2]), [1,3,5])
         self.assertEqual(sc.data, [(0,1),(1,3),(2,5)])
 
     def test_eval2(self):
         sa = DiscreteSeries([[0, 0], [1, 1], [2, 2]])
-        sb = FunctionBasedSeries(lambda x: x, '<0;2)')
+        sb = FunctionBasedSeries(lambda x: x, '<0;2>')
 
-        sc = sa.apply(sb, lambda a, b: a+b)
+        sc = sa.join_discrete(sb, lambda a, b: a+b)
         self.assertEqual(sc.eval_points([0,1,2]), [0,2,4])
 
         self.assertIsInstance(sc, DiscreteSeries)
         self.assertEqual(sc.data, [(0,0),(1,2),(2,4)])
 
+    def test_apply(self):
+        sa = DiscreteSeries([[0, 0], [1, 1], [2, 2]]).apply(lambda x: x+1)
+        self.assertEquals(sa.data, [(0,1),(1,2),(2,3)])
+
     def test_eval3(self):
         sa = FunctionBasedSeries(lambda x: x**2, '<-10;10)')
         sb = FunctionBasedSeries(lambda x: x, '<0;2)')
 
-        sc = sa.apply(sb, lambda a, b: a*b)
+        sc = sa.join(sb, lambda a, b: a*b)
 
         PTS = [0,1,1.9]
         EPTS = [x*x**2 for x in PTS]
@@ -75,6 +79,15 @@ class TestDiscreteSeries(unittest.TestCase):
         self.assertEqual(sc.eval_points(PTS), EPTS)
         self.assertTrue(Range('<0;2)') in sc.domain)
 
+    def test_discretize(self):
+        PTS = [0,1,2,3,4,5]
+        sa = FunctionBasedSeries(lambda x: x**2, '<-10;10)').discretize(PTS, '(-1;6)')
+        self.assertIsInstance(sa, DiscreteSeries)
+        self.assertEqual(sa.data, [(i, i**2) for i in PTS])
+
+        sa = FunctionBasedSeries(lambda x: x**2, '<-10;10)').discretize(PTS)
+        self.assertIsInstance(sa, DiscreteSeries)
+        self.assertEqual(sa.data, [(i, i**2) for i in PTS])
 
 class TestFunctionBasedSeries(unittest.TestCase):
     def test_slice(self):
@@ -89,6 +102,7 @@ class TestFunctionBasedSeries(unittest.TestCase):
         self.assertEqual(sp.domain.start, 0.5)
         self.assertEqual(sp.domain.stop, 1.5)
 
+
 class TestModuloSeries(unittest.TestCase):
     def test_base(self):
         series = ModuloSeries(DiscreteSeries([(0,1),(1,2),(2,3)], '<0;3)'))
@@ -102,5 +116,5 @@ class TestModuloSeries(unittest.TestCase):
         ser1 = ModuloSeries(FunctionBasedSeries(lambda x: x**2, '<0;3)'))
         ser2 = FunctionBasedSeries(lambda x: x, '<0;3)')
 
-        ser3 = ser1.apply(ser2, lambda x, y: x*y)
+        ser3 = ser1.join(ser2, lambda x, y: x*y)