diff --git a/+af_tools/nacaairfoil.m b/+af_tools/nacaairfoil.m
index 80791fb2f8ed73275e810e5f9139d856ca0077d1..674857c20a8a42a01d289e37be608b3aed60cdbc 100644
--- a/+af_tools/nacaairfoil.m
+++ b/+af_tools/nacaairfoil.m
@@ -40,12 +40,18 @@ function [x, y] = nacaairfoil(varargin)
 %     true   - (default) zero thickness trailing edge
 %     false  - finite trailing edge
 %
+%   [...] = NACAAIRFOIL(..., 'savedat', false) saves the output coordinates in a
+%   dat files. savedat can either be true or false (default).
+%
 % Inputs:
 %   digits  : Character vector representing the NACA airfoil (ex: '0012', '24120', '24012')
 %   nPoints : Number of output coordinates
 %
 % Output:
 %   [X, Y] : Airfoil coordinates, with 0 <= X <= 1.
+%            The output coordinates follow the most usual order. They start at
+%            the trailing edge, along the upper surface to the leading edge and
+%            back around the lower surface to trailing edge.
 %
 % Example:
 %   NACAAIRFOIL
@@ -69,11 +75,11 @@ narginchk(0,6);
 import af_tools.nacacamber;
 import af_tools.utils.parsenacainputs;
 
-OPTION_LIST = {'spacing', 'zerote'};
+OPTION_LIST = {'spacing', 'zerote', 'savedat'};
 
 % Extract and validate the inputs
 [digits, nPoints, idxOpts] = parsenacainputs(OPTION_LIST, varargin{:});
-[spacing, zerote] = parseoptioninputs(varargin{idxOpts:end});
+[spacing, zerote, savedat] = parseoptioninputs(varargin{idxOpts:end});
 
 %------------------------------------
 % Coordinates calculation
@@ -82,6 +88,10 @@ OPTION_LIST = {'spacing', 'zerote'};
 % number of points on the upper and lower surfaces.
 [xc, yc, gradY] = nacacamber(digits, ceil((nPoints+1)/2), 'spacing', spacing);
 
+% Flip camberline, so the first element is the trailing edge
+xc = fliplr(xc);
+yc = fliplr(yc);
+
 % Thickness value (checks were done in nacacamber)
 t = str2double(digits(end-1:end)) /100;
 
@@ -110,11 +120,16 @@ yl = yc - yt.*cos(theta);
 x = [xu, fliplr(xl(1:end-1))];
 y = [yu, fliplr(yl(1:end-1))];
 
+% Save coordinates to dat file
+if savedat
+    savetodat(digits, x, y, spacing);
+end
+
 end
 
 
 % ------------------------------------------------------------------------------
-function [spacing, zerote] = parseoptioninputs(varargin)
+function [spacing, zerote, savedat] = parseoptioninputs(varargin)
 % PARSEINPUTS Parses the input and checks their validity
 
 import af_tools.utils.*
@@ -122,6 +137,7 @@ import af_tools.utils.*
 % Constants and defaults
 DEFAULT_SPACING = 'halfcosine';
 DEFAULT_ZEROTE = true;
+DEFAULT_SAVEDAT = false;
 ALLOWED_SPACING = {'linear', 'cosine', 'halfcosine'};
 
 % Option validator
@@ -132,10 +148,21 @@ validSpacing = @(x) any(validatestring(x, ALLOWED_SPACING, mfilename()));
 p = inputParser;
 addParameter(p,'spacing',DEFAULT_SPACING, validSpacing);
 addParameter(p,'zerote',DEFAULT_ZEROTE, validLogical);
+addParameter(p,'savedat',DEFAULT_SAVEDAT, validLogical);
 parse(p,varargin{:});
 spacing = p.Results.spacing;
 zerote = p.Results.zerote;
+savedat = p.Results.savedat;
+
+end
 
+function savetodat(digits, x, y, spacing)
+% SAVETODAT Saves the airfoil coordinates to a dat file
 
+header = ['NACA',digits,' (',spacing,' spacing)'];
+data=num2str([x; y]');
+contents=char(header,data);
+fileName=['naca', digits, '.dat'];
+dlmwrite(fileName,contents,'delimiter','')
 
 end
diff --git a/.gitignore b/.gitignore
index 962e76b08d3e754b7cf9fb49b0825b0e100bffa8..b085236b0fffa791f30cbd87a5e8b88d7979c50e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@
 *.m~
 *.mat
 *.png
+*.dat
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 05ad16669418829589d95a0addd94b2ce4d0d797..ba71091c4a07751f81b0e6115ca13136448583a4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,8 +9,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added save to dat file for _nacaairfoil_
+
 ### Changed
 
+- Coordinates of nacaairfoil start at TE, go along upper side to LE, then TE
+
 ### Deprecated
 
 ### Removed
diff --git a/README.md b/README.md
index bb451a61f928cf79833ae315cdfa0e957922c88a..feda6b5e7a55033ab87f3401edbddc6ab3359893 100644
--- a/README.md
+++ b/README.md
@@ -35,7 +35,7 @@ Unfortunately, Matlab packages have a few limitations. Hence, you should **avoid
 changing the package name** if you do not know what you are doing.
 
 - The `+` in the folder's name is required for Matlab to treat it as a package.
-- Some functions of this package are interconnected (_e.g._, _nacairfoil_ calls
+- Some functions of this package are interconnected (_e.g._, _nacaairfoil_ calls
   _nacacamber_). In that situation, Matlab forces the main function to re-import
   explicitly all sub-functions before using them, even if they are part of the
   same package ([people have complained][matlabPackgeIssue] about that design
@@ -146,6 +146,7 @@ The following functions use the same two main input arguments.
 #### nacacamber
 
 Determines the camber line of a NACA 4 or 5 digits airfoil.
+The output is ordered from the leading edge to the trailing edge.
 
 ```matlab
 import af_tools.nacacamber
@@ -162,6 +163,9 @@ import af_tools.nacacamber
 #### nacaairfoil
 
 Generates the full coordinates of a NACA 4 or 5 digits airfoil.
+The output coordinates follow the most usual order. They start at the trailing
+edge, along the upper surface to the leading edge and back around the lower
+surface to trailing edge.
 
 ```matlab
 import af_tools.nacaairfoil
@@ -169,12 +173,14 @@ import af_tools.nacaairfoil
 [x, y] = nacaairfoil(digits)
 [x, y] = nacaairfoil(digits, nPoints)
 [x, y] = nacaairfoil(digits, nPoints, 'spacing', spacing, 'zerote', true)
+[x, y] = nacaairfoil(digits, nPoints, 'savedat', true)
 ```
 
 | Options   | Example                            | Default        |
 |---------- | ---------------------------------- | -------------- |
 | 'spacing' | `halfcosine`, `cosine`, `linear`   | `'halfcosine'` |
-| 'zerote'  | `true`, `false`                    | `true`        |
+| 'zerote'  | `true`, `false`                    | `true`         |
+| 'savedat' | `true`, `false`                    | `true`         |
 
 ### Polar manipulation
 
diff --git a/tests/test_nacaairfoil.m b/tests/test_nacaairfoil.m
index 93ead07c7316c3ec88c3af8942e3790186f8de68..def97d107228f4cc576f06e6164a48df5634cbf7 100644
--- a/tests/test_nacaairfoil.m
+++ b/tests/test_nacaairfoil.m
@@ -108,7 +108,7 @@ function test_invalidOptions(testCase)
 % Error if invalid option name or missing parameter value
 
 wrongName = 'testspacing';
- 
+
 verifyError(testCase, @() af_tools.nacaairfoil('0012', 100, wrongName), 'MATLAB:InputParser:ParamMissingValue')
 verifyError(testCase, @() af_tools.nacaairfoil('0012', 100, wrongName, false), 'MATLAB:InputParser:UnmatchedParameter')
 verifyError(testCase, @() af_tools.nacaairfoil('0012', 100, 'spacing', 'cosine', wrongName), 'MATLAB:InputParser:ParamMissingValue')
@@ -191,12 +191,12 @@ cosine = 1-cos(betacos);
 [xc_half, ~] = af_tools.nacaairfoil('0012', nPoints, 'spacing', 'halfcosine');
 
 % Test the two halves of the output vector (i.e. upper and lower surfaces)
-verifyEqual(testCase, xc_lin(1:(end+1)/2), linear)
-verifyEqual(testCase, xc_cos(1:(end+1)/2), cosine)
-verifyEqual(testCase, xc_half(1:(end+1)/2), halfcos)
-verifyEqual(testCase, xc_lin(end:-1:(end+1)/2), linear)
-verifyEqual(testCase, xc_cos(end:-1:(end+1)/2), cosine)
-verifyEqual(testCase, xc_half(end:-1:(end+1)/2), halfcos)
+verifyEqual(testCase, xc_lin(1:(end+1)/2), fliplr(linear))
+verifyEqual(testCase, xc_cos(1:(end+1)/2), fliplr(cosine))
+verifyEqual(testCase, xc_half(1:(end+1)/2), fliplr(halfcos))
+verifyEqual(testCase, xc_lin(end:-1:(end+1)/2), fliplr(linear))
+verifyEqual(testCase, xc_cos(end:-1:(end+1)/2), fliplr(cosine))
+verifyEqual(testCase, xc_half(end:-1:(end+1)/2), fliplr(halfcos))
 
 end