JonesTransform.m

function [rhoNB,h,theta] = JonesTransform(Z,Sig_s,Sig_e,periods,orient,xypairs,h,name)

% [rhoNB,h,theta] = JonesTransform(Z,Sig_s,Sig_e,periods,orient,xypairs,h,name)
%
% Alan Jones's 1D approximate transformation is an ad-hoc version of the
% (ad-hoc) traditional Niblett-Bostick that employs the phase:
%   rhoNB(h) = rhoa(T) * ((pi/(2*phia(T)) - 1)
%   h = sqrt( (rhoa(T) * T) / (2*pi*mu0) )
% Quoting:
% "I find that the most robust is the RhoMAX at a given depth, which I obtain
% by doing the N-B transportation at 0 deg and finding the RhoMAX at the
% depth I want, then rotating Z by 1 deg increments and doing the same. This
% gives me 91 estimates of RhoMAX(theta), and I take the maximum of them. I
% also have a lot of joy with the Anisotropy."

warning off;

T = periods;
nT = length(T);
mu0 = 4*pi*10^(-7);
rxy = zeros(91,nT);
ryx = zeros(91,nT);
pxy = zeros(91,nT);
pyx = zeros(91,nT);
hxy = zeros(91,nT);
hyx = zeros(91,nT);

rhoNBxy = zeros(91,nT);
rhoNByx = zeros(91,nT);

% use the given depths to interpolate (in meters)
if nargin <= 6
    h = 1000 * (10:10:200);
end

% if the site name if given, plot all profiles and print to file
if nargin > 7
    PLOT = 1;
else
    PLOT = 0;
end

for theta = 0:1:90
    
    i = theta+1;
    
    %  xform  takes TFs input from Z file, converts to fixed coordinate system
    %    (x points in direction 0 used as reference for channel orientations)
    %    rotates, and computes rh, phi + standard errors
    %   theta defines rotation angle, ixy defines pairs of channels to rotate
    [ryx(i,:),rxy(i,:),pyx(i,:),pxy(i,:)] = ...
        xform(Z,Sig_s,Sig_e,periods,orient,xypairs,theta);
    
    % compute two Bostick resistivity curves for XY and YX modes
    hxy(i,:) = sqrt(rxy(i,:).*T/(2*pi*mu0));
    rhoNBxy(i,:) = rxy(i,:) .* (180./(2*pxy(i,:))-1);
    
    hyx(i,:) = sqrt(ryx(i,:).*T/(2*pi*mu0));
    rhoNByx(i,:) = ryx(i,:) .* (180./(2*pyx(i,:))-1);
    
end

rhoNBxy(rhoNBxy<=0) = NaN;
rhoNByx(rhoNByx<=0) = NaN;

rhoNBmax = zeros(91,length(h));
for i = 1:91
    temp1 = interp1(hxy(i,:),rhoNBxy(i,:),h);
    temp2 = interp1(hyx(i,:),rhoNByx(i,:),h);
    rhoNBmax(i,:) = max(temp1,temp2);
end
[rhoNB,theta] = max(rhoNBmax,[],1);
theta = theta - 1;

if PLOT
    figure(101);
    for i=1:91
        semilogx(rhoNBxy(i,:),-hxy(i,:)/1000,'b-','linewidth',1.5); hold on;
        semilogx(rhoNByx(i,:),-hyx(i,:)/1000,'r-','linewidth',1.5); hold on;
    end
    ylim([-200 0]);
    ylabel('Depth [km]','FontSize',12,'FontWeight','demi');
    xlim([10^(-2) 10^4]);
    xlabel('Niblett-Bostick \rho_{NB} [\Omega m]',...
        'FontSize',12,'FontWeight','demi');
    title([name ' (XY=B; YX=R)'],'interpreter','none',...
        'FontSize',12,'FontWeight','demi');
    print('-r200','-dpng',[name '_JonesPlot']);
end