!
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2018 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS is distributed in the hope that it will be useful, but 
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
!
! File Contents:
! MODULE PUL_parameter_module
!   CONTAINS
!     SUBROUTINE allocate_and_reset_PUL_data
!     SUBROUTINE deallocate_PUL_data
!     FUNCTION point_is_inside
!     SUBROUTINE Z_Y_from_L_C
!     SUBROUTINE write_Dshape_gmsh
!
! plus the following in include files:
!Laplace.F90:        SUBROUTINE Laplace
!PUL_analytic.F90:   SUBROUTINE PUL_LC_calc_wide_separation_approximation
!PUL_analytic.F90:   SUBROUTINE PUL_LC_calc_overshield_wide_separation_approximation
!PUL_analytic.F90:   SUBROUTINE calculate_height_over_ground_plane
!PUL_LC_Laplace.F90: SUBROUTINE PUL_LC_Laplace
!Aprod.F90:          SUBROUTINE Aprod
!Aprod.F90:          SUBROUTINE ATprod
!Aprod.F90:          SUBROUTINE zAprod
!Aprod.F90:          SUBROUTINE zATprod
!CG_solve.F90:       SUBROUTINE solve_real_symm(n, b, x,tol)
!CG_solve.F90:       SUBROUTINE solve_complex_symm(n, b, x,tol)
!
!
!
! NAME
!     MODULE PUL_parameters
!
!     Data relating to the calculation of PUL_parameters
!     
! COMMENTS
!    The conjugate gradient solutions are included in this module rather than the maths module
!    as they are intimately linked with the matrix-vector product subroutine which in turn relies
!    on the sparse matrix storage structure in this module
!
! HISTORY
!    started 25/11/2015 CJS
!    30/3/2016 CJS fix factor of 2 in mutual inductance of wires over ground plane
!    18_4_2016 CJS include calculation for conductors within a cylindrical shield for oversheild domains
!    5_7_2016 CJS include numerical Laplace solver for L,C,G in inhomogeneous regions
!    13/3/2018CJS include iterative solver stuff based on Stanford University symmlq.f
!    16/3/2018 CJS add y offset for ML_flex_cable (PUL%ox, and PUL%oy)(PUL%rdox, and PUL%rdoy)
!    19/6/2018 CJS include iterative sparse matrix solver based on the biconjugate gradient method
!
MODULE PUL_parameter_module

USE type_specifications
USE filter_module
USE frequency_spec
USE general_module

IMPLICIT NONE

! Main structure used to hold data for the per-unit-length parameter calculation
TYPE::PUL_type

  integer                     :: n_conductors    ! number of conductors
  
! conductor based arrays:
  integer,allocatable         :: shape(:)        ! holds a nuber which indicates the shape of the conductor
  logical,allocatable         :: conductor_is_shield(:) ! flag to inducate a shield conductor for FastHenry2
  real(dp),allocatable        :: sigma(:)        ! electrical conductivity of a conductor
  real(dp),allocatable        :: r(:)            ! radius of a circular conductor
  real(dp),allocatable        :: x(:)            ! x coordinate of the centre of the cable to which this conductor belongs
  real(dp),allocatable        :: y(:)            ! y coordinate of the centre of the cable to which this conductor belongs
  real(dp),allocatable        :: ox(:)           ! offset in the x direction of this conductor from the cable centre (x():),y(:) above) 
  real(dp),allocatable        :: oy(:)           ! offset in the y direction of this conductor from the cable centre (x():),y(:) above) 
  real(dp),allocatable        :: rh(:)           ! height (y dimension) of rectangular conductor
  real(dp),allocatable        :: rw(:)           ! width1  (x dimension) of rectangular conductor/ Dshape
  real(dp),allocatable        :: rw2(:)          ! width2  (x dimension) of rectangular conductor/ Dshape
  real(dp),allocatable        :: rd(:)           ! radius of dielectric surrounding a circular conductor
  real(dp),allocatable        :: rdh(:)          ! height (y dimension) of rectangular dielectric around conductor
  real(dp),allocatable        :: rdw(:)          ! width  (x dimension) of rectangular dielectric around conductor
  real(dp),allocatable        :: rdox(:)          ! offset of dielectric in the x direction of this conductor from the cable centre
  real(dp),allocatable        :: rdoy(:)          ! offset of dielectric in the y direction of this conductor from the cable centre
  real(dp),allocatable        :: rtheta(:)       ! rotation angle of conductor
  type(Sfilter),allocatable   :: epsr(:)         ! relative permittivity of dielecrtric surrounding the conductor
  
  logical             :: ground_plane_present    ! flag to indicate the presence of a ground plane
  real(dp)            :: ground_plane_angle      ! angle of ground plane normal from the x axis
  real(dp)            :: ground_plane_offset     ! ground plane offset
  real(dp)            :: ground_plane_sigma      ! ground plane conductivity
  real(dp)            :: ground_plane_w          ! ground plane width
  real(dp)            :: ground_plane_h          ! ground plane height
  real(dp)            :: ground_plane_Rdc        ! ground plane dc resistance
  integer             :: ground_plane_nsegx         ! ground plane number of segments in x
  integer             :: ground_plane_nsegz         ! ground plane number of segments in z
  integer             :: ground_plane_nh         ! ground plane number of layers in h
    
  logical             :: overshield_present      ! flag to indicate the presence of an overshield
  integer             :: overshield_shape        ! holds a nuber which indicates the shape of the opvershield 
  real(dp)            :: overshield_x            ! x coordinate of the centre of the overshield
  real(dp)            :: overshield_y            ! y coordinate of the centre of the overshield 
  real(dp)            :: overshield_r            ! radius of a circular overshield
  real(dp)            :: overshield_h            ! height (y dimension) of Dshape overshield
  real(dp)            :: overshield_w            ! width1 (x) parameter of Dshape overshield
  real(dp)            :: overshield_w2           ! width2 (x) parameter of Dshape overshield
  
  type(Sfilter)       :: epsr_background         ! relative permittivity of background dielecrtric (unually air)
    
  type(matrix)        :: L     ! inductance matrix
  type(matrix)        :: C     ! capacitance matrix
  type(matrix)        :: G     ! conductance matrix
  
  type(Sfilter_matrix):: Zfilter   ! frequency dependent impedance matrix
  type(Sfilter_matrix):: Yfilter   ! frequency dependent impedance matrix
  
  integer                       :: filter_fit_order  ! order for filter fitting for frequency dependent models
  type(frequency_specification) :: filter_fit_freq   ! frequency range for filter fitting for frequency dependent models
  
END TYPE PUL_type 

integer,public :: n_entry

! 1D arrays used in the construction of the K matrix ( K(i_K(:),j_K(:))=K(i_K(:),j_K(:))+s_K(:) )
integer,allocatable,public         :: i_K(:)
integer,allocatable,public         :: j_K(:)
real(dp),allocatable,public        :: s_K_re(:)
complex(dp),allocatable,public     :: s_K(:)


CONTAINS 

include "PUL_analytic.F90"

include "PUL_LC_Laplace.F90"

include "Laplace.F90"

include "PUL_RL_FastHenry2.F90"

include "Aprod.F90"

include "CG_solve.F90"

! NAME
!     SUBROUTINE allocate_and_reset_PUL_data
!
!     allocate and reset per-unit-length calculation data
!     
! COMMENTS
!     
!
! HISTORY
!    started 10/10/2015 CJS
!

  SUBROUTINE allocate_and_reset_PUL_data(PUL,n_conductors)

  USE type_specifications
  USE general_module

  IMPLICIT NONE

! variables passed to subroutine

  type(PUL_type),intent(OUT)   :: PUL
  integer,intent(IN)           :: n_conductors
  
! local variables

  integer :: i

! START

! allocate memory for PUL structure
  PUL%n_conductors=n_conductors
  ALLOCATE( PUL%shape( PUL%n_conductors) ) 
  ALLOCATE( PUL%conductor_is_shield( PUL%n_conductors) ) 
  ALLOCATE( PUL%sigma( PUL%n_conductors) ) 
  ALLOCATE( PUL%x( PUL%n_conductors) ) 
  ALLOCATE( PUL%y( PUL%n_conductors) ) 
  ALLOCATE( PUL%r( PUL%n_conductors) ) 
  ALLOCATE( PUL%ox( PUL%n_conductors) ) 
  ALLOCATE( PUL%oy( PUL%n_conductors) ) 
  ALLOCATE( PUL%rd( PUL%n_conductors) ) 
  ALLOCATE( PUL%rdw( PUL%n_conductors) ) 
  ALLOCATE( PUL%rdh( PUL%n_conductors) ) 
  ALLOCATE( PUL%rdox( PUL%n_conductors) ) 
  ALLOCATE( PUL%rdoy( PUL%n_conductors) ) 
  ALLOCATE( PUL%epsr( PUL%n_conductors) )     
  ALLOCATE( PUL%rh( PUL%n_conductors) ) 
  ALLOCATE( PUL%rw( PUL%n_conductors) ) 
  ALLOCATE( PUL%rw2( PUL%n_conductors) ) 
  ALLOCATE( PUL%rtheta( PUL%n_conductors) ) 
  
! reset all the elements of the PUL structure
  do i=1,PUL%n_conductors
  
    PUL%shape(i)=0
    PUL%conductor_is_shield(i)=.FALSE.
    PUL%sigma(i)=0d0 
    PUL%x(i)=0d0 
    PUL%y(i)=0d0 
    PUL%r(i)=0d0 
    PUL%ox(i)=0d0 
    PUL%oy(i)=0d0 
    PUL%rh(i)=0d0 
    PUL%rw(i)=0d0 
    PUL%rw2(i)=0d0 
    PUL%rd(i)=0d0 
    PUL%rdw(i)=0d0 
    PUL%rdh(i)=0d0 
    PUL%rdox(i)=0d0 
    PUL%rdoy(i)=0d0 
    PUL%rtheta(i)=0d0
    PUL%epsr(i)=1d0

  end do
  
  PUL%ground_plane_present=.FALSE.
  PUL%ground_plane_angle=0d0
  PUL%ground_plane_offset=0d0
  PUL%ground_plane_sigma=0d0
  PUL%ground_plane_w=0d0
  PUL%ground_plane_h=0d0
  PUL%ground_plane_Rdc=0d0
  PUL%ground_plane_nsegx=0         
  PUL%ground_plane_nsegz=0       
  PUL%ground_plane_nh=1         

  PUL%overshield_present=.FALSE.
  PUL%overshield_shape=circle    ! circle by default
  PUL%overshield_x=0d0
  PUL%overshield_y=0d0
  PUL%overshield_r=0d0
  PUL%overshield_w=0d0
  PUL%overshield_w2=0d0
  PUL%overshield_h=0d0
  
  PUL%epsr_background=1d0

  RETURN

  END SUBROUTINE allocate_and_reset_PUL_data

! NAME
!     SUBROUTINE deallocate_PUL_data
!
!     deallocate per-unit-length calculation data
!     
! COMMENTS
!     
!
! HISTORY
!    started 2/12/2015 CJS
!    5/6/2016 CJS include additional variables for numerical Laplace solution.
!

  SUBROUTINE deallocate_PUL_data(PUL)

  USE type_specifications
  USE general_module

  IMPLICIT NONE

! variables passed to subroutine

  type(PUL_type),intent(INOUT)    :: PUL
  
! local variables

  integer :: i

! START

  if (allocated(PUL%shape )    ) DEALLOCATE( PUL%shape )
  if (allocated(PUL%conductor_is_shield )    ) DEALLOCATE( PUL%conductor_is_shield )
  if (allocated(PUL%sigma )    ) DEALLOCATE( PUL%sigma )
  if (allocated(PUL%x )    ) DEALLOCATE( PUL%x )
  if (allocated(PUL%y )    ) DEALLOCATE( PUL%y )
  if (allocated(PUL%r )    ) DEALLOCATE( PUL%r )
  if (allocated(PUL%ox )    ) DEALLOCATE( PUL%ox )
  if (allocated(PUL%oy )    ) DEALLOCATE( PUL%oy )
  if (allocated(PUL%rh )    ) DEALLOCATE( PUL%rh )
  if (allocated(PUL%rw )    ) DEALLOCATE( PUL%rw )
  if (allocated(PUL%rw2 )    ) DEALLOCATE( PUL%rw2 )
  if (allocated(PUL%rd )   ) DEALLOCATE( PUL%rd )
  if (allocated(PUL%rdw )   ) DEALLOCATE( PUL%rdw )
  if (allocated(PUL%rdh )   ) DEALLOCATE( PUL%rdh )
  if (allocated(PUL%rdox )   ) DEALLOCATE( PUL%rdox )
  if (allocated(PUL%rdoy )   ) DEALLOCATE( PUL%rdoy )
  if (allocated(PUL%rtheta )    ) DEALLOCATE( PUL%rtheta )
  
  if (allocated(PUL%epsr ) ) then
    do i=1,PUL%n_conductors
      CALL deallocate_Sfilter(PUL%epsr(i))
    end do
    DEALLOCATE( PUL%epsr )
  end if
  
  CALL deallocate_Sfilter(PUL%epsr_background)
  
  if (allocated(PUL%L%mat) ) DEALLOCATE( PUL%L%mat )
  if (allocated(PUL%C%mat) ) DEALLOCATE( PUL%C%mat )
  if (allocated(PUL%G%mat) ) DEALLOCATE( PUL%G%mat )
  
  CALL deallocate_Sfilter_matrix( PUL%Zfilter )
  CALL deallocate_Sfilter_matrix( PUL%Yfilter )

  RETURN

  END SUBROUTINE deallocate_PUL_data
 
! NAME
!     FUNCTION point_is_inside(xt,yt,shape_type,x,y,r,rh,rw,rtheta)
!
!     function to determine whether conductors are inside dielectrics - used to create region boundary
!     lists for mesh generation 
!     
! COMMENTS
!     
!
! HISTORY
!    started 6/10/2016 CJS
!    add rox, roy instead of only the x offset, ro 16/3/2018 CJS
!
  
  logical FUNCTION point_is_inside(xt,yt,shape_type,x,y,r,rh,rw,rox,roy,rtheta)
  
  real(dp),intent(IN) ::   xt           ! test point x
  real(dp),intent(IN) ::   yt           ! test point y
  integer,intent(IN)  ::   shape_type   ! holds a nuber which indicates the shape of the test object
  real(dp),intent(IN) ::   x            ! x coordinate of the centre of the test object
  real(dp),intent(IN) ::   y            ! y coordinate of the centre of the test object
  real(dp),intent(IN) ::   r            ! radius of cylindrical test object
  real(dp),intent(IN) ::   rh           ! height (y dimension) of rectangular test object
  real(dp),intent(IN) ::   rw           ! width  (x dimension) of  test object
  real(dp),intent(IN) ::   rox          ! offset in the x direction of this test object from the centre
  real(dp),intent(IN) ::   roy          ! offset in the y direction of this test object from the centre
  real(dp),intent(IN) ::   rtheta       ! rotation angle of test object
 
! local variables

  real(dp)::   d
  real(dp)::   xt_r,yt_r
  real(dp)::   xt_ro,yt_ro
  real(dp)::   xt_o,yt_o
  real(dp)::   dx,dy

! START
  
  if (shape_type.EQ.circle) then

! calculate the distance between the test point and the centre of the dielectric cylinder  
    d=sqrt((x-xt)**2+(y-yt)**2)
    
    if (d.LT.r) then
      point_is_inside=.TRUE.
    else
      point_is_inside=.FALSE.    
    end if
  
  else
! rectangular conductor

    write(*,*)'Test point',xt,yt
    
! apply the inverse of the offset to the shape i.e. move it into the shifted, un-rotated coordinate system of the rectangle
    xt_o=xt-x
    yt_o=yt-y
    
    write(*,*)'Offset on test point',xt_o,yt_o

! apply the inverse of the rotation to the test point i.e. move it into the un-rotated coordinate system of the rectangle
    xt_r= xt_o*cos(rtheta)+yt_o*sin(rtheta)
    yt_r=-xt_o*sin(rtheta)+yt_o*cos(rtheta)

    write(*,*)'Inverse rotation on test point',xt_r,yt_r
    
! apply the inverse of the offset to the shape i.e. move it into the shifted, un-rotated coordinate system of the rectangle
    xt_ro=xt_r-rox
    yt_ro=yt_r-roy

    write(*,*)'Offset Inverse rotation on test point',xt_ro,yt_ro
    
! xt_ro,yt_ro is the the distance from the transformed test point to the centre of the rectangle in x and y   
    
    dx=xt_ro
    dy=yt_ro
    
!    write(*,*)'dielectric centre point',x,y,' offset',rox,roy,' angle',rtheta
!    write(*,*)'dx',dx,' rw/2',rw/2d0
!    write(*,*)'dy',dy,' rh/2',rh/2d0
    
    if ( (abs(dx).LT.rw/2d0).AND.(abs(dy).LT.rh/2d0) ) then
      point_is_inside=.TRUE.
      write(*,*)'POINT IS INSIDE'
    else
      point_is_inside=.FALSE.    
      write(*,*)'POINT IS OUTSIDE'
    end if
   
  end if
  
  
  END FUNCTION point_is_inside
!
! NAME
!     SUBROUTINE Z_Y_from_L_C(unit)
!
!     calculate frequency dependent impedance and admittance (filter function) matrices 
!     for inductance and capacitance matrices with frequency independent properties
!     Z=jwL  Y=jwC
!
! COMMENTS
!     
!
! HISTORY
!    started 2/12/2015 CJS
!

  SUBROUTINE Z_Y_from_L_C(L,C,Z,Y)

  USE type_specifications
  USE general_module
  USE maths
  USE filter_module

  IMPLICIT NONE

! variables passed to subroutine

  type(matrix),intent(IN)    :: L
  type(matrix),intent(IN)    :: C
  
  type(Sfilter_matrix),intent(OUT)    :: Z
  type(Sfilter_matrix),intent(OUT)    :: Y

! local variables

  integer :: dim
  integer :: row,col
  type(Sfilter) :: jw

! START

  jw=jwA_filter(1d0)

! Impedance matrix
  dim=L%dim
  Z%dim=dim
  if (.NOT.ALLOCATED(Z%sfilter_mat)) ALLOCATE(Z%sfilter_mat(dim,dim))
  do row=1,dim
    do col=1,dim
      Z%sfilter_mat(row,col)=L%mat(row,col)*jw
    end do
  end do
  
! Capacitance matrix
  dim=C%dim
  Y%dim=dim
  if (.NOT.ALLOCATED(Y%sfilter_mat)) ALLOCATE(Y%sfilter_mat(dim,dim))
  do row=1,dim
    do col=1,dim
      Y%sfilter_mat(row,col)=C%mat(row,col)*jw
    end do
  end do
  
! Deallocate work filter
  CALL deallocate_Sfilter(jw)
  
  RETURN
  
  END SUBROUTINE Z_Y_from_L_C
!
! NAME
!     SUBROUTINE write_Dshape_gmsh
!
!     write the points required for a Dshape connector shell model in gmsh
!     
! COMMENTS
!     
!
! HISTORY
!    started 15/11/2016 CJS
!

  SUBROUTINE write_Dshape_gmsh(x,y,w_in,w2_in,h_in,r,ox,oy,theta,dl,point,number,unit)

  USE type_specifications
  USE general_module

  IMPLICIT NONE

! variables passed to subroutine
  
  real(dp),intent(IN) :: x      ! x coordinate of the centre of the Dshape
  real(dp),intent(IN) :: y      ! y coordinate of the centre of the Dshape
  real(dp),intent(IN) :: w_in   ! width1 (x dimension) of the Dshape
  real(dp),intent(IN) :: w2_in  ! width2 (x dimension) of the Dshape
  real(dp),intent(IN) :: h_in   ! height (x dimension) of the Dshape
  real(dp),intent(IN) :: r      ! radius of curves in Dshape
  real(dp),intent(IN) :: ox     ! x offset
  real(dp),intent(IN) :: oy     ! y offset
  real(dp),intent(IN) :: theta  ! rotation angle of Dshape
  real(dp),intent(IN) :: dl     ! edge length for the line segments making up the Dshape

  integer,intent(INOUT)        :: point  ! point counter
  integer,intent(IN)           :: number ! number of the Dshape - used as a label in the gmsh file
  integer,intent(IN)           :: unit   ! unit to write to
  
! local variables

  real(dp) :: w1,w2,h
  
  real(dp) :: vx,vy
  real(dp) :: voxl,voyl
  real(dp) :: norm

  real(dp) :: xt,yt
  real(dp) :: xp,yp

  integer :: i

! START

      w1=w_in/2d0-r
      w2=w2_in/2d0-r
      h=h_in/2d0-r

! vector from top left conductor to bottem left conductor      
      vx=w1-w2
      vy=-2d0*h

! perpendicular vector to -x edge
      norm=sqrt(vx*vx+vy*vy)
      voxl=vy*r/norm
      voyl=-vx*r/norm

      write(unit,*)' // Dshape ',number
      point=point+1

! POINT 1      ! top right
      xt=w1+ox
      yt=h+r+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 2      ! top left
      xt=-w1+ox
      yt=h+r+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 3      ! top left centre
      xt=-w1+ox
      yt=h+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 4      ! top left edge
      xt=-w1+ox+voxl
      yt=h+voyl+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 5      ! bottom left edge
      xt=-w2+ox+voxl
      yt=-h+voyl+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 6      ! bottom left centre
      xt=-w2+ox
      yt=-h+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
       
! POINT 7      ! bottom left
      xt=-w2+ox
      yt=-(h+r)+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
       
! POINT 8      ! bottom right
      xt=w2+ox
      yt=-(h+r)+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 9      ! bottom right centre
      xt=w2+ox
      yt=-h+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 10      ! bottom right edge
      xt=w2+ox-voxl
      yt=-h+voyl+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 11      ! top right edge
      xt=w1+ox-voxl
      yt=h+voyl+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
      point=point+1
      
! POINT 12      ! top right centre
      xt=w1+ox
      yt=h+oy
      xp=x+xt*cos(theta)-yt*sin(theta)
      yp=y+xt*sin(theta)+yt*cos(theta)
      write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'

  RETURN

  END SUBROUTINE write_Dshape_gmsh
  
END MODULE PUL_parameter_module