! get the normal vector of vector M
! C = M / SQRT(DOTPRODUCT(M))
SUBROUTINE VECTOR3DNORM(MX, MY, MZ, CX, CY, CZ)
REAL :: MAG
REAL :: MX,MY,MZ
MAG = SQRT(MX ** 2 + MY ** 2 + MZ ** 2)
CX = MX / MAG
CY = MY / MAG
CZ = MZ / MAG
END SUBROUTINE VECTOR3DNORM

! Calculate the normal vector of vector M and N by cross product
! C = M * N
SUBROUTINE VECTOR3DCROSS(MX, MY, MZ, NX, NY, NZ, CX, CY, CZ)
CX = MY * NZ - MZ * NY
    CY = MZ * NX - MX * NZ
    CZ = MX * NY - MY * NX
END SUBROUTINE VECTOR3DCROSS

SUBROUTINE TESTNORMALVECTOR(XN, YN, ZN)
REAL :: XA, YA, ZA, XB, YB, ZB, XC, YC, ZC
REAL :: XAB, YAB, ZAB, XBC, YBC, ZBC

PRINT*, '按顺时针方向输入平面A, B, C三点的坐标值：'
READ(*,*)XA, YA, ZA, XB, YB, ZB, XC, YC, ZC

XAB = XB - XA
YAB = YB - YA
ZAB = ZB - ZA

XBC = XC - XB
YBC = YC - YB
ZBC = ZC - ZB

! Get the normal vector by AB * BC
CALL VECTOR3DCROSS(XAB, YAB, ZAB, XBC, YBC, ZBC, XN, YN, ZN)
CALL VECTOR3DNORM(XN, YN, ZN, XN, YN, ZN)
END SUBROUTINE TESTNORMALVECTOR

PROGRAM TEST01
REAL :: XN, YN, ZN
CALL TESTNORMALVECTOR(XN, YN, ZN)
PRINT*, '法向量是:', XN, YN, ZN
END PROGRAM