pyramid_plans.f90

program plan
   implicit none
   
   real(8) :: pyramid(1:3,1:5)
   
   integer :: iSide
   integer :: nod(1:4)
   
   real(8) :: pt1(1:3)
   real(8) :: pt2(1:3)
   real(8) :: pt3(1:3)
   real(8) :: v1(1:3),v2(1:3),v3(1:3),d
 
 
!Vertices
!5
!-1 -1  0  1
!+1 -1  0  2
!+1 +1  0  3
!-1 +1  0  4
! 0  0 +1  5
!Quadrilaterals
!1
!1 4 3 2 1
!
!# Set of Triangles
!Triangles
!4
!1 2 5  2
!2 3 5  3
!3 4 5  4
!4 1 5  5
!Pyramids
!1
!1 2 3 4 5  0

   write(*,'(3x,"pyramid: ")',advance='no') ; read(*,*)pyramid(1:3,1:5)
   
   do iSide=1,5
     
     select case(iSide)
     case(1) ; nod(1:4)=[1,4,3,2] ; print '("# xyz=",4(3(i2,",")," ; ")$)',int(pyramid(1:3,nod(1:4)))
     case(2) ; nod(1:3)=[1,2,5]   ; print '("# xyz=",3(3(i2,",")," ; ")$)',int(pyramid(1:3,nod(1:3)))
     case(3) ; nod(1:3)=[2,3,5]   ; print '("# xyz=",3(3(i2,",")," ; ")$)',int(pyramid(1:3,nod(1:3)))
     case(4) ; nod(1:3)=[3,4,5]   ; print '("# xyz=",3(3(i2,",")," ; ")$)',int(pyramid(1:3,nod(1:3)))
     case(5) ; nod(1:3)=[4,1,5]   ; print '("# xyz=",3(3(i2,",")," ; ")$)',int(pyramid(1:3,nod(1:3)))
     end select
     
     pt1(1:3)=pyramid(1:3,nod(1))
     pt2(1:3)=pyramid(1:3,nod(2))
     pt3(1:3)=pyramid(1:3,nod(3))
     
     v2(1:3)=pt2(1:3)-pt1(1:3)
     v3(1:3)=pt3(1:3)-pt1(1:3)
     
     !> v1= v2 ^v3
     
     v1(1) =  v2(2)*v3(3)-v2(3)*v3(2)
     v1(2) =-(v2(1)*v3(3)-v2(3)*v3(1))
     v1(3) =  v2(1)*v3(2)-v2(2)*v3(1)
    !print '("vecteur normal au plan v=",3(f5.2,1x))',v1(1:3)
     
     d=-(v1(1)*pt1(1)+v1(2)*pt1(2)+v1(3)*pt1(3))
     print '( "Plan: "$)'
     if( .not.v1(1)==0d0 ) print '(    i2," u"$)',int(v1(1))
     if( .not.v1(2)==0d0 ) print '("+",i2," v"$)',int(v1(2))
     if( .not.v1(3)==0d0 ) print '("+",i2," w"$)',int(v1(3))
     if( .not.    d==0d0 ) print '("+",i2,    $)',int(d)
     print '("=0")'
     
    !print '( "Plan: ",i2," u + ",i2," v + ",i2," w + ",i2," = 0")',int(v1(1:3)),int(d)
     
    !d=-(v1(1)*pt2(1)+v1(2)*pt2(2)+v1(3)*pt2(3))
    !print '(/"plan: ",f5.2," x + ",f5.2," y + ",f5.2," z + ",f5.2," = 0")',v1(1:3),d
     
    !d=-(v1(1)*pt3(1)+v1(2)*pt3(2)+v1(3)*pt3(3))
    !print '(/"plan: ",f5.2," x + ",f5.2," y + ",f5.2," z + ",f5.2," = 0")',v1(1:3),d
   enddo
   
end program plan