diff --git a/DS3231/ds3231_pb.py b/DS3231/ds3231_pb.py
index 52fd292..1be6aa2 100644
--- a/DS3231/ds3231_pb.py
+++ b/DS3231/ds3231_pb.py
@@ -42,6 +42,7 @@ class DS3231:
         else:
             raise ValueError('Side must be "X" or "Y"')
         self.ds3231 = pyb.I2C(bus, mode=pyb.I2C.MASTER, baudrate=400000)
+        self.timebuf = bytearray(7)
         if DS3231_I2C_ADDR not in self.ds3231.scan():
             raise DS3231Exception("DS3231 not found on I2C bus at %d" % DS3231_I2C_ADDR)
 
@@ -49,7 +50,7 @@ class DS3231:
         if set_rtc:
             data = self.await_transition()      # For accuracy set RTC immediately after a seconds transition
         else:
-            data = self.ds3231.mem_read(7, DS3231_I2C_ADDR, 0) # don't wait
+            data = self.ds3231.mem_read(self.timebuf, DS3231_I2C_ADDR, 0) # don't wait
         ss = bcd2dec(data[0])
         mm = bcd2dec(data[1])
         if data[2] & 0x40:
@@ -91,10 +92,10 @@ class DS3231:
         return rtc_ms - 1000 * t_ds3231
 
     def await_transition(self):                 # Wait until DS3231 seconds value changes
-        data = self.ds3231.mem_read(7, DS3231_I2C_ADDR, 0)
+        data = self.ds3231.mem_read(self.timebuf, DS3231_I2C_ADDR, 0)
         ss = data[0]
         while ss == data[0]:
-            data = self.ds3231.mem_read(7, DS3231_I2C_ADDR, 0)
+            data = self.ds3231.mem_read(self.timebuf, DS3231_I2C_ADDR, 0)
         return data
 
 # Get calibration factor for Pyboard RTC. Note that the DS3231 doesn't have millisecond resolution so we
diff --git a/DS3231/ds3231_test.py b/DS3231/ds3231_test.py
new file mode 100644
index 0000000..690eabb
--- /dev/null
+++ b/DS3231/ds3231_test.py
@@ -0,0 +1,29 @@
+# Program to test/demonstrate consistency of results from getcal()
+from array import array
+from ds2321_pb import DS3231
+
+# This takes 12.5 hours to run: to sve you the trouble here are results from one sample of Pyboard.
+
+# Mean and standard deviation of RTC correction factors based on ten runs over different periods.
+# Conclusion: for the best possible accuracy, run for 20 minutes. However a ten minute run gave
+# a result within 2ppm (one minute/yr).
+# >>> test()
+# t =  5 -174 -175 -175 -172 -175 -175 -172 -174 -175 -172  avg -173.8 sd   1.3
+# t = 10 -175 -175 -175 -175 -173 -175 -176 -175 -174 -175  avg -174.8 sd   0.7
+# t = 20 -175 -175 -175 -175 -174 -175 -175 -175 -175 -174  avg -174.8 sd   0.4
+# t = 40 -175 -175 -175 -174 -174 -175 -174 -174 -175 -174  avg -174.4 sd   0.5
+
+def test():
+    NSAMPLES = 10
+    a = DS3231()
+    for t in (5, 10, 20, 40):
+        values = array('f', (0 for z in range(NSAMPLES)))
+        print('t = {:2d}'.format(t), end = '')
+        for x in range(NSAMPLES):
+            cal = a.getcal(t)
+            values[x] = cal
+            print('{:5d}'.format(cal), end = '')
+        avg = sum(values)/NSAMPLES
+        sd2 = sum([(v -avg)**2 for v in values])/NSAMPLES
+        sd = sd2 ** 0.5
+        print('  avg {:5.1f} sd {:5.1f}'.format(avg, sd))