如何在Fortran 90中创建圆形网格？

debugcn 发表于 Dev

凯克·奥利维拉（Kayke Oliveira）

我必须使用x² + y² = 1带有x和y从-1到的圆方程创建一个圆网格，并在每一步中1增加0.1。

我尝试了许多不同的尝试，但都失败了，所以如果有人可以帮助我，我将非常高兴。

我的尝试之一是：尝试之一

PROGRAM cir
    IMPLICIT real*8(a-h, o-z)

    real  :: x(20,20), py(20,20), ny(20,20)

    delx = -1.1
    dely = -1.1

    DO i=1,20
        delx = delx + 0.1

        DO j=1.20

            x(i,j) = delx

            ny(i,j)= -1.0 * SQRT ( 1.0 - ( x(i,j)**2 ) )
            py(i,j)= SQRT ( 1.0 - ( x(i,j)**2)  )

            PRINT*, x(i,j), ny(i,j), py(i,j)

        END DO


    END DO

    STOP

END PROGRAM cir

另一个：尝试两次

PROGRAM circle

    IMPLICIT real*8(a-h,o-z)

    real*8 :: x(20,20), y(20,20), z(20,20)

    delx = -1.0
    dely = -1.0

        DO i = 1, 20
        delx = delx + 0.1

        DO j = 1 , 20
            dely = dely + 0.1
            x(i,j) = delx
            y(i,j) = dely


            z(i,j) = x(i,j)**2 + y(i,j)**2 -1.0


        END DO
        delx = -1.0
        dely = -1.0

        PRINT*, x(i,j), y(i,j)


    END DO

    STOP

END PROGRAM circle

chw21

我们真的很难理解您想要实现的目标。目前，有两种可能性：

可能性1：你想对应的值的电网x和y中0.1之间的增量-1.0，并1.0告诉我们这个网格点是内部还是半径为1的圆外面。

这是一个这样的程序的示例，它的输出是：

program my_circle
    implicit none
    logical :: cir(-10:10, -10:10)
    real :: x, y
    integer :: i, j

    do i = -10, 10
        do j = -10, 10
            x = i / 10.0
            y = j / 10.0
            cir(j, i) = ((x**2 + y**2) <= 1.0)
        end do
    end do

    write(*, '(21A2)') print_cir(cir)

contains

    elemental function print_cir(c)
        ! returns " X" if true, " ." otherwise
        implicit none
        logical, intent(in) :: c
        character(len=2) :: print_cir

        if (c) then
            print_cir = " X"
        else
            print_cir = " ."
        end if
    end function print_cir

end program my_circle

输出：

 . . . . . . . . . . X . . . . . . . . . .
 . . . . . . X X X X X X X X X . . . . . .
 . . . . X X X X X X X X X X X X X . . . .
 . . . X X X X X X X X X X X X X X X . . .
 . . X X X X X X X X X X X X X X X X X . .
 . . X X X X X X X X X X X X X X X X X . .
 . X X X X X X X X X X X X X X X X X X X .
 . X X X X X X X X X X X X X X X X X X X .
 . X X X X X X X X X X X X X X X X X X X .
 . X X X X X X X X X X X X X X X X X X X .
 X X X X X X X X X X X X X X X X X X X X X
 . X X X X X X X X X X X X X X X X X X X .
 . X X X X X X X X X X X X X X X X X X X .
 . X X X X X X X X X X X X X X X X X X X .
 . X X X X X X X X X X X X X X X X X X X .
 . . X X X X X X X X X X X X X X X X X . .
 . . X X X X X X X X X X X X X X X X X . .
 . . . X X X X X X X X X X X X X X X . . .
 . . . . X X X X X X X X X X X X X . . . .
 . . . . . . X X X X X X X X X . . . . . .
 . . . . . . . . . . X . . . . . . . . . .

可能性2：您需要一个描述圆自身的x / y坐标列表。0.1但是，在那种情况下，增量的意义不大，因为每个象限在网格和圆上都只有3个点：（1.0 / 0.0），（0.8 / 0.6）和（0.6 / 0.8）。

因此，为此，我创建了一个小程序，在圆上打印出给定数量（当前为36个）的点，并使用三角函数计算x和y：

program my_circle2
    implicit none
    integer, parameter :: num_points = 36
    real, parameter :: pi = 4.0 * ATAN(1.0) ! Best way to calculate PI

    integer :: i
    real :: r
    real, dimension(2, num_points) :: cir

    do i = 1, num_points
        r = 2.0 * pi * real(i) / real(num_points)
        cir(1, i) = sin(r)
        cir(2, i) = cos(r)
    end do

    do i = 1, num_points
        write(*, '(I5, 3F8.4)') i, cir(:, i), my_dist(cir(:, i))
    end do

contains

    function my_dist(c)
        ! Calculates the distance of (c(1)/c(2)) from origin
        ! for verification
        implicit none
        real, dimension(2), intent(in) :: c
        real :: my_dist
        my_dist = sqrt(c(1)**2 + c(2)**2)
    end function my_dist

end program my_circle2

输出：

    1  0.1736  0.9848  1.0000
    2  0.3420  0.9397  1.0000
    3  0.5000  0.8660  1.0000
    4  0.6428  0.7660  1.0000
    5  0.7660  0.6428  1.0000
    6  0.8660  0.5000  1.0000
    7  0.9397  0.3420  1.0000
    8  0.9848  0.1736  1.0000
    9  1.0000 -0.0000  1.0000
   10  0.9848 -0.1736  1.0000
   11  0.9397 -0.3420  1.0000
   12  0.8660 -0.5000  1.0000
   13  0.7660 -0.6428  1.0000
   14  0.6428 -0.7660  1.0000
   15  0.5000 -0.8660  1.0000
   16  0.3420 -0.9397  1.0000
   17  0.1736 -0.9848  1.0000
   18 -0.0000 -1.0000  1.0000
   19 -0.1736 -0.9848  1.0000
   20 -0.3420 -0.9397  1.0000
   21 -0.5000 -0.8660  1.0000
   22 -0.6428 -0.7660  1.0000
   23 -0.7660 -0.6428  1.0000
   24 -0.8660 -0.5000  1.0000
   25 -0.9397 -0.3420  1.0000
   26 -0.9848 -0.1736  1.0000
   27 -1.0000  0.0000  1.0000
   28 -0.9848  0.1736  1.0000
   29 -0.9397  0.3420  1.0000
   30 -0.8660  0.5000  1.0000
   31 -0.7660  0.6428  1.0000
   32 -0.6428  0.7660  1.0000
   33 -0.5000  0.8660  1.0000
   34 -0.3420  0.9397  1.0000
   35 -0.1736  0.9848  1.0000
   36  0.0000  1.0000  1.0000