#! /usr/bin/env python
#
def naca4_cambered ( m, p, t, c, xc ):

#*****************************************************************************80
#
## NACA4_CAMBERED: (xu,yu), (xl,yl) for a NACA cambered 4-digit airfoil.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    13 July 2018
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Eastman Jacobs, Kenneth Ward, Robert Pinkerton,
#    "The characteristics of 78 related airfoil sections from tests in
#    the variable-density wind tunnel",
#    NACA Report 460, 1933.
#
#  Parameters:
#
#    Input, real M, the maximum camber.
#    0 < M.
#
#    Input, real P, the location of maximum camber.
#    0.0 < P < 1.0
#
#    Input, real T, the maximum relative thickness.
#    0.0 < T <= 1.0
#
#    Input, real C, the chord length.
#    0.0 < C.
#
#    Input, real XC(*), points along the chord length.  
#    0.0 <= XC(*) <= C.
#
#    Output, real XU(*), YU(*), XL(*), YL(*), for each value of XC, measured
#    along the camber line, the corresponding values (XU,YU) on the upper
#    airfoil surface and (XL,YL) on the lower airfoil surface.
#
  import numpy as np

  i = np.nonzero ( 0.0 <= xc / c ) and np.nonzero ( xc / c <= p   )
  j = np.nonzero ( p   <= xc / c ) and np.nonzero ( xc / c <= 1.0 )
  k = np.nonzero ( xc / c < 0.0 ) or np.nonzero ( 1.0 < xc / c )

  n = len ( xc )

  divisor = np.zeros ( n )
  divisor[i] =         p  ** 2
  divisor[j] = ( 1.0 - p ) ** 2
  divisor[k] = 1.0

  dycdx = 2.0 * m * ( p - xc / c ) / divisor

  theta = np.arctan ( dycdx )
   
  yt = 5.0 * t * c * ( \
     0.2969 * np.sqrt ( xc / c ) \
     + (((( \
       - 0.1015 ) * ( xc / c ) \
       + 0.2843 ) * ( xc / c ) \
       - 0.3516 ) * ( xc / c ) \
       - 0.1260 ) * ( xc / c ) )

  yc = np.zeros ( n )
  yc[i] = m * ( xc[i]     ) * ( 2.0 * p - xc[i] / c       ) /         p   ** 2
  yc[j] = m * ( xc[j] - c ) * ( 2.0 * p - xc[j] / c - 1.0 ) / ( 1.0 - p ) ** 2
  yc[k] = 0.0

  xu = xc - yt * np.sin ( theta )
  yu = yc + yt * np.cos ( theta )
  xl = xc + yt * np.sin ( theta )
  yl = yc - yt * np.cos ( theta )

  return xu, yu, xl, yl

def naca4_cambered_test ( m, p, t ):

#*****************************************************************************80
#
## NACA4_CAMBERED_TEST tests NACA4_CAMBERED.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    13 July 2018
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, real M, the maximum camber.
#    0 < M.
#
#    Input, real P, the location of maximum camber.
#    0.0 < P < 1.0.
#
#    Input, real T, the maximum relative thickness.
#    0.0 < T <= 1.0.
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'NACA4_CAMBERED_TEST' )
  print ( '  NACA4_CAMBERED evaluates (xu,yu) and (xl,yl) for a NACA' )
  print ( '  cambered airfoil defined by a 4-digit code.' )

  c = 10.0
  n = 51

  xc = np.linspace ( 0.0, c, n )

  xu, yu, xl, yl = naca4_cambered ( m, p, t, c, xc )

  x2 = np.append ( xu, np.flip ( xl, 0 ) )
  y2 = np.append ( yu, np.flip ( yl, 0 ) )
#
#  Plot the wing surface.
#
  plt.plot ( x2, y2, 'b-', linewidth = 3 )
  plt.axis ( 'equal' )
  plt.grid ( True )
  plt.xlabel ( '<---X--->', fontsize = 16 )
  plt.ylabel ( '<---Y--->', fontsize = 16 )
  plt.title ( 'NACA 4-digit cambered airfoil', fontsize = 24 )

  filename = 'naca4_cambered_test.png'
  plt.savefig ( filename )
  print ( '' )
  print ( '  Graphics saved in file "%s"' % ( filename ) )
  plt.show ( )
#
#  Save data to a file.
#
  filename = 'naca4_cambered_test.txt'
  output = open ( filename, 'w' )
  for i in range ( 0, 2 * n ):
    s = '  %g  %g\n' % ( x2[i], y2[i] )
    output.write ( s )
  output.close ( )
  print ( '  Data saved in file "%s"' % ( filename ) )

  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  from naca4_mpt import naca4_mpt
  timestamp ( )
  code = 2412
  m, p, t = naca4_mpt ( code )
  naca4_cambered_test ( m, p, t )
  timestamp ( )