
function input_EM

%
clear all;
close all;

ichoice = menu('Conductor','Very good (0.1 ohm)','Good (1 ohm)','Bad (3 ohm)','Very bad (10 ohm)')

if ichoice == 1; r=0.3; 	end
if ichoice == 2; r=1.0; 	end
if ichoice == 3; r=3.0; 	end
if ichoice == 4; r=10.0; 	end

l=0.02;
k_tr = 10;
k_tc = 1;
k_cr = 1;
dt =0.001;

t=[100,200,300,500,600,700,800]
npts = max(t);
ichoice = menu('TX waveform','Sawtooth','Pulsed sine')

if ichoice ==1 
% Define primary current (i_prim)
	for it=   1:t(1);     	i_prim(it)= 0;  								end
	for it= t(1)+1:t(2);   	i_prim(it)= (it-t(1))/(t(2)-t(1));  	end
	for it= t(2)+1:t(3);   	i_prim(it)= (t(3)-it)/(t(2)-t(1));  	end
	for it= t(3)+1:t(4);  	i_prim(it)=0;  								end
	for it= t(4)+1:t(5);  	i_prim(it)= -(it-t(4))/(t(2)-t(1));  	end
	for it= t(5)+1: t(6);  	i_prim(it)= -(t(6)-it)/(t(2)-t(1));  	end
	for it= t(6)+1: t(7);  	i_prim(it)=0; 									end
end

if ichoice ==2 
% Define primary current (i_prim)
	for it=   1:t(1);     	i_prim(it)= 0;  								end
	for it= t(1)+1:t(3);   	i_prim(it)= sin(pi*(it-t(1))/(t(3)-t(1)));  	end
	for it= t(3)+1:t(4);  	i_prim(it)=0;  								end
   for it= t(4)+1:t(6);   	i_prim(it)= -sin(pi*(it-t(4))/(t(6)-t(4)));  	end
	for it= t(6)+1: t(7);  	i_prim(it)=0; 									end
end



    
for it=1:npts
   i_sec(it)=0;
   time(it)=dt*it;
end
tmax=max(time);


% Compute primary magnetic field at conductor
h_prim_c = i_prim*k_tc;

% Compute primary magnetic field at receiver
h_prim_r = i_prim*k_tr;

% Compute secondary voltage in conductor by differentiating h_prim_r
for it=1:npts-1;   v_sec_c(it) = (h_prim_c(it+1)-h_prim_c(it))/dt;     end
v_sec_c(npts) = v_sec_c(npts-1);


% Integrate secondary voltage in conductor (v_sec_c) to get secondary current (IS)
for it = 1:npts
   sum=0;
   for ii=it:-1:1
      sum=sum+dt*(v_sec_c(ii)-r*i_sec(ii))/l;
   end
   i_sec(it)=sum;
end

% Compute secondary magnetic field at Receiver
h_sec_r = i_sec*k_cr;

% Compute total magnetic field at receiver
h_total_r = h_sec_r + h_prim_r;

% Compute secondary voltage in conductor by differentiating h_prim_r
for it=1:npts-1;   v_total_r(it) = (h_total_r(it+1)-h_total_r(it))/dt;     end
v_total_r(npts) = v_total_r(npts-1);


%===============================
% Plot results
%===============================

zero = [0.,0.]
tlimit = [min(time),max(time)];

figure(1);
subplot(5,1,1)
plot(time,i_prim,'k')
axis([0,tmax,1.3*min(i_prim),1.3*max(i_prim)])
title('Primary current in TX')
%xlabel('time(s)')
ylabel('Amps')
pause

subplot(5,1,2)
plot(time,v_sec_c,'r')
hold on
plot(tlimit,zero,'k')
hold off
axis([0,tmax,1.3*min(v_sec_c),1.3*max(v_sec_c)])
title('Secondary voltage in conductor')
%xlabel('time(s)')
ylabel('mV')
pause

subplot(5,1,3)
plot(time,i_sec,'r')
hold on
plot(tlimit,zero,'k')
hold off
axis([0,tmax,1.3*min(i_sec),1.3*max(i_sec)])
title('Secondary current in conductor')
%xlabel('time(s)')
ylabel('mA')
pause

subplot(5,1,4)
plot(time,h_prim_r)
hold on
plot(tlimit,zero,'k')
plot(time,h_sec_r,'r')
plot(time,h_prim_r,'k')
plot(time,h_total_r)
hold off
axis([0,tmax,1.3*min(h_total_r),1.3*max(h_total_r)])
title('Magnetic fields at RX')
%xlabel('time(s)')
ylabel('nT')
pause

subplot(5,1,5)
plot(time,v_total_r)
hold on
plot(tlimit,zero,'k')
hold off
axis([0,tmax,1.3*min(v_total_r),1.3*max(v_total_r)])
title('Total Voltage induced in RX')
xlabel('time(s)')
ylabel('mV')
pause



print (1,'-dpsc', 'input_EM.ps')


