! ! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. ! It was developed by the University of Nottingham and the Netherlands Aerospace ! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK. ! ! Copyright (C) 2016-2018 University of Nottingham ! ! SACAMOS is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! SACAMOS is distributed in the hope that it will be useful, but ! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY ! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License ! for more details. ! ! A copy of the GNU General Public License version 3 can be found in the ! file GNU_GPL_v3 in the root or at . ! ! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to ! the GNU Lesser General Public License. A copy of the GNU Lesser General Public ! License version can be found in the file GNU_LGPL in the root of EISPACK ! (/SRC/EISPACK ) or at . ! ! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk ! ! ! File Contents: ! SUBROUTINE time_domain_analysis ! ! NAME ! time_domain_analysis ! ! AUTHORS ! Chris Smartt ! ! DESCRIPTION ! This subroutine controls the analytic solution for the transient analysis of ! multi-conductor transmission lines. The solution is obtained using the Inverse ! Fourier Transform (IFT) method. ! The solution is obtained using the full dimension transmission line equations ! i.e. we are NOT using the weak form of transfer impedance coupling ! Note also that frequency dependent quantities are evaluated separately at ! each frequency of analysis, i.e. the frequency dependence of the solution is rigorous ! given only the frequency dependence of the dielectrics is modelled using impedance/ admittance ! matrices whose elements are rational frequency dependent filter functions. ! ! INPUTS: ! spice_bundle_model structure ! spice_validation_test structure ! ! OUTPUT ! analytic time domain termination voltage for the specified validation test case written to file ! ! Note that the analysis assumes that the excitation (and hence the response) are periodic ! with period equal to the maximimum time of the simulation. If the response of the circuit ! is longer than this then aliasing will occur and the comparison between Spice model results ! and analytic results will be in error so please ensure that all transients have reduced ! to insignificant levels within the simulation time. ! ! COMMENTS ! STAGE_1: frequency independent parameter solution ! STAGE_2: multi-conductor solution ! STAGE_3: shielded cable solution ! STAGE_4: frequency dependent model ! STAGE_5: transfer impedance model ! ! HISTORY ! ! started 7/12/2015 CJS: STAGE_1 developments ! 24/03/2016 CJS: STAGE_3 developments -shielded cables ! 22/04/2016 CJS: STAGE_4 developments -frequency dependent model ! 04/05/2016 CJS: Only write time domain output over the time period requested,not the full FFT time period. ! Include general conductor impedance model 12/05/2016 CJS ! Fix bug with conductor impedance contributions 12/05/2016 CJS ! 8/09/2016 CJS Correct the common mode/ differential mode loss terms for twisted pairs ! 13/10/2016 CJS Correct transfer impedance for multiple modes in external domain ! 7/3/2017 CJS: Add resistance and voltage source onto the reference coonductor ! 8/5/2017 CJS: Include references to Theory_Manual ! ! SUBROUTINE time_domain_analysis(spice_bundle_model,spice_validation_test) USE type_specifications USE general_module USE constants USE cable_module USE cable_bundle_module USE spice_cable_bundle_module USE maths USE frequency_spec IMPLICIT NONE ! variables passed to subroutine TYPE(spice_model_specification_type),intent(IN):: spice_bundle_model ! Spice cable bundle model structure TYPE(spice_validation_test_type),intent(IN) :: spice_validation_test ! Spice validation circuit structure ! local variables ! variables for time domain analytic solution using IFT integer :: n_timesteps ! number of timesteps requested integer :: n_FFT ! nummber of samples in the FFT complex(dp),allocatable :: Vs_time(:) ! Source voltage as a function of time complex(dp),allocatable :: Vs_freq(:) ! Source voltage as a function of frequency complex(dp),allocatable :: VTL_time(:) ! Output voltage as a function of time complex(dp),allocatable :: VTL_freq(:) ! Output voltage as a function of frequency real(dp),allocatable :: time(:) ! time values real(dp),allocatable :: freq(:) ! frequency values ! variables for the frequency domain MTL solution real(dp) :: f,w ! frequency and angular frequency integer :: frequency_loop ! frequency loop variable integer :: dim ! dimension of matrix system to solve ! domain based impedance and admittance matrices complex(dp),allocatable :: Z_domain(:,:) complex(dp),allocatable :: Y_domain(:,:) ! domain based conductor impedance terms complex(dp),allocatable ::Z_domain_conductor_impedance_correction(:,:) ! Vectors and matrices used in the frequency domain solution of the transmission line equations with termination conditions complex(dp),allocatable :: Vs1(:) complex(dp),allocatable :: Z1(:,:) complex(dp),allocatable :: Vs2(:) complex(dp),allocatable :: Z2(:,:) complex(dp) :: Vout ! complex output voltage value ! domain transformation matrices complex(dp),allocatable :: MI(:,:) complex(dp),allocatable :: MII(:,:) complex(dp),allocatable :: MV(:,:) complex(dp),allocatable :: MVI(:,:) ! temporary working matrices complex(dp),allocatable :: TM1(:,:) integer :: conductor,inner_domain,outer_domain integer :: domain1,inner_domain1,outer_domain1 integer :: conductor1,reference_conductor1 integer :: domain_conductor1,domain_reference_conductor1 logical :: is_shield1 integer :: domain2,inner_domain2,outer_domain2 integer :: conductor2,reference_conductor2 integer :: domain_conductor2,domain_reference_conductor2 logical :: is_shield2 ! conductor based impedance (loss) and transfer impedance model data complex(dp) :: Zint_c ! conductor surface impedance complex(dp) :: Zint_c_ref ! reference conductor surface impedance real(dp) :: Rdc_c ! d.c. resistance of conductor real(dp) :: Rdc_c_ref ! d.c. resistance of reference conductor complex(dp) :: Zint_t ! conductor transfer impedance complex(dp) :: Zint_t_ref ! reference conductor transfer impedance real(dp) :: Rdc_t ! d.c. resistance of conductor (from transfer impedance) real(dp) :: Rdc_t_ref ! d.c. resistance of reference conductor (from transfer impedance) ! complex amplitude of incident field complex(dp) :: Einc logical,allocatable :: is_shielded_flag(:) ! flags conductors which are not exposed to the incident field integer :: shield_conductor ! temporary variable, shield conductor number for shielded conductors real(dp),allocatable :: local_conductor_x_offset(:) ! x coordinate in bundle cross section of conductors real(dp),allocatable :: local_conductor_y_offset(:) ! y coordinate in bundle cross section of conductors integer :: n_conductors_outer_domain ! for shield conductors, the number of conductors in the domain outside the shield integer :: shield_conductor_number_in_outer_domain ! for shield conductors, the conductor number in the domain outside the shield ! loop variables integer :: i integer :: row,col integer :: ierr ! START ! initialise FFT stuff ! set the number of frequencies to be equal to the next power of 2 above the number of timesteps requested ! Theory_Manual_Section 2.3.3 n_timesteps=NINT(spice_validation_test%runtime/spice_validation_test%timestep)+1 if(verbose) write(*,*)'Number of timesteps in transient analysis:',n_timesteps i=1 10 CONTINUE i=i*2 if (i.ge.n_timesteps) then n_FFT=i else GOTO 10 end if if(verbose) write(*,*)'Number of frequency in transient analysis:',n_FFT ! allocate the arrays used in the FFT ALLOCATE( time(n_FFT) ) ALLOCATE( Vs_time(n_FFT) ) ALLOCATE( VTL_time(n_FFT) ) ALLOCATE( freq(n_FFT) ) ALLOCATE( Vs_freq(n_FFT) ) ALLOCATE( VTL_freq(n_FFT) ) ! Generate the time domain voltage source function array ! The time domain waveform is the Exponential Pulse, the same as that used in Spice ! See Ngspice manual, section 4.1.3 do i=1,n_FFT time(i)=(i-1)*spice_validation_test%timestep if (time(i).LT.spice_validation_test%width) then Vs_time(i)=(1d0-exp(-time(i)/spice_validation_test%risetime)) else if (time(i).GE.spice_validation_test%width) then Vs_time(i)=(exp(-(time(i)-spice_validation_test%width)/spice_validation_test%risetime)) end if end do ! FFT the time domain voltage source array to give the frequency domain excitation function ! Theory_Manual_Equation 2.48 CALL FFT_TIME_TO_FREQ(n_FFT,time,Vs_time,freq,Vs_freq) ! allocate memory dim=spice_bundle_model%bundle%system_dimension ALLOCATE( Z_domain(dim,dim) ) ALLOCATE( Y_domain(dim,dim) ) ALLOCATE( Z_domain_conductor_impedance_correction(dim,dim) ) ALLOCATE( Vs1(dim) ) ALLOCATE( Z1(dim,dim) ) ALLOCATE( Vs2(dim) ) ALLOCATE( Z2(dim,dim) ) ! domain transformation matrices ALLOCATE( MI(dim,dim) ) ALLOCATE( MII(dim,dim) ) ALLOCATE( MV(dim,dim) ) ALLOCATE( MVI(dim,dim) ) ! temporary working matrices ALLOCATE( TM1(dim,dim) ) ALLOCATE( is_shielded_flag(1:dim+1) ) ALLOCATE( local_conductor_x_offset(1:dim+1) ) ALLOCATE( local_conductor_y_offset(1:dim+1) ) ! loop over conductors to work out which are in shielded domains and which are in the external domain ! also get the position of the conductor in the bundle cross section for incident field excitation do i=1,dim+1 if (spice_bundle_model%bundle%terminal_conductor_to_outer_domain(i).EQ.spice_bundle_model%bundle%tot_n_domains) then is_shielded_flag(i)=.FALSE. local_conductor_x_offset(i)=spice_bundle_model%bundle%conductor_x_offset(i) local_conductor_y_offset(i)=spice_bundle_model%bundle%conductor_y_offset(i) else is_shielded_flag(i)=.TRUE. ! work out the conductor number of the shield shield_conductor=spice_bundle_model%bundle%terminal_conductor_to_reference_terminal_conductor(i) ! shielded conductors pick up the coordinate of the shield for the purposes on incident field excitation local_conductor_x_offset(i)=spice_bundle_model%bundle%conductor_x_offset(shield_conductor) local_conductor_y_offset(i)=spice_bundle_model%bundle%conductor_y_offset(shield_conductor) end if end do ! build the termination specifications and convert to complex Vs1(1:dim)=cmplx( spice_validation_test%Vs_end1(1:dim)-spice_validation_test%Vs_end1(dim+1) ) Vs2(1:dim)=cmplx( spice_validation_test%Vs_end2(1:dim)-spice_validation_test%Vs_end2(dim+1) ) Z1(:,:) =cmplx( spice_validation_test%R_end1(dim+1) ) Z2(:,:) =cmplx( spice_validation_test%R_end2(dim+1) ) do i=1,dim Z1(i,i) =Z1(i,i)+cmplx( spice_validation_test%R_end1(i) ) Z2(i,i) =Z2(i,i)+cmplx( spice_validation_test%R_end2(i) ) end do ! Copy the domain transformation matrices and calculate the inverses MI(:,:)=spice_bundle_model%bundle%global_MI%mat(:,:) MV(:,:)=spice_bundle_model%bundle%global_MV%mat(:,:) ierr=0 ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix CALL cinvert_Gauss_Jordan(MI,dim,MII,dim,ierr) ierr=0 ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix CALL cinvert_Gauss_Jordan(MV,dim,MVI,dim,ierr) ! loop over frequency to calculate the frequency domain transmission line transfer function write(6,8100,advance='no')'Frequency ',0,' of ',n_FFT/2+1 flush(6) do frequency_loop=1,n_FFT/2+1 write(6,'(A)',advance='no')char(13) write(6,8100,advance='no')'Frequency ',frequency_loop,' of ',n_FFT/2+1 flush(6) 8100 format(A10,I10,A4,I10) f=freq(frequency_loop) ! shift frequency from d.c. as solution fails for lossless systems at f=0 Hz if (f.eq.0d0) then f=1d0 end if w=2d0*pi*f ! Use the global domain based L and C matrices and the domain voltage and current ! domain transformation matrices to calculate the impedance [Z] and admittance [Y] matrices do row=1,dim do col=1,dim ! Evaluate the cable impedance filter function Z_domain(row,col)=evaluate_Sfilter_frequency_response(spice_bundle_model%bundle%global_Z%sfilter_mat(row,col),f) ! Evaluate the cable admittance filter function Y_domain(row,col)=evaluate_Sfilter_frequency_response(spice_bundle_model%bundle%global_Y%sfilter_mat(row,col),f) end do ! next col end do ! next row ! calculate the contribution to the matrices from the conductor based impedance models. ! See Theory_Manual_Section 2.2.3 !Initailly set to zero Z_domain_conductor_impedance_correction(1:dim,1:dim)=(0d0,0d0) ! new domain based model do conductor1=1,dim domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor1) reference_conductor1=spice_bundle_model%bundle%terminal_conductor_to_reference_terminal_conductor(conductor1) domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(conductor1) domain_reference_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(reference_conductor1) is_shield1=spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor1) ! evaluate the surface impedance for this conductor conductor CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(conductor1), & f,Zint_c,Rdc_c,Zint_t,Rdc_t) ! Apply multiplication factor to the conductor impedance to correct for common mode/ differential modes in twisted pairs ! See note at the end of Theory_Manual_Section 3.5.4 Zint_c=Zint_c*spice_bundle_model%bundle%conductor_impedance(conductor1)%Resistance_multiplication_factor ! evaluate the surface impedance for the reference conductor CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(reference_conductor1), & f,Zint_c_ref,Rdc_c_ref,Zint_t_ref,Rdc_t_ref) ! Apply multiplication factor to the conductor impedance to correct for common mode/ differential modes in twisted pairs Zint_c_ref=Zint_c_ref*spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%Resistance_multiplication_factor ! The surface impedance of the conductor and the reference conductor contribute to the diagonal element Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor1)=Zint_c+Zint_c_ref if (verbose) then write(*,*)'conductor:',conductor1,' reference conductor',reference_conductor1 write(*,*)'conductor loss model type',spice_bundle_model%bundle%conductor_impedance(conductor1)%impedance_model_type write(*,*)'refconductor loss model type', & spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%impedance_model_type write(*,*)'radius', & spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%radius write(*,*)'conductivity', & spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%conductivity write(*,*)'Resistance_multiplication_factor', & spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%Resistance_multiplication_factor write(*,*)'thickness', & spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%thickness write(*,*)'domain conductor:',domain_conductor1,' domain reference conductor',domain_reference_conductor1 write(*,*)'Contribution to Zc(',domain_conductor1,domain_conductor1,')' write(*,*)'Zc conductor:',Zint_c write(*,*)'Zc reference:',Zint_c_ref end if ! verbose ! conductor always gets the contribution from its own surface impedance do conductor2=conductor1+1,dim domain2=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor2) reference_conductor2=spice_bundle_model%bundle%terminal_conductor_to_reference_terminal_conductor(conductor2) domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(conductor2) domain_reference_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(reference_conductor2) is_shield2=spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor2) ! if the two conductors belong to the same domain then add the reference conductor impedance if (domain1.EQ.domain2) then Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)=Zint_c_ref Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)=Zint_c_ref if (verbose) then write(*,*)'Contribution to Zc(',domain_conductor1,domain_conductor2,')' write(*,*)'Contribution to Zc(',domain_conductor2,domain_conductor1,')' write(*,*)'Zc conductor:',Zint_c_ref end if ! verbose end if end do ! next conductor2 end do ! next conductor1 if (verbose) write(*,*)'Add transfer impedance contributions' ! add transfer impedance contributions ! loop over conductors looking for shields. Note include all conductors including the reference here do conductor=1,dim+1 is_shield1=spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor) if (is_shield1) then ! add transfer impedance contributions to inner and outer domain conductors inner_domain=spice_bundle_model%bundle%terminal_conductor_to_inner_domain(conductor) outer_domain=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor) CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(conductor), & f,Zint_c,Rdc_c,Zint_t,Rdc_t) ! Check whether the shield is the reference conductor in the outer domain - the contributions ! are different if this is the case. n_conductors_outer_domain=spice_bundle_model%bundle%n_conductors(outer_domain) shield_conductor_number_in_outer_domain=spice_bundle_model%bundle%terminal_conductor_to_local_domain_conductor(conductor) ! number of conductors in a domain is spice_bundle_model%bundle%n_conductors(domain) if (shield_conductor_number_in_outer_domain.NE.n_conductors_outer_domain) then ! loop over all conductors do row=1,dim ! get the domain of row conductor domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(row) if (domain1.EQ.inner_domain) then ! The row conductor is in the inner shield domain and so gets a transfer impedance contribution from the shield conductor ! the shield couples these two domains so add the transfer impedance term - also include term to make the matrix symmetric domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(row) domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(conductor) Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)= & Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2) -Zint_t Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)= & Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1) -Zint_t if (verbose) then write(*,*)'Shield conductor',conductor,' inner domain',inner_domain,' outer domain',outer_domain write(*,*)'row',row,' col',col,' row domain',domain1,' col domain',domain2 write(*,*)'Contribution to Zt(',domain_conductor1,domain_conductor2,')' write(*,*)'Contribution to Zt(',domain_conductor2,domain_conductor1,')' write(*,*)'Zt conductor:',-Zint_t end if ! verbose end if ! transfer impedance term required end do ! next row conductor else ! shield IS reference conductor in outer domain ! loop over all conductors do row=1,dim ! get the domain of row conductor domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(row) if (domain1.EQ.inner_domain) then ! The row conductor is in the inner shield domain and so gets a transfer impedance contribution from the shield conductor ! the shield couples these two domains so add the transfer impedance term - also include term to make the matrix symmetric domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(row) ! As the shield conductor is the reference we need to find all the conductors contributing to the shield current ! note that the contribution is then -ve of the normal transfer impedance contribution as the currents are in the ! opposite direction do col=1,dim domain2=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(col) ! Check the domain of the col conductor. If it is an outer domain conductor of the shield then it contributes if (domain2.EQ.outer_domain) then domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(col) Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)= & Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2) +Zint_t Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)= & Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1) +Zint_t if (verbose) then write(*,*)'Shield conductor',conductor,' inner domain',inner_domain,' outer domain',outer_domain write(*,*)'row',row,' col',col,' row domain',domain1,' col domain',domain2 write(*,*)'Contribution to Zt(',domain_conductor1,domain_conductor2,')' write(*,*)'Contribution to Zt(',domain_conductor2,domain_conductor1,')' write(*,*)'Zt conductor:',-Zint_t end if ! verbose end if ! transfer impedance term required for this col conductor end do ! next condutor to check end if ! transfer impedance term required for this row conductor end do ! next row conductor end if ! shield is/ is not reference conductor in outer domain end if ! conductor is a shield end do ! next conductor ! Add the conductor impedance contributions to the domain based impedance matrix Z_domain(:,:)=Z_domain(:,:)+Z_domain_conductor_impedance_correction(:,:) Einc=cmplx(spice_bundle_model%Eamplitude) ! Solve the frequency domain multi-conductor transmission line equations with the specified termination circuit and ! return the required conductor voltage in Vout. ! For the purposes of the frequency domain solution, the excitations are assumed to be constant in frequency ! The frequency response of the excitation is included afterwards if (.NOT.run_validation_test_Vbased) then CALL frequency_domain_MTL_solution(dim,Z_domain,Y_domain,MV,MVI,MI,MII, & Einc,spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez, & spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz, & spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz, & local_conductor_x_offset, & local_conductor_y_offset, & spice_bundle_model%bundle%ground_plane_present, & spice_bundle_model%bundle%ground_plane_x, & spice_bundle_model%bundle%ground_plane_y, & spice_bundle_model%bundle%ground_plane_theta, & spice_bundle_model%length,Vs1,Z1,Vs2,Z2, & is_shielded_flag, & f,spice_validation_test%output_end,spice_validation_test%output_conductor, & spice_validation_test%output_conductor_ref,Vout) else CALL frequency_domain_MTL_solution_V(dim,Z_domain,Y_domain,MV,MVI,MI,MII, & Einc,spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez, & spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz, & spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz, & local_conductor_x_offset, & local_conductor_y_offset, & spice_bundle_model%bundle%ground_plane_present, & spice_bundle_model%bundle%ground_plane_x, & spice_bundle_model%bundle%ground_plane_y, & spice_bundle_model%bundle%ground_plane_theta, & spice_bundle_model%length,Vs1,Z1,Vs2,Z2, & is_shielded_flag, & f,spice_validation_test%output_end,spice_validation_test%output_conductor, & spice_validation_test%output_conductor_ref,Vout) end if ! multiply the Frequency domain source and transmission line transfer function ! Save transfer function multiplied by frequency domain source voltage function i.e. include the time response of ! the excitation here. Theory_Manual_Section 2.3.3 VTL_freq(frequency_loop)=Vout*Vs_freq(frequency_loop) if ( (frequency_loop.NE.1).AND.(frequency_loop.NE.N_FFT/2+1) ) then VTL_freq(N_FFT-frequency_loop+2)=conjg(VTL_freq(frequency_loop)) end if end do ! next frequency in frequency loop ! Inverse FFT to give the time domain transmission line response ! Theory_Manual_Equation 2.50 CALL FFT_FREQ_TO_TIME(n_FFT,time,VTL_time,freq,VTL_freq) ! write the time domain response to file ! Open output file open(unit=analytic_soln_file_unit,file=trim(analytic_soln_filename)) ! write the file header line write(analytic_soln_file_unit,'(A)')time_header do i=1,n_timesteps ! changed from n_FFT 04/5/2016 CJS write(analytic_soln_file_unit,*)time(i),real(VTL_time(i)) end do ! Close output file Close(unit=analytic_soln_file_unit) ! deallocate memory DEALLOCATE( Z_domain ) DEALLOCATE( Y_domain ) DEALLOCATE( Z_domain_conductor_impedance_correction ) DEALLOCATE( Vs_time ) DEALLOCATE( Vs_freq ) DEALLOCATE( VTL_time ) DEALLOCATE( VTL_freq ) DEALLOCATE( time ) DEALLOCATE( freq ) ! domain transformation matrices DEALLOCATE( MI ) DEALLOCATE( MII ) DEALLOCATE( MV ) DEALLOCATE( MVI ) DEALLOCATE( is_shielded_flag ) DEALLOCATE( local_conductor_x_offset ) DEALLOCATE( local_conductor_y_offset ) ! temporary working matrices DEALLOCATE( TM1 ) write(6,*) RETURN END SUBROUTINE time_domain_analysis