clear all
%   ECE539 - wrapper for Ansoft (time domain)   
%   HW problem to lecture #5 - 4.2mm spacing

%---------------------------------------------------------------
%   Input parameters - ramp
V0      = 1.0;          % ramp height in V
tau     = 62.5e-12;     % ramp duration in sec
t_min   = -10*tau;      % minimum observable time
t_max   = +40*tau;      % maximum observable time

%---------------------------------------------------------------
%   Import Ansoft data
string1  = strcat('S41.tab');
fid      = fopen(string1, 'r'); DUMMY  = fscanf(fid, '%g', [1 inf]); fclose(fid);
%   Extract data (tested)
f       = DUMMY(1:3:end)*1e9; delta_f = f(2)-f(1)
S41re   = DUMMY(2:3:end)*1.0;
S41im   = DUMMY(3:3:end)*1.0;
%   Continue to DC
dummy   =   f(1):-delta_f:0; dummy = dummy(end:-1:1);
f       = [dummy f];
S41re   = [zeros(size(dummy)) S41re];
S41im   = [zeros(size(dummy)) S41im];

%---------------------------------------------------------------
%   Plot the data
f1      = figure; hold on; grid on;
plot(f, S41re, '-k', 'LineWidth', 2)
plot(f, S41im, '--k','LineWidth', 2)
xlabel('frequency, Hz');
ylabel('S41, \Omega');
title('real(S41)-solid line; imag(S41)-dashed line')

%---------------------------------------------------------------
%   Determine parameters for IDFT
T       = max(abs(t_min),abs(t_max));
N       = round(1/(delta_f*T));
if N<6  %   Less than 6 points per period 
    warning('Accuracy of IDFT is lost! Please decrease the frequency step in Ansoft')
end
t               =   [t_min:0.01*tau:t_max];    % times
omega           =   2*pi*f;
delta_omega     =   2*pi*delta_f; 

%---------------------------------------------------------------
%   Ramp response in frequency domain (no delta-function)
F       = f + eps;              % to eliminate singularity
OMEGA   = 2*pi*F;
V       = (1./(OMEGA.^2*tau)).*(exp(-j*OMEGA*tau)-1);
V       = V0*V;

%---------------------------------------------------------------
%   Test the initial ramp (how good is the Ansoft frequency data)
H       = ones(size(f));        % transfer function of 1 Ohm resistance
V       = H.*V;                 % units are V*sec
%   IDFT
for m =1:length(t)    
    %   positive frequencies
    temp_plus        =   1/(2*pi)*delta_omega*sum(V.*exp(+j*omega*t(m)));  
    %   negative frequencies
    temp_minus       =   1/(2*pi)*delta_omega*sum(conj(V).*exp(-j*omega*t(m)));  
    % zero frequencies - delta function (only if H(1)~=0 - resistance)
    temp_zero        =   +1/(2*pi)*V0*pi*H(1);
    %   integral
    v(m)             = temp_plus + temp_minus + temp_zero;
end
f2      = figure; hold on; grid on;
plot(t*1e9, v, '-k', 'LineWidth', 2)
xlabel('time, ns');
ylabel('v(t), Volts');
title('Test of the ramp response with the availble frequency data')

%---------------------------------------------------------------
%   Coupled voltage on the second trace
H       = S41re + j*S41im;
V       = H.*V;                 % units are V*sec
%   IDFT
for m =1:length(t)    
    %   positive frequencies
    temp_plus        =   1/(2*pi)*delta_omega*sum(V.*exp(+j*omega*t(m)));  
    %   negative frequencies
    temp_minus       =   1/(2*pi)*delta_omega*sum(conj(V).*exp(-j*omega*t(m)));  
    % zero frequencies - delta function (only if H(1)~=0 - resistance)
    temp_zero        =   +1/(2*pi)*V0*pi*H(1);
    %   integral
    v(m)             = temp_plus + temp_minus + temp_zero;
end
f3      = figure; hold on; grid on;
plot(t*1e9, v, '-k', 'LineWidth', 2)
xlabel('time, ns');
ylabel('v(t), Volts');
title('Coupled voltage at port 4 - the second trace')
return;