取一定飞机翼的参数返回系列点的坐标绘出其外形

!********************************************
!********************************************
!取一定飞机翼的参数返回系列点的坐标绘出其外形
!********************************************
! SINOMACH
!
! 取一定飞机翼的参数返回系列点的坐标绘出其外形
!
! 周勇 20221125
!*********************************************
!c$$Run Output:
!*********************************************
!*********************************************
Program main

!*********************************************************************72
!
!c MAIN is the main program for NACA_PRB.
!
! Discussion:
!
! NACA_PRB tests the NACA library.
!
Implicit None

Call timestamp()
Write (, '(a)') ''
Write (, '(a)') 'NACA_PRB:'
Write (, '(a)') ' FORTRAN95,2013&18 version'
Write (, '(a)') ' Test the NACA library.'

Call test01()
Call test02()
!
! Terminate.
!
Write (, '(a)') ''
Write (, '(a)') 'NACA_PRB:'
Write (, '(a)') ' Normal end of execution.'
Write (, '(a)') ''
Call timestamp()

Return
End Program main
Subroutine test01()

!*********************************************************************72
!
!c TEST01 tests NACA4_SYMMETRIC.
!
Implicit None

Integer n
Parameter (n=51)

Double Precision c
Character *(255) command_filename
Integer command_unit
Character (255) data_filename
Integer i
Double Precision r8vec_max
Double Precision r8vec_min
Double Precision ratio
Double Precision t
Double Precision x(n)
Double Precision x_max
Double Precision x_min
Double Precision xy(2, 2n)
Double Precision y(n)
Double Precision y_max
Double Precision y_min

Write (, '(a)') ''
Write (, '(a)') 'TEST01'
Write (, '(a)') ' NACA4_SYMMETRIC evaluates y(x) for a NACA'
Write (, '(a)') ' symmetric airfoil defined by a 4-digit code.'

c = 10.0D+00
t = 0.15D+00
Call r8vec_linspace(n, 0.0D+00, c, x)
Call naca4_symmetric(t, c, n, x, y)
!
! Reorganize data into a single object.
!
Do i = 1, n
xy(1, i) = x(i)
xy(2, i) = -y(i)
End Do

Do i = 1, n
xy(1, n+i) = x(n+1-i)
xy(2, n+i) = y(n+1-i)
End Do
!
! Determine size ratio.
!
x_min = r8vec_min(n, x)
x_max = r8vec_max(n, x)
y_max = r8vec_max(n, y)
y_min = -y_max
ratio = (y_max-y_min)/(x_max-x_min)
!
! Save data to a file.
!
data_filename = 'symmetric_data.txt'
Call r8mat_write(data_filename, 2, 2n, xy)
Write (, '(a)') ' Data saved in file "' // trim(data_filename) // '"'
!
! Create the command file.
!
command_filename = 'symmetric_commands.txt'
Call get_unit(command_unit)
Open (Unit=command_unit, File=command_filename, Status='replace')
Write (command_unit, '(a)') 'set term png'
Write (command_unit, '(a)') 'set grid'
Write (command_unit, '(a)') 'set style data lines'
Write (command_unit, '(a,g14.6)') 'set size ratio ', ratio
Write (command_unit, '(a)') 'set timestamp'
Write (command_unit, '(a)') 'unset key'
Write (command_unit, '(a)') 'set output "symmetric.png"'
Write (command_unit, '(a)') 'set xlabel "<---X--->"'
Write (command_unit, '(a)') 'set ylabel "<---Y--->"'
Write (command_unit, '(a)') 'set title "NACA Symmetric Airfoil"'
Write (command_unit, '(a)') 'plot "' // trim(data_filename) // '" using 1:2 with lines lw 3'
Write (command_unit, '(a)') 'quit'

Close (Unit=command_unit)

Write (*, '(a)') ' Created command file "' // trim(command_filename) // '".'

Return
End Subroutine test01
Subroutine test02()

!*********************************************************************72
!
!c TEST02 tests NACA4_CAMBERED.
!
! 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.
!
Implicit None

Integer n
Parameter (n=51)

Double Precision c
Character *(255) command_filename
Integer command_unit
Character (255) data_filename
Integer i
Double Precision m
Double Precision p
Double Precision r8vec_max
Double Precision r8vec_min
Double Precision ratio
Double Precision t
Double Precision x_max
Double Precision x_min
Double Precision xc(n)
Double Precision xl(n)
Double Precision xu(n)
Double Precision xy(2, 2n)
Double Precision y_max
Double Precision y_min
Double Precision yl(n)
Double Precision yu(n)

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

m = 0.02D+00
p = 0.4D+00
t = 0.12D+00
c = 10.0D+00

Call r8vec_linspace(n, 0.0D+00, c, xc)
Call naca4_cambered(m, p, t, c, n, xc, xu, yu, xl, yl)
!
! Reorganize data into a single object.
!
Do i = 1, n
xy(1, i) = xl(i)
xy(2, i) = yl(i)
End Do

Do i = 1, n
xy(1, n+i) = xu(n+1-i)
xy(2, n+i) = yu(n+1-i)
End Do
!
! Determine size ratio.
!
x_min = min(r8vec_min(n,xl), r8vec_min(n,xu))
x_max = max(r8vec_max(n,xl), r8vec_max(n,xu))
y_min = min(r8vec_min(n,yl), r8vec_min(n,yu))
y_max = max(r8vec_max(n,yl), r8vec_max(n,yu))
ratio = (y_max-y_min)/(x_max-x_min)
!
! Save data to a file.
!
data_filename = 'cambered_data.txt'
Call r8mat_write(data_filename, 2, 2n, xy)
Write (, '(a)') ' Data saved in file "' // trim(data_filename) // '"'
!
! Create the command file.
!
command_filename = 'cambered_commands.txt'
Call get_unit(command_unit)
Open (Unit=command_unit, File=command_filename, Status='replace')
Write (command_unit, '(a)') 'set term png'
Write (command_unit, '(a)') 'set grid'
Write (command_unit, '(a)') 'set style data lines'
Write (command_unit, '(a,g14.6)') 'set size ratio ', ratio
Write (command_unit, '(a)') 'set timestamp'
Write (command_unit, '(a)') 'unset key'
Write (command_unit, '(a)') 'set output "cambered.png"'
Write (command_unit, '(a)') 'set xlabel "<---X--->"'
Write (command_unit, '(a)') 'set ylabel "<---Y--->"'
Write (command_unit, '(a)') 'set title "NACA Cambered Airfoil"'
Write (command_unit, '(a)') 'plot "' // trim(data_filename) // '" using 1:2 with lines lw 3'
Write (command_unit, '(a)') 'quit'

Close (Unit=command_unit)

Write (*, '(a)') ' Created command file "' // trim(command_filename) // '".'

Return
End Subroutine test02
!**************************************
Subroutine get_unit(iunit)

!*********************************************************************72
!
!c GET_UNIT returns a free FORTRAN unit number.
!
! Discussion:
!
! A "free" FORTRAN unit number is a value between 1 and 99 which
! is not currently associated with an I/O device. A free FORTRAN unit
! number is needed in order to open a file with the OPEN command.
!
! If IUNIT = 0, then no free FORTRAN unit could be found, although
! all 99 units were checked (except for units 5, 6 and 9, which
! are commonly reserved for console I/O).
!
! Otherwise, IUNIT is a value between 1 and 99, representing a
! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6
! are special, and will never return those values.
!
! Parameters:
!
! Output, integer IUNIT, the free unit number.
!
Implicit None

Integer i
Integer iunit
Logical value

iunit = 0

Do i = 1, 99

If (i/=5 .And. i/=6 .And. i/=9) Then

  Inquire (Unit=i, Opened=value, Err=10)

  If (.Not. value) Then
    iunit = i
    Return
  End If

End If

10 Continue

End Do

Return
End Subroutine get_unit
Subroutine naca4_cambered(m, p, t, c, n, xc, xu, yu, xl, yl)

!********************************************************************72
!
!c NACA4_CAMBERED: (xu,yu), (xl,yl) for a NACA cambered 4-digit airfoil.
!
! Parameters:
!
! Input, double precision M, the maximum camber.
! 0.0 < M.
!
! Input, double precision P, the location of maximum camber.
! 0.0 < P < 1.0
!
! Input, double precision T, the maximum relative thickness.
! 0.0 < T <= 1.0
!
! Input, double precision C, the chord length.
! 0.0 < C.
!
! Input, integer N, the number of sample points.
!
! Input, double precision XC(N), points along the chord length.
! 0.0 <= XC() <= C.
!
! Output, double precision XU(N), YU(N), XL(N), YL(N), 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.
!
Implicit None

Integer n

Double Precision c
Double Precision divisor
Double Precision dycdx
Integer i
Double Precision m
Double Precision p
Double Precision t
Double Precision theta
Double Precision xc(n)
Double Precision xl(n)
Double Precision xu(n)
Double Precision yc
Double Precision yl(n)
Double Precision yt
Double Precision yu(n)

Do i = 1, n

If (0.0D+00<=xc(i)/c .And. xc(i)/c<=p) Then
  divisor = p**2
Else If (p<=xc(i)/c .And. xc(i)/c<=1.0D+00) Then
  divisor = (1.0D+00-p)**2
Else
  divisor = 1.0D+00
End If

dycdx = 2.0D+00*m*(p-xc(i)/c)/divisor

theta = atan(dycdx)

yt = 5.0D+00*t*c*&
(0.2969D+00*sqrt(xc(i)/c)+((((-0.1015D+00)*(xc(i)/c)+0.2843D+00)*&
(xc(i)/c)-0.3516D+00)*(xc(i)/c)-0.1260D+00)*(xc(i)/c))

If (0.0D+00<=xc(i)/c .And. xc(i)/c<=p) Then
  yc = m*xc(i)*(2.0D+00*p-xc(i)/c)/p**2
Else If (p<=xc(i)/c .And. xc(i)/c<=1.0D+00) Then
  yc = m*(xc(i)-c)*(2.0D+00*p-xc(i)/c-1.0D+00)/(1.0D+00-p)**2
Else
  yc = 0.0D+00
End If

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

End Do

Return
End Subroutine naca4_cambered
Subroutine naca4_symmetric(t, c, n, x, y)

!********************************************************************72
!
!c NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
!
! Parameters:
!
! Input, double precision T, the maximum relative thickness.
!
! Input, double precision C, the chord length.
!
! Input, integer N, the number of sample points.
!
! Input, double precision X(N), points along the chord length.
! 0.0 <= X() <= C.
!
! Output, double precision Y(N), for each value of X, the corresponding
! value of Y so that (X,Y) is on the upper wing surface, and (X,-Y) is on the
! lower wing surface.
!
Implicit None

Integer n

Double Precision c
Integer i
Double Precision t
Double Precision x(n)
Double Precision y(n)

Do i = 1, n
y(i) = 5.0D+00*&
tc(0.2969D+00sqrt(x(i)/c)+((((-0.1015D+00)(x(i)/c)+0.2843D+00)(x(i)/c)&
-0.3516D+00)(x(i)/c)-0.1260D+00)*(x(i)/c))
End Do

Return
End Subroutine naca4_symmetric
Subroutine r8mat_write(output_filename, m, n, table)

!*********************************************************************72
!
!c R8MAT_WRITE writes a R8MAT file.
!
! Discussion:
!
! An R8MAT is an array of R8's.
!
! Parameters:
!
! Input, character * ( * ) OUTPUT_FILENAME, the output file name.
!
! Input, integer M, the spatial dimension.
!
! Input, integer N, the number of points.
!
! Input, double precision TABLE(M,N), the data.
!
Implicit None

Integer m
Integer n

Integer j
Character () output_filename
Integer output_unit
Character *(30) string
Double Precision table(m, n)
!
! Open the file.
!
Call get_unit(output_unit)

Open (Unit=output_unit, File=output_filename, Status='replace')
!
! Create the format string.
!
If (0<m .And. 0<n) Then

Write (string, '(a1,i8,a1,i8,a1,i8,a1)') '(', m, 'g', 24, '.', 16, ')'

!
! Write the data.
!
Do j = 1, n
Write (output_unit, string) table(1:m, j)
End Do

End If
!
! Close the file.
!
Close (Unit=output_unit)

Return
End Subroutine r8mat_write
Subroutine r8vec_linspace(n, a, b, x)

!*********************************************************************72
!
!c R8VEC_LINSPACE creates a vector of linearly spaced values.
!
! Discussion:
!
! An R8VEC is a vector of R8's.
!
! 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
!
! In other words, the interval is divided into N-1 even subintervals,
! and the endpoints of intervals are used as the points.
!
! Parameters:
!
! Input, integer N, the number of entries in the vector.
!
! Input, double precision A, B, the first and last entries.
!
! Output, double precision X(N), a vector of linearly spaced data.
!
Implicit None

Integer n

Double Precision a
Double Precision b
Integer i
Double Precision x(n)

If (n==1) Then

x(1) = (a+b)/2.0D+00

Else

Do i = 1, n
  x(i) = (dble(n-i)*a+dble(i-1)*b)/dble(n-1)
End Do

End If

Return
End Subroutine r8vec_linspace
Function r8vec_max(n, a)

!*********************************************************************72
!
!c R8VEC_MAX returns the maximum value in an R8VEC.
!
! Discussion:
!
! An R8VEC is a vector of R8's.
!
! Parameters:
!
! Input, integer N, the number of entries in the array.
!
! Input, double precision A(N), the array.
!
! Output, double precision R8VEC_MAX, the value of the largest entry.
!
Implicit None

Integer n

Double Precision a(n)
Integer i
Double Precision r8vec_max
Double Precision value

value = a(1)
Do i = 2, n
value = max(value, a(i))
End Do

r8vec_max = value

Return
End Function r8vec_max
Function r8vec_min(n, a)

!*********************************************************************72
!
!c R8VEC_MIN returns the minimum value in an R8VEC.
!
! Discussion:
!
! An R8VEC is a vector of R8's.
!
! Parameters:
!
! Input, integer N, the number of entries in the array.
!
! Input, double precision A(N), the array.
!
! Output, double precision R8VEC_MIN, the value of the smallest entry.
!
Implicit None

Integer n

Double Precision a(n)
Integer i
Double Precision r8vec_min
Double Precision value

value = a(1)
Do i = 2, n
value = min(value, a(i))
End Do

r8vec_min = value

Return
End Function r8vec_min
Subroutine timestamp()

!*********************************************************************72
!
!c TIMESTAMP prints out the current YMDHMS date as a timestamp.
!
! Parameters:
!
! None
!
Implicit None

Character *(8) ampm
Integer d
Character *(8) date
Integer h
Integer m
Integer mm
Character *(9) month(12)
Integer n
Integer s
Character *(10) time
Integer y

Save month

Data month/'January ', 'February ', 'March ', 'April ',&
'May ', 'June ', 'July ', 'August ', 'September', &
'October ', 'November ', 'December '/

Call date_and_time(date, time)

Read (date, '(i4,i2,i2)') y, m, d
Read (time, '(i2,i2,i2,1x,i3)') h, n, s, mm

If (h<12) Then
ampm = 'AM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Noon'
Else
ampm = 'PM'
End If
Else
h = h - 12
If (h<12) Then
ampm = 'PM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Midnight'
Else
ampm = 'AM'
End If
End If
End If

Write (*, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') d,&
month(m), y, h, ':', n, ':', s, '.', mm, ampm

Return
End Subroutine timestamp
!*********************************************
!c$$Run Output:
!*********************************************
!*********************************************

set term png

0.000000000000000 0.000000000000000

0.000000000000000 0.000000000000000

set term png
................

0.000000000000000 -0.000000000000000
......................................
0.000000000000000 0.000000000000000

NACA_PRB:
FORTRAN90,95,&2018 version
Test the NACA library.

TEST01
NACA4_SYMMETRIC evaluates y(x) for a NACA
symmetric airfoil defined by a 4-digit code.
Data saved in file "symmetric_data.txt"
Created command file "symmetric_commands.txt".

TEST02
NACA4_CAMBERED evaluates (xu,yu) and (xl,yl) for a NACA
cambered airfoil defined by a 4-digit code.
Data saved in file "cambered_data.txt"
Created command file "cambered_commands.txt".

NACA_PRB:
Normal end of execution.

26 April 2022 4:33:47.472 AM
...Program finished with exit code 0
Press ENTER to exit console.

取一定飞机翼的参数返回系列点的坐标绘出其外形,
示意如对称和曲面的.

广州