Как создать нормали сферы, используя перекрестное произведение в декартовой координате? - PullRequest
Новая
       38

Как создать нормали сферы, используя перекрестное произведение в декартовой координате?

0 голосов
/ 22 апреля 2020

Ссылка на изображение Я написал код, который вычисляет нормали полушария. И изображение показывает цветную карту, показывающую (нормаль) * (вектор скорости) от z-направления до полушария. Итак, я должен получить самый высокий результат на северном полюсе. Но я не мог получить результаты. Я не уверен, но я думаю, что перекрестные продукты не верны. Что не так с моим кодом? Генерация вектора не так? Тогда как преодолеть трудности? Есть ли способ вычислить нормали в декартовой координате?

Может быть, мой английский sh был ужасен. Спасибо за чтение предложений. 

do 100 i=0,i_max
        do 101 j=0,j_max
                x(i,j) = r*cos(phi(i))*sin(theta(j)) 
                y(i,j) = r*sin(phi(i))*sin(theta(j))
                z(i,j) = r*cos(theta(j))
        101 continue
    100 continue

subroutine normal_and_Cp
    ! calculate cross product
    ! normal is n(i,j)
    ! First, creating 4 vectors using surrounded 4 points
    ! 

    !           (U)
    !           o
    !           |
    !           |
    !(L)o-------o-------o (R)
    !           |
    !           |
    !           o
    !           (B)

    double precision :: vec_R(3),vec_L(3), vec_U(3), vec_B(3), &
        & cross_RU(3), cross_UL(3), cross_LB(3), cross_BR(3), &
        & cross_avg(3)
    do 200 i=0,i_max-1
        do 201 j=0,j_max-1


                ! phi direction's change ?
                vec_R(1) = x(i+1,j) - x(i,j) ! x
                vec_R(2) = y(i+1,j) - y(i,j) ! y
                vec_R(3) = z(i+1,j) - z(i,j) ! z
                vec_L(1) = x(i-1,j) - x(i,j) ! x
                vec_L(2) = y(i-1,j) - y(i,j) ! y
                vec_L(3) = z(i-1,j) - z(i,j) ! z
                ! theta direction's change ?
                vec_U(1) = x(i,j+1) - x(i,j) ! x
                vec_U(2) = y(i,j+1) - y(i,j) ! y
                vec_U(3) = z(i,j+1) - z(i,j) ! z
                vec_B(1) = x(i,j-1) - x(i,j) ! x
                vec_B(2) = y(i,j-1) - y(i,j) ! y
                vec_B(3) = z(i,j-1) - z(i,j) ! z

                ! cross product using Right and upper
                cross_RU(1) = vec_R(2)*vec_U(3) - vec_R(3)*vec_U(2)
                cross_RU(2) = vec_R(3)*vec_U(1)- vec_R(1)*vec_U(3)
                cross_RU(3) = vec_R(1)*vec_U(2) - vec_R(2)*vec_U(1)

              ! cross product using upper and Left
                cross_UL(1) = vec_U(2)*vec_L(3) - vec_U(3)*vec_L(2)
                cross_UL(2) = vec_U(3)*vec_L(1) - vec_U(1)*vec_L(3) 
                cross_UL(3) = vec_U(1)*vec_L(2) - vec_U(2)*vec_L(1)

                ! cross product using Left and Bottom
                cross_LB(1) = vec_L(2)*vec_B(3) - vec_L(3)*vec_B(2)
                cross_LB(2) = vec_L(3)*vec_B(1) - vec_L(1)*vec_B(3)
                cross_LB(3) = vec_L(1)*vec_B(2) - vec_L(2)*vec_B(1)

                ! cross product using Bottom and Right
                cross_BR(1) = vec_B(2)*vec_R(3) - vec_B(3)*vec_R(2)
                cross_BR(2) = vec_B(3)*vec_R(1) - vec_B(1)*vec_R(3)
                cross_BR(3) = vec_B(1)*vec_R(2) - vec_B(2)*vec_R(1)

                cross_avg(1) = (cross_RU(1) + cross_UL(1) + cross_LB(1) + cross_BR(1))/4.0
                cross_avg(2) = (cross_RU(2) + cross_UL(2) + cross_LB(2) + cross_BR(2))/4.0
                cross_avg(3) = (cross_RU(3) + cross_UL(3) + cross_LB(3) + cross_BR(3))/4.0
                cross_avg(1) = cross_avg(1)/sqrt(cross_avg(1)*cross_avg(1)+1.0d-5)
                cross_avg(2) = cross_avg(2)/sqrt(cross_avg(2)*cross_avg(2)+1.0d-5)
                cross_avg(3) = cross_avg(3)/sqrt(cross_avg(3)*cross_avg(3)+1.0d-5)

                v(1)=0
                v(2)=0
                v(3)=-100*7
                if (dot_product(v,cross_avg) > 0) then
                    cp(i,j) = 2*((dot_product(v,cross_avg)) &
                        & /sqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)))**2
                else
                    cp(i,j)=0
                endif
        201 continue
    200 continue
end subroutine

'' '

...