!------------------------------------------------------------------------------
!        IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE         : Mohid Model
! PROJECT       : Digital Terrain Creator
! MODULE        : Drawing
! URL           : http://www.mohid.com
! AFFILIATION   : IST/MARETEC, Marine Modelling Group
! DATE          : May 2003
! REVISION      : Luis Fernandes / Frank Braunschweig - v4.0
! DESCRIPTION   : Module to handle Geometric objects
!
!------------------------------------------------------------------------------
Module ModuleDrawing

    use ModuleGlobalData
    use ModuleEnterData
    use ModuleFunctions

    implicit none
    
    private
    
    !Subroutines--------------------------------------------------------------
    public  ::    New
    private ::    Add
    public  ::    SetLimits
    public  ::    IsVisible
    public  ::    IsPointInsidePolygon
    public  ::    IsPointInsideCircle
    public  ::    WriteItem

   
    private ::    NewPolygon
    private ::    NewXYZPoint
    private ::    NewXYZPoint_V2
    private ::    NewXYZPoint_V3
    private ::    NewXYZPoint_V4
    interface     New
        module procedure NewPolygon
        module procedure NewXYZPoint
        module procedure NewXYZPoint_V2 
        module procedure NewXYZPoint_V3 
        module procedure NewXYZPoint_V4 
        module procedure NewLine
    end interface New
    
    private ::    AddPolygon
    private ::    AddXYZPoint
    private ::    AddLine
    interface     Add
        module procedure AddPolygon
        module procedure AddXYZPoint
        module procedure AddLine
    end interface Add

    private ::    SetLimitsPolygon
    private ::    SetLimitsXYZ
    interface     SetLimits
        module procedure SetLimitsPolygon
        module procedure SetLimitsXYZ
    end interface SetLimits

    private ::    WriteItemPolygon
    private ::    WriteItemXYZ
    private ::    WriteItemXYZ_v2
    interface     WriteItem
        module procedure WriteItemPolygon
        module procedure WriteItemXYZ
        module procedure WriteItemXYZ_v2
    end interface WriteItem

    !Parameter-----------------------------------------------------------------
    integer(4), parameter  :: TypeX_Y_Z     =  1
    integer(4), parameter  :: TypeX_Y_Z_P   =  2    
    integer(4), parameter  :: TypeT_X_Y_Z_P =  3    
    
    
    !Types--------------------------------------------------------------
    public   T_Limits
    type     T_Limits
        real                                    :: Left         = null_real
        real                                    :: Right        = null_real
        real                                    :: Top          = null_real
        real                                    :: Bottom       = null_real
        real                                    :: MaximumValue = 999999.
        real                                    :: MinimumValue =-999999.
    end type T_Limits
    
    public   T_Point                                                
    type     T_Point                                                
        integer                                 :: I            = null_int
        integer                                 :: J            = null_int
    end type T_Point                                                                                          

    public   T_PointF
    type     T_PointF
        real                                    :: X            = null_real
        real                                    :: Y            = null_real
    end type T_PointF
   
    public   T_XYZPoints
    type     T_XYZPoints
        integer                                 :: ID           = null_int
        real,           dimension(:), pointer   :: X
        real,           dimension(:), pointer   :: Y
        real,           dimension(:), pointer   :: Z
        real,           dimension(:), pointer   :: T
        real,           dimension(:), pointer   :: P                
        logical,        dimension(:), pointer   :: Inside
        integer                                 :: FileFormat
        integer                                 :: Count
        type(T_Limits)                          :: Limits
        type (T_XYZPoints),           pointer   :: Next
    end type T_XYZPoints
    
    public   T_Segment
    type     T_Segment
        type(T_PointF)                          :: StartAt
        type(T_PointF)                          :: EndAt
    end type T_Segment

    public   T_Polygon
    type     T_Polygon
        integer                                 :: ID           = null_int
        type(T_Point ),  dimension(:), pointer  :: Vertices
        type(T_PointF),  dimension(:), pointer  :: VerticesF
        type(T_Limits)                          :: Limits
        integer                                 :: Count
        type(T_Polygon),              pointer   :: Next 
    end type T_Polygon


    public   T_Lines
    type     T_Lines
        integer                                 :: ID           = null_int
        real, dimension(:), pointer                 :: X
        real, dimension(:), pointer                 :: Y
        integer                                     :: nNodes
        type(T_Lines) , pointer                     :: Next
    end type T_Lines

    contains
    
    
    !--------------------------------------------------------------------------
    
    subroutine NewPolygon(Polygons, PolygonsFileName)
        
        !Arguments-------------------------------------------------------------
        type(T_Polygon),                 pointer    :: Polygons
        character(len=*),  intent(IN)               :: PolygonsFileName

        !Local-----------------------------------------------------------------
        type(T_Polygon),                 pointer    :: Polygon
        integer                                     :: PolygonsFile = 0
        integer                                     :: ClientNumber, StartLine, EndLine
        integer                                     :: FromBlock, CurrentLineNumber, CurrentVertix
        character(len = line_length)                :: FullBufferLine
        logical                                     :: BlockFound = OFF
        real, dimension(:),              pointer    :: PointCoordinates
        integer                                     :: STAT_CALL, iflag

        !Begin-----------------------------------------------------------------

        call ConstructEnterData(PolygonsFile, PolygonsFileName, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'NewPolygon - ModuleDrawing - ERR01'
    
do1 :       do
                call ExtractBlockFromBuffer(PolygonsFile,                                &
                                            ClientNumber,                                &
                                            '<beginpolygon>',                           &
                                            '<endpolygon>',                             &
                                            BlockFound,                                  &
                                            STAT = STAT_CALL)

if1 :           if      (STAT_CALL .EQ. SUCCESS_ ) then    

if2 :               if (BlockFound) then

                        call Add(Polygons, Polygon)

                        call GetExtractType(FromBlock = FromBlock)

                        call GetBlockSize(PolygonsFile,                                  &
                                          ClientNumber,                                  &
                                          StartLine,                                     &
                                          EndLine,                                       &
                                          FromBlock,                                     &
                                          STAT = STAT_CALL)

                        !ATTENTION!!! there is a bug in module EnterData regarding
                        !             the STAT variable in sub GetBlockSize, see also
                        !             modules Lagrangian and TimeSerie

                        Polygon%Count = EndLine - StartLine 
                    
                        allocate(Polygon%VerticesF(1:Polygon%Count))
                        allocate(PointCoordinates (1:2))
                    
                        CurrentVertix = 1

                        do CurrentLineNumber = StartLine + 1 , EndLine - 1 
                        
                            call GetFullBufferLine(PolygonsFile,                         &
                                                   CurrentLineNumber,                    &
                                                   FullBufferLine,                       &
                                                   STAT = STAT_CALL)
                            if(STAT_CALL .ne. SUCCESS_)stop 'NewPolygon - ModuleDrawing - ERR02'
                        
                            call GetExtractType(FromBlock = FromBlock)


                            call GetData(PointCoordinates,                               &
                                         PolygonsFile,                                   &
                                         flag = iflag,                                   &
                                         SearchType = FromBlock,                         &
                                         Buffer_Line = CurrentLineNumber,                &
                                         STAT        = STAT_CALL)
                            if(STAT_CALL .ne. SUCCESS_ .and.                             &
                               STAT_CALL .ne. SIZE_ERR_) stop 'NewPolygon - ModuleDrawing - ERR03'
                            
                            Polygon%VerticesF(CurrentVertix)%X  = PointCoordinates(1)
                            Polygon%VerticesF(CurrentVertix)%Y  = PointCoordinates(2)
                       
                            CurrentVertix = CurrentVertix + 1
                    
                        end do
                    
                        !Close polygon
                        Polygon%VerticesF(CurrentVertix)%X  = Polygon%VerticesF(1)%X
                        Polygon%VerticesF(CurrentVertix)%Y  = Polygon%VerticesF(1)%Y

                        call SetLimits(Polygon)

                        deallocate(PointCoordinates)
                        nullify   (Polygon)
                    else

                        call Block_Unlock(PolygonsFile, ClientNumber, STAT = STAT_CALL) 
                        if(STAT_CALL .ne. SUCCESS_)stop 'NewPolygon - ModuleDrawing - ERR04'
                       

                        exit do1    !No more blocks

                    end if if2

                else if (STAT_CALL .EQ. BLOCK_END_ERR_) then if1

                    if(STAT_CALL .ne. SUCCESS_)stop 'NewPolygon - ModuleDrawing - ERR05'

                end if if1
            end do do1
        
            call KillEnterData(PolygonsFile, STAT = STAT_CALL)
            if(STAT_CALL .ne. SUCCESS_)stop 'NewPolygon - ModuleDrawing - ERR06'

    end subroutine NewPolygon
    !--------------------------------------------------------------------------

    
    !--------------------------------------------------------------------------
    subroutine AddPolygon (Polygons, ObjPolygon)

        !Arguments-------------------------------------------------------------
        type (T_Polygon), pointer                   :: Polygons
        type (T_Polygon), pointer                   :: ObjPolygon

        !Local-----------------------------------------------------------------
        type (T_Polygon), pointer                   :: NewPolygon
        type (T_Polygon), pointer                   :: PreviousPolygon
        integer, save                               :: NextID = 1
        
        !Allocates new instance
        allocate (NewPolygon)
        nullify  (NewPolygon%Next)

        !Insert New Instance into list and makes Current point to it
        if (.not. associated(Polygons)) then
            Polygons             => NewPolygon
            ObjPolygon           => NewPolygon
        else
            PreviousPolygon      => Polygons
            ObjPolygon           => Polygons%Next
            do while (associated(ObjPolygon))
                PreviousPolygon  => ObjPolygon
                ObjPolygon       => ObjPolygon%Next
            enddo
            ObjPolygon           => NewPolygon
            PreviousPolygon%Next => NewPolygon
        endif

        !Attributes ID
        ObjPolygon%ID            = NextID
        NextID                   = NextID + 1

    end subroutine AddPolygon

    !--------------------------------------------------------------------------
   
    subroutine NewLine(Lines, LinesFileName)
        
        !Arguments-------------------------------------------------------------
        type(T_Lines),                 pointer      :: Lines
        character(len=*),  intent(IN)               :: LinesFileName

        !Local-----------------------------------------------------------------
        integer                                     :: ObjEnterData = 0
        logical                                     :: BlockFound
        integer                                     :: ClientNumber, STAT_CALL
        integer                                     :: iflag, FromBlock
        integer                                     :: ret, StartLine, EndLine
        integer                                     :: iPoint, CurrentLineNumber
        type(T_Lines), pointer                      :: CurrLine
        real, dimension(:), pointer                 :: PointCoordinates

        nullify(CurrLine)

        call ConstructEnterData(ObjEnterData, LinesFileName, STAT = ret) 
        if (ret == 0) then

            !ColorStuff
            BlockFound = .true.
            do while (BlockFound)

                call ExtractBlockFromBuffer(ObjEnterData, ClientNumber,              &
                                            "<begin_line>", "<end_line>",            &
                                            BlockFound,StartLine, EndLine,           &
                                            STAT = STAT_CALL)
                if (BlockFound) then
                    
                    call AddLine (Lines, CurrLine)

                    CurrLine%nNodes = EndLine - StartLine - 1
                
                    iPoint = 1

                    allocate(CurrLine%X(1:CurrLine%nNodes))
                    allocate(CurrLine%Y(1:CurrLine%nNodes))
                    allocate(PointCoordinates(1:2))

                    do CurrentLineNumber = StartLine + 1 , EndLine - 1 
                    
                        call GetData(PointCoordinates,                               &
                                     ObjEnterData,                                   &
                                     flag = iflag,                                   &
                                     SearchType = FromBlock,                         &
                                     Buffer_Line = CurrentLineNumber,                &
                                     STAT        = STAT_CALL)
                        if(STAT_CALL .ne. SUCCESS_ .and.                             &
                           STAT_CALL .ne. SIZE_ERR_) stop 'NewLine - ModuleDrawing - ERR03'

                        CurrLine%X(iPoint) = PointCoordinates(1)
                        CurrLine%Y(iPoint) = PointCoordinates(2)

                        iPoint = iPoint + 1

                    end do

                endif
            enddo
        end if

        call KillEnterData(ObjEnterData, ret)

    end subroutine NewLine   
   
    !--------------------------------------------------------------------------

    !--------------------------------------------------------------------------
    subroutine AddLine (Lines, ObjLine)

        !Arguments-------------------------------------------------------------
        type (T_Lines), pointer                   :: Lines
        type (T_Lines), pointer                   :: ObjLine

        !Local-----------------------------------------------------------------
        type (T_Lines), pointer                   :: NewLine
        type (T_Lines), pointer                   :: PreviousLine
        integer, save                            :: NextID = 1
        
        !Allocates new instance
        allocate (NewLine)
        nullify  (NewLine%Next)

        !Insert New Instance into list and makes Current point to it
        if (.not. associated(Lines)) then
            Lines             => NewLine
            ObjLine           => NewLine
        else
            PreviousLine      => Lines
            ObjLine           => Lines%Next
            do while (associated(ObjLine))
                PreviousLine  => ObjLine
                ObjLine       => ObjLine%Next
            enddo
            ObjLine           => NewLine
            PreviousLine%Next => NewLine
        endif

        !Attributes ID
        ObjLine%ID            = NextID
        NextID                = NextID + 1

    end subroutine AddLine

    !--------------------------------------------------------------------------
    subroutine NewXYZPoint(XYZPoints, XYZPointsFileName, DepthReferential, AddFactor, FileFormat)
        
        !Arguments-------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoints
        character(len=*),  intent(IN)               :: XYZPointsFileName
        real,    optional, intent(IN)               :: DepthReferential, AddFactor
        integer, optional, intent(IN)               :: FileFormat

        !Local-----------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoint
        integer                                     :: XYZPointFile = 0
        integer                                     :: ClientNumber
        integer                                     :: StartLine, EndLine
        integer                                     :: FromBlock, CurrentLineNumber
        character(len = line_length)                :: FullBufferLine
        logical                                     :: BlockFound = OFF
        real, dimension(:),              pointer    :: PointCoordinates
        integer                                     :: i, STAT_CALL
        real                                        :: DepthReferential_, AddFactor_

        !Begin-----------------------------------------------------------------

        call ConstructEnterData(XYZPointFile, XYZPointsFileName, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'NewXYZPoint - ModuleDrawing - ERR01'

        if (present(DepthReferential)) then
            DepthReferential_ = DepthReferential
        else
            DepthReferential_ = 1.
        endif

        if (present(AddFactor)) then
            AddFactor_ = AddFactor
        else
            AddFactor_ = 0.
        endif
        
do1 :       do
                call ExtractBlockFromBuffer(XYZPointFile,                                &
                                            ClientNumber,                                &
                                            '<begin_xyz>',                               &
                                            '<end_xyz>',                                 &
                                            BlockFound,                                  &
                                            STAT = STAT_CALL)

if1 :           if      (STAT_CALL .EQ. SUCCESS_ ) then    

if2 :               if (BlockFound) then
                        
                        call GetExtractType(FromBlock = FromBlock)
                        
                        call GetBlockSize(XYZPointFile,                                  &
                                          ClientNumber,                                  &
                                          StartLine,                                     &
                                          EndLine,                                       &
                                          FromBlock,                                     &
                                          STAT = STAT_CALL)

                        !ATTENTION!!! there is a bug in module EnterData regarding
                        !             the STAT variable in sub GetBlockSize, see also
                        !             modules Lagrangian and TimeSerie

                        
                        call Add(XYZPoints, XYZPoint)
                        
                        if (present(FileFormat)) then
                            XYZPoint%FileFormat = FileFormat
                        else
                            XYZPoint%FileFormat = TypeX_Y_Z
                        endif

                        if      (XYZPoint%FileFormat == TypeX_Y_Z)     then
                            allocate(PointCoordinates(1:3))
                        else if (XYZPoint%FileFormat == TypeX_Y_Z_P) then
                            allocate(PointCoordinates(1:4))
                        else if (XYZPoint%FileFormat == TypeT_X_Y_Z_P) then
                            allocate(PointCoordinates(1:5))
                        endif

                        XYZPoint%Count = EndLine - StartLine - 1
                        
                        allocate(XYZPoint%X     (1:(XYZPoint%Count))) 
                        allocate(XYZPoint%Y     (1:(XYZPoint%Count)))
                        allocate(XYZPoint%Z     (1:(XYZPoint%Count))) 
                        allocate(XYZPoint%Inside(1:(XYZPoint%Count))) 

                        if      (XYZPoint%FileFormat == TypeX_Y_Z_P) then
                            allocate(XYZPoint%P(1:(XYZPoint%Count)))                            
                        else if (XYZPoint%FileFormat == TypeT_X_Y_Z_P) then
                            allocate(XYZPoint%T(1:(XYZPoint%Count)))
                            allocate(XYZPoint%P(1:(XYZPoint%Count)))                            
                        endif
                        
                        XYZPoint%Inside = .false.

                        i = 0

                        do CurrentLineNumber = StartLine + 1 , EndLine - 1 
                            
                            i = i + 1

                            call GetFullBufferLine(XYZPointFile,                         &
                                                   CurrentLineNumber,                    &
                                                   FullBufferLine,                       &
                                                   STAT = STAT_CALL)
                            if(STAT_CALL .ne. SUCCESS_)stop 'NewXYZPoint - ModuleDrawing - ERR02'
                        
                            call GetExtractType(FromBlock = FromBlock)

                            !Read XYZ even if there is a comment after Z
                            call GetData(PointCoordinates,                               &
                                         XYZPointFile,                                   &
                                         FromBlock,                                      &
                                         Buffer_Line = CurrentLineNumber,                &
                                         STAT        = STAT_CALL)
                            if(STAT_CALL .ne. SUCCESS_ .and.                             &
                               STAT_CALL .ne. SIZE_ERR_) then
                                
                                write(*,*)'Error in Line Number = ', CurrentLineNumber
                                write(*,*)'File name            = ', trim(XYZPointsFileName)
                                stop 'NewXYZPoint - ModuleDrawing - ERR03'

                            end if
                            
                            if (XYZPoint%FileFormat == TypeX_Y_Z) then
                                XYZPoint%X(i)  = PointCoordinates(1)
                                XYZPoint%Y(i)  = PointCoordinates(2)
                                XYZPoint%Z(i)  = PointCoordinates(3)                            
                            else if (XYZPoint%FileFormat == TypeX_Y_Z_P) then
                                XYZPoint%X(i)  = PointCoordinates(1)
                                XYZPoint%Y(i)  = PointCoordinates(2)
                                XYZPoint%Z(i)  = PointCoordinates(3)                            
                                XYZPoint%P(i)  = PointCoordinates(4)                            
                            else if (XYZPoint%FileFormat == TypeT_X_Y_Z_P) then
                                XYZPoint%T(i)  = PointCoordinates(1)                            
                                XYZPoint%X(i)  = PointCoordinates(2)
                                XYZPoint%Y(i)  = PointCoordinates(3)
                                XYZPoint%Z(i)  = PointCoordinates(4)                            
                                XYZPoint%P(i)  = PointCoordinates(5)                            
                            endif
                            
                            XYZPoint%Z(i)  = XYZPoint%Z(i) * DepthReferential_ + AddFactor_

                        end do
                        
                        call SetLimits(XYZPoint)
                        
                        deallocate(PointCoordinates)
                        nullify   (XYZPoint)
                    else

                        call Block_Unlock(XYZPointFile, ClientNumber, STAT = STAT_CALL) 
                        if(STAT_CALL .ne. SUCCESS_)stop 'NewXYZPoint - ModuleDrawing - ERR04'
                       

                        exit do1    !No more blocks

                    end if if2

                else if (STAT_CALL .EQ. BLOCK_END_ERR_) then if1

                    if(STAT_CALL .ne. SUCCESS_)stop 'NewXYZPoint - ModuleDrawing - ERR05'

                end if if1
            end do do1
        
            call KillEnterData(XYZPointFile, STAT = STAT_CALL)
            if(STAT_CALL .ne. SUCCESS_)stop 'NewXYZPoint - Polygon - ModuleDrawing'

    end subroutine NewXYZPoint
    !--------------------------------------------------------------------------
    
   
    !--------------------------------------------------------------------------
    subroutine NewXYZPoint_V2(XYZPoints, XVector, YVector, ZVector)
        
        !Arguments-------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoints
        real,   dimension(:),            pointer    :: XVector, YVector, ZVector

        !Local-----------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoint
        integer                                     :: i
 
        !Begin-----------------------------------------------------------------

 
        
        call Add(XYZPoints, XYZPoint)

        XYZPoint%Count = size(XVector)
        
        allocate(XYZPoint%X     (1:(XYZPoint%Count))) 
        allocate(XYZPoint%Y     (1:(XYZPoint%Count)))
        allocate(XYZPoint%Z     (1:(XYZPoint%Count))) 
        allocate(XYZPoint%Inside(1:(XYZPoint%Count))) 
        
        XYZPoint%FileFormat = TypeX_Y_Z
       
        XYZPoint%Inside = .false.

         do i = 1 , XYZPoint%Count
           
            XYZPoint%X(i)  = XVector(i)
            XYZPoint%Y(i)  = YVector(i)
            XYZPoint%Z(i)  = ZVector(i)

        end do
        
        call SetLimits(XYZPoint)
        

        nullify   (XYZPoint)

    end subroutine NewXYZPoint_V2
    !--------------------------------------------------------------------------
    
    !--------------------------------------------------------------------------
    subroutine NewXYZPoint_V3(XYZPoints, XVector, YVector, ZVector, PVector)
        
        !Arguments-------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoints
        real,   dimension(:),            pointer    :: XVector, YVector, ZVector, PVector

        !Local-----------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoint
        integer                                     :: i
 
        !Begin-----------------------------------------------------------------

 
        
        call Add(XYZPoints, XYZPoint)

        XYZPoint%Count = size(XVector)
        
        allocate(XYZPoint%X     (1:(XYZPoint%Count))) 
        allocate(XYZPoint%Y     (1:(XYZPoint%Count)))
        allocate(XYZPoint%Z     (1:(XYZPoint%Count))) 
        allocate(XYZPoint%Inside(1:(XYZPoint%Count))) 
        allocate(XYZPoint%P     (1:(XYZPoint%Count)))                 
        
        XYZPoint%FileFormat = TypeX_Y_Z_P
       
        XYZPoint%Inside = .false.

         do i = 1 , XYZPoint%Count
            XYZPoint%X(i)  = XVector(i)
            XYZPoint%Y(i)  = YVector(i)
            XYZPoint%Z(i)  = ZVector(i)
            XYZPoint%P(i)  = PVector(i)
        end do
        
        call SetLimits(XYZPoint)
        
        nullify   (XYZPoint)

    end subroutine NewXYZPoint_V3
    !--------------------------------------------------------------------------
    !--------------------------------------------------------------------------
    subroutine NewXYZPoint_V4(XYZPoints, TVector, XVector, YVector, ZVector, PVector)
        
        !Arguments-------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoints
        real,   dimension(:),            pointer    :: TVector, XVector, YVector, ZVector, PVector

        !Local-----------------------------------------------------------------
        type(T_XYZPoints),               pointer    :: XYZPoint
        integer                                     :: i
 
        !Begin-----------------------------------------------------------------

 
        
        call Add(XYZPoints, XYZPoint)

        XYZPoint%Count = size(XVector)
        
        allocate(XYZPoint%X     (1:(XYZPoint%Count))) 
        allocate(XYZPoint%Y     (1:(XYZPoint%Count)))
        allocate(XYZPoint%Z     (1:(XYZPoint%Count))) 
        allocate(XYZPoint%Inside(1:(XYZPoint%Count))) 
        allocate(XYZPoint%T     (1:(XYZPoint%Count)))         
        allocate(XYZPoint%P     (1:(XYZPoint%Count)))                 
        
        XYZPoint%FileFormat = TypeT_X_Y_Z_P
       
        XYZPoint%Inside = .false.

         do i = 1 , XYZPoint%Count
            XYZPoint%T(i)  = TVector(i)
            XYZPoint%X(i)  = XVector(i)
            XYZPoint%Y(i)  = YVector(i)
            XYZPoint%Z(i)  = ZVector(i)
            XYZPoint%P(i)  = PVector(i)
        end do
        
        call SetLimits(XYZPoint)
        
        nullify   (XYZPoint)

    end subroutine NewXYZPoint_V4
    !--------------------------------------------------------------------------    
    !--------------------------------------------------------------------------
    subroutine AddXYZPoint (XYZPoints, ObjXYZPoint)

        !Arguments-------------------------------------------------------------
        type (T_XYZPoints), pointer                 :: XYZPoints
        type (T_XYZPoints), pointer                 :: ObjXYZPoint

        !Local-----------------------------------------------------------------
        type (T_XYZPoints), pointer                 :: NewXYZPoint
        type (T_XYZPoints), pointer                 :: PreviousXYZPoint
        integer, save                               :: NextID = 1
        
        !Allocates new instance
        allocate (NewXYZPoint)
        nullify  (NewXYZPoint%Next)

        if (.not. associated(XYZPoints)) then
            XYZPoints            => NewXYZPoint
            ObjXYZPoint          => NewXYZPoint
        else
            PreviousXYZPoint     => XYZPoints
            ObjXYZPoint          => XYZPoints%Next
            do while (associated(ObjXYZPoint))
                PreviousXYZPoint => ObjXYZPoint
                ObjXYZPoint      => ObjXYZPoint%Next
            enddo
            ObjXYZPoint          => NewXYZPoint
            PreviousXYZPoint%Next=> NewXYZPoint
        endif


        !Attributes ID
        ObjXYZPoint%ID           = NextID
        NextID                   = NextID + 1

    end subroutine AddXYZPoint


    !--------------------------------------------------------------------------
    
    subroutine SetLimitsPolygon(Polygon)

        !Arguments-------------------------------------------------------------
        type (T_Polygon), pointer                   :: Polygon
        
        !Begin-----------------------------------------------------------------

        Polygon%Limits%Left   = minval(Polygon%VerticesF%X)
        Polygon%Limits%Right  = maxval(Polygon%VerticesF%X)
        Polygon%Limits%Bottom = minval(Polygon%VerticesF%Y)
        Polygon%Limits%Top    = maxval(Polygon%VerticesF%Y)

    end subroutine SetLimitsPolygon
    
    !--------------------------------------------------------------------------

    subroutine SetLimitsXYZ(XYZ)

        !Arguments-------------------------------------------------------------
        type (T_XYZPoints), pointer                 :: XYZ
        
        !Begin-----------------------------------------------------------------

        XYZ%Limits%Left     = minval(XYZ%X)
        XYZ%Limits%Right    = maxval(XYZ%X)
        XYZ%Limits%Bottom   = minval(XYZ%Y)
        XYZ%Limits%Top      = maxval(XYZ%Y)

    end subroutine SetLimitsXYZ

    !--------------------------------------------------------------------------
    
    logical function IsVisible(Polygons, Point)
        
        !Arguments-------------------------------------------------------------
        type (T_Polygon), pointer                   :: Polygons
        type (T_PointF),   pointer                  :: Point
        
        !Local-----------------------------------------------------------------
        type (T_Polygon), pointer                   :: CurrentPolygon

        !Begin-----------------------------------------------------------------

        IsVisible = .false.

        CurrentPolygon => Polygons

        do while(associated(CurrentPolygon))

            if(IsPointInsidePolygon(Point, CurrentPolygon))then
                IsVisible = .true.
                return
            end if

            CurrentPolygon => CurrentPolygon%Next
        enddo

    end function IsVisible
    !--------------------------------------------------------------------------
    
    
    
    !--------------------------------------------------------------------------
    !Description    :   The first action performed is an acceleration test. if 
    !                   point is not inside maximum limits of the polygon
    !                   then it is not inside the polygon. If point is inside
    !                   the maximum bounds of the polygon then there is a 
    !                   possibility of being inside the polygon. Thus, it is 
    !                   drawn a "semi-recta" on the X axis to the right of the
    !                   point. If the "semi-recta" intersects the polygon an 
    !                   odd ("ímpar") number of times then the point is 
    !                   inside the polygon. The intersection test is performed 
    !                   for every segment of the polygon. If the point belongs 
    !                   is a vertix of the polygon or belongs to one of the 
    !                   segments that defines the polygon then it is considered 
    !                   to be inside the polygon. There are a few cases in which
    !                   special care needs to be taken, regarding the intersection
    !                   counting. If the intersection point coincides with one of 
    !                   the vertices then one intersection is counted only if it is
    !                   the vertice with the higher Y value. This is the way to 
    !                   count only one intersection once one polygon vertix belongs
    !                   to two segments and otherwise it would be tested twice.
    !                   Another case is when the segment belongs to the "semi-recta".
    !                   In this case no intersection is counted. This in the next 
    !                   segment to be tested the "semi-recta" will intersect the first
    !                   vertice of the segment. If the segments slope is negative then
    !                   the vertice has the higher Y value and one intersection is 
    !                   counted. 
     
    logical function IsPointInsidePolygon(Point, Polygon)

        !Arguments-------------------------------------------------------------
        type (T_PointF),   pointer                  :: Point
        type (T_Polygon), pointer                   :: Polygon

        !Local-----------------------------------------------------------------
        integer                                     :: CurrentVertix
        integer                                     :: NumberOfIntersections
        type(T_Segment),  pointer                   :: Segment
        real                                        :: Slope, IntersectionPointX
        real                                        :: HigherY

        !Begin-----------------------------------------------------------------

        IsPointInsidePolygon  = .false.
        NumberOfIntersections = 0

        !Acceleration test. If point is not inside maximum limits of the polygon
        !then it is not inside the polygon
        if(Point%X < Polygon%Limits%Left .or. Point%X > Polygon%Limits%Right .or.        &
           Point%Y > Polygon%Limits%Top  .or. Point%Y < Polygon%Limits%Bottom)return
        
        allocate(Segment)

do1:    do CurrentVertix = 1, Polygon%Count - 1
            
            !construct segment
            Segment%StartAt%X = Polygon%VerticesF(CurrentVertix)%X
            Segment%StartAt%Y = Polygon%VerticesF(CurrentVertix)%Y
            Segment%EndAt%X   = Polygon%VerticesF(CurrentVertix + 1)%X
            Segment%EndAt%Y   = Polygon%VerticesF(CurrentVertix + 1)%Y
            
            !if point coincides with one of the segments vertices then
            !it in inside the polygon
cd1:        if((Point%X .eq. Segment%StartAt%X .and. Point%Y .eq. Segment%StartAt%Y).or. &
               (Point%X .eq. Segment%EndAt%X   .and. Point%Y .eq. Segment%EndAt%Y  ))then

                IsPointInsidePolygon = .true.
                return
            
            !if point is placed between the segment vertices in the Y axis
            elseif((Point%Y <= Segment%StartAt%Y .and. Point%Y >= Segment%EndAt%Y) .or.  &
                   (Point%Y >= Segment%StartAt%Y .and. Point%Y <= Segment%EndAt%Y))then cd1

                !if segment has a slope
cd3:            if(Segment%StartAt%Y .ne. Segment%EndAt%Y)then
                    
                    !compute slope
                    Slope = (Segment%StartAt%X - Segment%EndAt%X)  /                     &
                            (Segment%StartAt%Y - Segment%EndAt%Y)
                    
                    !compute intersection point X coordinate of the "semi-recta"
                    !with the segment 
                    IntersectionPointX = Segment%StartAt%X - Slope *                     &
                                        (Segment%StartAt%Y - Point%Y)
                    
                    !if point belongs to the segment then it is inside the polygon
cd4:                if(IntersectionPointX .eq. Point%X)then

                        IsPointInsidePolygon = .true.
                        return

                    elseif(IntersectionPointX > Point%X)then cd4
                        
                        !if the intersection point coincides with one of the
                        !vertices then one intersection is counted only if it is
                        !the vertice with the higher Y value 
cd5:                    if(Point%Y .eq. Segment%StartAt%Y .or.                           &
                           Point%Y .eq. Segment%EndAt%Y)then
                            
                            !find higher segment ends higher Y coordinate
                            HigherY = max(Segment%StartAt%Y, Segment%EndAt%Y)

cd6:                        if(Point%Y .eq. HigherY)then
                                NumberOfIntersections = NumberOfIntersections + 1
                            end if cd6

                        else cd5
                            !if the intersection point is placed to the right 
                            !of the point or if point belongs to the segment 
                            !then one intersection is counted
                            NumberOfIntersections = NumberOfIntersections + 1
                        end if cd5
                        
                    end if cd4
                
                elseif(Segment%StartAt%Y .eq. Segment%EndAt%Y)then cd3

cd7:                if(Point%Y .eq. Segment%StartAt%Y)then
                        
                        if(Segment%StartAt%X > Segment%EndAt%X)then

                            if(Point%X <= Segment%StartAt%X .and. Point%X >= Segment%EndAt%X)then
                                !if point belongs to the segment then it is inside the polygon
                                IsPointInsidePolygon = .true.
                                return
                            end if

                        elseif(Segment%StartAt%X < Segment%EndAt%X)then
                            
                            if(Point%X >= Segment%StartAt%X .and. Point%X <= Segment%EndAt%X)then
                                !if point belongs to the segment then it is inside the polygon
                                IsPointInsidePolygon = .true.
                                return
                            end if


                        elseif(Segment%StartAt%X == Segment%EndAt%X)then

                            if(Point%X == Segment%StartAt%X)then
                                !if point belongs to the segment then it is inside the polygon
                                IsPointInsidePolygon = .true.
                                return
                            end if

                        end if    
                                           
                        end if cd7

                end if cd3
                
            end if cd1
        enddo do1

        deallocate(Segment)

        !if number of intersections is odd (odd = ímpar) then
        !point is inside the polygon
        IsPointInsidePolygon = IsOdd(NumberOfIntersections)

    end function IsPointInsidePolygon
    
    !--------------------------------------------------------------------------

     


    logical function IsPointInsideCircle(Point, CircleOrigin, Radius)

        !Arguments-------------------------------------------------------------
        type (T_PointF),   pointer                  :: Point, CircleOrigin
        real                                        :: Radius

        !Local-----------------------------------------------------------------
        !type(T_Segment),  pointer                   :: Segment
        real                                        :: Distance

        !Begin-----------------------------------------------------------------

        IsPointInsideCircle  = .false.

        !allocate(Segment)

        !Segment%StartAt%X = CircleOrigin%X
        !Segment%StartAt%Y = CircleOrigin%Y
        !Segment%EndAt%X   = Point%X
        !Segment%EndAt%Y   = Point%Y

        Distance = sqrt((CircleOrigin%X-Point%X)**2.+(CircleOrigin%Y-Point%Y)**2.)

        if(Distance <= Radius)IsPointInsideCircle = .true.

        !deallocate(Segment)



    end function IsPointInsideCircle

    
    !--------------------------------------------------------------------------
    
    
    subroutine WriteItemPolygon(Polygon, FilePath)

        !Arguments-------------------------------------------------------------
        type (T_Polygon), pointer           :: Polygon
        character(len=*), intent(IN)        :: FilePath
        
        !Local-----------------------------------------------------------------
        integer                             :: UnitID, STAT_CALL, CurrentVertix

        !Begin-----------------------------------------------------------------
        
        call UnitsManager(UnitID, OPEN_FILE, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'WriteItemPolygon - ModuleDrawing - ERR01'

    

        open(unit=UnitID, status = "unknown", file = FilePath)
        rewind(UnitID)


        write(UnitID, *)"<begin_polygon>"

        do CurrentVertix = 1, Polygon%Count
            write(UnitID, *)Polygon%VerticesF(CurrentVertix)%X, Polygon%VerticesF(CurrentVertix)%Y
        enddo
        write(UnitID, *)"<end_polygon>"


        call UnitsManager(UnitID, CLOSE_FILE, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'WriteItemPolygon - ModuleDrawing - ERR02'


    end subroutine WriteItemPolygon
    !--------------------------------------------------------------------------
    

    !--------------------------------------------------------------------------
    subroutine WriteItemXYZ(XYZ, FilePath)

        !Arguments-------------------------------------------------------------
        type (T_XYZPoints), pointer         :: XYZ
        character(len=*),   intent(IN)      :: FilePath
        
        !Local-----------------------------------------------------------------
        integer                             :: UnitID, STAT_CALL, CurrentPoint

        !Begin-----------------------------------------------------------------
        
        call UnitsManager(UnitID, OPEN_FILE, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'WriteItemXYZ - ModuleDrawing - ERR01'

        open(unit=UnitID, status = "unknown", file = FilePath)
        rewind(UnitID)

        write(UnitID, *)"<begin_xyz>"

        do CurrentPoint = 1, XYZ%Count
            write(UnitID, *)XYZ%X(CurrentPoint), XYZ%Y(CurrentPoint), XYZ%Z(CurrentPoint)
        enddo
        
        write(UnitID, *)"<end_xyz>"


        call UnitsManager(UnitID, CLOSE_FILE, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'WriteItemXYZ - ModuleDrawing - ERR02'


    end subroutine WriteItemXYZ

    !--------------------------------------------------------------------------

    !--------------------------------------------------------------------------
    subroutine WriteItemXYZ_v2(XYZArray, TotalXYZ, FilePath)

        !Arguments-------------------------------------------------------------
        real, dimension(:,:), pointer                :: XYZArray
        integer,                         intent(IN)  :: TotalXYZ
        character(len=*),                intent(IN)  :: FilePath

        
        !Local-----------------------------------------------------------------
        integer                             :: UnitID, STAT_CALL, CurrentPoint

        !Begin-----------------------------------------------------------------
        
        call UnitsManager(UnitID, OPEN_FILE, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'WriteItemXYZ - ModuleDrawing - ERR01'

        open(unit=UnitID, status = "unknown", file = FilePath)
        rewind(UnitID)

        write(UnitID, *)"<begin_xyz>"

        do CurrentPoint = 1, TotalXYZ
            write(UnitID, *)XYZArray(CurrentPoint,1), XYZArray(CurrentPoint,2), XYZArray(CurrentPoint,3)
        enddo
        
        write(UnitID, *)"<end_xyz>"


        call UnitsManager(UnitID, CLOSE_FILE, STAT = STAT_CALL)
        if(STAT_CALL .ne. SUCCESS_)stop 'WriteItemXYZ - ModuleDrawing - ERR02'


    end subroutine WriteItemXYZ_v2

    !--------------------------------------------------------------------------


        

end module ModuleDrawing
