% radioloc1 plots SNR for uplink and downlink radiolocation beacon

% based on PHDmagloop3

%clear
filetype = '-dtiff';

mu = 4E-7 * pi;
KT4 = 1.38E-23 * 290 * 4;
N = 50;
fig = 1;

%figure(1); clf;

D = linspace(10, 500, N);

%from spreadsheet noise-5.xls...
freq = 4; %kHz
FaMaxdB = 172 + 23*LOG10(10/freq);  % ref point is at 10kHz
FaMindB = 148 + 31*LOG10(10/freq);
Fa_dB = [FaMaxdB; FaMindB];
Fk_dB = Fa_dB + 40*LOG10(freq*1000)-294.147 + 3;  %3 is fiddle - see thesis

Fk_dB = [34; 13]

sigma = 10;  % mS/m
skinDepth = sqrt(1/(pi*freq*mu*sigma))

TxSA = 30;
TxPower = 0.5;
Txmoment = TxSA * sqrt(TxPower);

RxSA = 100;
NEB	= 0.35;

T = D/skinDepth;
T_attn = exp(-T) .* sqrt(T.^2 + (1+T).^2) ;
T_attn_dB = 20*log10(T_attn);
T_attn_dB = [T_attn_dB; T_attn_dB];

Skin_attn_dB = 20*log10( exp(1) ) * T;

% make skin_attn have two rows
Skin_attn_dB = [Skin_attn_dB; Skin_attn_dB];

% Make Fk have N columns
Fk_dB = Fk_dB * ones(1,N);

DownlinkFk_dB = Fk_dB - Skin_attn_dB;
DownlinkFk = 10.^(DownlinkFk_dB./20);

SurfaceFk = 10.^(Fk_dB./20);

DownlinkNF = sqrt(1 + (DownlinkFk.*RxSA).^2);
DownlinkNF_dB = 20*log10(DownlinkNF);
SurfaceNF = sqrt(1 + (SurfaceFk.*RxSA).^2);
SurfaceNF_dB = 20*log10(SurfaceNF);

% H near field at receiver

Hfield = Txmoment ./ (2*pi * D.^3) ;

% now work out overall SNR
SNR = 2*pi * freq*1000 * mu * Hfield * RxSA / sqrt(KT4 * NEB);
SNR_dB = 20*log10(SNR);

% make SNR have two rows
SNR_dB = [SNR_dB; SNR_dB];

Surface_SNR_dB = SNR_dB - SurfaceNF_dB;
Uplink_SNR_dB = Surface_SNR_dB + T_attn_dB;
Downlink_SNR_dB = SNR_dB -  DownlinkNF_dB + T_attn_dB;

    
plot(flipud(rot90(Surface_SNR_dB)),'k'); hold on
plot(flipud(rot90(Uplink_SNR_dB)),'r')
plot(flipud(rot90(Downlink_SNR_dB)),'b')

legend('surface - noisy', 'surface - quiet', ...
    'uplink - noisy', 'uplink - quiet', ...
    'downlink - noisy', 'downlink - quiet')   
break

%axis([0.1 1000 -60 120]);
%set(gca, 'YTick', linspace(-60, 120,10));
%set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
%parms = ['D=', num2str(D), ' depth=50:50:500 & 600:100:1000, Z=',num2str(Z), ', mode=',num2str(mode) ];
%xlabel('Frequency (kHz)','fontsize',12)
%ylabel('Signal/Noise ratio (dB)','fontsize',12)
%title(['Opt Freq: ', parms],'fontsize',12)
%legend([H1 H5], ...
%['no noise, ', num2str(sigma), ' S/m'],...
%'maximum noise',2)

%print(filetype, '', ['opt-' num2str(fig)] )


