Update to proximity effects: working for simple cases

Chris Smartt
1 parent ca596cf4
Showing 34 changed files with 1548 additions and 81 deletions Show diff stats
PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE/single_wire.cable_spec
PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/single_wire.cable_spec
PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire.bundle_spec
PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire_ac.spice_model_spec
PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_LAPLACE/single_wire.cable_spec
PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/plot_comparison.plt
PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/status
PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/single_wire.cable_spec
PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground.bundle_spec
PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground_ac.spice_model_spec
PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/status
SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90
SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90
SRC/CABLE_MODULES/cable_module.F90
SRC/CABLE_MODULES/coax.F90
SRC/CABLE_MODULES/conductor_impedance_model.F90
SRC/CABLE_MODULES/cylindrical.F90
SRC/CREATE_SPICE_CIRCUIT_MODEL/create_spice_subcircuit_model.F90
SRC/GENERAL_MODULES/frequency_spec.F90
SRC/GENERAL_MODULES/general_module.F90
@@ -5,7 +5,7 @@ Cylindrical
 3                # number of parameters
     1.0000E-04   # parameter 1: conductor radius
     1.0000E-04   # parameter 2: dielectric radius
-    0.0000E+00   # parameter 3: conductivity
+    4.461E+07   # parameter 3: conductivity
 1                # number of frequency dependent parameters
 # Dielectric relative permittivity model follows
    1.0000000000000000        # w normalisation constant
@@ -5,7 +5,7 @@ Cylindrical
 3                # number of parameters
     1.0000E-04   # parameter 1: conductor radius
     1.0000E-04   # parameter 2: dielectric radius
-    0.0000E+00   # parameter 3: conductivity
+    4.461E+07   # parameter 3: conductivity
 1                # number of frequency dependent parameters
 # Dielectric relative permittivity model follows
    1.0000000000000000        # w normalisation constant
@@ -10,7 +10,7 @@ single_wire
 no_ground_plane
 -10     # order for vector fit model
 log          # frequency scale (log or lin)
-1e3 1e9 60 # fmin fmax number_of_frequencies
+1e2 1e8 60 # fmin fmax number_of_frequencies
 verbose
 use_laplace
 use_fasthenry
@@ -31,3 +31,8 @@ log          # frequency scale (log or lin)
 # Output conductor number and end number
 1    1
 lin   # output type (lin or dB)
+-10     # order for vector fit model
+log          # frequency scale (log or lin)
+1e2 1e8 60 # fmin fmax number_of_frequencies
+verbose
+plot_propagation_correction_filter_fit_data
@@ -5,7 +5,7 @@ Cylindrical
 3                # number of parameters
     1.0000E-04   # parameter 1: conductor radius
     1.0000E-04   # parameter 2: dielectric radius
-    0.0000E+00   # parameter 3: conductivity
+    4.461E+07   # parameter 3: conductivity
 1                # number of frequency dependent parameters
 # Dielectric relative permittivity model follows
    1.0000000000000000        # w normalisation constant
@@ -5,6 +5,7 @@ set xlabel &#39;Frequency(Hz)&#39;
 set ylabel 'V(f)'
 set logscale x
 plot './TWO_WIRE/RUN_DIRECTORY/analytic_solution.dat' u 1:2 title 'Analytic solution' w p lc 3,\
+     './TWO_WIRE_FASTHENRY/RUN_DIRECTORY/analytic_solution.dat' u 1:2 title 'Analytic solution2' w p lc 3,\
      './TWO_WIRE/RUN_DIRECTORY/spice_solution.dat' u 2:3 title 'TWO_WIRE' w l lc 1,\
      './TWO_WIRE_LAPLACE/RUN_DIRECTORY/spice_solution.dat' u 2:3 title 'TWO_WIRE_LAPLACE' w l lc 2,\
      './TWO_WIRE_FASTHENRY/RUN_DIRECTORY/spice_solution.dat' u 2:3 title 'TWO_WIRE_FASTHENRY' w l lc 4
-generate_spice_cable_bundle_model clean
+generate_spice_cable_bundle_model run
-TWO_WIRE: Finished correctly  
-TWO_WIRE_LAPLACE: Finished correctly  
-TWO_WIRE_FASTHENRY: Finished correctly  
+TWO_WIRE_FASTHENRY: Finished correctly Result comparison:    0.00000000
@@ -5,7 +5,7 @@ Cylindrical
 3                # number of parameters
     1.0000E-04   # parameter 1: conductor radius
     1.0000E-04   # parameter 2: dielectric radius
-    0.0000E+00   # parameter 3: conductivity
+    4.461E+07    # parameter 3: conductivity
 1                # number of frequency dependent parameters
 # Dielectric relative permittivity model follows
    1.0000000000000000        # w normalisation constant
@@ -5,11 +5,10 @@
 1    # Number of cables in bundle, cable list follows: name then x,y of cable
 single_wire
 0.0  4.e-3   ! x and y coordinates of the cable centre
-ground_plane
+ground_plane  4.461E+07 25e-3 0.1e-3 10 10 1     ! sigma w h nx nz nh
 -10     # order for vector fit model
 log          # frequency scale (log or lin)
-1e3 1e9 60 # fmin fmax number_of_frequencies
+1e2 1e8 60 # fmin fmax number_of_frequencies
 verbose
 use_laplace
-use_fasthenry
-
+use_fasthenry 4 3  3  2.0 2.0    ! condcutor meshing parameters: nlayers_radius nw   nh   rw   rh
@@ -27,7 +27,12 @@ wire_over_ground
 # Type of analysis
 AC
 log          # frequency scale (log or lin)
-1e3 1e8 1000 # fmin fmax number_of_frequencies
+1e3 1e8 1000 # fmin fmax number_of_frequencies   1e3 1e8 1000 # fmin fmax number_of_frequencies
 # Output conductor number and end number
 1    1
 lin   # output type (lin or dB)
+-10     # order for vector fit model
+log          # frequency scale (log or lin)
+1e2 1e8 60 # fmin fmax number_of_frequencies
+verbose
+plot_propagation_correction_filter_fit_data
-generate_spice_cable_bundle_model clean
+generate_spice_cable_bundle_model plot
-WIRE_OVER_GROUND_PLANE: Finished correctly  
-WIRE_OVER_GROUND_PLANE_LAPLACE: Finished correctly  
-WIRE_OVER_GROUND_PLANE_FASTHENRY: Finished correctly  
+WIRE_OVER_GROUND_PLANE_FASTHENRY/: Finished correctly  
@@ -1035,6 +1035,8 @@ USE maths
             PUL%y(conductor)=bundle%cable_y_offset(cable)   
             PUL%rtheta(conductor)=bundle%cable_angle(cable)
+            PUL%sigma(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_sigma
+	    
             PUL%ox(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_ox
             PUL%oy(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_oy
             PUL%r(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_radius
@@ -1064,12 +1066,19 @@ USE maths
       overshield=overshield+1
-      if (verbose) write(*,*)'PUL parameter calculation for oversheild domain'
+      if (verbose) write(*,*)'PUL parameter calculation for overshield domain'
 ! no ground plane
       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     =0
+           
 ! add overshield information   
       if (verbose) write(*,*)'Add overshield information'
       is_overshield_domain=.TRUE.
@@ -1080,7 +1089,7 @@ USE maths
       PUL%overshield_y = overshield_y(overshield)
       PUL%overshield_r = overshield_r(overshield)
       PUL%overshield_w = overshield_w(overshield)
-      PUL%overshield_w2 = overshield_w2(overshield)
+      PUL%overshield_w2= overshield_w2(overshield)
       PUL%overshield_h = overshield_h(overshield)
       PUL%epsr_background = 1d0 ! background permittivity =1.0 within an overshield i.e. cables are in air
@@ -1109,14 +1118,21 @@ USE maths
       PUL%ground_plane_present=bundle%ground_plane_present
       PUL%ground_plane_angle  =bundle%ground_plane_angle
       PUL%ground_plane_offset =bundle%ground_plane_offset
-     
+      PUL%ground_plane_sigma  =bundle%ground_plane_sigma   ! ground plane conductivity
+      PUL%ground_plane_w      =bundle%ground_plane_w       ! ground plane width
+      PUL%ground_plane_h      =bundle%ground_plane_h       ! ground plane height
+      PUL%ground_plane_Rdc    =bundle%ground_plane_Rdc     ! ground plane dc resistance
+      PUL%ground_plane_nsegx  =bundle%ground_plane_nsegx      ! ground plane x discretisation
+      PUL%ground_plane_nsegz  =bundle%ground_plane_nsegz      ! ground plane z discretisation
+      PUL%ground_plane_nh     =bundle%ground_plane_nh      ! ground plane h discretisation
+
 ! No overshield information   
       is_overshield_domain=.FALSE.
       PUL%overshield_present=.FALSE.
       PUL%overshield_shape=0
       PUL%overshield_r= 0d0
       PUL%overshield_w= 0d0
-      PUL%overshield_w2= 0d0
+      PUL%overshield_w2=0d0
       PUL%overshield_h= 0d0
       PUL%epsr_background = 1d0 ! background permittivity =1.0 i.e. cables are in air
@@ -1134,6 +1150,10 @@ USE maths
       if (use_FastHenry) then
          CALL PUL_RL_FastHenry2(PUL,bundle%bundle_name,bundle%Y_fit_model_order,bundle%Y_fit_freq_spec,domain)
+
+!! Change the conductor model type as we now let FastHenry2 deal with this. 	 
+!	 bundle%cable(cable)%conductor_impedance(local_conductor)%impedance_model_type=impedance_model_type_FH2
+!	 bundle%cable(cable)%conductor_impedance(local_conductor)%Rdc,
       end if
@@ -1278,11 +1298,12 @@ USE maths
         conductor=conductor+1
 ! copy the conductor impedance model for this cable conductor to the bundle structure        
+
+        write(*,*)'cable=',cable,' conductor=',conductor,' local_conductor=',local_conductor
         bundle%conductor_impedance(conductor)%impedance_model_type=  &
                     bundle%cable(cable)%conductor_impedance(local_conductor)%impedance_model_type
-
         bundle%conductor_impedance(conductor)%radius=  &
                     bundle%cable(cable)%conductor_impedance(local_conductor)%radius
@@ -1300,6 +1321,9 @@ USE maths
         bundle%conductor_impedance(conductor)%Resistance_multiplication_factor=  &
                     bundle%cable(cable)%conductor_impedance(local_conductor)%Resistance_multiplication_factor
+        
+        bundle%conductor_impedance(conductor)%Rdc=  &
+                    bundle%cable(cable)%conductor_impedance(local_conductor)%Rdc
         if((bundle%conductor_impedance(conductor)%impedance_model_type.EQ.impedance_model_type_filter).OR. &  
            (bundle%conductor_impedance(conductor)%impedance_model_type.EQ.impedance_model_type_cylindrical_shield)) then
@@ -101,6 +101,14 @@ TYPE::bundle_specification_type
   real(dp)            :: ground_plane_angle
   real(dp)            :: ground_plane_offset
+  real(dp)            :: ground_plane_sigma      ! ground plane conductivity for FastHenry2
+  real(dp)            :: ground_plane_w          ! ground plane width for FastHenry2
+  real(dp)            :: ground_plane_h          ! ground plane height for FastHenry2
+  real(dp)            :: ground_plane_Rdc        ! ground plane dc resistance for FastHenry2
+  integer             :: ground_plane_nsegx      ! ground plane number of segments in x for FastHenry2
+  integer             :: ground_plane_nsegz      ! ground plane number of segments in z for FastHenry2
+  integer             :: ground_plane_nh         ! ground plane number of layers in h for FastHenry2
+  
   real(dp)            :: ground_plane_x,ground_plane_y,ground_plane_theta
   real(dp)            :: ground_plane_nx,ground_plane_ny
@@ -270,6 +278,20 @@ CONTAINS
 ! calculate offset
     bundle%ground_plane_offset=-bundle%ground_plane_x*sin(bundle%ground_plane_angle)+  &
                                 bundle%ground_plane_y*cos(bundle%ground_plane_angle)
+
+    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_sigma
+    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_w  
+    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_h      
+    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_Rdc
+    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nsegx
+    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nsegz
+    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nh
+    
+    if (ierr.NE.0) then 
+      run_status='ERROR reading ground plane sigma, w h and Rdc'
+      CALL write_program_status()
+      STOP 1
+    end if 
 ! set the geometric data for the last cable in the list i.e. the ground plane
@@ -508,6 +530,15 @@ CONTAINS
                       ' x y coordinates and angle of ground plane '
     write(file_unit,*)bundle%ground_plane_nx,bundle%ground_plane_ny,' ground plane normal direction'
     write(file_unit,*)bundle%ground_plane_cable_side,' orientation of cables wrt ground plane'
+
+    write(file_unit,*)bundle%ground_plane_sigma,' ground plane conductivity'
+    write(file_unit,*)bundle%ground_plane_w,' ground plane width '  
+    write(file_unit,*)bundle%ground_plane_h,' ground plane height (thickness) '      
+    write(file_unit,*)bundle%ground_plane_Rdc,' ground plane Rdc '
+    write(file_unit,*)bundle%ground_plane_nsegx,' ground plane nsegx '
+    write(file_unit,*)bundle%ground_plane_nsegz,' ground plane nsegz '
+    write(file_unit,*)bundle%ground_plane_nh,' ground plane nh '
+
   else
     write(file_unit,'(A)')'no_ground_plane'
   end if  
@@ -126,6 +126,7 @@
 !    put the external condcutor model information into its own structure 6/10/2016 CJS
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
 !     16/3/2018 CJS Add ML_flex_cable
+!     24/10/2023 CJS Add impedance_model_type_FH2
 !
 MODULE cable_module
@@ -146,6 +147,7 @@ TYPE:: conductor_impedance_model
 ! impedance_model_type_cylidrical_shell_with_conductivity
 ! impedance_model_type_cylindrical_shield                
 ! impedance_model_type_rectangular_with_conductivity    
+! impedance_model_type_FH2   
   integer  :: impedance_model_type
@@ -173,6 +175,7 @@ TYPE:: external_conductor_model
   integer  :: conductor_type
   real(dp) :: conductor_radius
+  real(dp) :: conductor_sigma
   real(dp) :: conductor_width
   real(dp) :: conductor_width2
   real(dp) :: conductor_height
@@ -277,6 +280,7 @@ integer,parameter    :: impedance_model_type_filter                             
 integer,parameter    :: impedance_model_type_cylidrical_shell_with_conductivity =3
 integer,parameter    :: impedance_model_type_cylindrical_shield                 =4
 integer,parameter    :: impedance_model_type_rectangular_with_conductivity      =5
+integer,parameter    :: impedance_model_type_FH2                                =6
 CONTAINS
@@ -335,6 +339,7 @@ include &#39;conductor_impedance_model.F90&#39;
   ext_model%conductor_type=0d0
   ext_model%conductor_radius=0d0
+  ext_model%conductor_sigma=0d0
   ext_model%conductor_width=0d0
   ext_model%conductor_width2=0d0
   ext_model%conductor_height=0d0
@@ -512,6 +517,7 @@ END SUBROUTINE reset_external_conductor_model
   do i=1,cable%n_external_conductors
     read(file_unit,*) cable%external_model(i)%conductor_type
     read(file_unit,*) cable%external_model(i)%conductor_radius
+    read(file_unit,*) cable%external_model(i)%conductor_sigma
     read(file_unit,*) cable%external_model(i)%conductor_width
     read(file_unit,*) cable%external_model(i)%conductor_width2
     read(file_unit,*) cable%external_model(i)%conductor_height
@@ -703,6 +709,7 @@ END SUBROUTINE reset_external_conductor_model
   do i=1,cable%n_external_conductors
     write(file_unit,*) cable%external_model(i)%conductor_type,'  conductor type'
     write(file_unit,*) cable%external_model(i)%conductor_radius  ,' conductor_radius '
+    write(file_unit,*) cable%external_model(i)%conductor_sigma  ,' conductor_sigma '
     write(file_unit,*) cable%external_model(i)%conductor_width   ,' conductor_width  '
     write(file_unit,*) cable%external_model(i)%conductor_width2   ,' conductor_width2  '
     write(file_unit,*) cable%external_model(i)%conductor_height  ,' conductor_height '
@@ -245,6 +245,7 @@ IMPLICIT NONE
   CALL reset_external_conductor_model(cable%external_model(1))
   cable%external_model(1)%conductor_type=circle
   cable%external_model(1)%conductor_radius=rs
+  cable%external_model(1)%conductor_sigma=sigma_s
   cable%external_model(1)%dielectric_radius=rd
   cable%external_model(1)%dielectric_epsr=epsr2
@@ -39,6 +39,7 @@
 ! impedance_model_type_cylidrical_shell_with_conductivity
 ! impedance_model_type_cylindrical_shield                
 ! impedance_model_type_rectangular_with_conductivity    
+! impedance_model_type_FH2      
 ! File Contents:
 ! SUBROUTINE read_conductor_impedance_model
@@ -133,6 +134,10 @@ IMPLICIT NONE
     read(unit,*,ERR=9000)conductor_impedance%height
     read(unit,*,ERR=9000)conductor_impedance%conductivity
     read(unit,*,ERR=9000)conductor_impedance%Resistance_multiplication_factor
+  
+  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_FH2) then
+  
+    read(unit,*,ERR=9000)conductor_impedance%Rdc
   else
@@ -224,6 +229,10 @@ IMPLICIT NONE
     write(unit,*)conductor_impedance%height,' # conductor height'
     write(unit,*)conductor_impedance%conductivity,' # conductivity'
     write(unit,*)conductor_impedance%Resistance_multiplication_factor,' # Resistance_multiplication_factor'
+  
+  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_FH2) then
+  
+    write(unit,*)conductor_impedance%Rdc,' # conductor dc resistance'
   end if
@@ -344,6 +353,13 @@ IMPLICIT NONE
     t=conductor_impedance%height
     CALL calculate_internal_impedance_rectangular(sigma,w,t,f,Z_c,Rdc_c)
+ 
+  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_FH2) then
+  
+    Z_c=conductor_impedance%Rdc
+    Rdc_c=conductor_impedance%Rdc
+    Z_t=(0d0,0d0)
+    Rdc_t=0d0
   end if
@@ -174,6 +174,7 @@ IMPLICIT NONE
   CALL reset_external_conductor_model(cable%external_model(1))
   cable%external_model(1)%conductor_type=circle
   cable%external_model(1)%conductor_radius=r
+  cable%external_model(1)%conductor_sigma=sigma
   cable%external_model(1)%dielectric_radius=rd 
   cable%external_model(1)%dielectric_epsr=epsr
@@ -57,7 +57,7 @@
 !
 ! STAGE_9, generate the information required for and then write the transfer impedance coupling model(s) as required
 ! STAGE_10, generate the information required for and then write the incident field excitation model as required
-! STAGE 11: Dallocate the domain based information
+! STAGE 11: Deallocate the domain based information
 !     
 ! COMMENTS
 !     1. The notation needs to be organised a bit better.
@@ -69,6 +69,7 @@ real(dp)        :: fmin            ! minimum frequency
 real(dp)        :: fmax            ! maximum frequency
 integer         :: n_frequencies   ! number of frequencies
 real(dp),allocatable ::freq_list(:)! list of frequencies
+real(dp)        :: ndec            ! number of frequencies per decade: used in FastHenry2
 END TYPE frequency_specification
@@ -247,12 +248,16 @@ integer         :: ierr
       log_fstep=(log_fmax-log_fmin)/dble(freq_data%n_frequencies-1)
     end if
+    freq_data%ndec=dble(freq_data%n_frequencies-1)/log10(freq_data%fmax/freq_data%fmin)
+
   else if (freq_data%freq_range_type.EQ.'lin') then
     fstep=0d0                      ! this is the value used if freq_data%n_frequencies=1
     if (freq_data%n_frequencies.ne.1) then
       fstep=(freq_data%fmax-freq_data%fmin)/dble(freq_data%n_frequencies-1)
     end if
+    
+    freq_data%ndec=0d0
   else 
@@ -267,6 +267,14 @@ integer,parameter     :: local_file_unit=90
 integer,parameter     :: temp_file_unit=99
+! FastHenry2 default conductor meshing parameters
+
+integer :: FH2_nlayers_radius=4
+integer :: FH2_nw=3
+integer :: FH2_nh=3
+real(dp):: FH2_rw=2.0
+real(dp):: FH2_rh=2.0
+
 CONTAINS
 include 'add_integer_to_string.F90'
@@ -278,7 +278,7 @@ integer :: ierr    ! error code for matrix inverse calls
       end do  ! next col      
     end do ! next row
-! calculate the contribution to the matrices from the conductor based impedance models. Initailly set to zero
+! calculate the contribution to the matrices from the conductor based impedance models. Initially set to zero
 ! See Theory_Manual_Section 2.2.3 
     Z_domain_conductor_impedance_correction(1:dim,1:dim)=(0d0,0d0)
@@ -479,7 +479,7 @@ integer :: ierr    ! error code for matrix inverse calls
 ! Add the conductor impedance contributions to the domain based impedance matrix
     Z_domain(:,:)=Z_domain(:,:)+Z_domain_conductor_impedance_correction(:,:)
-    
+        
     if (verbose) then
       write(*,*)'[R_domain]=Re[Z_domain]'
@@ -173,7 +173,7 @@ integer  :: ierr
   ALLOCATE( Correction(1:dim,1:freq_spec%n_frequencies) )
   ALLOCATE( Hp(1:freq_spec%n_frequencies) )
-!  if (plot_propagation_correction_filter_fit_data) then
+  if (plot_propagation_correction_filter_fit_data) then
 ! open files for plotting the propagation correction functions
@@ -206,7 +206,7 @@ integer  :: ierr
     close(unit=80)
-!  end if ! plot_propagation_correction_filter_fit_data
+  end if ! plot_propagation_correction_filter_fit_data
 ! loop over frequency working out the correction required for each mode at each frequency
 ! Note we work downwards in frequency as we use high frequency values of L and C to define the
@@ -243,7 +243,7 @@ integer  :: ierr
     end do ! next row
-! Use the pre_calulcated delay (input to the subroutine) to obtain GAMMA_C_hf 
+! Use the pre_calculated delay (input to the subroutine) to obtain GAMMA_C_hf 
 ! Theory_Manual_Eqns 3.81, 3.82, 3.83
     GAMMA_C_hf(:)=j*w*delay(:)/length
@@ -422,9 +422,9 @@ integer  :: ierr
        if (plot_propagation_correction_filter_fit_data) then
 ! write propagation correction data for plotting
-!         write(80+row,8000)f,real(Correction(row,frequency_loop)), &
-!                           aimag(Correction(row,frequency_loop)),abs(Correction(row,frequency_loop))
-         write(80+row,8000)f,real(GAMMA_C(local_mode)),aimag(GAMMA_C(local_mode)),real(GAMMA_C_hf(row)),aimag(GAMMA_C_hf(row))
+         write(80+row,8000)f,real(Correction(row,frequency_loop)), &
+                           aimag(Correction(row,frequency_loop)),abs(Correction(row,frequency_loop))
+!         write(80+row,8000)f,real(GAMMA_C(local_mode)),aimag(GAMMA_C(local_mode)),real(GAMMA_C_hf(row)),aimag(GAMMA_C_hf(row))
        end if
 8000 format(5E16.6)
@@ -452,6 +452,7 @@ integer  :: ierr
 ! Call the filter fitting process with this mode propagation correction data 
 ! with border=aorder and fit_type=1 i.e. Hfilter(f=0)=1 is imposed
+
     CALL Calculate_Sfilter(Hp,freq_spec%freq_list,freq_spec%n_frequencies,Hfilter(i),model_order,0,1) 
     if (plot_propagation_correction_filter_fit_data) then
@@ -459,7 +460,10 @@ integer  :: ierr
       open(unit=80,file=trim(filename(i))//'_fit')
       do frequency_loop=1,freq_spec%n_frequencies
-        
+
+        f=freq_spec%freq_list(frequency_loop)  
+! shift the frequency from d.c. as jwL, jwC are zero at d.c. and this messes up the process
+        if (f.EQ.0d0) f=1d0  
         write(80,8010)f,evaluate_Sfilter_frequency_response(Hfilter(i),f)
 8010 format(3E16.6)
@@ -199,7 +199,8 @@ cable_bundle_model_builder \
 spice_cable_bundle_model_builder \
 shield_conductor_and_transfer_impedance_model_builder \
 rational_function_fit \
-compare_results 
+compare_results \
+write_FH_input_file
 MAKE_COMPILATION_DATE:
 	echo "SPICE_CABLE_MODEL_BUILDER_compilation_date='$(COMPILATION_DATE)' " > compilation_date.inc 
@@ -264,6 +265,11 @@ rational_function_fit: rational_function_fit.F90 $(MODS)
 compare_results: compare_results.F90  $(TYPE_SPEC_MODULE)
 	$(FC) $(FLAGS) -o $(EXECUTABLE_DIR)/compare_results compare_results.F90 $(OBJS) $(LIBS)
+write_FH_input_file: write_FH_input_file.F90 WRITE_FH2_IPFILE/create_grid.F90 \
+                                             WRITE_FH2_IPFILE/get_grid_type.F90  \
+                                             WRITE_FH2_IPFILE/set_segments_from_grid.F90
+	$(FC) $(FLAGS) -o $(EXECUTABLE_DIR)/write_FH_input_file write_FH_input_file.F90 $(TYPE_SPEC_OBJS)
+
 clean:
 	( rm -f $(EXECUTABLE_DIR)/cable_model_builder* )
 	( rm -f $(EXECUTABLE_DIR)/cable_bundle_model_builder* )
@@ -271,6 +277,7 @@ clean:
 	( rm -f $(EXECUTABLE_DIR)/shield_conductor_and_transfer_impedance_model_builder* )
 	( rm -f $(EXECUTABLE_DIR)/rational_function_fit* )
 	( rm -f $(EXECUTABLE_DIR)/compare_results* )
+	( rm -f $(EXECUTABLE_DIR)/write_FH_input_file* )
 	( rm -f *.o )
 	( rm -f $(OBJ_MOD_DIR)/*.mod )
 	( rm -f $(OBJ_MOD_DIR)/*.o )
@@ -1237,10 +1237,9 @@ IMPLICIT NONE
     else
-! There are frequency dependent dielectrics or proximity effects present             
-! Secondly we calculate the capacitance matrix at a number of frequencies before fitting filter functions
+! There are frequency dependent dielectrics  present             
+! We calculate the capacitance matrix at a number of frequencies before fitting filter functions
 ! to each of the frequency dependent admittance matrix entries.
-! 23/10/2023: Call FastHenry2 for proximity effects if required
 ! allocate the frequency dependent matrices
       ALLOCATE( C_freq(1:freq_spec%n_frequencies) )
@@ -1305,7 +1304,6 @@ IMPLICIT NONE
 ! make the matrix symmetrical
             PUL%Yfilter%sfilter_mat(col,row)=PUL%Yfilter%sfilter_mat(row,col)
             PUL%Zfilter%sfilter_mat(col,row)=PUL%Zfilter%sfilter_mat(row,col)
-
           end if
         end do ! next col
@@ -1341,7 +1339,6 @@ IMPLICIT NONE
     CALL dwrite_matrix(PUL%C%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
     write(*,*)'Conductance matrix, G'
     CALL dwrite_matrix(PUL%G%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
-
   end if
   if (allocated( dielectric_region_epsr )) DEALLOCATE( dielectric_region_epsr )
@@ -75,8 +75,6 @@ IMPLICIT NONE
   integer,parameter :: outside=2
 ! local variables
-
-  integer :: ndec
   integer :: first_conductor,last_conductor,conductor
@@ -85,6 +83,13 @@ IMPLICIT NONE
 ! Process_Zc variables:
+  integer  :: matrix_dimension
+    
+  complex(dp),allocatable :: function_to_fit(:)    ! complex function to be fitted using Sfilter_fit
+  
+  real(dp) :: ferr
+  
+  real(dp) :: gp_half_width
 character(LEN=256) :: line
 character(LEN=256) :: freq_and_dim_string
@@ -96,11 +101,13 @@ integer :: i
 integer :: n_freq_in
 real(dp) :: freq,w
-integer :: nrows,ncols,r,c
+integer :: nrows,ncols,r,c,row,col
 complex(dp),allocatable :: Zc_in(:,:,:)
 real(dp),allocatable    :: freq_in(:)
+real(dp),allocatable    :: Rdc_domain(:,:)
+
 real(dp),allocatable    :: values(:),re,im
 logical :: subtract_dc_resistance
@@ -111,16 +118,25 @@ integer :: ierr
 ! START
-  if (verbose) write(*,*)'CALLED: PUL_LC_Laplace'
+  if (verbose) write(*,*)'CALLED: PUL_RL_FastHenry2'
-! STAGE 1: work out the configuration for the calculation i.e. is there a ground plane, is the outer boundary free space or a conductor (overshield)
+! work out the local frequency sampling: FastHenry uses a specific form based on log frequency
+! and number of samples per decade
   write(*,*)'Filter fit order=',fit_order
   write(*,*)'fmin=',freq_spec%fmin
   write(*,*)'fmax=',freq_spec%fmax
-  write(*,*)'n_frequencies=',freq_spec%n_frequencies
+  write(*,*)'n_frequencies=',freq_spec%n_frequencies  
+  write(*,*)'ndec=',freq_spec%ndec 
-  ndec=NINT(real(freq_spec%n_frequencies)/log10(freq_spec%fmax/freq_spec%fmin))
+  if (freq_spec%freq_range_type.EQ.'lin') then
+    write(run_status,*)'ERROR in PUL_RL_FastHenry2: We must have a logarithmic frequency range specified'
+    CALL write_program_status()
+    STOP 1
+
+  end if
+
+! STAGE 1: work out the configuration for the calculation i.e. is there a ground plane, is the outer boundary free space or a conductor (overshield)
 ! write the file which is used to generate to FastHenry2 input file
   OPEN(unit=fh2_input_file_unit,file='fh2.txt')
@@ -136,21 +152,42 @@ integer :: ierr
     write(fh2_input_file_unit,'(A)')'ground_plane'
-! ********* HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************    
-    write(fh2_input_file_unit,'(A)')'-25e-3 0.0e-3  0.0e-3	! x y z of ground plane point 1	              '
-    write(fh2_input_file_unit,'(A)')' 25e-3 0.0e-3  0.0e-3	! x y z of ground plane point 2 	      '
-    write(fh2_input_file_unit,'(A)')' 25e-3 0.0e-3 1000.0e-3	! x y z of ground plane point 3 	      '
-    write(fh2_input_file_unit,'(A)')'0.01e-3			! thickness				      '
-    write(fh2_input_file_unit,'(A)')'10 10			! discretisation in p1-p2 and p2=p3 directions'
-    write(fh2_input_file_unit,'(A)')'4.461E+07  		! conductivity, sigma			      '
-    write(fh2_input_file_unit,'(A)')'1  			! gp discretisation in thickness,  nhinc      '
+! ********* STILL SOME HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************  
+
+! ALSO NEED TO CHECK THAT EVERYTHING WE NEED IS SET...
+
+    write(*,*)'Ground plane parameters:'
+    write(*,*)'PUL%ground_plane_w=',PUL%ground_plane_w
+    write(*,*)'PUL%ground_plane_h=',PUL%ground_plane_h
+    write(*,*)'PUL%ground_plane_sigma=',PUL%ground_plane_sigma
+    write(*,*)'PUL%ground_plane_nsegx=',PUL%ground_plane_nsegx
+    write(*,*)'PUL%ground_plane_nsegz=',PUL%ground_plane_nsegz
+    write(*,*)'PUL%ground_plane_nh=',PUL%ground_plane_nh
+
+    if ( (PUL%ground_plane_w.EQ.0d0).OR.(PUL%ground_plane_h.EQ.0d0).OR.(PUL%ground_plane_sigma.EQ.0d0) &
+     .OR.(PUL%ground_plane_nsegx.EQ.0).OR.(PUL%ground_plane_nsegz.EQ.0).OR.(PUL%ground_plane_nh.EQ.0) )then
+      write(run_status,*)'ERROR in PUL_RL_FastHenry2: ground plane parameters not defined'
+      CALL write_program_status()
+      STOP 1
+    end if
+
+    gp_half_width=PUL%ground_plane_w/2.0
+    
+    write(fh2_input_file_unit,'(ES12.4,A)')-gp_half_width,' 0.0e-3  0.0e-3	! x y z of ground plane point 1'
+    write(fh2_input_file_unit,'(ES12.4,A)') gp_half_width,' 0.0e-3  0.0e-3	! x y z of ground plane point 2'
+    write(fh2_input_file_unit,'(ES12.4,A)') gp_half_width,' 0.0e-3 1000.0e-3	! x y z of ground plane point 3'
+    write(fh2_input_file_unit,'(ES12.4,A)')PUL%ground_plane_h,'	 ! thickness				      '
+    write(fh2_input_file_unit,'(2I6,A)')PUL%ground_plane_nsegx,PUL%ground_plane_nsegz,  &
+                                        '  ! discretisation in p1-p2 and p2=p3 directions'
+    write(fh2_input_file_unit,'(ES12.4,A)')PUL%ground_plane_sigma,' ! conductivity, sigma			      '
+    write(fh2_input_file_unit,'(I4,A)')PUL%ground_plane_nh,'	! gp discretisation in thickness,  nhinc      '
     write(fh2_input_file_unit,'(A)')'2.0			! gp discretisation in thickness ratio, rh    '
     write(fh2_input_file_unit,'(A)')'0e-3 0.0e-3 0.0e-3 	! x y z of ground plane node at end 1	      '
     write(fh2_input_file_unit,'(A)')'0e-3 0.0e-3 1000.0e-3	! x y z of ground plane node at end 2	      '
-! ********* HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************    
+! ********* STILL SOME HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************    
   else if ((.NOT.PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then 
-    
+
     first_conductor=1
     last_conductor=PUL%n_conductors
@@ -174,10 +211,10 @@ integer :: ierr
       write(fh2_input_file_unit,'(A)')'cylindrical	! conductor 1 shape (cylinder, rectangle, annulus)'
       write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor),PUL%y(conductor),' ! centre xc yc'
       write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor),' ! radius rc'
-      write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor)/8,' ! cylindrical conductor discretisation size, dl'
-      write(fh2_input_file_unit,'(A)')'4.461E+07	  ! conductivity, sigma '
-      write(fh2_input_file_unit,'(A)')'3  3		  ! segment discretisation in x and y, nwinc nhinc '
-      write(fh2_input_file_unit,'(A)')'2.0 2.0  	  ! segment discretisation ratio in x and y, rw rh '
+      write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor)/FH2_nlayers_radius,' ! cylindrical conductor discretisation size, dl'
+      write(fh2_input_file_unit,'(ES12.4,A)')PUL%sigma(conductor),'  ! conductivity, sigma '
+      write(fh2_input_file_unit,'(2I4,A)')FH2_nw,FH2_nh,' ! segment discretisation in x and y, nwinc nhinc '
+      write(fh2_input_file_unit,'(2ES12.4,A)')FH2_rw,FH2_rh,'  ! segment discretisation ratio in x and y, rw rh '
     else if (PUL%shape(conductor).EQ.rectangle) then
@@ -195,12 +232,12 @@ integer :: ierr
   write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%fmin,'   ! Minimum frequency, fmin'
   write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%fmax,'   ! Maximum frequency, fmax'
-  write(fh2_input_file_unit,'(I8,A)')ndec,'   ! Number of points per decade, ndec'
+  write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%ndec,'   ! Number of points per decade, ndec'
   CLOSE(unit=fh2_input_file_unit)
   write(*,*)'CALLING FastHenry2'
-  command='/home/chris/SACAMOS_FASTHENRY2/FASTHENRY2_WRITE_INPUT_FILE/bin/write_FH_input_file < fh2.txt'
+  command='write_FH_input_file < fh2.txt'
   CALL execute_command_line(command,EXITSTAT=exit_stat)
 ! Check that the FastHenry2 input file generated correctly
@@ -306,6 +343,7 @@ integer :: ierr
       if (verbose) write(*,*)'Number of frequencies=',n_freq_in
       ALLOCATE( freq_in(n_freq_in) )
       ALLOCATE( Zc_in(n_freq_in,nrows,ncols) )
+      ALLOCATE( Rdc_domain(nrows,ncols) )
       ALLOCATE( values(2*ncols) )
       rewind(unit=fh2_output_file_unit)
     end if
@@ -313,24 +351,102 @@ integer :: ierr
   end do ! next file read loop
   CLOSE(unit=fh2_output_file_unit)
+  
+  matrix_dimension=PUL%n_conductors-1
+  
+  if (matrix_dimension.NE.nrows) then
+    write(run_status,*)'ERROR in PUL_RL_FastHenry2 matrix dimension',matrix_dimension,nrows
+    CALL write_program_status()
+    STOP 1
+  end if
-! For now use the inductance from the highest frequency calculated 
+! For the modal decomposition, use the inductance from the highest frequency calculated 
   do r=1,nrows
     do c=1,ncols
       freq=freq_in(n_freq_in)
       w=2.0*pi*freq
-      PUL%L%mat(r,c)=real(Zc_in(n_freq_in,r,c)/(j*w))
-      PUL%Zfilter%sfilter_mat(r,c)=jwA_filter( PUL%L%mat(r,c) )
+      PUL%L%mat(r,c)=real(Zc_in(n_freq_in,r,c)/(j*w))      
-      write(*,*)r,c,PUL%L%mat(r,c)
+    end do
+  end do
+
+! Save the resistance matrix from the lowest frequency calculated 
+  do r=1,nrows
+    do c=1,ncols
+    
+      Rdc_domain(r,c)=real(Zc_in(1,r,c))      
     end do
   end do
+  
+! We calculate the inductance and resistance matrix at a number of frequencies before fitting filter functions
+! to each of the frequency dependent impedance matrix entries.
+  
+  write(*,*)'**************************'
+  write(*,*)'Process FastHenry2 results'
+  write(*,*)'**************************'
+  write(*,*)'Number of frequencies (freq_spec) =',freq_spec%n_frequencies
+  write(*,*)'Number of frequencies (FastHenry2)=',n_freq_in
+
+! Frequency check  
+  if (n_freq_in.NE.freq_spec%n_frequencies) then
+    write(run_status,*)'ERROR in PUL_RL_FastHenry2 n_frequencies specification',n_freq_in,freq_spec%n_frequencies
+    CALL write_program_status()
+    STOP 1
+  end if
+  
+  do floop=1,n_freq_in
+    ferr=abs(freq_spec%freq_list(floop)-freq_in(floop))/(freq_spec%freq_list(floop)+freq_in(floop))
+    if (ferr.GT.0.001) then
+      write(run_status,*)'ERROR in PUL_RL_FastHenry2 frequency samples, floop=',floop, &
+                         ' freq_spec=',freq_spec%freq_list(floop),' freq_FH2=',freq_in(floop)
+      CALL write_program_status()
+      STOP 1   
+    end if
+  end do
+ 
+! loop over frequency    
+  if (verbose) then
+    do floop=1,n_freq_in
+      write(*,*)freq_spec%freq_list(floop),freq_in(floop)
+    end do
+  end if
+  
+! we have frequency domain impedance matrix, Zc
+
+! loop over the elements of Zc and the Zfilter matrix by filter fitting 
+  ALLOCATE( function_to_fit(1:freq_spec%n_frequencies) )
+    
+  do row=1,matrix_dimension
+    do col=row,matrix_dimension
+
+! get the function values for this matrix element function_to_fit=R-Rdc+jwL
+      do floop=1,freq_spec%n_frequencies
+        function_to_fit(floop)=Zc_in(floop,row,col)-Rdc_domain(row,col)
+      end do
+        
+! calculate the Zfilter matrix using the filter fitting process
+! note aorder=border and no restrictions are put on the function
+
+      CALL Calculate_Sfilter(function_to_fit,freq_spec%freq_list,freq_spec%n_frequencies,  &
+                             PUL%Zfilter%sfilter_mat(row,col),fit_order+1,-1,2)    ! note type 2 filter fit=> a0=0.0
+
+      if (col.NE.row) then       
+! make the matrix symmetrical
+        PUL%Zfilter%sfilter_mat(col,row)=PUL%Zfilter%sfilter_mat(row,col)
+      end if
+        
+    end do ! next col
+      
+  end do ! next row
+    
+  DEALLOCATE( function_to_fit )  
   DEALLOCATE( freq_in )
   DEALLOCATE( Zc_in )
+  DEALLOCATE( Rdc_domain )
   DEALLOCATE( values )
   if (verbose) then
@@ -85,6 +85,7 @@ TYPE::PUL_type
 ! conductor based arrays:
   integer,allocatable         :: shape(:)        ! holds a nuber which indicates the shape of the conductor
+  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
@@ -104,7 +105,14 @@ TYPE::PUL_type
   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
@@ -184,6 +192,7 @@ include &quot;CG_solve.F90&quot;
 ! allocate memory for PUL structure
   PUL%n_conductors=n_conductors
   ALLOCATE( PUL%shape( 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) ) 
@@ -204,6 +213,7 @@ include &quot;CG_solve.F90&quot;
   do i=1,PUL%n_conductors
     PUL%shape(i)=0
+    PUL%sigma(i)=0d0 
     PUL%x(i)=0d0 
     PUL%y(i)=0d0 
     PUL%r(i)=0d0 
@@ -225,7 +235,14 @@ include &quot;CG_solve.F90&quot;
   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
@@ -272,6 +289,7 @@ include &quot;CG_solve.F90&quot;
 ! START
   if (allocated(PUL%shape )    ) DEALLOCATE( PUL%shape )
+  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 )
@@ -39,9 +39,11 @@
 ! If the model order is negative on input then the 'optimium' model  is returned
 ! between model order 0 and abs(model_order)
 !
-! Two types of model are considered - s general model in which no constraints are placed on the filter function
-! and a model for the propagation correction in which the constraint that the value of the filter function at
+! Three types of model are considered - 
+! 1: a general model in which no constraints are placed on the filter function
+! 2: a model for the propagation correction in which the constraint that the value of the filter function at
 ! zero frequency is equal to 1.0
+! 3: a model for the Z-Rdc matrices where the value of the filter function at zero frequency is equal to 0.0
 !     
 ! COMMENTS
 !     The stopping criterion for calculation of the optimum model is an inflection in the function log_error(model_order)
@@ -53,6 +55,7 @@
 !                    return the order model specified.
 !     14/10/2017 CJS include a check for unstable filters and add an absolute convergence check
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
+!     25/10/2023 CJS Add a mode fit type where a0=0.0 for fitting to Z-Rdc matrices 
 !
 !
 SUBROUTINE Calculate_Sfilter(Hp,f,n_frequencies,Hfilter,model_order,relative_border,fit_type)
@@ -75,6 +78,7 @@ integer,intent(IN)                 :: model_order             ! Numerator model 
 integer,intent(IN)                 :: relative_border         ! Denominator model order to calculate relative to numerator order, border=aorder+relative_border
 integer,intent(IN)                 :: fit_type                ! Type of model fit. fit_type=0 has no restrictions on f=0 value
                                                               !                    fit_type=1 imposes Hfilter(f=0)=1.0
+                                                              !                    fit_type=2 imposes Hfilter(f=0)=0.0
 ! local variables
   integer  :: order
@@ -38,9 +38,11 @@
 ! as described in Dawson paper [**** REFERENCE REQUIRED ****]
 !     
 ! COMMENTS
+!     If fit_type=0 then we fit directly to the given function values.
 !     If fit_type=1 then we impose Hfilter(f=0)=1.0. In this case we fit to a finction Hp-1 to
 !     calculate a filter function with a(0)=0, then add 1 to this filter to give the result.
-!     If fit_type=0 then we fit directly to the given function values.
+!     If fit_type=2 then we impose Hfilter(f=0)=0.0. In this case we fit to a finction Hp to
+!     calculate a filter function with a(0)=0.
 !
 ! HISTORY
 !
@@ -121,6 +123,12 @@ integer,intent(IN)                 :: fit_type                ! Type of model fi
       Hfilter%a%coeff(0)=1d0
       RETURN
+   
+    else if (fit_type.EQ.2) then   ! this type has the requirement that Hfilter(f=0)=1.0 so return  Hfilter=1 
+
+      Hfilter%a%coeff(0)=0d0
+      
+      RETURN
     else
 ! Hfilter is effectively a real constant and can be calculated as the average of the real part of the input function values
@@ -146,10 +154,17 @@ integer,intent(IN)                 :: fit_type                ! Type of model fi
     coeff_dim=aorder+1+border 
     first_a=0
     N_a=aorder+1
-  else                     ! Hfilter(f=0)=1.0 thus a(0)=1.0 and b(0)=1.0
+  else if (fit_type.EQ.1) then  ! Hfilter(f=0)=1.0 thus a(0)=1.0 and b(0)=1.0
     coeff_dim=aorder+border     
     first_a=1
     N_a=aorder
+  else if (fit_type.EQ.2) then  ! Hfilter(f=0)=0.0 thus a(0)=0.0 and b(0)=1.0
+    coeff_dim=aorder+border     
+    first_a=1
+    N_a=aorder
+  else
+    write(*,*)'ERROR in Weiner_Hopf: fit_type should be 0, 1 or 2, here fit_type=',fit_type
+    STOP 1
   end if
 ! Allocate memory for the matrices and vectors used in the calculations
@@ -255,6 +270,12 @@ integer,intent(IN)                 :: fit_type                ! Type of model fi
       row=i
       Hfilter%a%coeff(i)=VC(row)+Hfilter%b%coeff(i)    
     end do
+  else if (fit_type.eq.2) then            ! we have a fit to Hp(f) so calculate the coefficients of H=A(jomega)/B(jomega)
+    Hfilter%a%coeff(0)=0d0
+    do i=1,aorder
+      row=i
+      Hfilter%a%coeff(i)=VC(row)   
+    end do
   else
     do i=0,aorder
       row=i+1
@@ -0,0 +1,22 @@
+    if (conductor_data(conductor)%auto_grid_density) then
+      required_dl=conductor_data(conductor)%skin_depth         ! could add reduction factor here e.g. /2.0          
+      conductor_data(conductor)%dl=required_dl
+    end if
+
+    conductor_data(conductor)%nxmax=NINT(conductor_data(conductor)%grid_dim/conductor_data(conductor)%dl)+1
+    conductor_data(conductor)%nxmin=-conductor_data(conductor)%nxmax
+
+    conductor_data(conductor)%nymax=conductor_data(conductor)%nxmax
+    conductor_data(conductor)%nymin=conductor_data(conductor)%nxmin
+
+    nxmin=conductor_data(conductor)%nxmin
+    nxmax=conductor_data(conductor)%nxmax
+    nymin=conductor_data(conductor)%nymin
+    nymax=conductor_data(conductor)%nymax
+        
+    ALLOCATE( conductor_data(conductor)%grid(nxmin:nxmax,nymin:nymax) )
+    ALLOCATE( conductor_data(conductor)%depth(nxmin:nxmax,nymin:nymax) )
+
+! Reset the grid   
+    conductor_data(conductor)%grid(nxmin:nxmax,nymin:nymax)=0
+    conductor_data(conductor)%depth(nxmin:nxmax,nymin:nymax)=0.0
@@ -0,0 +1,16 @@
+    conductor_data(conductor)%mesh_type=mesh_type_layer  ! default 
+    conductor_data(conductor)%mesh_to_layer_type=1
+
+    if ( INDEX(line_string,'auto').NE.0) conductor_data(conductor)%auto_grid_density=.TRUE.
+          
+    if ( INDEX(line_string,'grid_layer').NE.0) then
+    
+      conductor_data(conductor)%mesh_type=mesh_type_grid   
+      conductor_data(conductor)%mesh_to_layer_type=2
+      
+    else if ( INDEX(line_string,'grid').NE.0) then
+    
+      conductor_data(conductor)%mesh_type=mesh_type_grid   
+      conductor_data(conductor)%mesh_to_layer_type=1
+            
+    end if
@@ -0,0 +1,90 @@
+
+    nxmin=conductor_data(conductor)%nxmin
+    nxmax=conductor_data(conductor)%nxmax
+    nymin=conductor_data(conductor)%nymin
+    nymax=conductor_data(conductor)%nymax
+    
+    layer_number=0
+      
+    if (conductor_data(conductor)%mesh_to_layer_type.EQ.1) then
+! each marked cell in the grid becomes a layer and hence a segment
+    
+      do ix=nxmin,nxmax
+        do iy=nymin,nymax
+      if ( conductor_data(conductor)%grid(ix,iy).EQ.1 ) then
+    
+        layer_number=layer_number+1
+      
+        apx=conductor_data(conductor)%xc+real(ix)*conductor_data(conductor)%dl
+        apy=conductor_data(conductor)%yc+real(iy)*conductor_data(conductor)%dl
+    
+            conductor_data(conductor)%x(layer_number)=apx
+            conductor_data(conductor)%y(layer_number)=apy
+            conductor_data(conductor)%w(layer_number)=conductor_data(conductor)%dl
+            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 if
+        end do
+      end do
+
+    else if (conductor_data(conductor)%mesh_to_layer_type.EQ.2) then
+    
+! each row of cells in the grid becomes a segment
+    do iy=nymin,nymax
+    
+      in_conductor=.FALSE.
+      cx1=0
+      cx2=0
+      
+      do ix=nxmin,nxmax
+      
+        if ( (.NOT.in_conductor).AND.(conductor_data(conductor)%grid(ix,iy).EQ.1) ) then
+    
+! we have moved from air to conductor so save this x cell number
+           cx1=ix
+       in_conductor=.TRUE.
+       
+    else if ( (in_conductor).AND.(conductor_data(conductor)%grid(ix,iy).EQ.0) ) then
+! we have moved from air to conductor so save this x cell number
+
+           cx2=ix-1
+       in_conductor=.FALSE.
+       
+       layer_number=layer_number+1
+       
+       cw=real(cx2-cx1+1)*conductor_data(conductor)%dl
+       ch=conductor_data(conductor)%dl
+       
+       apx=conductor_data(conductor)%xc+real(cx2+cx1)/2.0
+       apy=conductor_data(conductor)%yc+real(iy)*conductor_data(conductor)%dl
+    
+           conductor_data(conductor)%x(layer_number)=apx
+           conductor_data(conductor)%y(layer_number)=apy
+           conductor_data(conductor)%w(layer_number)=cw
+           conductor_data(conductor)%h(layer_number)=ch
+               
+           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
+
+           cx1=0
+           cx2=0
+    
+    end if
+    
+      end do ! next col
+
+    end do ! next row
+
+    else
+      write(*,*)'Unknown mesh_to_layer_type:',conductor_data(conductor)%mesh_to_layer_type, &
+                'in conductor',conductor
+      STOP 1
+    end if
+    
+! the number of layers may have changed if we have combined grid cells      
+    conductor_data(conductor)%tot_n_layers=layer_number
@@ -111,6 +111,7 @@ type(bundle_specification_type)    ::bundle_spec
 logical    :: file_exists
 character(len=line_length)    :: line
+character(len=line_length)    :: stripped_line
 integer        :: ierr,ierr2   ! integers to return error codes from file reads
@@ -246,7 +247,8 @@ logical    :: must_use_laplace
   bundle_spec%n_cables_without_ground_plane=bundle_spec%n_cables
-  if (line.eq.'ground_plane') then
+!  if (line.eq.'ground_plane') then
+  if (index(line,'ground_plane').EQ.1) then
     bundle_spec%ground_plane_present=.TRUE.
@@ -272,13 +274,49 @@ logical    :: must_use_laplace
     bundle_spec%cable_x_offset(cable)=bundle_spec%ground_plane_x
     bundle_spec%cable_y_offset(cable)=bundle_spec%ground_plane_y
     bundle_spec%cable_angle(cable)=bundle_spec%ground_plane_theta
+    
+! Attempt to read additional parameters for FastHenry2 calculation, sigma, width, thickness
+    stripped_line=line(13:line_length)
+    read(stripped_line,*,ERR=100,END=100)bundle_spec%ground_plane_sigma &
+                                        ,bundle_spec%ground_plane_w     &
+                                        ,bundle_spec%ground_plane_h     &     
+                                        ,bundle_spec%ground_plane_nsegx     &     
+                                        ,bundle_spec%ground_plane_nsegz     &     
+                                        ,bundle_spec%ground_plane_nh        
+  
+    bundle_spec%ground_plane_Rdc=1d0/(bundle_spec%ground_plane_sigma*bundle_spec%ground_plane_w*bundle_spec%ground_plane_h)
 ! create conductor impedance model for the new cable    
     ALLOCATE(bundle_spec%cable(cable)%conductor_impedance(1:1))
+    bundle_spec%cable(cable)%conductor_impedance(1)%impedance_model_type=impedance_model_type_FH2
+    bundle_spec%cable(cable)%conductor_impedance(1)%Resistance_multiplication_factor=1d0
+    bundle_spec%cable(cable)%conductor_impedance(1)%radius=0d0
+    bundle_spec%cable(cable)%conductor_impedance(1)%conductivity=bundle_spec%ground_plane_sigma
+    bundle_spec%cable(cable)%conductor_impedance(1)%width=bundle_spec%ground_plane_w
+    bundle_spec%cable(cable)%conductor_impedance(1)%height=bundle_spec%ground_plane_h
+    bundle_spec%cable(cable)%conductor_impedance(1)%Rdc=bundle_spec%ground_plane_Rdc
+  
+  GOTO 200
+
+100 CONTINUE ! jump here if we have failed to read FastHenry2 parameters
+! create conductor impedance model for the new cable    
+    ALLOCATE(bundle_spec%cable(cable)%conductor_impedance(1:1))
     bundle_spec%cable(cable)%conductor_impedance(1)%impedance_model_type=impedance_model_type_PEC
     bundle_spec%cable(cable)%conductor_impedance(1)%Resistance_multiplication_factor=1d0
     bundle_spec%cable(cable)%conductor_impedance(1)%radius=0d0
     bundle_spec%cable(cable)%conductor_impedance(1)%conductivity=0d0
+    bundle_spec%cable(cable)%conductor_impedance(1)%width=0d0
+    bundle_spec%cable(cable)%conductor_impedance(1)%height=0d0
+    bundle_spec%cable(cable)%conductor_impedance(1)%Rdc=0d0
+    bundle_spec%ground_plane_sigma=0d0
+    bundle_spec%ground_plane_w=0d0
+    bundle_spec%ground_plane_h=0d0
+    bundle_spec%ground_plane_Rdc=0d0
+    bundle_spec%ground_plane_nsegx=0
+    bundle_spec%ground_plane_nsegz=0
+    bundle_spec%ground_plane_nh=0
+
+200 CONTINUE ! jump here if we do have FastHenry2 parameters
   else if (line.eq.'no_ground_plane') then
@@ -311,14 +349,14 @@ logical    :: must_use_laplace
   if (ierr.NE.0) then
 ! Assume there is no filter fit information specified so move on to the next stage
     backspace(bundle_spec_file_unit)
-    goto 100
+    goto 300
   end if
   write(*,*)'Reading the filter fit frequency range'
   CALL read_and_set_up_frequency_specification(bundle_spec%Y_fit_freq_spec,bundle_spec_file_unit)
-100 continue
+300 continue
 ! the file can contain flags to control the running of the software and the output
   rewind(bundle_spec_file_unit)
@@ -326,7 +364,7 @@ logical    :: must_use_laplace
   write(*,*)'Processing flags'
   do
-    read(bundle_spec_file_unit,*,END=110,ERR=110)line
+    read(bundle_spec_file_unit,'(A)',END=310,ERR=310)line
     CALL convert_to_lower_case(line,line_length)
 ! Set flags according to the information at the end of the .spice_model_spec file
@@ -336,8 +374,7 @@ logical    :: must_use_laplace
     if (INDEX(line,'no_s_xfer').NE.0) use_s_xfer=.FALSE.
     if (INDEX(line,'use_laplace').NE.0) use_Laplace=.TRUE.
     if (INDEX(line,'no_laplace').NE.0) use_Laplace=.FALSE.
-    if (INDEX(line,'use_fasthenry').NE.0) use_FastHenry=.TRUE.
-    if (INDEX(line,'no_fasthenry').NE.0) use_FastHenry=.FALSE.
+        
     if (INDEX(line,'plot_potential').NE.0) plot_potential=.TRUE.
     if (INDEX(line,'no_plot_potential').NE.0) plot_potential=.FALSE.
     if (INDEX(line,'plot_mesh').NE.0) plot_mesh=.TRUE.
@@ -376,10 +413,24 @@ logical    :: must_use_laplace
     if (INDEX(line,'cg_tol').NE.0) then
       read(bundle_spec_file_unit,*,END=9000,ERR=9000)cg_tol
     end if
+    
+    if (INDEX(line,'no_fasthenry').NE.0) use_FastHenry=.FALSE.
+    
+    if (INDEX(line,'use_fasthenry').EQ.1) then
+      use_FastHenry=.TRUE.
+    
+      stripped_line=line(14:line_length)
+      read(stripped_line,*,ERR=305,END=305)FH2_nlayers_radius,FH2_nw,FH2_nh,FH2_rw,FH2_rh
+    
+      write(*,*)'use_FastHenry=.TRUE.'
+      write(*,*)'Parameters:',FH2_nlayers_radius,FH2_nw,FH2_nh,FH2_rw,FH2_rh
+305   CONTINUE
+    end if
+    
   end do  ! continue until all flags are read - indicated by an end of file.
-110 CONTINUE
+310 CONTINUE
 ! close the bundle_spec file
@@ -88,7 +88,7 @@
 !     24/2/2017 CJS Allow the input name to include a path i.e. the _spec file does not need to be local.
 !     28/2/2017 CJS Make the validation test circuit model optional
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
-!
+!     24/11/2023 CJS Start to embed FastHenry2 frequency dependent impedance matrix 
 !
 PROGRAM spice_cable_bundle_model_builder
@@ -574,6 +574,8 @@ integer    :: i
     if (INDEX(line,'use_analytic_i').NE.0) run_validation_test_Vbased=.FALSE.
     if (INDEX(line,'use_analytic_v').NE.0) run_validation_test_Vbased=.TRUE.
+
+    if (INDEX(line,'plot_propagation_correction_filter_fit_data').NE.0) plot_propagation_correction_filter_fit_data=.TRUE.
     if (INDEX(line,'min_delay').NE.0) then
@@ -0,0 +1,1020 @@
+PROGRAM write_FH_input_file
+
+USE type_specifications
+
+IMPLICIT NONE
+
+integer :: n_conductors
+
+logical :: ground_plane
+
+integer :: mesh_type
+
+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
+
+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
+
+TYPE::conductor_type
+
+  integer :: type
+  integer :: mesh_type
+  integer :: mesh_to_layer_type
+
+  real(dp)    :: rc
+  real(dp)    :: rci
+  real(dp)    :: rco
+  real(dp)    :: xc
+  real(dp)    :: yc
+  real(dp)    :: width
+  real(dp)    :: height
+  integer :: n_layers2
+  integer :: tot_n_layers
+  
+  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
+  real(dp)    :: dl
+  integer :: nwinc 
+  integer :: nhinc 
+  real(dp)    :: rw
+  real(dp)    :: rh
+
+  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 :: mesh_type_layer=1
+integer,parameter :: mesh_type_grid=2
+
+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
+
+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.'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
+  
+  if (conductor_data(conductor)%auto_grid_density) then ! uniform grid in segments
+    conductor_data(conductor)%rw=1.0
+    conductor_data(conductor)%rh=1.0
+  end if
+
+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
+
+open(unit=20,file=trim(FH2_filename))
+
+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
+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
+
+! 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
+    
+      conductor_data(conductor)%n_layers2=0
+
+    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)),')'
+                       
+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"
+     
+    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)/2.0
+    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)%nhinc*conductor_data(conductor)%nwinc
+    
+  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 conductor 1 to all other conductors in turn'
+  
+  do conductor=2,n_conductors
+
+    write(20,'(A)')'*'
+    write(20,'(A,I4)')'* Conductor 1 to conductor',conductor
+  
+    write(conductor_string,'(I4)')conductor-1
+    loop_string='loop_'//trim(adjustl(conductor_string))   
+  
+    write(20,'(A,A,A,A,A,A)')'.external ',trim(conductor_data(1)%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))
+  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
+  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
+