!
! 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) 2023 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:
! PROGRAM write_FH_input_file
!
! NAME
!     write_FH_input_file
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!    This program creates an input file for FastHenry2 for the calculation
!    of the per-unit-length impedance matrix over a specified range of
!    frequencies.
!    
!
! COMMENTS
!     This software is still at the experimental stage. There are a number of
!     issues related the most efficient manner in which to run fasthenry to 
!     obtain accurate results over the specified frequency range.  These issues 
!     are mostly related to the manner in which ciruclar conductors are represented 
!     as a set of rectangles and subsequently decomposed into filaments to 
!     represent the current distribution. 
!
! HISTORY
!
!     started October 2023 CJS in proximity_effects branch
!
PROGRAM write_FH_input_file

USE type_specifications

IMPLICIT NONE

integer :: n_conductors

logical :: ground_plane

integer :: mesh_type

! Ground plane data

real(dp) :: gp_x(3),gp_y(3),gp_z(3)
  
real(dp)    :: gp_t,gp_sigma,gp_rh,gp_ex1,gp_ey1,gp_ez1,gp_ex2,gp_ey2,gp_ez2
real(dp)    :: gp_skin_depth
integer :: gp_nhinc,gp_seg1,gp_seg2

real(dp) :: gp_xmin,gp_xmax
real(dp) :: gp_xc,gp_yc,gp_w,gp_h

character(len=80),allocatable    :: gp_node1
character(len=80),allocatable    :: gp_node2

character(len=80),allocatable    :: gp_node1_external
character(len=80),allocatable    :: gp_node2_external

! Conductor data

TYPE::conductor_type

  integer :: type                    ! conductor shape
  integer :: mesh_type               ! way to decompose circular conductors into rectangles
  integer :: mesh_to_layer_type

  real(dp)    :: rc                  ! circular conductor radius
  real(dp)    :: rci                 ! inner conductor radius (for annular conductor)
  real(dp)    :: rco                 ! outer conductor radius (for annular conductor)
  real(dp)    :: xc                  ! x position of conductor centre
  real(dp)    :: yc                  ! y position of conductor centre
  real(dp)    :: width               ! rectangular conductor width
  real(dp)    :: height              ! rectangular conductor height
  
  real(dp)    :: rss                 ! seven strand conductor radius
  real(dp)    :: rss_outer           ! seven strand conductor radius
  real(dp)    :: xcss(7)             ! x centre coordinate of the seven conductors
  real(dp)    :: ycss(7)             ! y centre coordinate of the seven conductors
  real(dp)    :: rot_angle
  integer     :: n_layers2ss
  
  integer :: n_layers2               ! number of rectangular layers in a circular conductor radius 
  integer :: tot_n_layers            ! total number of layers in this conductor 
  
  real(dp),allocatable    :: x(:)
  real(dp),allocatable    :: y(:)
  real(dp),allocatable    :: w(:)
  real(dp),allocatable    :: h(:)
  real(dp),allocatable    :: d(:)
  integer,allocatable :: anwinc(:)
  integer,allocatable :: anhinc(:)
  
  real(dp),allocatable    :: wx(:)
  real(dp),allocatable    :: wy(:)
  real(dp),allocatable    :: wz(:)
  
  real(dp)    :: sigma              ! electrical conductivity (Siemens/metre)
  real(dp)    :: dl
  integer :: nwinc                  ! fastherny parameter nwinc: number of filaments in rectangle width
  integer :: nhinc                  ! fastherny parameter nhinc: number of filaments in rectangle height
  real(dp)    :: rw                 ! fastherny parameter rw: filament size ratio in width
  real(dp)    :: rh                 ! fastherny parameter rh: filament size ratio in height

  character(len=80),allocatable    :: node1_list(:)
  character(len=80),allocatable    :: node2_list(:)
  
  real(dp)    :: grid_dim
  integer :: nxmin,nxmax,nymin,nymax
  integer,allocatable :: grid(:,:)
  real(dp),allocatable    :: depth(:,:)
  
  real(dp)    :: skin_depth
  logical :: auto_grid_density

END TYPE conductor_type

integer,parameter :: type_cyl=1
integer,parameter :: type_rect=2
integer,parameter :: type_annulus=3
integer,parameter :: type_gnd=4
integer,parameter :: type_seven_strand=5

integer,parameter :: mesh_type_layer=1   ! divide circles into uniform thickness layers 
integer,parameter :: mesh_type_grid=2    ! divide circles into a uniform grid of squares
integer,parameter :: mesh_type_shell=3

type(conductor_type),allocatable :: conductor_data(:)

real(dp)    :: fmin,fmax
real(dp)    :: ndec

real(dp)    :: rh,rw
integer :: nhinc,nwinc

real(dp) :: x,y,ymin,ymax,w,h,wx,hy
integer :: conductor,layer,layer_number,nc

real(dp) :: angle
real(dp) :: lrc,lxc,lyc,lw,lh
real(dp) :: vx,vy

real(dp) :: dout,din

integer :: nxmin,nxmax,nymin,nymax,ix,iy
real(dp)    :: dpt,apx,apy

integer :: cx1,cx2
real(dp)    :: cw,ch
logical :: in_conductor

real(dp)    :: required_dl,nsd
integer :: nfilaments

character(LEN=6) :: conductor_string
character(LEN=6) :: layer_string
character(LEN=80) :: node_string
character(LEN=80) :: segment_string
character(LEN=80) :: loop_string
character(LEN=80) :: line_string

character(LEN=12) :: x_string,y_string,z_string

integer :: line
character(LEN=80) :: FH2_filename
character :: type_ch

integer :: tot_n_segments,tot_n_filaments
integer :: gp_n_segments,gp_n_filaments

integer :: i

integer :: nsegments_recommended
integer :: nh_auto,nw_auto

logical :: found_square

real(dp) :: xpt,ypt

real(dp),parameter :: pi=3.1415926535

! START

! READ THE PROBLEM SPECIFICATION

line=1
write(*,*)'Enter the name of the FastHenry2 input file to write:'
read(*,'(A80)',ERR=9000)FH2_filename

write(*,*)"Enter the number of conductors or 'ground_plane'"
line=line+1
read(*,'(A80)',ERR=9000)line_string
type_ch=line_string(1:1)

if ( (type_ch.EQ.'g').OR.(type_ch.EQ.'G') ) then

  write(*,*)'Reading ground plane specification...'
  ground_plane=.TRUE.
  
  write(*,*)'Enter the ground plane point 1 coordinates, x y z in metres'
  line=line+1
  read(*,*,ERR=9000)gp_x(1),gp_y(1),gp_z(1)
  write(*,*)'Enter the ground plane point 2 coordinates, x y z in metres'
  line=line+1
  read(*,*,ERR=9000)gp_x(2),gp_y(2),gp_z(2)
  write(*,*)'Enter the ground plane point 3 coordinates, x y z in metres'
  line=line+1
  read(*,*,ERR=9000)gp_x(3),gp_y(3),gp_z(3)
  
  write(*,*)'Enter the ground plane thickness in metres'
  line=line+1
  read(*,*,ERR=9000)gp_t
  
  write(*,*)'Enter the ground plane discretisation in p1-p2 and p2=p3 directions'
  line=line+1
  read(*,*,ERR=9000)gp_seg1,gp_seg2
  
  write(*,*)'Enter the ground conductivity in Siemens/metre'
  line=line+1
  read(*,*,ERR=9000)gp_sigma
  
  write(*,*)'Enter the ground plane discretisation in thickness,  nhinc'
  line=line+1
  read(*,*,ERR=9000)gp_nhinc
  
  write(*,*)'Enter the ground plane discretisation in thickness ratio,  rh'
  line=line+1
  read(*,*,ERR=9000)gp_rh

  write(*,*)'Enter the ground plane end 1 node coordinates, x y z in metres'
  line=line+1
  read(*,*,ERR=9000)gp_ex1,gp_ey1,gp_ez1
  write(*,*)'Enter the ground plane end 2 node coordinates, x y z in metres'
  line=line+1
  read(*,*,ERR=9000)gp_ex2,gp_ey2,gp_ez2
  
  gp_node1='Ngp_e1'
  gp_node2='Ngp_e2'
  
  gp_node1_external='Ngp_e1_ext'
  gp_node2_external='Ngp_e2_ext'

  write(*,*)"Enter the number of conductors"
  line=line+1
  read(*,*,ERR=9000)n_conductors

else

  write(*,*)'No ground plane'
  ground_plane=.FALSE.
  read(line_string,*,ERR=9000)n_conductors
  
end if

write(*,*)'Number of conductors (excluding ground plane)=',n_conductors

ALLOCATE( conductor_data(n_conductors) )

do conductor=1,n_conductors

  write(*,*)'Enter the conductor number'
  line=line+1
  read(*,*,ERR=9000)nc
  
  if (nc.NE.conductor) GOTO 9010

  write(*,*)'Enter the conductor type (cylindrical rectangular or annulus)'
  line=line+1
  read(*,'(A80)',ERR=9000)line_string
  type_ch=line_string(1:1)
  
  conductor_data(conductor)%auto_grid_density=.FALSE.
  
  if ( (type_ch.EQ.'c').OR.(type_ch.EQ.'C') ) then
    conductor_data(conductor)%type=type_cyl    

! work out the grid type    

INCLUDE "WRITE_FH2_IPFILE/get_grid_type.F90"

    write(*,*)Conductor,conductor,' mesh type=',conductor_data(conductor)%mesh_type

    write(*,*)'Enter the cylindrical conductor centre coordinates, xc yc in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%xc,conductor_data(conductor)%yc
    
    write(*,*)'Enter the cylindrical conductor radius, rc in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%rc
    
    write(*,*)'Enter the cylindrical conductor discretisation, dl in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%dl
    
  else if ( (type_ch.EQ.'s').OR.(type_ch.EQ.'S') ) then
    conductor_data(conductor)%type=type_seven_strand   

! work out the grid type    

INCLUDE "WRITE_FH2_IPFILE/get_grid_type.F90"

    write(*,*)Conductor,conductor,' mesh type=',conductor_data(conductor)%mesh_type

    write(*,*)'Enter the seven strand conductor centre coordinates, xc yc in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%xc,conductor_data(conductor)%yc
    
    write(*,*)'Enter the seven strand conductor equivalent radius, rc in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%rc
    
    write(*,*)'Enter the seven strand conductor rotation angle in degrees'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%rot_angle
    conductor_data(conductor)%rot_angle=conductor_data(conductor)%rot_angle*pi/180d0
    
    write(*,*)'Enter the seven strand conductor discretisation, dl in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%dl
       
  else if ( (type_ch.EQ.'r').OR.(type_ch.EQ.'R') ) then
    conductor_data(conductor)%type=type_rect

    write(*,*)'Enter the rectangular conductor centre coordinates, xc yc in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%xc,conductor_data(conductor)%yc
    
    write(*,*)'Enter the rectangular conductor width (x dimension) height (y dimension) in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%width,conductor_data(conductor)%height
 
    conductor_data(conductor)%n_layers2=1
    conductor_data(conductor)%tot_n_layers=1
  
  else if ( (type_ch.EQ.'a').OR.(type_ch.EQ.'A') ) then
    conductor_data(conductor)%type=type_annulus
    
INCLUDE "WRITE_FH2_IPFILE/get_grid_type.F90"

    write(*,*)Conductor,conductor,' mesh type=',conductor_data(conductor)%mesh_type
     
    write(*,*)'Enter the cylindrical conductor centre coordinates, xc yc in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%xc,conductor_data(conductor)%yc
    
    write(*,*)'Enter the cylindrical inner conductor radius, rci in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%rci
    
    write(*,*)'Enter the cylindrical outer conductor radius, rco in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%rco
    
    write(*,*)'Enter the cylindrical conductor discretisation, dl in metres'
    line=line+1
    read(*,*,ERR=9000)conductor_data(conductor)%dl

  else
    GOTO 9020
  end if   ! Conductor type
      
  write(*,*)'Enter the conductor conductivity, sigma in Siemens/metre'
  line=line+1
  read(*,*,ERR=9000)conductor_data(conductor)%sigma
            
  write(*,*)'Enter the conductor discretisations, nwinc nhinc'
  line=line+1
  read(*,*,ERR=9000)conductor_data(conductor)%nwinc,conductor_data(conductor)%nhinc
            
  write(*,*)'Enter the conductor discretisation ratios, rw rh'
  line=line+1
  read(*,*,ERR=9000)conductor_data(conductor)%rw,conductor_data(conductor)%rh

end do ! read next conductor

! Read the frequency range data
     
write(*,*)'Enter the minimum frequency, fmin (Hz)'
line=line+1
read(*,*,ERR=9000)fmin
     
write(*,*)'Enter the maximum frequency, fmax (Hz)'
line=line+1
read(*,*,ERR=9000)fmax
     
write(*,*)'Enter the number of frequency samples per decade, ndec'
line=line+1
read(*,*,ERR=9000)ndec

! END OF PROBLEM SPECIFICATION

! FastHenry input file
open(unit=20,file=trim(FH2_filename))

! files for plotting the fasthenry 2 input geometry and decomposition into filaments
open(unit=10,file='cross_section.dat')
open(unit=12,file='grid.dat')

! write header
write(20,'(A)')'* conductors in free space'
write(20,'(A)')'*'
write(20,'(A)')'.units m'
write(20,'(A)')'*'

! WORK OUT THE SKIN DEPTH IN EACH CONDUCTOR

if (ground_plane) then
  gp_skin_depth=sqrt(1.0/(pi*fmax*4.0*pi*1e-7*gp_sigma))
  write(*,*)'Minimum ground plane skin depth =',gp_skin_depth
  CALL calc_nmin(gp_skin_depth,gp_t,gp_rh,nsegments_recommended)
end if

do conductor=1,n_conductors
  
  conductor_data(conductor)%skin_depth=sqrt(1.0/(pi*fmax*4.0*pi*1e-7*conductor_data(conductor)%sigma))
  write(*,*)'Conductor',conductor,' Minimum skin depth =',conductor_data(conductor)%skin_depth

end do ! next conductor

! DECOMPOSE CIRCULAR CONDUCTORS INTO RECTANGULAR LAYERS
! ALLOCATE GRIDS IN CONDUCTORS IF REQUIRED...

do conductor=1,n_conductors

  if (conductor_data(conductor)%type.EQ.type_cyl) then
  
     if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then
    
      conductor_data(conductor)%n_layers2=NINT(conductor_data(conductor)%rc/conductor_data(conductor)%dl)
      conductor_data(conductor)%tot_n_layers=2*conductor_data(conductor)%n_layers2
    
    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
    
! allocate grid for the circular geometry
      conductor_data(conductor)%grid_dim=conductor_data(conductor)%rc
     
INCLUDE "WRITE_FH2_IPFILE/create_grid.F90"
         
! Loop over the grid and set segments within the conductor
    
      conductor_data(conductor)%tot_n_layers=0
    
      do ix=nxmin,nxmax
        do iy=nymin,nymax
          dpt=conductor_data(conductor)%dl*sqrt(real(ix)**2+real(iy)**2)
          if ( dpt.LE.conductor_data(conductor)%rc ) then
            conductor_data(conductor)%grid(ix,iy)=1
            conductor_data(conductor)%depth(ix,iy)=conductor_data(conductor)%rc-dpt
            conductor_data(conductor)%tot_n_layers=conductor_data(conductor)%tot_n_layers+1
          end if
        end do
      end do
      
      write(*,*)'Set cylindrical grid, ncells=',conductor_data(conductor)%tot_n_layers
    
      conductor_data(conductor)%n_layers2=0

    else
      write(*,*)'Unknown grid type'
      STOP 1
    end if
  
  else if (conductor_data(conductor)%type.EQ.type_seven_strand) then

! radius of each strand for the same conductor area  
    conductor_data(conductor)%rss=conductor_data(conductor)%rc/sqrt(7.0)
! maximum radius of the combined conductors
    conductor_data(conductor)%rss_outer=conductor_data(conductor)%rss*3.0

! seven conductor centre coordinates
! central conductor
    conductor_data(conductor)%xcss(1)=conductor_data(conductor)%xc
    conductor_data(conductor)%ycss(1)=conductor_data(conductor)%yc
    do i=1,6
      angle=real(i-1)*2.0*pi/6.0+conductor_data(conductor)%rot_angle
      conductor_data(conductor)%xcss(i+1)=conductor_data(conductor)%xc+2d0*conductor_data(conductor)%rss*cos(angle)
      conductor_data(conductor)%ycss(i+1)=conductor_data(conductor)%yc+2d0*conductor_data(conductor)%rss*sin(angle)    
    end do
    
    if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then
    
      conductor_data(conductor)%n_layers2ss=NINT(conductor_data(conductor)%rss/conductor_data(conductor)%dl)
      conductor_data(conductor)%tot_n_layers=7*2*conductor_data(conductor)%n_layers2ss
    
    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
    
! allocate grid for the circular geometry
      conductor_data(conductor)%grid_dim=conductor_data(conductor)%rss_outer
     
INCLUDE "WRITE_FH2_IPFILE/create_grid.F90"
         
! Loop over the grid and set segments within the conductor
    
      conductor_data(conductor)%tot_n_layers=0

! loop over 7 strands      
      do i=1,7

! loop over grid    
        do ix=nxmin,nxmax
          do iy=nymin,nymax
          
! calculate dpt, the distance to the centre of strand i
            xpt=conductor_data(conductor)%dl*real(ix)- &
                (conductor_data(conductor)%xcss(i)-conductor_data(conductor)%xcss(1))
            ypt=conductor_data(conductor)%dl*real(iy)- &
                (conductor_data(conductor)%ycss(i)-conductor_data(conductor)%ycss(1))
            
            dpt=sqrt(xpt**2+ypt*2)
            
            if ( dpt.LE.conductor_data(conductor)%rss ) then
              conductor_data(conductor)%grid(ix,iy)=1
              conductor_data(conductor)%depth(ix,iy)=conductor_data(conductor)%rss-dpt
              conductor_data(conductor)%tot_n_layers=conductor_data(conductor)%tot_n_layers+1
            end if
          end do
        end do
    
        conductor_data(conductor)%n_layers2=0

      end do ! next strand

    else
      write(*,*)'Unknown grid type'
      STOP 1
    end if
 
  else if (conductor_data(conductor)%type.EQ.type_annulus) then
      
      
    if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then
    
      conductor_data(conductor)%n_layers2=0
      conductor_data(conductor)%tot_n_layers=20  ! number of segments in the loop
    
    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
    
! allocate grid for the annular geometry

      conductor_data(conductor)%grid_dim=conductor_data(conductor)%rco
    
INCLUDE "WRITE_FH2_IPFILE/create_grid.F90"
    
! Loop over the grid and set segments within the annular conductor
    
      conductor_data(conductor)%tot_n_layers=0
    
      do ix=nxmin,nxmax
        do iy=nymin,nymax
          dpt=conductor_data(conductor)%dl*sqrt(real(ix)**2+real(iy)**2)
          
          if ( (dpt.GE.conductor_data(conductor)%rci).AND.(dpt.LE.conductor_data(conductor)%rco) ) then
            conductor_data(conductor)%grid(ix,iy)=1
            dout=conductor_data(conductor)%rco-dpt
            din=dpt-conductor_data(conductor)%rci
            conductor_data(conductor)%depth(ix,iy)=min(dout,din)
            conductor_data(conductor)%tot_n_layers=conductor_data(conductor)%tot_n_layers+1
          end if
        end do
      end do
    
      conductor_data(conductor)%n_layers2=0

    else
      
      write(*,*)'We can only use the grid mesh type for an annulus'
      STOP 1
      
    end if ! mesh type grid

  end if  ! annular
  
end do ! next conductor

! WRITE GROUND PLANE INFORMATION IF REQUIRED

gp_n_segments=0
gp_n_filaments=0

if (ground_plane) then

  write(20,'(A)')'*'
  write(20,'(A)')'* Ground plane'
  write(20,'(A,ES12.4,A,ES12.4,A,ES12.4)')'ground_plane x1=',gp_x(1),' y1=',gp_y(1),' z1=',gp_z(1)
  write(20,'(A,ES12.4,A,ES12.4,A,ES12.4)')'+            x2=',gp_x(2),' y2=',gp_y(2),' z2=',gp_z(2)
  write(20,'(A,ES12.4,A,ES12.4,A,ES12.4)')'+            x3=',gp_x(3),' y3=',gp_y(3),' z3=',gp_z(3)
  write(20,'(A,ES12.4)')'+  thick=',gp_t
  write(20,'(A,ES12.4)')'+  sigma=',gp_sigma
  write(20,'(A,I4)')'+  nhinc=',gp_nhinc
  write(20,'(A,ES12.4)')'+  rh=',gp_rh
  write(20,'(A,I4,A,I4)')'+  seg1=',gp_seg1,' seg2=',gp_seg2
 
  gp_n_segments =(gp_seg1+1)*gp_seg2+(gp_seg2+1)*gp_seg1
  gp_n_filaments=gp_nhinc*gp_n_segments
 
  write(x_string,'(ES12.4)')gp_ex1
  write(y_string,'(ES12.4)')gp_ey1
  write(z_string,'(ES12.4)')gp_ez1  
  write(20,'(9A)')'+ ',trim(gp_node1),' (',trim(adjustl(x_string)),  &
                                       ',',trim(adjustl(y_string)),  &
                       ',',trim(adjustl(z_string)),')'
  
  write(x_string,'(ES12.4)')gp_ex2
  write(y_string,'(ES12.4)')gp_ey2
  write(z_string,'(ES12.4)')gp_ez2  
  write(20,'(9A)')'+ ',trim(gp_node2),' (',trim(adjustl(x_string)),  &
                                       ',',trim(adjustl(y_string)),  &
                       ',',trim(adjustl(z_string)),')'
  
  gp_xmin=min(gp_x(1),gp_x(2),gp_x(3))
  gp_xmax=max(gp_x(1),gp_x(2),gp_x(3))
  gp_xc=(gp_xmin+gp_xmax)/2d0
  gp_yc=gp_y(1)
  gp_w=(gp_xmax-gp_xmin)
  gp_h=gp_t
  
  CALL plot_layer(gp_xc,gp_yc,gp_w,gp_h,1d0,0d0,10)
    
  CALL plot_grid(gp_xc,gp_yc,gp_w,gp_h,1d0,0d0,1d0,gp_rh,gp_seg1,gp_nhinc,12)
                       
end if

! LOOP OVER THE CONDUCTORS AND DEFINE THEN WRITE THE NODES FOR EACH LAYER ON THE CONDUCTOR

write(20,'(A)')'*'
write(20,'(A)')'* Specify the conductor nodes'
do conductor=1,n_conductors

  write(20,'(A)')'*'
  write(20,'(A,I4)')'* Conductor',conductor
  
  write(conductor_string,'(I4)')conductor
  
  ALLOCATE( conductor_data(conductor)%x(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%y(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%w(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%h(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%d(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%anwinc(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%anhinc(1:conductor_data(conductor)%tot_n_layers) )

  ALLOCATE( conductor_data(conductor)%wx(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%wy(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%wz(1:conductor_data(conductor)%tot_n_layers) )

  conductor_data(conductor)%wx(1:conductor_data(conductor)%tot_n_layers)=1.0
  conductor_data(conductor)%wy(1:conductor_data(conductor)%tot_n_layers)=0.0
  conductor_data(conductor)%wz(1:conductor_data(conductor)%tot_n_layers)=0.0  

  ALLOCATE( conductor_data(conductor)%node1_list(1:conductor_data(conductor)%tot_n_layers) )
  ALLOCATE( conductor_data(conductor)%node2_list(1:conductor_data(conductor)%tot_n_layers) )
  
  if (conductor_data(conductor)%type.EQ.type_cyl) then 
    
    if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then

      do layer=-conductor_data(conductor)%n_layers2+1,conductor_data(conductor)%n_layers2
  
        layer_number=conductor_data(conductor)%n_layers2+layer
    
        ymin=conductor_data(conductor)%dl*real(layer)
        ymax=conductor_data(conductor)%dl*real(layer-1)
        y=(ymin+ymax)/2.0
        x=sqrt(conductor_data(conductor)%rc**2-y**2)
    
        conductor_data(conductor)%x(layer_number)=conductor_data(conductor)%xc
        conductor_data(conductor)%y(layer_number)=conductor_data(conductor)%yc+y
        conductor_data(conductor)%w(layer_number)=2.0*x
        conductor_data(conductor)%h(layer_number)=conductor_data(conductor)%dl
        conductor_data(conductor)%d(layer_number)=0.0
        conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
        conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
   
      end do  ! next layer
        
    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
    
! Loop over the grid and set segments within the circular conductor

INCLUDE "WRITE_FH2_IPFILE/set_segments_from_grid.F90"
        
    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_shell) then
       
      write(*,*)'SET SHELL DATA FOR CYLINDER'
      STOP 1
     
    end if   ! mesh_type_grid
  
  else if (conductor_data(conductor)%type.EQ.type_seven_strand) then 
    
    if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then
    
      layer_number=0

! loop over 7 strands      
      do i=1,7

        do layer=-conductor_data(conductor)%n_layers2ss+1,conductor_data(conductor)%n_layers2ss
  
          layer_number=layer_number+1
    
          ymin=conductor_data(conductor)%dl*real(layer)
          ymax=conductor_data(conductor)%dl*real(layer-1)
          y=(ymin+ymax)/2.0
          x=sqrt(conductor_data(conductor)%rss**2-y**2)
    
          conductor_data(conductor)%x(layer_number)=conductor_data(conductor)%xcss(i)
          conductor_data(conductor)%y(layer_number)=conductor_data(conductor)%ycss(i)+y
          conductor_data(conductor)%w(layer_number)=2.0*x
          conductor_data(conductor)%h(layer_number)=conductor_data(conductor)%dl
          conductor_data(conductor)%d(layer_number)=0.0
          conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
          conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
   
        end do  ! next layer
   
      end do  ! next strand
        
    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
    
! Loop over the grid and set segments within the circular conductor

INCLUDE "WRITE_FH2_IPFILE/set_segments_from_grid.F90"
     
    end if   ! mesh_type_grid
    
  else if (conductor_data(conductor)%type.EQ.type_rect) then 
  
    layer_number=1
        
    conductor_data(conductor)%x(layer_number)=conductor_data(conductor)%xc
    conductor_data(conductor)%y(layer_number)=conductor_data(conductor)%yc
    conductor_data(conductor)%w(layer_number)=conductor_data(conductor)%width
    conductor_data(conductor)%h(layer_number)=conductor_data(conductor)%height
    conductor_data(conductor)%d(layer_number)=0.0
    conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
    conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
     
  else if (conductor_data(conductor)%type.EQ.type_annulus) then 
  
    if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then

      do layer_number=1,conductor_data(conductor)%tot_n_layers
      
        angle=real(layer_number)*2.0*pi/real(conductor_data(conductor)%tot_n_layers)
        conductor_data(conductor)%wx(layer_number)=-sin(angle)
        conductor_data(conductor)%wy(layer_number)=cos(angle)
        conductor_data(conductor)%wz(layer_number)=0.0
    
    lrc=(conductor_data(conductor)%rco+conductor_data(conductor)%rci)/2.0
    lh=(conductor_data(conductor)%rco-conductor_data(conductor)%rci)
    lw=lrc*2.0*pi/real(conductor_data(conductor)%tot_n_layers)
    lxc=lrc*cos(angle)
    lyc=lrc*sin(angle)
          
        conductor_data(conductor)%x(layer_number)=conductor_data(conductor)%xc+lxc
        conductor_data(conductor)%y(layer_number)=conductor_data(conductor)%yc+lyc
        conductor_data(conductor)%w(layer_number)=lw
        conductor_data(conductor)%h(layer_number)=lh
        conductor_data(conductor)%d(layer_number)=0.0
        conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
        conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
   
      end do  ! next layer
        
    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
    
! Loop over the grid and set segments within the circular conductor

INCLUDE "WRITE_FH2_IPFILE/set_segments_from_grid.F90"
     
    end if   ! mesh_type_grid
    
  else
    write(*,*)'Cant deal with this conductor type yet:',conductor_data(conductor)%type
    STOP 1
  end if

! WRITE_THE NODES TO THE INPUT FILE AND SAVE IN THE NODE LIST FOR THIS CONDUCTOR

  do layer_number=1,conductor_data(conductor)%tot_n_layers
  
    write(layer_string,'(I4)')layer_number
    
    node_string='n_c'//trim(adjustl(conductor_string))//'_e1_l'//trim(adjustl(layer_string))   
    write(20,'(A,A,ES12.4,A,ES12.4,A)')trim(node_string),                      &
                         ' x=',conductor_data(conductor)%x(layer_number),          &
                         ' y=',conductor_data(conductor)%y(layer_number),' z=0.0'
    conductor_data(conductor)%node1_list(layer_number)=node_string
    
    node_string='n_c'//trim(adjustl(conductor_string))//'_e2_l'//trim(adjustl(layer_string))   
    write(20,'(A,A,ES12.4,A,ES12.4,A)')trim(node_string),                      &
                         ' x=',conductor_data(conductor)%x(layer_number),          &
                         ' y=',conductor_data(conductor)%y(layer_number),' z=1.0'
    conductor_data(conductor)%node2_list(layer_number)=node_string
    
    wx=conductor_data(conductor)%w(layer_number)/2.0
    hy=conductor_data(conductor)%h(layer_number)/2.0
    rw=conductor_data(conductor)%rw
    rh=conductor_data(conductor)%rh              
    nwinc=conductor_data(conductor)%anwinc(layer_number)
    nhinc=conductor_data(conductor)%anhinc(layer_number)
    
    vx=conductor_data(conductor)%wx(layer_number)
    vy=conductor_data(conductor)%wy(layer_number)
     
    CALL plot_layer(conductor_data(conductor)%x(layer_number),conductor_data(conductor)%y(layer_number), &
                    wx,hy,vx,vy,10)
    
    CALL plot_grid(conductor_data(conductor)%x(layer_number),conductor_data(conductor)%y(layer_number), &
                    wx,hy,vx,vy,rw,rh,nwinc,nhinc,12)
  
  end do  ! next layer

end do  ! next conductor

! LOOP OVER THE CONDUCTORS AND WRITE THE WIRE SEGMENTS FOR EACH LAYER ON THE CONDUCTOR

tot_n_segments=0
tot_n_filaments=0

write(20,'(A)')'*'
write(20,'(A)')'* conductor segments'
do conductor=1,n_conductors

  write(20,'(A)')'*'
  write(20,'(A,I4)')'* Conductor',conductor

  write(conductor_string,'(I4)')conductor
  
  do layer_number=1,conductor_data(conductor)%tot_n_layers
  
    write(layer_string,'(I4)')layer_number
    
    segment_string='E_c'//trim(adjustl(conductor_string))//'_l'//trim(adjustl(layer_string)) 
    
    write(20,'(A,A,A,A,A,A,ES12.4,A,ES12.4,A,ES12.4,A,ES12.4,A,ES12.4,A,ES12.4,A,I4,A,I4,A,ES12.4,A,ES12.4)')         &
                           trim(segment_string),' ',                                       &
                           trim(conductor_data(conductor)%node1_list(layer_number)),' ',   &
               trim(conductor_data(conductor)%node2_list(layer_number)),       &  
                           ' h=',conductor_data(conductor)%h(layer_number),                &
                           ' w=',conductor_data(conductor)%w(layer_number),                &
               ' sigma=',conductor_data(conductor)%sigma,                      &
                           ' wx=',conductor_data(conductor)%wx(layer_number),              &
                           ' wy=',conductor_data(conductor)%wy(layer_number),              &
                           ' wz=',conductor_data(conductor)%wz(layer_number),              &
               ' nhinc=',conductor_data(conductor)%anhinc(layer_number),                      &
               ' nwinc=',conductor_data(conductor)%anwinc(layer_number),                      &
               ' rh=',conductor_data(conductor)%rh,                            &
               ' rw=',conductor_data(conductor)%rw
               
    tot_n_segments=tot_n_segments+1
    tot_n_filaments=tot_n_filaments+conductor_data(conductor)%anhinc(layer_number)   &
                                   *conductor_data(conductor)%anwinc(layer_number)
    
  end do
  
end do

! EQUIVALENCE CONDUCTOR NODES AT NEAR END

write(20,'(A)')'*'
write(20,'(A)')'* Near end equivalent nodes on conductors'
do conductor=1,n_conductors

  if (conductor_data(conductor)%tot_n_layers.GT.1) then

    write(20,'(A)')'*'
    write(20,'(A,I4)')'* Conductor',conductor
  
    write(20,'(A)')'.equiv'  
    do layer_number=1,conductor_data(conductor)%tot_n_layers 
      write(20,'(A,A)')'+ ',trim(conductor_data(conductor)%node1_list(layer_number))  
    end do 
    
  end if
  
end do

! EQUIVALENCE CONDUCTOR NODES AT FAR END

write(20,'(A)')'*'
write(20,'(A)')'* Far end equivalent nodes on conductors'
do conductor=1,n_conductors

  if (conductor_data(conductor)%tot_n_layers.GT.1) then

    write(20,'(A)')'*'
    write(20,'(A,I4)')'* Conductor',conductor
  
    write(20,'(A)')'.equiv'  
    do layer_number=1,conductor_data(conductor)%tot_n_layers
      write(20,'(A,A)')'+ ',trim(conductor_data(conductor)%node2_list(layer_number))   
    end do
  
  end if
  
end do

! DEFINE LOOPS

if (ground_plane) then

  write(20,'(A)')'*'
  write(20,'(A)')'* Make loops from ground plane to all other conductors in turn'
  
  do conductor=1,n_conductors

    write(20,'(A)')'*'
    write(20,'(A,I4)')'* Ground plane to conductor',conductor
  
    write(conductor_string,'(I4)')conductor
    loop_string='loop_'//trim(adjustl(conductor_string))   
  
    write(20,'(A,A,A,A,A,A)')'.external ',trim(gp_node1),   &
                                      ' ',trim(conductor_data(conductor)%node1_list(1)), &
                        ' ',trim(loop_string)    
  end do

! JOIN ALL FAR END CONDUCTORS

  write(20,'(A)')'*'
  write(20,'(A)')'* Join all Far end conductors'
  
  write(20,'(A)',ADVANCE='NO')'.equiv'
  write(20,'(A,A)',ADVANCE='NO')' ',trim(gp_node2) 
  do conductor=1,n_conductors    

    layer_number=1  
    write(20,'(A,A)',ADVANCE='NO')' ',trim(conductor_data(conductor)%node2_list(layer_number)) 
  
  end do
    
  write(20,*)

else
! No ground plane 

  write(20,'(A)')'*'
  write(20,'(A)')'* Make loops from the last conductor to all other conductors in turn'
  
  do conductor=1,n_conductors-1

    write(20,'(A)')'*'
    write(20,'(A,I4)')'* last Conductor to conductor',conductor
  
    write(conductor_string,'(I4)')conductor
    loop_string='loop_'//trim(adjustl(conductor_string))   
  
    write(20,'(A,A,A,A,A,A)')'.external ',trim(conductor_data(n_conductors)%node1_list(1)),   &
                                      ' ',trim(conductor_data(conductor)%node1_list(1)), &
                        ' ',trim(loop_string)    
  end do

! JOIN ALL FAR END CONDUCTORS

  write(20,'(A)')'*'
  write(20,'(A)')'* Join all Far end conductors'
  
  write(20,'(A)',ADVANCE='NO')'.equiv'
  do conductor=1,n_conductors    

    layer_number=1  
    write(20,'(A,A)',ADVANCE='NO')' ',trim(conductor_data(conductor)%node2_list(layer_number)) 
  
  end do
    
  write(20,*)
  
end if ! ground plane

! WRITE THE FREQUENCY RANGE AND END

write(20,'(A)')'*'
write(20,'(A,ES12.4,A,ES12.4,A,ES12.4)')'.freq fmin=',fmin,' fmax=',fmax,' ndec=',ndec
write(20,'(A)')'*'
write(20,'(A)')'.end'

! CLOSE FILES

close(unit=10)
close(unit=12)
close(unit=20)

! DEALLOCATE MEMORY

do conductor=1,n_conductors

  if (ALLOCATED( conductor_data(conductor)%x )) DEALLOCATE( conductor_data(conductor)%x )
  if (ALLOCATED( conductor_data(conductor)%y )) DEALLOCATE( conductor_data(conductor)%y )
  if (ALLOCATED( conductor_data(conductor)%w )) DEALLOCATE( conductor_data(conductor)%w )
  if (ALLOCATED( conductor_data(conductor)%h )) DEALLOCATE( conductor_data(conductor)%h )

  if (ALLOCATED( conductor_data(conductor)%wx )) DEALLOCATE( conductor_data(conductor)%wx )
  if (ALLOCATED( conductor_data(conductor)%wy )) DEALLOCATE( conductor_data(conductor)%wy )
  if (ALLOCATED( conductor_data(conductor)%wz )) DEALLOCATE( conductor_data(conductor)%wz )

  if (ALLOCATED( conductor_data(conductor)%d )) DEALLOCATE( conductor_data(conductor)%d )
  if (ALLOCATED( conductor_data(conductor)%anwinc )) DEALLOCATE( conductor_data(conductor)%anwinc )
  if (ALLOCATED( conductor_data(conductor)%anhinc )) DEALLOCATE( conductor_data(conductor)%anhinc )

  if (ALLOCATED( conductor_data(conductor)%node1_list )) DEALLOCATE( conductor_data(conductor)%node1_list )
  if (ALLOCATED( conductor_data(conductor)%node2_list )) DEALLOCATE( conductor_data(conductor)%node2_list )

  if (ALLOCATED( conductor_data(conductor)%grid )) DEALLOCATE( conductor_data(conductor)%grid )
  if (ALLOCATED( conductor_data(conductor)%depth )) DEALLOCATE( conductor_data(conductor)%depth )
  
end do

if ( ALLOCATED(conductor_data) ) DEALLOCATE( conductor_data )

write(*,*)
write(*,*)'Total number of conductor segments =',tot_n_segments
write(*,*)'Total number of conductor filaments=',tot_n_filaments
write(*,*)

if (ground_plane) then
  write(*,*)'Total number of ground plane segments =',gp_n_segments
  write(*,*)'Total number of ground plane filaments=',gp_n_filaments
  write(*,*)
end if

write(*,*)'Total number of segments =',gp_n_segments+tot_n_segments
write(*,*)'Total number of filaments=',gp_n_filaments+tot_n_filaments
write(*,*)

STOP 0

9000 write(*,*)'ERROR reading the Cross Section Specification data, line',line
STOP 1

9010 write(*,*)'ERROR conductors should be numbered in order, line',line
STOP 1

9020 write(*,*)'ERROR Unknown conductor type:',type_ch,', line',line
STOP 1


END PROGRAM write_FH_input_file
!
! ___________________________________________________
!
!
SUBROUTINE plot_layer(xc,yc,wx,wy,vxx,vxy,file_unit)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

real(dp)    :: xc,yc,wx,wy,vxx,vxy
integer :: file_unit

! local variables

real(dp) :: rotx,roty

! START

 CALL rotate(-wx,+wy,vxx,vxy,rotx,roty)
 write(file_unit,*)xc+rotx,yc+roty
 CALL rotate(+wx,+wy,vxx,vxy,rotx,roty)
 write(file_unit,*)xc+rotx,yc+roty
 CALL rotate(+wx,-wy,vxx,vxy,rotx,roty)
 write(file_unit,*)xc+rotx,yc+roty
 CALL rotate(-wx,-wy,vxx,vxy,rotx,roty)
 write(file_unit,*)xc+rotx,yc+roty
 CALL rotate(-wx,+wy,vxx,vxy,rotx,roty)
 write(file_unit,*)xc+rotx,yc+roty
 write(file_unit,*)

RETURN

END SUBROUTINE plot_layer
!
! _______________________________________________________
!
!
SUBROUTINE plot_grid(xc,yc,wx,wy,vxx,vxy,rx,ry,nx,ny,file_unit)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

real(dp)    :: xc,yc,wx,wy,vxx,vxy,rx,ry
integer :: nx,ny,file_unit

! local variables

real(dp)    :: dx,dy,den,ox,oy
integer :: nl2
integer :: i
logical :: odd
real(dp) :: rotx,roty

! START

! lines in the x direction
if ( mod(nx,2).EQ.0 ) then
  odd=.FALSE.
else
  odd=.TRUE.
end if

if (odd) then
  nl2=(nx-1)/2
else
  nl2=nx/2 
end if

den=0.0
do i=1,nl2
  den=den+2.0*rx**(i-1)
end do
if (odd) den=den+rx**(nl2)

dx=wx*2.0/den

!write(*,*)'xc=',xc
!write(*,*)'nx=',nx
!write(*,*)'odd=',odd
!write(*,*)'nl2=',nl2
!write(*,*)'wx=',wx
!write(*,*)'den=',den
!write(*,*)'dx=',dx

ox=wx
do i=1,nl2
  ox=ox-dx*(rx**(i-1))
!  write(*,*)'i=',i,' ox=',ox
  CALL rotate(-ox,+wy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  CALL rotate(-ox,-wy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  write(file_unit,*)
  CALL rotate(+ox,+wy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  CALL rotate(+ox,-wy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  write(file_unit,*)
end do

if (.NOT.odd) then  ! write centre line
  CALL rotate(0d0,+wy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  CALL rotate(0d0,-wy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  write(file_unit,*)
end if

! lines in the y direction
if ( mod(ny,2).EQ.0 ) then
  odd=.FALSE.
else
  odd=.TRUE.
end if

if (odd) then
  nl2=(ny-1)/2
else
  nl2=ny/2 
end if

den=0.0
do i=1,nl2
  den=den+2.0*ry**(i-1)
end do
if (odd) den=den+ry**(nl2)

dy=wy*2.0/den

oy=wy
do i=1,nl2
  oy=oy-dy*(ry**(i-1))
  
  CALL rotate(+wx,-oy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  CALL rotate(-wx,-oy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  write(file_unit,*)
  CALL rotate(+wx,+oy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  CALL rotate(-wx,+oy,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  write(file_unit,*)
  
!  write(file_unit,*)xc+wx,yc-oy
!  write(file_unit,*)xc-wx,yc-oy
!  write(file_unit,*)
!  write(file_unit,*)xc+wx,yc+oy
!  write(file_unit,*)xc-wx,yc+oy
!  write(file_unit,*)

end do

if (.NOT.odd) then  ! write centre line

  CALL rotate(+wx,0d0,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  CALL rotate(-wx,0d0,vxx,vxy,rotx,roty)
  write(file_unit,*)xc+rotx,yc+roty
  write(file_unit,*)

!  write(file_unit,*)xc+wx,yc
!  write(file_unit,*)xc-wx,yc
!  write(file_unit,*)
end if

RETURN

END SUBROUTINE plot_grid
!
! _______________________________________________________
!
!
SUBROUTINE rotate(x,y,vxx,vxy,rx,ry)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

real(dp)    :: x,y,vxx,vxy,rx,ry

! local variables

real(dp)    :: vyx,vyy

! START

vyx=-vxy
vyy=vxx

rx=x*vxx+y*vyx
ry=x*vxy+y*vyy

RETURN

END SUBROUTINE rotate
!
! _______________________________________________________
!
!
SUBROUTINE calc_nmin(delta,t,r,n)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

real(dp)    :: delta,t,r
integer     :: n

! local variables

integer,parameter :: imax=21
integer :: i
real(dp) :: ndelta
real(dp) :: nequiv,nadd

! START

write(*,*)'CALLED: calc_nmin'

write(*,*)'delta=',delta
write(*,*)'t=',t
write(*,*)'r=',r

ndelta=NINT(t/delta)

write(*,*)'Number of skin depths in thickness=',ndelta

nadd=1d0/r
nequiv=0d0

do i=1,imax
  
  if(mod(i,2).GT.0) then
! i is odd 
    nadd=nadd*r
    nequiv=nequiv+nadd
  else
! i is even
    nequiv=nequiv+nadd
  end if
  
  write(*,*)'i=',i,' nadd=',nadd,' nequiv=',nequiv
  if (nequiv.GT.ndelta) then
    n=i
    write(*,*)'Number of layers required=',n
    RETURN
  end if

end do

n=i
RETURN

RETURN

END SUBROUTINE calc_nmin