
function fluxgate
clear all;
close all;
% Illustrates operation of a fluxgate magnetometer

% Sinusoidal H field generated with frequency f Hz
% External field added (H_ext); Can be shifted up and down in steps on inc
h_ext=0.0;

% Option 1  ince =0.001    output +/-  20 mv shows spike that is prop to H_ext
% Option 2  ince =0.01     output +/-  50 mv shows earlier and later saturation

ichoice = MENU('Options','Large H','Small H','Exit')
if ichoice ==1
   inc=0.01; hlim = 50;
end

if ichoice ==2
   inc=0.001;
end

% B computed with assumption that saturation occurs at B=sat and B=-sat
mu =1;       sat=0.75;
% V is voltage induced in secondary;


irun=1;


t=[1:1:1000];
dt=0.001;
n=length(t);
f=5;

% Start main loop
while irun == 1

  for i=1:n
    h1(i)=h_ext+sin(f*i*2*pi/n);    h2(i)=h_ext-sin(f*i*2*pi/n);   
  end
  b1 = saturation(h1,mu,sat);       b2 = saturation(h2,mu,sat);

  v1 =-time_der(b1,dt);             v2 = -time_der(b2,dt);

  % Plot primary H-fields
  figure(1);
  subplot(4,1,1);  plot(t,h1,'k');  hold on;   plot(t,h2,'r');  hold off;  axis([0,max(t),-1.5,1.5])
  title ('Operation of a fluxgate magnetometer')
  ylabel('H_P (A/m)');
  
  subplot(4,1,2);  plot(t,b1,'k');  hold on;   plot(t,b2,'r');  hold off;  axis([0,max(t),-1.2,1.2])
  ylabel('B(nT)');
  subplot(4,1,3);  plot(t,v1,'k');  hold on;   plot(t,v2,'r');  hold off;  axis([0,max(t),-50,50])
  ylabel('V(mv)');
  
  subplot(4,1,4);  plot(t,v1+v2);              axis([0,max(t),-hlim,hlim])
  ylabel('Output V(mv)');
  xlabel('time(ms)');
   

  ichoice = MENU('Options','Increase H_Ext','Decrease H_ext','Exit')

  if ichoice ==1; h_ext=h_ext+inc; end
  if ichoice ==2; h_ext=h_ext-inc; end
  if ichoice ==3; irun = -1 ; end
  h_ext
end

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

return

%====================================
function b=saturation (h,mu,sat)
n=length(h);
for i=1:n
  b(i)=mu*h(i);
  if b(i) >= sat; b(i)=sat; end
  if b(i) <= -sat; b(i)=-sat; end
end
return
%====================================


%====================================
function v=time_der(b,dt)
n=length(b);
for i=1:n-1
   v(i)=(b(i+1)-b(i))/dt;
end
   v(n)=v(n-1);
return
%====================================

