From ca8ab9e90b29cfc11d3d1b8358076c5512aa66f5 Mon Sep 17 00:00:00 2001 From: chris.smartt Date: Wed, 25 Oct 2023 16:50:59 +0100 Subject: [PATCH] Update to proximity effects: working for simple cases --- PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE/single_wire.cable_spec | 2 +- PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/single_wire.cable_spec | 2 +- PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire.bundle_spec | 2 +- PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire_ac.spice_model_spec | 5 +++++ PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_LAPLACE/single_wire.cable_spec | 2 +- PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/plot_comparison.plt | 1 + PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/status | 6 ++---- PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/single_wire.cable_spec | 2 +- PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground.bundle_spec | 7 +++---- PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground_ac.spice_model_spec | 7 ++++++- PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/status | 6 ++---- SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 | 36 ++++++++++++++++++++++++++++++------ SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 | 31 +++++++++++++++++++++++++++++++ SRC/CABLE_MODULES/cable_module.F90 | 7 +++++++ SRC/CABLE_MODULES/coax.F90 | 1 + SRC/CABLE_MODULES/conductor_impedance_model.F90 | 16 ++++++++++++++++ SRC/CABLE_MODULES/cylindrical.F90 | 1 + SRC/CREATE_SPICE_CIRCUIT_MODEL/create_spice_subcircuit_model.F90 | 2 +- SRC/GENERAL_MODULES/frequency_spec.F90 | 5 +++++ SRC/GENERAL_MODULES/general_module.F90 | 8 ++++++++ SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis.F90 | 4 ++-- SRC/MTL_ANALYTIC_SOLUTION/propagation_correction_filters.F90 | 18 +++++++++++------- SRC/Makefile | 9 ++++++++- SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 | 7 ++----- SRC/PUL_PARAMETER_CALCULATION/PUL_RL_FastHenry2.F90 | 170 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------- SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 | 22 ++++++++++++++++++++-- SRC/SFILTER_FIT/Calculate_Sfilter.F90 | 8 ++++++-- SRC/SFILTER_FIT/Weiner_Hopf.F90 | 25 +++++++++++++++++++++++-- SRC/WRITE_FH2_IPFILE/create_grid.F90 | 22 ++++++++++++++++++++++ SRC/WRITE_FH2_IPFILE/get_grid_type.F90 | 16 ++++++++++++++++ SRC/WRITE_FH2_IPFILE/set_segments_from_grid.F90 | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/cable_bundle_model_builder.F90 | 65 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------- SRC/spice_cable_bundle_model_builder.F90 | 4 +++- SRC/write_FH_input_file.F90 | 1020 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 34 files changed, 1548 insertions(+), 81 deletions(-) create mode 100644 SRC/WRITE_FH2_IPFILE/create_grid.F90 create mode 100644 SRC/WRITE_FH2_IPFILE/get_grid_type.F90 create mode 100644 SRC/WRITE_FH2_IPFILE/set_segments_from_grid.F90 create mode 100644 SRC/write_FH_input_file.F90 diff --git a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE/single_wire.cable_spec b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE/single_wire.cable_spec index 238da3e..072de6e 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE/single_wire.cable_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE/single_wire.cable_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/single_wire.cable_spec b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/single_wire.cable_spec index 238da3e..072de6e 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/single_wire.cable_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/single_wire.cable_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire.bundle_spec b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire.bundle_spec index dc8a1a2..349628d 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire.bundle_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire.bundle_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire_ac.spice_model_spec b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire_ac.spice_model_spec index 2859f3a..31aa331 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire_ac.spice_model_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_FASTHENRY/two_wire_ac.spice_model_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_LAPLACE/single_wire.cable_spec b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_LAPLACE/single_wire.cable_spec index 238da3e..072de6e 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_LAPLACE/single_wire.cable_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/TWO_WIRE_LAPLACE/single_wire.cable_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/plot_comparison.plt b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/plot_comparison.plt index b275957..61b0469 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/plot_comparison.plt +++ b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/plot_comparison.plt @@ -5,6 +5,7 @@ set xlabel 'Frequency(Hz)' 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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/status b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/status index 3a85b9c..61ce6a6 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/status +++ b/PROXIMITY_EFFECT_TEST_CASES/TWO_WIRE/status @@ -1,5 +1,3 @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/single_wire.cable_spec b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/single_wire.cable_spec index 238da3e..b711849 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/single_wire.cable_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/single_wire.cable_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground.bundle_spec b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground.bundle_spec index 016ff44..4096a26 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground.bundle_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground.bundle_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground_ac.spice_model_spec b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground_ac.spice_model_spec index 396329f..d50e661 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground_ac.spice_model_spec +++ b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/WIRE_OVER_GROUND_PLANE_FASTHENRY/wire_over_ground_ac.spice_model_spec @@ -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 diff --git a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/status b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/status index fd314e8..9682e1e 100644 --- a/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/status +++ b/PROXIMITY_EFFECT_TEST_CASES/WIRE_OVER_GP/status @@ -1,5 +1,3 @@ -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 diff --git a/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 b/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 index dc20a81..965f499 100644 --- a/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 +++ b/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 @@ -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 diff --git a/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 b/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 index 9312735..754423e 100644 --- a/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 +++ b/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 @@ -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 diff --git a/SRC/CABLE_MODULES/cable_module.F90 b/SRC/CABLE_MODULES/cable_module.F90 index db7c2e7..62f09e7 100644 --- a/SRC/CABLE_MODULES/cable_module.F90 +++ b/SRC/CABLE_MODULES/cable_module.F90 @@ -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 'conductor_impedance_model.F90' 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 ' diff --git a/SRC/CABLE_MODULES/coax.F90 b/SRC/CABLE_MODULES/coax.F90 index 8128d74..adacf20 100644 --- a/SRC/CABLE_MODULES/coax.F90 +++ b/SRC/CABLE_MODULES/coax.F90 @@ -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 diff --git a/SRC/CABLE_MODULES/conductor_impedance_model.F90 b/SRC/CABLE_MODULES/conductor_impedance_model.F90 index 628a8df..7926d18 100644 --- a/SRC/CABLE_MODULES/conductor_impedance_model.F90 +++ b/SRC/CABLE_MODULES/conductor_impedance_model.F90 @@ -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 diff --git a/SRC/CABLE_MODULES/cylindrical.F90 b/SRC/CABLE_MODULES/cylindrical.F90 index ce2ec88..5a58201 100644 --- a/SRC/CABLE_MODULES/cylindrical.F90 +++ b/SRC/CABLE_MODULES/cylindrical.F90 @@ -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 diff --git a/SRC/CREATE_SPICE_CIRCUIT_MODEL/create_spice_subcircuit_model.F90 b/SRC/CREATE_SPICE_CIRCUIT_MODEL/create_spice_subcircuit_model.F90 index 287e840..5b99cde 100644 --- a/SRC/CREATE_SPICE_CIRCUIT_MODEL/create_spice_subcircuit_model.F90 +++ b/SRC/CREATE_SPICE_CIRCUIT_MODEL/create_spice_subcircuit_model.F90 @@ -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. diff --git a/SRC/GENERAL_MODULES/frequency_spec.F90 b/SRC/GENERAL_MODULES/frequency_spec.F90 index fa9f025..c800f3a 100644 --- a/SRC/GENERAL_MODULES/frequency_spec.F90 +++ b/SRC/GENERAL_MODULES/frequency_spec.F90 @@ -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 diff --git a/SRC/GENERAL_MODULES/general_module.F90 b/SRC/GENERAL_MODULES/general_module.F90 index 74faeaa..4dc781a 100644 --- a/SRC/GENERAL_MODULES/general_module.F90 +++ b/SRC/GENERAL_MODULES/general_module.F90 @@ -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' diff --git a/SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis.F90 b/SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis.F90 index 5d11e68..c99e897 100644 --- a/SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis.F90 +++ b/SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis.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]' diff --git a/SRC/MTL_ANALYTIC_SOLUTION/propagation_correction_filters.F90 b/SRC/MTL_ANALYTIC_SOLUTION/propagation_correction_filters.F90 index 22dc9ce..ed1ef46 100644 --- a/SRC/MTL_ANALYTIC_SOLUTION/propagation_correction_filters.F90 +++ b/SRC/MTL_ANALYTIC_SOLUTION/propagation_correction_filters.F90 @@ -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) diff --git a/SRC/Makefile b/SRC/Makefile index 4bc788e..f23c582 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -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 ) diff --git a/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 b/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 index 5e7caa0..04a86b2 100644 --- a/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 +++ b/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 @@ -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 ) diff --git a/SRC/PUL_PARAMETER_CALCULATION/PUL_RL_FastHenry2.F90 b/SRC/PUL_PARAMETER_CALCULATION/PUL_RL_FastHenry2.F90 index 3f4b502..ea26757 100644 --- a/SRC/PUL_PARAMETER_CALCULATION/PUL_RL_FastHenry2.F90 +++ b/SRC/PUL_PARAMETER_CALCULATION/PUL_RL_FastHenry2.F90 @@ -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 diff --git a/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 b/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 index 314a754..adc6b04 100644 --- a/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 +++ b/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 @@ -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 "CG_solve.F90" ! 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 "CG_solve.F90" 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 "CG_solve.F90" 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 "CG_solve.F90" ! 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 ) diff --git a/SRC/SFILTER_FIT/Calculate_Sfilter.F90 b/SRC/SFILTER_FIT/Calculate_Sfilter.F90 index 811aefc..374bef8 100644 --- a/SRC/SFILTER_FIT/Calculate_Sfilter.F90 +++ b/SRC/SFILTER_FIT/Calculate_Sfilter.F90 @@ -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 diff --git a/SRC/SFILTER_FIT/Weiner_Hopf.F90 b/SRC/SFILTER_FIT/Weiner_Hopf.F90 index 031129a..0347ce3 100644 --- a/SRC/SFILTER_FIT/Weiner_Hopf.F90 +++ b/SRC/SFILTER_FIT/Weiner_Hopf.F90 @@ -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 diff --git a/SRC/WRITE_FH2_IPFILE/create_grid.F90 b/SRC/WRITE_FH2_IPFILE/create_grid.F90 new file mode 100644 index 0000000..fedcd71 --- /dev/null +++ b/SRC/WRITE_FH2_IPFILE/create_grid.F90 @@ -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 diff --git a/SRC/WRITE_FH2_IPFILE/get_grid_type.F90 b/SRC/WRITE_FH2_IPFILE/get_grid_type.F90 new file mode 100644 index 0000000..f44ed19 --- /dev/null +++ b/SRC/WRITE_FH2_IPFILE/get_grid_type.F90 @@ -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 diff --git a/SRC/WRITE_FH2_IPFILE/set_segments_from_grid.F90 b/SRC/WRITE_FH2_IPFILE/set_segments_from_grid.F90 new file mode 100644 index 0000000..960bb0f --- /dev/null +++ b/SRC/WRITE_FH2_IPFILE/set_segments_from_grid.F90 @@ -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 diff --git a/SRC/cable_bundle_model_builder.F90 b/SRC/cable_bundle_model_builder.F90 index d3311ea..217dff1 100644 --- a/SRC/cable_bundle_model_builder.F90 +++ b/SRC/cable_bundle_model_builder.F90 @@ -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 diff --git a/SRC/spice_cable_bundle_model_builder.F90 b/SRC/spice_cable_bundle_model_builder.F90 index 9fd49b0..46f58dd 100644 --- a/SRC/spice_cable_bundle_model_builder.F90 +++ b/SRC/spice_cable_bundle_model_builder.F90 @@ -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 diff --git a/SRC/write_FH_input_file.F90 b/SRC/write_FH_input_file.F90 new file mode 100644 index 0000000..1ec161d --- /dev/null +++ b/SRC/write_FH_input_file.F90 @@ -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 + -- libgit2 0.21.2