
function MT1Dfit

clear all
close all

rho_floor = 0.2;
pha_floor = 5.;
fmin = 0.00003;
fmax = 300;

%load data  (In NLCG6 format)
load MTdata.txt
T_data = MTdata(:,1);
ndata=length (T_data);
for i=1:ndata
   f_data(i)=1./T_data(i);
end

f=f_data; nf = length(f);

rho_data =  MTdata(:,2);  rho_err  =  MTdata(:,4);
pha_data = -MTdata(:,3);  pha_err  =  180.*MTdata(:,5)/pi;


% APPLY ERROR FLOORS
i_error=menu('Data errors:', 'Use defined errors','Apply error floor','Set to constant')
    if i_error ==2  
      for i=1:nf    
        if rho_err(i) <= rho_floor
           rho_err(i) = rho_floor;
        end
    
        if pha_err(i) <= pha_floor
          pha_err(i) = pha_floor;
        end
      end
    end  
    
    if i_error ==3  
      for i=1:nf    
        rho_err(i) = rho_floor;
        pha_err(i) = pha_floor;
      end
    end  
    

% Start loop requesting model change
irun =1
while irun ~= -1

%load model
  load MTmodel.txt
  z = MTmodel(:,1);  r = MTmodel(:,2); nl=length(r);
  for i=1:nl
    s(i)= 1./r(i);   
  end
  
  
 % Define frequency set with logarithmic spacing
 % ifreq=0;
 % for logfreq = -5:0.1:3 
 %   ifreq=ifreq+1;
 %   f(ifreq)=10^logfreq;
 %end
 


% Compute MT responses
%==========================================================================
muo=pi*4.0e-7;
ci = sqrt(-1.);
degree= 180.0/pi;
sqrt2 = sqrt(2.);
c3 =  pi*muo;

nf=length(f);

% Start frequency loop
for i=1:nf
  e(nl) = 1.;
  if s(nl) >= 0.
     epr(nl) = -(ci+1)*sqrt(c3 * f(i) *s(nl) );
     d   = 1./epr(nl);
  else
     epr(nl) = 0.;
  end
  jp1   = nl;
     
% Start layer loop  
  for j= nl-1:-1:1
     dz = z(j) - z(jp1);
     if s(j)>= 0.
       epr(j) = -(ci+1)*sqrt(c3 * f(i) *s(j) );
       d   = 1./epr(j);        
       if abs(dz/d)>= 30
% field completely attennuated in this layer
         e(j) = 1.;
         epr(j) = 1./d;
         for k=jp1:nl
           e(k) = 0.;
           epr(k) = 0.;
         end
       else
         ex= exp( dz/d );
         eprd = epr(jp1) * d;
         a = ( e(jp1) + eprd ) * ex / 2.;
         b = ( e(jp1) - eprd ) / (2. * ex );
         e(j)=a+b;
         epr(j) = (a-b)/d;
       end
     else
  % Zero conductivity option   
       e(j)   = dz *epr(jp1) + e(jp1);
       epr(j) = epr(jp1);
     end
        
     et = e(j);
     for k=j:nl
        e(k)   = e(k)/et;
        epr(k) = epr(k)/et;
     end
     jp1  = j;
% End layer loop
end

c =  -e(1) / epr(1);
rho(i)= abs(c*c) * muo * 2.0*pi* f(i);
phase(i)= 90.0 + atan(imag(c)/real(c))*degree;
   
% end of frequency loop 
end         
%===========================================================================
  
% Compute misfit
misfit_res=0.0; misfit_pha=0.0;
for i=1:nf
   rho_res(i)=(rho(i)-rho_data(i))/(rho(i)*rho_err(i));
   misfit_res = misfit_res+rho_res(i)^2;
   pha_res(i)=(phase(i)-pha_data(i))/(pha_err(i));
   misfit_pha = misfit_pha+pha_res(i)^2;
end

misfit = misfit_res+misfit_pha;
misfit_res = sqrt(misfit_res/nf);
misfit_pha = sqrt(misfit_pha/nf) ;
misfit = sqrt(misfit/(2*nf));
 
for i=1:nf
    rho_upper(i) = rho_data(i)*(1.+rho_err(i));
    rho_lower(i) = rho_data(i)*(1.-rho_err(i));
    pha_upper(i) = pha_data(i)+pha_err(i);
    pha_lower(i) = pha_data(i)-pha_err(i);
end

%===========================================
% PLOT FIT TO RHO DATA
%===========================================
subplot (3,2,1)

  loglog(f_data,rho_data,'b.')
  hold on
  loglog (f,rho,'r-')
  loglog (f,rho_upper,'b:');    loglog (f,rho_lower,'b:')
  set (gca,'YTICK',[0.01 0.1 1 10 100 1000 10000])
  set (gca,'XTICK',[0.0001 0.01  1 100])
  axis([fmin fmax 0.3 1000])
  title(['r.m.s. = ',num2str(misfit_res), '  Floor = ' num2str(100*rho_floor),' %']);
  %xlabel ('frequency')
  ylabel ('rho_a (ohm-m)')
 
%===========================================  
% PLOT FIT TO PHASE DATA
%===========================================
  subplot (3,2,3)
  semilogx(f,phase,'r-')
  hold on
  semilogx(f_data,pha_data,'b+')
  semilogx(f,pha_upper,'b:'); semilogx(f,pha_lower,'b:')
  set (gca,'YTICK',[0 45 90])
  set (gca,'XTICK',[0.0001 0.01 1 100])
  axis([fmin fmax 0 120])
  title(['r.m.s. = ',num2str(misfit_pha), '  Floor = ' num2str(pha_floor),' deg']);
  
  %xlabel ('frequency')
  ylabel ('phase (degrees)')
  
%=========================================== 
% RESIDUALS
%===========================================
  subplot (3,2,5)
  semilogx(f,rho_res,'b.')
  hold on
  semilogx(f_data,pha_res,'b+')

  set (gca,'XTICK',[0.0001 0.01  1 100])
  axis([fmin fmax -10 10])

    
  lp = [-4,-2,0,2,4]
  for i=1:length(lp)
      xp = [fmin,fmax]; yp = [lp(i),lp(i)];
      plot(xp,yp,':')
  end

  
  xlabel ('frequency')
  ylabel ('residuals')
  if i_error == 1; title([' Residuals : Original errors']); end
  if i_error == 2; title([' Residuals : Error floor used']); end
  if i_error == 3; title([' Residuals : Errors constant']); end
  
%===========================================
% PLOT RESISTIVITY MODEL
%===========================================
  for ilayer =1:nl-1
     rhoplot(2*ilayer-1)=r(ilayer);
     rhoplot(2*ilayer)  =r(ilayer);
     dplot(2*ilayer-1)  =z(ilayer);
     dplot(2*ilayer)    =z(ilayer+1);
  end
  
  rhoplot(2*nl-1) = r(nl); rhoplot(2*nl)   = r(nl);
  dplot (2*nl-1)  = z(nl); dplot (2*nl)    = 20000.;
  
  subplot(1,2,2)
    semilogx(rhoplot,-dplot*0.001,'r')
    hold on
    semilogx(rhoplot,-dplot*0.001,'r')
    set (gca,'XTICK',[0.01 0.1 1 10 100 1000 10000])
    axis([0.01 10000 -20 0])
    xlabel('resistivity (ohm.m)')
    ylabel('depth(km)')
    title(['Total r.m.s. misfit ',num2str(misfit)]);
    

        
  choice=menu('Options', 'Change model','Print result and quit')
    if choice ==1  
      irun=1
    end

    if choice ==2
      irun = -1
      print -dpsc  MT1Dfit.ps
    end   
  end   
end
