% radiometer channel frequencies (lower limit), GHz
channels = [2.35 3.2 4.05 4.9 5.75 6.6 7.45 8.3 9.15 10 10.85 11.7 12.6 13.4 14.25 15.1];
broken_ch_num = [1, 9, 10, 12, 13, 14, 15, 16];
n_max = 4.6e18/1e6; % max electron density, cm-3
B_max = 0.5; % central magnetic foeld, T
% heterodyne frequency, GHz
f_mix = 25;
channel_width = 0.85; %GHz
a = 8.5; % minor radius, cm
R0 = 40; % major radius, cm
% observation line
x = linspace(-a, a, 500);
y = zeros(size(x));
% density profile
n = n_max*(1 - (x.^2 + y.^2)./(a*a));
n(n < 0) = 0;
% magnetic field profile
B = B_max*R0./(R0 + x);
k1 = (4*pi*(4.8032e-10)*(4.8032e-10))/(9.1094e-28); % square electron plasma frequency coeff
k2 = ((4.8032e-10))/((9.1094e-28)*3e10); % electron cyclotron frequency coeff
f_cutoffX = (sqrt(k1*n + power(k2*B*1e4,2)/4) + k2*B*1e4/2)/(1e9*2*pi); % X-mode cut-off frequency, GHz
f_cutoffO = sqrt(k1*n)/(1e9*2*pi); % O-mode cut-off frequency, GHz
f_uhr = sqrt(k1*n + power(k2*B*1e4,2))/(1e9*2*pi); % UHR frequency, GHz
f_ce = k2*B*1e4/(1e9*2*pi); % electron cyclotron frequency, GHz
% ECE 1st harmonic channel location
x_ch1 = zeros(size(channels));
for i = 1:length(channels)
[~, num] = min(abs((f_ce-f_mix) - channels(i)+channel_width/2));
x_ch1(i) = x(num);
end
% ECE 2nd harmonic channel location
x_ch2 = zeros(size(channels));
for i = 1:length(channels)
[~, num] = min(abs((2*f_ce-f_mix) - channels(i)+channel_width/2));
x_ch2(i) = x(num);
end
%% X-mode radiometery
final_intermediate_chanels = channels-channel_width/2+f_mix;
final_intermediate_chanels(abs(x_ch2./a) == 1) = NaN;
figure;
title('X-mode radiometry');
hold on;
plot(x./a, f_ce); %1st ECE harmonic
plot(x./a, 2*f_ce); %2nd ECE harmonic
plot(x./a, 3*f_ce); % %3nd ECE harmonic
plot(x./a, f_cutoffX, 'b'); % X-mode cut-off
plot(x_ch2./a, final_intermediate_chanels, 'ok'); % channels
plot(x_ch2(broken_ch_num)./a, final_intermediate_chanels(broken_ch_num ), 'xr'); % broken channels
plot(x./a, 3*f_ce(end)*ones(size(x)), '--');
grid on;
xlabel('r/a, cm');
ylabel('f, GHz');
legend('f_c_e', '2f_c_e', '3f_c_e', 'f_c_u_t_-_o_f_f_X', 'channels', 'broken channels');
%% O-mode radiometery
final_intermediate_chanels = channels-channel_width/2+f_mix;
final_intermediate_chanels(abs(x_ch1./a) == 1) = NaN;
figure;
title('O-mode radiometry');
hold on;
plot(x./a, f_ce); %1st ECE harmonic
plot(x./a, 2*f_ce); %2nd ECE harmonic
plot(x./a, f_cutoffO, 'b'); % O-mode cut-off
plot(x_ch1./a, final_intermediate_chanels, 'ok'); % channels
plot(x_ch1(broken_ch_num)./a, final_intermediate_chanels(broken_ch_num ), 'xr'); % broken channels
plot(x./a, 2*f_ce(end)*ones(size(x)), '--');
grid on;
xlabel('r/a, cm');
ylabel('f, GHz');
legend('f_c_e', '2f_c_e', 'f_c_u_t_-_o_f_f_O', 'channels', 'broken channels');
% figure;
% %������ ����, ��������������� ������� ������
% g = linspace(-a,a, 1000);
% g1 = sqrt(a*a-g.*g);
% hold on;
% plot(g,g1, 'k');
% plot(g, -g1, 'k');
% ylim([-a, a]);
% xlim([-a, a]);
% hold on;
% contour(x, y, f_ce*2, 20);
% contour(x, y, f_cutoff, 20);
% grid on;
