Merge branch 'master' into master-ninfinity

# Conflicts:
#	CMakeLists.txt
#	README.md

CMakeLists.txt

@@ -8,37 +8,36 @@ set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} \
 -finit-real=inf \
 -ffpe-trap=invalid,zero \
 -fdefault-real-8")
-set(SRC_DIR "src/fortran")
 
 add_library(xfoil SHARED
         # Include modules
-        ${SRC_DIR}/i_pindex.f90
-        ${SRC_DIR}/i_xfoil.f90
-        ${SRC_DIR}/i_blpar.f90
-        ${SRC_DIR}/i_circle.f90
-        ${SRC_DIR}/i_xbl.f90
+        src/i_pindex.f90
+        src/i_xfoil.f90
+        src/i_blpar.f90
+        src/i_circle.f90
+        src/i_xbl.f90
         # Shared modules to avoid circular dependencies
-        ${SRC_DIR}/s_xbl.f90
-        ${SRC_DIR}/s_xoper.f90
-        ${SRC_DIR}/s_xfoil.f90
+        src/s_xbl.f90
+        src/s_xoper.f90
+        src/s_xfoil.f90
         # Main modules
-        ${SRC_DIR}/m_aread.f90
-        ${SRC_DIR}/m_iopol.f90
-        ${SRC_DIR}/m_naca.f90
-        ${SRC_DIR}/m_sort.f90
-        ${SRC_DIR}/m_spline.f90
-        ${SRC_DIR}/m_userio.f90
-        ${SRC_DIR}/m_xbl.f90
-        ${SRC_DIR}/m_xblsys.f90
-        ${SRC_DIR}/m_xfoil.f90
-        ${SRC_DIR}/m_xgdes.f90
-        ${SRC_DIR}/m_xgeom.f90
-        ${SRC_DIR}/m_xmdes.f90
-        ${SRC_DIR}/m_xoper.f90
-        ${SRC_DIR}/m_xpanel.f90
-        ${SRC_DIR}/m_xpol.f90
-        ${SRC_DIR}/m_xqdes.f90
-        ${SRC_DIR}/m_xsolve.f90
-        ${SRC_DIR}/m_xutils.f90
+        src/m_aread.f90
+        src/m_iopol.f90
+        src/m_naca.f90
+        src/m_sort.f90
+        src/m_spline.f90
+        src/m_userio.f90
+        src/m_xbl.f90
+        src/m_xblsys.f90
+        src/m_xfoil.f90
+        src/m_xgdes.f90
+        src/m_xgeom.f90
+        src/m_xmdes.f90
+        src/m_xoper.f90
+        src/m_xpanel.f90
+        src/m_xpol.f90
+        src/m_xqdes.f90
+        src/m_xsolve.f90
+        src/m_xutils.f90
         # API
-        ${SRC_DIR}/api.f90)
+        src/api.f90)

README.md

@@ -112,4 +112,4 @@ to analyze for an angle of attack of 10 degrees, a lift coefficient of 1.0, and
 -0.5 to 0.5 with steps of 0.05.
 
 For other features and specifics, see the documentation in the Python source files.
-
+

setup.py

@@ -17,14 +17,18 @@
 #   along with XFoil.  If not, see <https://www.gnu.org/licenses/>.
 import os
 import platform
+import re
 import subprocess
 import sys
 
 from setuptools import setup
 from setuptools.extension import Extension
 from setuptools.command.build_ext import build_ext
 
-version = '0.0.16'
+__version__ = re.findall(
+    r"""__version__ = ["']+([0-9\.]*)["']+""",
+    open('xfoil/__init__.py').read(),
+)[0]
 
 options = {k: 'OFF' for k in ['--opt', '--debug', '--cuda']}
 for flag in options.keys():
@@ -115,9 +119,10 @@ def readme():
 
 setup(
     name='xfoil',
-    version=version,
+    version=__version__,
     description='Stripped down version of XFOIL as compiled python module ',
     long_description=readme(),
+    long_description_content_type='text/markdown',
     classifiers=[
         'Development Status :: 3 - Alpha',
         'Intended Audience :: Science/Research',
@@ -132,12 +137,12 @@ def readme():
     ],
     keywords='xfoil airfoil aerodynamic analysis',
     url='https://github.com/daniel-de-vries/xfoil-python',
-    download_url='https://github.com/daniel-de-vries/xfoil-python/tarball/' + version,
+    download_url='https://github.com/daniel-de-vries/xfoil-python/tarball/' + __version__,
     author='Daniël de Vries',
     author_email='[email protected]',
     license='GNU General Public License v3 or later (GPLv3+)',
     packages=['xfoil'],
-    package_dir={'': 'src'},
+    # package_dir={'': 'src'},
     ext_modules=[CMakeExtension('xfoil.xfoil')],
     cmdclass={'build_ext': CMakeBuild},
     install_requires=['numpy'],

src/fortran/api.f90 → src/api.f90

@@ -266,13 +266,9 @@ end function get_max_iter
 
     subroutine reset_bls() bind(c, name='reset_bls')
         use i_xfoil, only: LBLini, LIPan, show_output
-        LBLini = .not.LBLini
-        if (LBLini) then
-            if (show_output) write (*, *) 'BLs are assumed to be initialized'
-        else
-            if (show_output) write (*, *) 'BLs will be initialized on next point'
-            LIPan = .false.
-        endif
+        LBLini = .false.
+        LIPan = .false.
+        if (show_output) write (*, *) 'BLs will be initialized on next point'
     end subroutine reset_bls
 
     subroutine repanel(n_panel, cv_par, cte_ratio, ctr_ratio, &
@@ -483,4 +479,35 @@ subroutine cseq(cl_start, cl_end, n_step, &
         end do
     end subroutine cseq
 
+
+    function get_n_cp() bind(c, name='get_n_cp')
+        use i_xfoil, only: N
+        integer(c_int) :: get_n_cp
+        get_n_cp = N
+    end function get_n_cp
+
+    subroutine get_cp(x_out, cp_out, n_points) bind(c, name='get_cp')
+        use s_xfoil, only: comset
+        use m_xoper
+        use i_xfoil
+        implicit none
+
+        integer(c_int), intent(in) :: n_points
+        real(c_float), dimension(n_points), intent(inout) :: x_out, cp_out
+
+        real :: beta, bfac, cpcom, cpinc, den
+        integer :: i
+
+        call comset
+        beta = sqrt(1.0 - MINf**2)
+        bfac = 0.5 * MINf**2 / (1.0 + beta)
+
+        do i = 1, N
+            cpinc = 1.0 - (GAM(i) / QINf)**2
+            den = beta + bfac * cpinc
+            cp_out(i) = cpinc / den
+            x_out(i) = X(i)
+        enddo
+    end subroutine get_cp
+
 end module api

src/fortran/i_xfoil.f90 → src/i_xfoil.f90

@@ -26,8 +26,8 @@ module i_xfoil
     !*** Start of declarations rewritten by SPAG
     !
     ! PARAMETER definitions
-    !
-    integer, parameter :: IQX = 370, IPX = 5, ISX = 2, IWX = IQX / 8 + 2, IBX = 4 * IQX, IZX = IQX + IWX, &
+    ! IQX = 370 originally
+    integer, parameter :: IQX = 640, IPX = 5, ISX = 2, IWX = IQX / 8 + 2, IBX = 4 * IQX, IZX = IQX + IWX, &
             & IVX = IQX / 2 + IWX + 50, NAX = 800, NPX = 12, NFX = 128, NTX = 2 * IBX
     !
     ! Local variables

src/fortran/m_sort.f90 → src/m_sort.f90

@@ -1,4 +1,4 @@
-!***********************************************************************
+! ***********************************************************************
 !   Copyright (c) 2019 D. de Vries
 !   Original Copyright (c) 2000 Mark Drela
 !

src/fortran/m_xoper.f90 → src/m_xoper.f90

@@ -2514,10 +2514,13 @@ function viscal(Niter1)
         !
         ! Local variables
         !
-        real :: cdpdif, ddef, delp, dpds, duds, edef, fnum, hfrac, hk, hki, hkm, hkmax, patt, pdef, pdefi, &
-                & pdefm, psep
+        real :: cdpdif, delp, hfrac, hki, hkm, hkmax, patt, pdefi, pdefm, psep
         real, save :: eps1
-        integer :: i, ibl, iblm, iblp, idec, ione, is, iten, iter, izero, lu, niter
+        integer :: i, ibl, is, iter, niter
+
+! These are used by the weird writing to a file code
+!        real :: ddef, dpds, duds, edef, fnum, hk, pdef
+!        integer :: i, iblm, iblp, idec, ione, iten, izero, lu
         !
         !*** End of declarations rewritten by SPAG
         !
@@ -2698,35 +2701,36 @@ function viscal(Niter1)
         if (show_output) write (*, 99004) ACRit(is), hkmax, CD, 2.0 * psep, 2.0 * patt, 2.0 * delp, XOCtr(is)
         99004 format (1x, f10.3, f10.4, f11.6, 3F11.6, f10.4, '     #')
 
-        izero = ichar('0')
-
-        !c      fnum = acrit(is)
-        fnum = XSTrip(is) * 100.0
-
-        iten = int(fnum / 9.99999)
-        ione = int((fnum - float(10 * iten)) / 0.99999)
-        idec = int((fnum - float(10 * iten) - float(ione)) / 0.09999)
-
-        FNAme = char(iten + izero) // char(ione + izero) // char(idec + izero) // '.bl'
-        lu = 44
-        open (lu, file = FNAme, status = 'unknown')
-        rewind (lu)
-        write (lu, '(a,a)') '#       s         ue          H          P         K ', '        x    -m du/dx'
-        !       1234567890 1234567890 1234567890 1234567890 1234567890 1234567890
-        do ibl = 2, IBLte(is)
-            iblm = max(ibl - 1, 2)
-            iblp = min(ibl + 1, IBLte(is))
-            i = IPAn(ibl, is)
-            hk = DSTr(ibl, is) / THEt(ibl, is)
-            ddef = DSTr(ibl, is) * UEDg(ibl, is)
-            pdef = THEt(ibl, is) * UEDg(ibl, is)**2
-            edef = TSTr(ibl, is) * UEDg(ibl, is)**3 * 0.5
-            duds = (UEDg(iblp, is) - UEDg(iblm, is)) / (XSSi(iblp, is) - XSSi(iblm, is))
-            dpds = -ddef * duds
-            write (lu, 99005) XSSi(ibl, is), UEDg(ibl, is), hk, pdef, edef, X(i), dpds
-            99005 format (1x, 3F11.4, 2F11.6, f11.3, e14.6)
-        enddo
-        close (lu)
+! This stuff seems to only write stuff to a file, which seems unecessary to do each time viscal is called...
+!        izero = ichar('0')
+!
+!        !c      fnum = acrit(is)
+!        fnum = XSTrip(is) * 100.0
+!
+!        iten = int(fnum / 9.99999)
+!        ione = int((fnum - float(10 * iten)) / 0.99999)
+!        idec = int((fnum - float(10 * iten) - float(ione)) / 0.09999)
+!
+!        FNAme = char(iten + izero) // char(ione + izero) // char(idec + izero) // '.bl'
+!        lu = 44
+!        open (lu, file = FNAme, status = 'unknown')
+!        rewind (lu)
+!        write (lu, '(a,a)') '#       s         ue          H          P         K ', '        x    -m du/dx'
+!        !       1234567890 1234567890 1234567890 1234567890 1234567890 1234567890
+!        do ibl = 2, IBLte(is)
+!            iblm = max(ibl - 1, 2)
+!            iblp = min(ibl + 1, IBLte(is))
+!            i = IPAn(ibl, is)
+!            hk = DSTr(ibl, is) / THEt(ibl, is)
+!            ddef = DSTr(ibl, is) * UEDg(ibl, is)
+!            pdef = THEt(ibl, is) * UEDg(ibl, is)**2
+!            edef = TSTr(ibl, is) * UEDg(ibl, is)**3 * 0.5
+!            duds = (UEDg(iblp, is) - UEDg(iblm, is)) / (XSSi(iblp, is) - XSSi(iblm, is))
+!            dpds = -ddef * duds
+!            write (lu, 99005) XSSi(ibl, is), UEDg(ibl, is), hk, pdef, edef, X(i), dpds
+!            99005 format (1x, 3F11.4, 2F11.6, f11.3, e14.6)
+!        enddo
+!        close (lu)
 
     end function viscal
     !*==DCPOUT.f90  processed by SPAG 7.21DC at 11:25 on 11 Jan 2019

src/fortran/m_xpanel.f90 → src/m_xpanel.f90

@@ -1350,6 +1350,9 @@ subroutine qdcalc
         do i = N + 1, N + NW
             !
             iw = i - N
+            if (iw > iwx .or. iw <= 0) then
+                exit
+            end if
             !
             !------ airfoil contribution at wake panel node
             call psilin(i, X(i), Y(i), NX(i), NY(i), psi, psi_n, .false., .true.)

src/xfoil/__init__.py → xfoil/__init__.py

@@ -16,3 +16,5 @@
 #   You should have received a copy of the GNU General Public License
 #   along with XFoil.  If not, see <https://www.gnu.org/licenses/>.
 from .xfoil import XFoil
+
+__version__ = "1.1.1"

src/xfoil/model.py → xfoil/model.py

@@ -47,7 +47,7 @@ def x(self):
     @x.setter
     def x(self, value):
         v = value.flatten()
-        self.coords.resize((v.size, 2))
+        self.coords = np.resize(self.coords, (v.size, 2))
         self.coords[:, 0] = v[:]
 
     @property
@@ -58,5 +58,5 @@ def y(self):
     @y.setter
     def y(self, value):
         v = value.flatten()
-        self.coords.resize((v.size, 2))
+        self.coords = np.resize(self.coords, (v.size, 2))
         self.coords[:, 1] = v[:]

src/xfoil/xfoil.py → xfoil/xfoil.py

@@ -313,3 +313,20 @@ def cseq(self, cl_start, cl_end, cl_step):
         cm[isnan] = np.nan
         cp[isnan] = np.nan
         return a.astype(float), cl.astype(float), cd.astype(float), cm.astype(float), cp.astype(float)
+
+    def get_cp_distribution(self):
+        """Get the Cp distribution from the last converged point.
+
+        Returns
+        -------
+        x : np.array
+            X-coordinates
+        cp : np.ndarray
+            Pressure coefficients at the corresponding x-coordinates
+        """
+        n = self._lib.get_n_cp()
+        x = np.zeros(n, dtype=c_float)
+        cp = np.zeros(n, dtype=c_float)
+
+        self._lib.get_cp(x.ctypes.data_as(fptr), cp.ctypes.data_as(fptr), byref(c_int(n)))
+        return x.astype(float), cp.astype(float)