Contents

Shot numbers

 ShotNo = 39393;
 vacuum_shot = 39390;


 shotNum_vac = vacuum_shot;
 shotNum = ShotNo;

Download data

data = data_golem_coils(shotNum);
data_vac = data_golem_coils(shotNum_vac);

% Coil's normalization parameters
Aeff = 1e-4*[-68.9 -140.70 138.9 -140.4 -68.6 134.5 -134.3 142.5 -67.6 -142.8 -140.4 -138 -76.3 -142.2 -139.8 -139.3]';

% Normalization and integration to get B_theta

data.Bcoils = integ(data.coils, 1e-6)./repelem((Aeff)', size(data.t_coils,1), size(data.t_coils,2));
%data.Bcoils(:,12) = data.Bcoils(:,12)*10;
data_vac.Bcoils = integ(data_vac.coils, 1e-6)./repelem((Aeff)', size(data.t_coils,1), size(data.t_coils,2));
% remove external fields
data.Bnorm = data.Bcoils-data_vac.Bcoils;

Check offset removal by vacuum shot

f = figure('visible', 'off');
    hold on

for i = 1:16
    subplot(4,4,i);
    plot(data.t_coils*1000, data.Bcoils(:,i));
    plot(data.t_coils*1000, data_vac.Bcoils(:,i));
    plot(data.t_coils*1000, data.Bnorm(:,i), 'k');
    title(['coil ' num2str(i)]);
    xlabel('t, ms');
    grid on;
end
    hold off

set(gcf, 'Position', [0 0 1000 800]);
saveas(f, 'offset.png');
close(f)

safety factor calculation

B_theta =(1.25663706e-6)*data.Ip*1000/(2*pi*0.093);

q = (8.5/40)*data.Bt./B_theta;

% discharge time
t1 = find(data.t_basic > data.t_plasma_start, 1, 'first');
t2 = find(data.t_basic < data.t_plasma_end, 1, 'last');

f = figure('visible','off');

hold on
plot(data.t_basic(t1+1000:t2),q(t1+1000:t2));
grid on;
xlabel('t, ms');
ylabel('q');
title('Safety factor');
set(gcf, 'Position', [0 0 1000 800]);

hold off
saveas(f, 'q.png');
close(f)

Spectrogram

step = 128*2;   % window

fh = figure('visible','off');
hold on

%figure;
for i = 1:1
[sp, t, f] = Spectrs(data.coils(t1:t2,i)-data_vac.coils(t1:t2,i), step, 1, step/2, 1);

%subplot(4,4,i);
n = find(f*1000 < 200, 1,'last'); % frequency limit
h = pcolor(t+data.t_coils(t1-1)*1000,f(1:n)*1000, sp(1:n,:));
set(h,'EdgeColor','None');
ylabel('f, kHz');
xlabel('t, ms');

title('Spectrogram');
set(gcf, 'Position', [0 0 1000 800]);

end

hold off
saveas(fh, 'Spectrogram.png');
close(fh)

crosscorrelation

step = 128*2; %window

[r, lag, t] = corr_cross(diff(data.Bnorm(t1:t2,:)), step, 1, step/2, 1);

polarplot gif animation

plot a croscorrelation coefficients for the base coil as an animation

lag_num = step; % lag = 0



radius = 1;
for i = 1:size(r,2) % set the mean radius of the polar plot
    r(lag_num,i,:) = (r(lag_num,i,:) + radius);
end

f = figure('visible','off');

for k = 1:length(t)
    clf;
    y = zeros(1,size(data.Bnorm,2)+1);

    hold on;

    for i = 1:size(data.Bnorm,2)
        x2(i,:) = r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k);
    end
    % normalize, average and prepare data from different coils for plotting
    for i = 1:1 %  display only the first coils's correlation
        x = [r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k) r(lag_num,(i-1)*size(data.Bnorm,2)+1,k)];

        x1= x;
        x1(end) = [];

        x1 = x1 - mean(x1);
        st(k,i) = std(x1); % activity level

        y = x;
        x = (0:16)/8*pi;

        % spline approximation
        [x1,fx, phasespline]=phaseapprox(x, y, 0.9999);

       %subplot(4,4,i);
        hold on;
        polar(x, y, 'or'); %plot points
        polar(x1, fx, 'b');  %plot approximation
        polar(x1, ones(size(fx)), 'k--');  %plot medium level
        text(1.2, 1.5, [num2str(t(k) + data.t_plasma_start) ' ms']) %plot time
        [~, t_num] = min(abs((t(k) + data.t_plasma_start)-data.t_basic));
        text(-1.2, 1.5, ['q = ' num2str(q(t_num))]); %plot q
        xlim([-2 2]);
        ylim([-2 2]);
        grid on;
        set(gcf, 'Position', [0 0 1000 800]);
    end

    save frame
    frame(k) = getframe(gcf);
   % pause(0.3);

end

Save animation as .gif

for n = 1:size(frame,2)
    [A,map] = rgb2ind(frame2im(frame(n)),256);
    if n == 1
        imwrite(A, map, 'animation.gif','gif','LoopCount',Inf,'DelayTime',0.3);
    else
        imwrite(A, map, 'animation.gif','gif','WriteMode','append','DelayTime',0.3);
    end
end

Calculate and plot MHD activity

MHD activity is estimated as standart deviation of coils signals. Points for all coils are ploted. The Black line is a base coil

for k = 1:length(t)
    y = zeros(1,size(data.Bnorm,2)+1);

    hold on;

    for i = 1:size(data.Bnorm,2)
        x2(i,:) = r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k);
    end

    for i = 1:size(data.Bnorm,2)

        x1= r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k);

        x1 = x1 - mean(x1);
        st(k,i) = std(x1); % activity level

    end
end

f = figure('visible','off');

hold on;
plot(t + data.t_plasma_start, st, 'ro');
plot(t + data.t_plasma_start, st(:,1), 'ko-', 'LineWidth', 2);
grid on;
title('MHD activity level');
xlabel('t, ms');
set(gcf, 'Position', [0 0 1000 800]);

saveas(f, 'MHD_activity.png');
close(f)

polarplot

plot timeslices of crosscorrelation coefficients for the base coil

time_step = ceil(length(t)/16);
s = 0;

f = figure('visible','off');

%figure;
for k = 1:time_step:length(t)
    y = zeros(1,size(data.Bnorm,2)+1);

    hold on;

    s = s + 1;
    % normalize, average and prepare data from different coils for plotting
    for i = 1:1%size(data.Bnorm,2)  %  display only the first coils's correlation
        x = [r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k) r(lag_num,(i-1)*size(data.Bnorm,2)+1,k)];

        y = x;
        x = (0:16)/8*pi;

        % spline approximation
        [x1,fx, phasespline]=phaseapprox(x, y, 0.9999);

        subplot(4,4,s);
        hold on;
        polar(x, y, 'or'); %plot points
        polar(x1, fx, 'b');  %plot approximation
        polar(x1, ones(size(fx)), 'k--');  %plot medium level
        text(1, 1.7, [num2str(t(k) + data.t_plasma_start,3) 'ms']) %plot time
        [~, t_num] = min(abs((t(k) + data.t_plasma_start)-data.t_basic));
        text(-1.9, 1.7, ['q =', num2str(q(t_num), 2)]); %plot q
        xlim([-2 2]);
        ylim([-2 2]);
        grid on;
        set(gcf, 'Position', [0 0 1200 1000]);
    end

end


saveas(f, 'pictures.png');
close(f)

Published with MATLAB® R2020b