Contents
Shot numbers
ShotNo = 39420; 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)