#!/usr/bin/env python # -*- coding: utf-8 -*- # Main script of the plasma position algorithm. In this script are data prepared # (removed false signals and integration drifts) and the resultes from Plasmposition.py are finaly plotted # # Autor: Tomas Odstrcil import matplotlib matplotlib.rcParams['backend'] = 'Agg' matplotlib.rc('font', size='10') matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!! import time cas = time.time() from numpy import * import os import sys from pygolem_lite.modules import list2array,deconvolveExp,save_adv,saveconst from scipy.signal import fftconvolve,gaussian from scipy.interpolate import interpolate from pygolem_lite import Shot print 'including time: ', time.time()-cas #constants a = 0.085 #[m] R_0 = 0.4 #[m] shot = Shot()['shotno'] if shot > 12079: calib =-array([ 261, 261, -261, 261]) #T/(V*s) else: calib =-array([ 261, 261, -261, -261]) #T/(V*s) def LoadData(): Data = Shot() t0 = time.time() Bt = array(Data['toroidal_field'], copy = False).T Uloop = array(Data['loop_voltage'], copy = False).T Ip = array(Data['plasma_current'], copy = False).T t0 = time.time() gd = Shot().get_data das, m1 = gd('any', 'mirnov_1', return_channel = True) das, m5 = gd('any', 'mirnov_5', return_channel = True) das, m9 = gd('any', 'mirnov_9', return_channel = True) das, m13 = gd('any', 'mirnov_13', return_channel = True) t1 = time.time() Papouch = list2array( Data[das, [m1, m5, m9, m13]] ) #FIXME odkazuje to přímo na kanál, mělo by to odkazovat na mirnonovy cívky CD_trigger = Data['Tcd'] BD_trigger = Data['Tbd'] Bt_trigger = Data['Tb'] CD_voltage = Data['Ucd'] BD_voltage = Data['Ubd'] Bt_voltage = Data['Ub'] plasma_start = Data['plasma_start'] plasma_end = Data['plasma_end'] plasma = Data['plasma'] if CD_voltage == 0: CD_trigger = nan if BD_voltage == 0: BD_trigger = nan if Bt_voltage == 0: Bt_trigger = nan shot = Shot()['shotno'] return Bt,Uloop,Ip,Papouch,CD_trigger,BD_trigger, Bt_trigger,plasma_end,plasma_start, plasma,shot def Resample(tvec,vec,t_min,t_max,n): std = (t_max-t_min)/(tvec[1]-tvec[0])/n/2 gauss_win = gaussian(std*3,std) gauss_win/= sum(gauss_win) vec_smooth = fftconvolve(vec,gauss_win, mode = 'same') tvec_new = linspace(t_min,t_max,n) vec_new = interp(tvec_new,tvec, vec_smooth, left=0, right=None) #f = interpolate.interp1d(tvec,vec_smooth,bounds_error = False, fill_value = 0, copy = False) #vec_new = f(tvec_new) #plot(tvec_new,vec_new, '--') #plot(tvec,vec_smooth) #plot(tvec,vec) #show() return tvec_new,vec_new def CorrectRCcircuit(tvec,sig, tau): N = len(sig) dt = (tvec[-1]-tvec[0])/N sig-= mean(sig) f = fft.fftfreq(N, d=dt) fsig = fft.fft(sig) fsig[0] = 0 #=> mean(x) = 0 ifac = 1-exp(-2*pi*1j*f*dt) #invert integ factor q = exp(-(1./tau+2*pi*1j*f)*dt) fexp = (1-q**(N/2-1))/(1-q) fexp/= fexp[0] #cca dt/tau filter = ones_like(fsig) filter/= fexp#deconvolution filter[1:]/= ifac[1:] #integration integ_sig = real(fft.ifft(fsig*filter)) integ_sig -= integ_sig[0]+sig[0] integ_sig*= dt ##O(n) calculation vfor realtime aplication #integ = zeros_like(sig) #xn = sig[0] #for j in xrange(step,N,step): #xn*= exp(-dt*step/tau) #deriv = sig[j]-xn #integ[j] = deriv+integ[j-step] #xn = sig[j] #integ*= dt/tau return integ_sig def get_position(): from PlasmaPosition import CalcPlasmaPosition,RemoveDriftsAuto cas = time.time() Bt,Uloop,Ip,signal,CD_trigger,BD_trigger, Bt_trigger,plasma_end,plasma_start, plasma,shotno = LoadData() tvec = signal[:,0] signal = signal[:,1:] t0 = tvec[0] dt = tvec[1]-tvec[0] #BUG 1. chanell has RC integrator if shotno >= 11381 and shotno <= 11837: signal[:,0] = CorrectRCcircuit(tvec,signal[:,0], 0.1e-3) signal[:,1:] = cumsum(signal[:,1:],axis = 0, out=signal[:,1:])*dt #možná by šla udělat nějaká lepší integrace else: signal = cumsum(signal,axis = 0, out=signal)*dt #možná by šla udělat nějaká lepší integrace print 'load data: ', time.time()-cas reduce = len(signal)/400 # calculate in 400 points if mean(abs(Uloop[:,1])) == 0 or mean(abs(Bt[:,1])) == 0 : #selhání diagnostik print 'all diagnostics failured ' return None #correct effects of the chamber currents (projection constants) Uproj = array([-0.5e-07, -1.0e-07, -0.2e-07,-0.8e-07]) U_loop_ds = interp(tvec, Uloop[:,0],Uloop[:,1], left=0, right=None) #f_uloop = interpolate.interp1d(Uloop[:,0],Uloop[:,1],bounds_error = False, fill_value = 0, copy = False) #U_loop_ds = f_uloop(tvec) for i in range(4): signal[:,i]-= Uproj[i]*U_loop_ds Bt_crosstalk = True Uloop_crosstalk = False if abs(median(Bt[:,1])) < 0.001: Bt_crosstalk = False if abs(median(Uloop[:,1])) < 0.05: Uloop_crosstalk = False #prepare signals E_trigger = nanmin([CD_trigger, BD_trigger]) #remove crosstalks+drifts signal,err, AutoRemoveGraph = RemoveDriftsAuto(vstack((tvec,signal.T)).T, Bt,Uloop,plasma_start,plasma_end,Bt_trigger,E_trigger) #now is a plasma presence is necassary if not plasma: return [AutoRemoveGraph,Bt_trigger] #create detectors rho = empty((4,1)) zeta = empty((4,1)) r_det = 0.093 phi_det = linspace(0,1.5*pi,4 ) rho[:,0] = r_det*cos(phi_det)+R_0 zeta[:,0] = r_det*sin(phi_det) detectorPos = hstack((rho,zeta)) tvec_ds, Ip_ds = Resample(Ip[:,0],Ip[:,1],plasma_start,plasma_end,reduce) detectorSignal = signal[:,1:] tvec = signal[:,0] detectorSignal_ds = list() for i in range(size(detectorSignal,1)): #print shape(tvec), shape(detectorSignal[:,i]) tvec_ds,data_ds = Resample(tvec,detectorSignal[:,i],plasma_start,plasma_end,reduce) detectorSignal_ds.append(data_ds) detectorSignal_ds = array(detectorSignal_ds, copy = False).T # estimated errors #TODO here can by applied information about lower precision of some detector IpDriftError =7 detectorDriftError = ones(4)* 1.5e-7 detectorDriftError*= abs(array(calib)) detectorSignal_ds*= array(calib) pos_tvec, position, radius, residuum, retrofit, data,chi2 = CalcPlasmaPosition(single(tvec_ds),single(detectorPos), single(detectorSignal_ds), single(detectorDriftError), single(Ip_ds), single(IpDriftError)) save_adv('results/R_position', pos_tvec, position[:,0]) save_adv('results/Z_position', pos_tvec, position[:,1]) save_adv('results/plasma_radius', pos_tvec, radius) savetxt('results/R_position.txt', vstack((pos_tvec, position[:,0])).T, fmt='%.4e %.3e' ) savetxt('results/Z_position.txt', vstack((pos_tvec, position[:,1])).T, fmt='%.4e %.2e' ) savetxt('results/plasma_radius', vstack((pos_tvec, radius)).T, fmt='%.4e %.2e' ) #savetxt('results/residuum', mean(residuum), fmt='%.1e') saveconst('results/residuum', median(residuum/ data[:,0]**2*1e6)) #f_Ip = interpolate.interp1d( Ip[:,0],Ip[:,1] , bounds_error = False, fill_value = 0,copy = False) Ip = interp(tvec, Ip[:,0],Ip[:,1]) return pos_tvec, position, residuum, retrofit, single(data),single(tvec),single(Ip), IpDriftError, single(detectorSignal), detectorDriftError, plasma_start, plasma_end,Bt_trigger,chi2,AutoRemoveGraph def graphs( inputs , filetype = 'png' ): if len(inputs) == 2: #vacuum discharge [AutoRemoveGraph,Bt_trigger] = inputs else: [pos_tvec,position,residuum,retrofit,data,tvec,Ip,IpDriftError, detectorSignal,detectorDriftError,plasma_start,plasma_end,Bt_trigger,chi2,AutoRemoveGraph] = inputs #=============== Plot results====================== import matplotlib matplotlib.rcParams['backend'] = 'Agg' matplotlib.rc('font', size='10') matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!! import matplotlib.pyplot as plt class MyFormatter(plt.ScalarFormatter): def __call__(self, x, pos=None): if pos==0: return '' else: return plt.ScalarFormatter.__call__(self, x, pos) #plot graph from auto removing of crosstalks t1 = time.time() (frames, vlines,hlines,rectangles) = AutoRemoveGraph fig = plt.figure(num=None, figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k') plt.subplots_adjust(hspace=0, wspace = 0) for i,frame in enumerate(frames): (t, curves, ytext,text) = frame ax = fig.add_subplot(len(frames),1,i+1) ax.xaxis.set_major_formatter( plt.NullFormatter() ) ax.yaxis.set_major_formatter( MyFormatter() ) ax.text(.05,.1,text,horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes,backgroundcolor = 'w') for curve in curves: (y,label, linestyle) = curve plt.plot(t,y, linestyle , label = label) for hline in hlines: (y0,linestyle) = hline plt.axhline(y = y0,linestyle = linestyle) for vline in vlines: (x0,linestyle) = vline plt.axvline(x = x0,linestyle = linestyle) plt.ylabel(ytext) plt.axis('tight') plt.xlim(Bt_trigger*1000, None) (y_min,y_max) = plt.ylim() for (x_min,x_max) in rectangles: r = plt.Rectangle((x_min, y_min),x_max-x_min,y_max-y_min) r.set_clip_box(ax.bbox) r.set_alpha(0.05) ax.add_artist(r) handles, labels = ax.get_legend_handles_labels() handles.append(r) labels.append('Minimized interval') ax.xaxis.set_major_formatter( plt.ScalarFormatter() ) plt.xlabel('t [ms]') leg = plt.legend(handles, labels,loc='best', fancybox=True) leg.get_frame().set_alpha(0.5) plt.savefig('./graphs/signal_correction.'+filetype,bbox_inches='tight') plt.close() print ' plot graph from auto removing of crosstalks ',time.time()-t1 if len(inputs) == 2: #no plasma return n_det = size(detectorSignal,1) t0 = time.time() t = time.time() #plot position fig = plt.figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k') plt.plot(pos_tvec*1e3, (position[:,0]-0.4)*100,'b', label = 'R-R$_0$', linewidth=0.5) plt.plot(pos_tvec*1e3, position[:,1]*100,'r', label = 'Z', linewidth=0.5) plt.axhline(y = a*100, linestyle = '--' ,label = 'limiter') plt.axhline(y = -a*100, linestyle = '--' ) leg = plt.legend(loc='lower left', fancybox=True) leg.get_frame().set_alpha(0.7) plt.xlim(pos_tvec[0]*1e3, pos_tvec[-1]*1e3) plt.ylim(-10,10) plt.ylabel('R,Z [cm]') plt.xlabel('t [ms]') plt.subplots_adjust( bottom=0.2) plt.savefig('graphs/plasma_position.'+filetype,bbox_inches='tight') fig.set_figwidth(4) fig.set_figheight(3) fig.dpi = 40 leg.get_frame().set_alpha(0) plt.savefig('icon.png', bbox_inches='tight', dpi= 40) plt.close() #plot position polar r = hypot(position[:,0]-0.4,position[:,1]) phi = arctan2(position[:,1],-(position[:,0]-0.4)) fig = plt.figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k') axL = fig.add_subplot(1,1,1) plt.plot(pos_tvec*1e3, r*100,'b', label='radial coordinate $\\rho$', linewidth=0.5) axL.set_ylabel('$\\rho$ [cm]') axL.set_ylim(0,10) axL.set_xlim(pos_tvec[0]*1e3, pos_tvec[-1]*1e3) axL.set_xlabel('t [ms]') axL.yaxis.tick_left() handlesL, labelsL = axL.get_legend_handles_labels() axR = fig.add_subplot(1,1,1, sharex=axL, frameon=False) axR.yaxis.tick_right() axR.yaxis.set_label_position("right") plt.plot(pos_tvec*1e3, phi,'r', label='angular coordinate $\phi$', linewidth=0.5) axR.set_ylabel('$\phi$ [rad]') axR.set_xlim(pos_tvec[0]*1e3, pos_tvec[-1]*1e3) axR.set_ylim(-pi,pi) handlesR, labelsR = axR.get_legend_handles_labels() handles = hstack((handlesL,handlesR)) labels = hstack((labelsL,labelsR)) leg = axL.legend(handles,labels,loc='lower left', fancybox=True) #leg = plt.legend(loc='lower left', fancybox=True) leg.get_frame().set_alpha(0.7) #plt.ylim(-10,10) #plt.xlabel('t [ms]') plt.subplots_adjust( bottom=0.2) fig.savefig('graphs/plasma_position_polar.'+filetype,bbox_inches='tight') #fig.set_figwidth(4) #fig.set_figheight(3) #fig.dpi = 40 #leg.get_frame().set_alpha(0) #plt.savefig('icon.png', bbox_inches='tight', dpi= 40) plt.close() print 'plot position ',time.time()-t t = time.time() #plot residuum fig = plt.figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k') plt.semilogy(pos_tvec*1e3, residuum/ data[:,0]**2*1e6, 'r', label = 'residuum of the fit/$I_p^2$') plt.axhline(y = 10) plt.axhline(y = 1e2) plt.text(mean(pos_tvec*1e3),3 , r'reliable fit', color='k') plt.text(mean(pos_tvec*1e3),50 , r'considerable fit', color='k') plt.text(mean(pos_tvec*1e3),1e3 , r'poor fit', color='k') plt.xlim(pos_tvec[0]*1e3, pos_tvec[-1]*1e3) plt.ylim(1,1e4) plt.xlabel('t [ms]') plt.ylabel('SSE/$I^2_p$ [kA$^{-2}$]') plt.subplots_adjust( bottom=0.2) leg = plt.legend(loc='upper left', fancybox=True) leg.get_frame().set_alpha(0.5) plt.savefig('graphs/residuum.'+filetype, bbox_inches='tight') plt.close() print 'plot residuum ',time.time()-t t = time.time() cas = time.time() # plot plasma radius t = time.time() # plot plasma radius radius = a-sqrt((position[:,0]-0.4)**2+position[:,1]**2) fig = plt.figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k') plt.plot(pos_tvec*1e3, radius*100,'b', label = 'Plasma radius', linewidth=0.5) plt.xlim(pos_tvec[0]*1e3, pos_tvec[-1]*1e3) plt.ylim(0,10) plt.axhline(y = a*100, linestyle = '--' ,label = 'limiter') plt.ylabel('r [cm]') plt.xlabel('t [ms]') leg = plt.legend(loc='upper right', fancybox=True) leg.get_frame().set_alpha(0.7) plt.subplots_adjust( bottom=0.2) plt.savefig('graphs/plasma_radius.'+filetype,bbox_inches='tight') plt.close() print ' plasma radius ',time.time()-t t = time.time() #plot raw data fig = plt.figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k') plt.plot(tvec[:-50]*1000,Ip[:-50]/IpDriftError,'-.', label = '$I_p$') for i in range(4): plt.plot(tvec*1e3,detectorSignal[:,i]/detectorDriftError[i]*calib[i], label = 'mc'+str(i*4+1), linewidth=0.5) plt.axvline(x = plasma_start*1000, linestyle= '--') plt.axvline(x = plasma_end*1000, linestyle = '--') plt.axhline(y = 0, linestyle = '-.') plt.xlabel('t [ms]') plt.ylabel('signal/estimated error [-]') leg = plt.legend(loc='upper left', fancybox=True) leg.get_frame().set_alpha(0.7) plt.subplots_adjust( bottom=0.2) plt.axis('tight') plt.xlim(Bt_trigger*1000, min(1000*plasma_end+5, 40)) plt.ylim(-10,None) plt.savefig('graphs/preprocesed_signal.'+filetype,bbox_inches='tight') print ' plot raw data',time.time()-t plt.close('all') t = time.time() ## plot plasma profile R,Z #time_horizont_profile = empty((50,size(data,0))) #time_vertical_profile = empty((50,size(data,0))) #x = linspace(-a,a,50) #for i in range(size(data,0)): #circle = 1-((x-position[i,0]+0.4)/radius[i])**2 #circle*= (1+sign(circle))/2 #time_horizont_profile[:,i] = sqrt(circle) #for i in range(size(data,0)): #circle = 1-((x-position[i,1])/radius[i])**2 #circle*= (1+sign(circle))/2 #time_vertical_profile[:,i] = sqrt(circle) #fig = figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k') #pcolor( pos_tvec*1e3, (x+0.4)*100, time_horizont_profile) #xlabel('t [ms]') #ylabel( 'R [cm]') #axis('tight') #subplots_adjust( bottom=0.2) #savefig('graphs/plasma_profile_r.png') #close() #fig = figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k') #pcolor( pos_tvec*1e3, x*100, time_vertical_profile) #xlabel('t [ms]') #ylabel( 'Z [cm]') #axis('tight') #subplots_adjust( bottom=0.2) #savefig('graphs/plasma_profile_z.png') #close() #plot retrofit from scipy.stats.mstats import mquantiles fig = plt.figure(num=None, figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k') plt.subplots_adjust(hspace=0, wspace = 0) c = ('r', 'g', 'b', 'y') y_min_lim = amin(data[:,1:]) y_max_lim = amax(data[:,1:]) for i in range(1,n_det+1): ax = fig.add_subplot(n_det+1,1,i) ax.xaxis.set_major_formatter( plt.NullFormatter() ) ax.yaxis.set_major_formatter( MyFormatter() ) plt.plot(pos_tvec*1e3, data[:,i]/data[:,0]*1e6, c[i%4], label = 'mc'+str((i-1)*4+1)+'/I$_p$') plt.plot(pos_tvec*1e3, retrofit[:,i]/data[:,0]*1e6, c[i%4]+'-.', label = 'retrofit') ax.text(.05,.1,'$\\chi^2$ = %2.1f'% chi2[i],horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes) leg = plt.legend(loc='upper left', fancybox=True) leg.get_frame().set_alpha(0.5) data_max = mquantiles(data[:,i]/data[:,0]*1e6, 0.95)*2 plt.ylim(0, data_max) plt.ylabel('mc'+str((i-1)*4+1)+'/$I_p$ [mT/kA]') ax = fig.add_subplot(n_det+1,1,n_det+1) plt.plot(pos_tvec*1e3, ones(size(data,0)), label = 'I$_p$/I$_p$') plt.plot(pos_tvec*1e3, retrofit[:,0]/data[:,0], '--', label = 'retrofit') ax.text(.05,.1,'$\\chi^2$ = %2.1f'% chi2[0],horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes) leg = plt.legend(loc='upper left', fancybox=True) leg.get_frame().set_alpha(0.5) plt.ylim(0,2) plt.ylabel('$I_p/I_p$ [-]') plt.xlabel('t [ms]') plt.savefig('./graphs/retrofit.'+filetype,bbox_inches='tight') plt.close() print 'retrofit ',time.time()-t t = time.time() #plot animation -- slow!!! #print "plotting animation" #import MagFieldProfileGen #MagFieldProfileGen.PlotProfile(pos_tvec,data[:,0],position, 'animation') print 'graphs plotted in %g s' % (time.time()-t0) def main(): for path in ['graphs', 'results','constants' ]: if not os.path.exists(path): os.mkdir(path) if sys.argv[1] == "analysis": out = get_position() save('out', out) if sys.argv[1] == "plots": out = load('out.npy') os.remove('out.npy') graphs(out, 'png') saveconst('status', 0) if __name__ == "__main__": main()