diff --git a/matlab/binData.m b/matlab/binData.m
new file mode 100644
index 00000000..5d789eda
--- /dev/null
+++ b/matlab/binData.m
@@ -0,0 +1,21 @@
+function [g6B,bin] = binData(g6C,binMax,bins)
+%Bin the data returned from g6_struct into equi-separated partitions
+% g6C:       Matrix containing the g6 values output from g6_struct
+% binMax:    The maximum binning value
+% bins:      The number of bins to take
+%Returns
+% g6B:       The binned g6 data
+% bin:       The bin values of g6B
+
+bin = linspace(0,binMax,bins);
+g6B = zeros(length(bin)-1,1);
+for kk=1:(length(bin)-1)
+    idx = find( g6C(1,:) >= bin(kk) & g6C(1,:) < bin(kk+1) );
+    if length(idx) >0
+        g6B(kk) = sum(g6C(2,idx))./length(idx); %Average the values in a bin
+    else
+        g6B(kk) = 0;
+    end
+end
+
+end
\ No newline at end of file
diff --git a/matlab/findNN.m b/matlab/findNN.m
new file mode 100644
index 00000000..86351847
--- /dev/null
+++ b/matlab/findNN.m
@@ -0,0 +1,8 @@
+function [neighbours,locations] = findNN(pos,X,Y,radius)
+% Returns values and indices of nearest neighbours to location pos in
+% positions X,Y within radius.
+
+    r = sqrt(bsxfun(@minus,pos(1),X).^2 +  bsxfun(@minus,pos(2),Y).^2);
+    locations = find(r < radius & r~=0);
+    neighbours = [X(locations) Y(locations)];
+end
\ No newline at end of file
diff --git a/matlab/g6_struct.m b/matlab/g6_struct.m
new file mode 100644
index 00000000..4fb0da2d
--- /dev/null
+++ b/matlab/g6_struct.m
@@ -0,0 +1,31 @@
+function [g6C] = g6_struct(X,Y,rad)
+%Determine the orientational correlation function between every pairing of
+%points, and sorts the results based on distance between points. This gives
+%g6(r).
+% X,Y:   Vectors of x,y values for points.
+% rad:   Radius in which to examine for neighbouring points.
+%Return:
+% g6C:   Matrix of g6 values and parameters
+
+%Determine all orientational order values for points
+psi6p = zeros(length(X),1);
+for ii=1:length(X)
+    psi6p(ii) = psi6([X(ii),Y(ii)],X,Y,rad);
+end
+
+%Calculate all n choose k pairings of points, and give distance between
+%them, values for the respective orientational orders, and g6 values. Sort
+%the values based on separated distance
+S = uniqPairIdx_precalc([X,Y],psi6p);
+[~,order] = sort([S(:).rabs],'ascend');
+SSorted = S(order);
+
+%Create output structure for next stage of calculation
+for ii=1:size(SSorted,2)
+    g6S(ii).r = SSorted(ii).rabs;
+    psi6_1 = SSorted(ii).cor0;
+    psi6_2 = SSorted(ii).cor1;
+    g6S(ii).val = abs(conj(psi6_1)*conj(psi6_1));
+end
+g6C = cell2mat(squeeze(struct2cell(g6S)));
+end
\ No newline at end of file
diff --git a/matlab/psi6.m b/matlab/psi6.m
new file mode 100644
index 00000000..049b8a2e
--- /dev/null
+++ b/matlab/psi6.m
@@ -0,0 +1,12 @@
+function [psi6_pos,nn] = psi6(pos,X,Y,radius)
+    [nn,idx] = findNN(pos,X,Y,radius);
+    psi6_pos = 0;
+    
+    if size(nn,1) > 0
+        for ii=1:size(nn,1)
+            psi6_pos = psi6_pos + exp(6*1i*getAngle2(pos, [X(idx(ii)), Y(idx(ii))] ));
+        end
+        psi6_pos = psi6_pos./length(nn);
+    end
+    nn=length(nn);
+end
diff --git a/matlab/uniqPairIdx_precalc.m b/matlab/uniqPairIdx_precalc.m
new file mode 100644
index 00000000..e3ff7f92
--- /dev/null
+++ b/matlab/uniqPairIdx_precalc.m
@@ -0,0 +1,19 @@
+function [S] = uniqPairIdx_precalc(R,psi6p)
+%Calculate all unique pairings of the orientational correlations, and
+%create a struct to return that holds the distance between paired elements,
+%and the respective orientational order values.
+% R:        Vectors of R=[X Y] values for points.
+% psi6p:    Orientational order values defined over the range of points
+%Return:
+% S:        Struct of orientational corder values cor0,cor1 and distance
+%            between elements
+    count = 0;
+    for ii=1:(size(R,1)-1)
+        for jj=ii+1:size(R,1)
+            count = count +1;
+            S(count).cor0 = psi6p(ii);
+            S(count).cor1 = psi6p(jj);
+            S(count).rabs = sqrt( sum( (R(ii,:) - R(jj,:) ).^2 ) );
+        end
+    end
+end
\ No newline at end of file