!
! 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 .
!
! 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 .
!
! 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