!!-----------------------------------------------------------------------
!!--
!!-- Copyright (c) 2010 VINAS
!!-- All rights reserved.
!!--
!!-- This sample fortran script is not supported
!!--                                by VINAS and VINAS provides
!!-- no warranties or assurances about its fitness or merchantability.
!!-- It is provided at no cost and is for demonstration purposes only.
!!--
!!--
!!--                                          SurfExport2STL.f
!!--
!!--                                          T.Y.   Thu. June 10, 2010 
!!--
!!----------------------------------------------------------------------

!     /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
      program SurfExport2STL            ! 

      parameter(md1=1071155, md2=100385)

      character(len=256) FileInput      ! input Export file
      character(len=256) FileOutput     ! output stl file

      integer  nPlot                    ! number of plot points
      integer  nPatch                   ! number of patches

      character(len=100) charLine       ! one-line character variable
      integer nVert(0:10)               ! vertex list
      real*4  Xp(3), Yp(3), Zp(3)       ! vertex of triangle 
      real*4  Scalar(3), Threshold(3)   ! 
      real*4  U(3), V(3), W(3)          ! 

      real*4   Xnorm, Ynorm, Znorm      ! normal vector

      character(len=100) fmtVertex      ! 
      character(len=100) fmtNormal      ! 


      real*4   Xarray(md1)              ! array of X
      real*4   Yarray(md1)              ! array of Y
      real*4   Zarray(md1)              ! array of Z
      real*4   Sarray(md1)              ! array of Scalar
      real*4   Uarray(md1)              ! array of U
      real*4   Varray(md1)              ! array of V
      real*4   Warray(md1)              ! array of W
      real*4   Tarray(md1)              ! array of Z


      integer  nTriangle(md2,3)         ! vertex list


!!---
      fmtNormal="(a12,3(1x,e12.5))"
      fmtVertex="(a6,3(1x,e12.5))"
!!---

      call ReadCtrlFile(                ! read control file
     &     FileInput, FileOutput)

      call ReadPlotData(                ! read plot data
     &     FileInput, nPlot,
     &     Xarray, Yarray, Zarray, Sarray,
     &     Uarray, Varray, Warray, Tarray, md1)

      call ReadPatchData(               ! read patch data
     &     FileInput, nPatch, nTriangle, md2)


      write(*,*) "------ main loop start ----------"

      open(31, file=FileOutput, status="unknown")
      write(31,"(a)") "solid test"

      do loop=1, nPatch

         do ii=1, 3                     ! for every vertex
            Xp(ii)=Xarray(nTriangle(loop, ii))
            Yp(ii)=Yarray(nTriangle(loop, ii))
            Zp(ii)=Zarray(nTriangle(loop, ii))
         end do

         call CreateNormal(          ! normal vector detection
     &   Xnorm, Ynorm, Znorm, Xp, Yp, Zp)

          !! write(*,*) Xnorm, Ynorm, Znorm

         write(31,fmtNormal)
     &   "facet normal", Xnorm, Ynorm, Znorm
         write(31, "(a10)")"outer loop"
         do ii=1, 3
            write(31,fmtVertex)
     &      "vertex", Xp(ii), Yp(ii), Zp(ii)
         end do
         write(31, "(a7)")"endloop"
         write(31, "(a8)")"endfacet"

      end do
      close(11)


      write(31,"(a)") "end solid test"
      close(31)

      write(*,*) "------ main loop end ----------"


      stop
      end program


!     /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
      subroutine ReadCtrlFile(          ! control file の読み込み
     &     FileInput, FileOutput)

      character(len=*) FileInput        ! input Export file
      character(len=*) FileOutput       ! output stl file

!!---
      write(*,*) "sub : ReadCtrlFile"
!!---
      open(99,file='./SurfExport2STL.ctl',status='old')

      read(99,*)                        ! read through
      read(99,'(a)') FileInput
      read(99,'(a)') FileOutput
      close(99)

      write(*,*)           '-------------------------------------------'
      write(*,*)           '-----        SurfExport2STL.f     --------'
      write(*,*)           '-------------------------------------------'
      write(*,"(a20,a40)") ' input Export file :',FileInput
      write(*,"(a20,a40)") 'output    STL file :',FileOutput
      write(*,*)           '--------------------'
      write(*,*)           ''

      return
      end subroutine

!     /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
      subroutine ReadPlotData(          ! read plot data
     &     FileInput, nPlot,
     &     Xarray, Yarray, Zarray, Sarray,
     &     Uarray, Varray, Warray, Tarray, md1)

      character(len=*) FileInput        ! input Export file
      integer  nPlot                    ! number of plot points
      real*4   Xarray(md1)              ! array of X
      real*4   Yarray(md1)              ! array of Y
      real*4   Zarray(md1)              ! array of Z
      real*4   Sarray(md1)              ! array of Scalar
      real*4   Uarray(md1)              ! array of U
      real*4   Varray(md1)              ! array of V
      real*4   Warray(md1)              ! array of W
      real*4   Tarray(md1)              ! array of Z

      character(len=50) charCheck

!!---
      write(*,*) "sub : ReadPlotData"
!!---
      open(11, file=FileInput, status="old")
      do loop=1, 3
         read(11,"(a)") charCheck       ! skip header
         write(*,*) charCheck
      end do
      read(11,*) nPlot
      write(*,*) "nPlot:", nPlot
      read(11,*)                        ! skip header
      do loop=1, nPlot
         read(11, *) 
     &   Xarray(loop), Yarray(loop), Zarray(loop), Sarray(loop),
     &   Uarray(loop), Varray(loop), Warray(loop), Tarray(loop)
      end do
      close(11)

      return
      end subroutine


!     /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
      subroutine ReadPatchData(         ! read patch data
     &     FileInput, nPatch, nTriangle, md2)

      character(len=*) FileInput        ! input Export file
      integer  nPatch                   ! number of patches
      integer  nTriangle(md2,3)         ! vertex list
      character(len=100) charLine       ! one-line character variable
      integer nVert(0:10)               ! vertex list

!!---
      write(*,*) "sub : ReadPatchData"
!!---
      open(11, file=FileInput, status="old")
      do loop=1, 3
         read(11,"(a)") charCheck       ! skip header
         write(*,*) charCheck
      end do
      read(11,*) nPlot
      read(11,*)                        ! skip header
      do loop=1, nPlot
         read(11, *) 
      end do
      read(11,*)                        ! skip "GEOMETRY"
!!---
      read(11,*) nPatch
      write(*,*) "nPatch", nPatch
      nCount=0
      do loop=1, nPatch
         read(11,"(a)") charLine

         nCount=nCount+1

         if(charLine(1:1) == "3" .and.
     &      charLine(2:2) == " " )then
            read(charLine,*) (nVert(i),i=0, 3)
            do ii=1, 3
               nTriangle(nCount, ii)=nVert(ii)
            end do
         else if(charLine(1:1) == "4" .and.
     &      charLine(2:2) == " " )then
            read(charLine,*) (nVert(i),i=0, 4)

            do ii=1, 3
               nTriangle(nCount, ii)=nVert(ii)
            end do

            nCount=nCount+1
            nTriangle(nCount, 1)=nVert(1)
            nTriangle(nCount, 2)=nVert(3)
            nTriangle(nCount, 3)=nVert(4)

         else if(charLine(1:1) == "5" .and.
     &      charLine(2:2) == " " )then
            read(charLine,*) (nVert(i),i=0, 5)

            do ii=1, 3
               nTriangle(nCount, ii)=nVert(ii)
            end do

            nCount=nCount+1
            nTriangle(nCount, 1)=nVert(1)
            nTriangle(nCount, 2)=nVert(3)
            nTriangle(nCount, 3)=nVert(4)

            nCount=nCount+1
            nTriangle(nCount, 1)=nVert(1)
            nTriangle(nCount, 2)=nVert(4)
            nTriangle(nCount, 3)=nVert(5)

         end if
      end do
      nPatch=nCount
      write(*,*) "---> nPatch", nPatch
      close(11)


      return
      end subroutine




!     /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
       subroutine CreateNormal(         ! create normal vector
     &          Xnorm, Ynorm, Znorm, Xp, Yp, Zp)

      real*4   Xnorm, Ynorm, Znorm      ! normal vector
      real*4  Xp(3), Yp(3), Zp(3)       ! vertex of triangle 

      real*4   x1, y1, z1               ! 第１ベクトル成分
      real*4   x2, y2, z2               ! 第２ベクトル成分

      real*4   Dist                     ! ベクトル絶対値
!!---
      x1=Xp(3)-Xp(2)
      y1=Yp(3)-Yp(2)
      z1=Zp(3)-Zp(2)

      x2=Xp(1)-Xp(2)
      y2=Yp(1)-Yp(2)
      z2=Zp(1)-Zp(2)
!!---                                   ! outer product
      Xnorm=y1*z2-y2*z1
      Ynorm=z1*x2-z2*x1
      Znorm=x1*y2-x2*y1

!!---
      Dist=(Xnorm**2.0)+(Ynorm**2.0)+(Znorm**2.0)
      Dist=sqrt(Dist)

      Xnorm=Xnorm/Dist
      Ynorm=Ynorm/Dist
      Znorm=Znorm/Dist


      return
      end subroutine

