螺旋线数据计算一个速度向量场其满足连续方程

!******************************************
!******************************************
!螺旋线数据计算一个速度向量场其满足连续方程
!******************************************
! SINOMACH
!
!
! 螺旋线数据计算一个速度向量场其
! 满足连续方程,数据可用于图形示意表达
!
! 周勇 20221125
!*****************************************
!c$$Run, Output:
!*****************************************
Program main

!*********************************************************************72
!
!c SPIRAL_DATA_PRB tests the SPIRAL_DATA library.
!
Implicit None

Call timestamp()
Write (, '(a)') ''
Write (, '(a)') 'SPIRAL_DATA_PRB'
Write (, '(a)') ' FORTRAN95,2013,&2018 version'
Write (, '(a)') ' Test the SPIRAL_DATA library.'

Call test01()
Call test02()
Call test03()
!
! Terminate.
!
Write (, '(a)') ''
Write (, '(a)') 'SPIRAL_DATA_TEST'
Write (*, '(a)') ' Normal end of execution.'
Call timestamp()

Return
End Program main
Subroutine test01()

!*********************************************************************72
!
!c TEST01 generates a field and estimates its range.
!
Implicit None

Integer n
Parameter (n=1000)

Double Precision c
Double Precision r8vec_max
Double Precision r8vec_min
Integer seed
Double Precision u(n)
Double Precision v(n)
Double Precision x(n)
Double Precision xy_hi
Double Precision xy_lo
Double Precision y(n)

Write (, '(a)') ''
Write (, '(a)') 'TEST01'
Write (, '(a)') ' Sample a spiral velocity field and estimate'
Write (, '(a)') ' the range of the solution values.'

xy_lo = +0.0D+00
xy_hi = +1.0D+00
seed = 123456789

Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x)
Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y)
c = 1.0D+00

Call uv_spiral(n, x, y, c, u, v)

Write (, '(a)') ''
Write (, '(a)') ' Minimum Maximum'
Write (, '(a)') ''
Write (, '(a,2x,g15.7,2x,g15.7)') ' U: ', r8vec_min(n, u), r8vec_max(n, u)
Write (*, '(a,2x,g15.7,2x,g15.7)') ' V: ', r8vec_min(n, v), r8vec_max(n, v)

Return
End Subroutine test01
Subroutine test02()

!*********************************************************************72
!
!c TEST02 generates a field and samples its residuals.
!
Implicit None

Integer n
Parameter (n=1000)

Double Precision c
Double Precision pr(n)
Double Precision r8vec_amax
Double Precision r8vec_amin
Integer seed
Double Precision x(n)
Double Precision xy_hi
Double Precision xy_lo
Double Precision y(n)

Write (, '(a)') ''
Write (, '(a)') 'TEST02:'
Write (, '(a)') ' Sample a spiral velocity field and estimate the'
Write (, '(a)') ' range of residuals in the continuity equation.'

xy_lo = +0.0D+00
xy_hi = +1.0D+00
seed = 123456789

Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x)
Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y)
c = 1.0D+00

Call resid_spiral(n, x, y, c, pr)

Write (, '(a)') ''
Write (, '(a)') ' Minimum Maximum'
Write (, '(a)') ''
Write (, '(a,2x,g15.7,2x,g15.7)') ' Pr: ', r8vec_amin(n, pr), &
r8vec_amax(n, pr)

Return
End Subroutine test02
Subroutine test03()

!*********************************************************************72
!
!c TEST03 generates a field on a regular grid and plots it.
!
Implicit None

Integer x_num
Parameter (x_num=21)
Integer y_num
Parameter (y_num=21)

Double Precision c
Character *(255) header
Integer n
Double Precision s
Integer seed
Double Precision u(x_num, y_num)
Double Precision v(x_num, y_num)
Double Precision x(x_num, y_num)
Double Precision x_hi
Double Precision x_lo
Double Precision y(x_num, y_num)
Double Precision y_hi
Double Precision y_lo

Write (, '(a)') ''
Write (, '(a)') 'TEST03:'
Write (, '(a)') ' Generate a spiral velocity field on a regular grid.'
Write (, '(a)') ' Store in GNUPLOT data and command files.'

x_lo = 0.0D+00
x_hi = 1.0D+00

y_lo = 0.0D+00
y_hi = 1.0D+00

Call grid_2d(x_num, x_lo, x_hi, y_num, y_lo, y_hi, x, y)

n = x_num*y_num
c = 1.0D+00

Call uv_spiral(n, x, y, c, u, v)

header = 'spiral'
s = 0.05D+00
Call spiral_gnuplot(header, n, x, y, u, v, s)

Return
End Subroutine test03
!*****************************
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 grid_2d(x_num, x_lo, x_hi, y_num, y_lo, y_hi, x, y)

!*********************************************************************72
!
!c GRID_2D returns a regular 2D grid.
!
! Parameters:
!
! Input, integer X_NUM, the number of X values to use.
!
! Input, double precision X_LO, X_HI, the range of X values.
!
! Input, integer Y_NUM, the number of Y values to use.
!
! Input, double precision Y_LO, Y_HI, the range of Y values.
!
! Output, double precision X(X_NUM,Y_NUM), Y(X_NUM,Y_NUM),
! the coordinates of the grid.
!
Implicit None

Integer x_num
Integer y_num

Integer i
Integer j
Double Precision x(x_num, y_num)
Double Precision x_hi
Double Precision x_lo
Double Precision xi
Double Precision y(x_num, y_num)
Double Precision y_hi
Double Precision y_lo
Double Precision yj

If (x_num==1) Then
  Do j = 1, y_num
    Do i = 1, x_num
      x(i, j) = (x_lo+x_hi)/2.0D+00
    End Do
  End Do
Else
  Do i = 1, x_num
    xi = (dble(x_num-i)*x_lo+dble(i-1)*x_hi)/dble(x_num-1)
    Do j = 1, y_num
      x(i, j) = xi
    End Do
  End Do
End If

If (y_num==1) Then
  Do j = 1, y_num
    Do i = 1, x_num
      y(i, j) = (y_lo+y_hi)/2.0D+00
    End Do
  End Do
Else
  Do j = 1, y_num
    yj = (dble(y_num-j)*y_lo+dble(j-1)*y_hi)/dble(y_num-1)
    Do i = 1, x_num
      y(i, j) = yj
    End Do
  End Do
End If

Return

End Subroutine grid_2d
Function r8vec_amax(n, a)

!*********************************************************************72
!
!c R8VEC_AMAX returns the maximum absolute 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_AMAX, the value of the entry
! of largest magnitude.
!
Implicit None

Integer n

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

value = 0.0D+00
Do i = 1, n
  value = max(value, abs(a(i)))
End Do

r8vec_amax = value

Return

End Function r8vec_amax
Function r8vec_amin(n, a)

!*********************************************************************72
!
!c R8VEC_AMIN returns the minimum absolute 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 precisionA(N), the array.
!
! Output, double precision R8VEC_AMIN, the value of the entry
! of smallest magnitude.
!
Implicit None

Integer n

Double Precision a(n)
Integer i
Double Precision r8_huge
Parameter (r8_huge=1.79769313486231571D+308)
Double Precision r8vec_amin
Double Precision value

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

r8vec_amin = value

Return

End Function r8vec_amin
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 r8_huge
Parameter (r8_huge=1.79769313486231571D+308)
Double Precision r8vec_max
Double Precision value

value = -r8_huge
Do i = 1, 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 r8_huge
Parameter (r8_huge=1.79769313486231571D+308)
Double Precision r8vec_min
Double Precision value

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

r8vec_min = value

Return

End Function r8vec_min
Subroutine r8vec_uniform_ab(n, a, b, seed, r)

!*********************************************************************72
!
!c R8VEC_UNIFORM_AB returns a scaled pseudorandom R8VEC.
!
! Discussion:
!
! Each dimension ranges from A to B.
!
! Parameters:
!
! Input, integer N, the number of entries in the vector.
!
! Input, double precision A, B, the lower and upper limits.
!
! Input/output, integer SEED, the "seed" value, which should NOT be 0.
! On output, SEED has been updated.
!
! Output, double precision R(N), the vector of pseudorandom values.
!
Implicit None

Integer n

Double Precision a
Double Precision b
Integer i
Integer i4_huge
Parameter (i4_huge=2147483647)
Integer k
Integer seed
Double Precision r(n)

If (seed==0) Then
  Write (*, '(a)') ' '
  Write (*, '(a)') 'R8VEC_UNIFORM_AB - Fatal error!'
  Write (*, '(a)') '  Input value of SEED = 0.'
  Stop 1
End If

Do i = 1, n

  k = seed/127773

  seed = 16807*(seed-k*127773) - k*2836

  If (seed<0) Then
    seed = seed + i4_huge
  End If

  r(i) = a + (b-a)*dble(seed)*4.656612875D-10

End Do

Return

End Subroutine r8vec_uniform_ab
Subroutine resid_spiral(n, x, y, c, pr)

!*********************************************************************72
!
!c RESID_SPIRAL computes the residual for a spiral velocity vector field.
!
! Discussion:
!
! Note that the continuous velocity field (U,V)(X,Y) that is discretely
! sampled here satisfies the homogeneous continuity equation, that is,
! it has zero divergence. In other words:
!
! dU/dX + dV/dY = 0.
!
! This is by construction, since we have
!
! U(X,Y) = 10 * d/dY ( PHI(X) * PHI(Y) )
! V(X,Y) = -10 * d/dX ( PHI(X) * PHI(Y) )
!
! which guarantees zero divergence.
!
! The underlying function PHI is defined by
!
! PHI(Z) = ( 1 - cos ( C * pi * Z ) ) * ( 1 - Z )^2
!
! where C is a parameter.
!
! Parameters:
!
! Input, integer N, the number of evaluation points.
!
! Input, double precision X(N), Y(N), the coordinates of the
! evaluation points.
!
! Input, double precision C, a parameter, typically between 0 and 2 * PI.
!
! Output, double precision PR(N), the residual in the continuity equation.
!
Implicit None

Integer n

Double Precision c
Integer i
Double Precision pr(n)
Double Precision r8_pi
Parameter (r8_pi=3.141592653589793D+00)
Double Precision u
Double Precision ux
Double Precision v
Double Precision vy
Double Precision x(n)
Double Precision y(n)

Do i = 1, n

  u = 10.0D+00*(1.0D+00-cos(c*r8_pi*x(i)))*(1.0D+00-x(i))**2*&
  (c*r8_pi*sin(c*r8_pi*y(i))*(1.0D+00-y(i))**2-&
  (1.0D+00-cos(c*r8_pi*y(i)))*2.0D+00*(1.0D+00-y(i)))

  ux = 10.0D+00*(c*r8_pi*sin(c*r8_pi*x(i))*(1.0D+00-x(i))**2-&
  (1.0D+00-cos(c*r8_pi*x(i)))*2.0D+00*(1.0D+00-x(i)))*&
  (c*r8_pi*sin(c*r8_pi*y(i))*(1.0D+00-y(i))**2-&
  (1.0D+00-cos(c*r8_pi*y(i)))*2.0D+00*(1.0D+00-y(i)))

  v = -10.0D+00*(1.0D+00-cos(c*r8_pi*y(i)))*(1.0D+00-y(i))**2*&
  (c*r8_pi*sin(c*r8_pi*x(i))*(1.0D+00-x(i))**2-&
  (1.0D+00-cos(c*r8_pi*x(i)))*2.0D+00*(1.0D+00-x(i)))

  vy = -10.0D+00*(c*r8_pi*sin(c*r8_pi*x(i))*(1.0D+00-x(i))**2-&
  (1.0D+00-cos(c*r8_pi*x(i)))*2.0D+00*(1.0D+00-x(i)))*&
  (c*r8_pi*sin(c*r8_pi*y(i))*(1.0D+00-y(i))**2-&
  (1.0D+00-cos(c*r8_pi*y(i)))*2.0D+00*(1.0D+00-y(i)))

  pr(i) = ux + vy

End Do

Return

End Subroutine resid_spiral
Subroutine spiral_gnuplot(header, n, x, y, u, v, s)

!*********************************************************************72
!
!c SPIRAL_GNUPLOT writes the spiral vector field to files for GNUPLOT.
!
! Parameters:
!
! Input, character * ( * ) HEADER, a header to be used to name the files.
!
! Input, integer N, the number of evaluation points.
!
! Input, real X(N), Y(N), the coordinates of the evaluation points.
!
! Input, real U(N), V(N), the velocity components.
!
! Input, real S, a scale factor for the velocity vectors.
!
Implicit None

Integer n

Character *(255) command_filename
Integer command_unit
Character *(255) data_filename
Integer data_unit
Character *(*) header
Integer i
Character *(255) plot_filename
Double Precision s
Double Precision u(n)
Double Precision v(n)
Double Precision x(n)
Double Precision y(n)

!
! Write the data file.
!
data_filename = trim(header) // '_data.txt'

Call get_unit(data_unit)

Open (Unit=data_unit, File=data_filename, Status='replace')

Do i = 1, n
  Write (data_unit, '(2x,g15.7,2x,g15.7,2x,g15.7,2x,g15.7)')&
  x(i), y(i), s*u(i), s*v(i)
End Do

Close (Unit=data_unit)

Write (*, '(a)') ''
Write (*, '(a)') '  Data written to "' // trim(data_filename) // '".'

!
! Write the command file.
!
command_filename = trim(header) // '_commands.txt'
Call get_unit(command_unit)

plot_filename = trim(header) // '.png'

Open (Unit=command_unit, File=command_filename, Status='replace')

Write (command_unit, '(a)') '#  ' // trim(command_filename)
Write (command_unit, '(a)') '#'
Write (command_unit, '(a)') 'set term png'
Write (command_unit, '(a)') 'set output "' // trim(plot_filename) // '"'
Write (command_unit, '(a)') '#'
Write (command_unit, '(a)') '#  Add s and labels.'
Write (command_unit, '(a)') '#'
Write (command_unit, '(a)') 'set xlabel "<--- X --->"'
Write (command_unit, '(a)') 'set ylabel "<--- Y --->"'
Write (command_unit, '(a)') 'set "Spiral velocity flow"'
Write (command_unit, '(a)') 'unset key'
Write (command_unit, '(a)') '#'
Write (command_unit, '(a)') '#  Add grid lines.'
Write (command_unit, '(a)') '#'
Write (command_unit, '(a)') 'set grid'
Write (command_unit, '(a)') 'set size ratio -1'
Write (command_unit, '(a)') '#'
Write (command_unit, '(a)') '#  Timestamp the plot.'
Write (command_unit, '(a)') '#'
Write (command_unit, '(a)') 'set timestamp'
Write (command_unit, '(a)') 'plot "' // trim(data_filename) //&
'" using 1:2:3:4 with vectors \'
Write (command_unit, '(a)') '  head filled lt 2 linecolor rgb "blue"'
Write (command_unit, '(a)') 'quit'

Close (Unit=command_unit)

Write (*, '(a)') '  Commands written to "' // trim(command_filename) // '".'

Return

End Subroutine spiral_gnuplot
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
Subroutine uv_spiral(n, x, y, c, u, v)

!*********************************************************************72
!
!c UV_SPIRAL computes a spiral velocity vector field.
!
! Discussion:
!
! Note that the continuous velocity field (U,V)(X,Y) that is discretely
! sampled here satisfies the homogeneous continuity equation, that is,
! it has zero divergence. In other words:
!
! dU/dX + dV/dY = 0.
!
! This is by construction, since we have
!
! U(X,Y) = 10 * d/dY ( PHI(X) * PHI(Y) )
! V(X,Y) = -10 * d/dX ( PHI(X) * PHI(Y) )
!
! which guarantees zero divergence.
!
! The underlying function PHI is defined by
!
! PHI(Z) = ( 1 - cos ( C * pi * Z ) ) * ( 1 - Z )^2
!
! where C is a parameter.
!
! Parameters:
!
! Input, integer N, the number of evaluation points.
!
! Input, double precision X(N), Y(N), the coordinates of the
! evaluation points.
!
! Input, double precision C, a parameter, typically between 0 and 2 * PI.
!
! Output, double precision U(N), V(N), the velocity components.
!
Implicit None

Integer n

Double Precision c
Integer i
Double Precision r8_pi
Parameter (r8_pi=3.141592653589793D+00)
Double Precision u(n)
Double Precision v(n)
Double Precision x(n)
Double Precision y(n)

Do i = 1, n

  u(i) = 10.0D+00*(1.0D+00-cos(c*r8_pi*x(i)))*(1.0D+00-x(i))&
  **2*(c*r8_pi*sin(c*r8_pi*y(i))*(1.0D+00-y(i))**2-&
  (1.0D+00-cos(c*r8_pi*y(i)))*2.0D+00*(1.0D+00-y(i)))

  v(i) = -10.0D+00*(1.0D+00-cos(c*r8_pi*y(i)))*&
  (1.0D+00-y(i))**2*(c*r8_pi*sin(c*r8_pi*x(i))*&
  (1.0D+00-x(i))**2-(1.0D+00-cos(c*r8_pi*x(i)))*2.0D+00*(1.0D+00-x(i)))

End Do

Return

End Subroutine uv_spiral

SPIRAL_DATA_PRB
FORTRAN90,95,&2013 version
Test the SPIRAL_DATA library.

TEST01
Sample a spiral velocity field and estimate
the range of the solution values.

Minimum       Maximum

U: -1.854039 2.224943
V: -2.231540 1.852338

TEST02:
Sample a spiral velocity field and estimate the
range of residuals in the continuity equation.

Minimum       Maximum

Pr: 0.000000 0.000000

TEST03:
Generate a spiral velocity field on a regular grid.
Store in GNUPLOT data and command files.

Data written to "spiral_data.txt".
Commands written to "spiral_commands.txt".

SPIRAL_DATA_TEST
Normal end of execution.
March 2022
...Program finished with exit code 0
Press ENTER to exit console.
!****************************
0.000000 0.000000 0.000000 -0.000000
0.5000000E-01 0.000000 0.000000 -0.000000
0.1000000 0.000000 0.000000 -0.000000
.......................................
1.000000 1.000000 0.000000 -0.000000
!****************************
螺旋线数据计算一个速度向量场其满足连续方程,数据可用于图形示意表达.
!****************************
程序运行验证数值如上, 图形处理输出示意.
!****************************
广州