Source code :: VAcharSweeped

# coding: utf-8

# # Initials

# In[1]:

import matplotlib
matplotlib.use('Agg') 

import matplotlib.pyplot as plt
import numpy as np
from urllib import urlopen
import os
from IPython import get_ipython
import string
from scipy.signal import argrelextrema

def is_interactive(): # are we in jupyter ??
    import __main__ as main
    return not hasattr(main, '__file__')



if is_interactive():
    get_ipython().magic(u'matplotlib inline')
    
   


baseURL = "http://golem.fjfi.cvut.cz/utils/data/" #global
#baseURL = "/golem/database/operation/shots/" #local

        
# ack http://stackoverflow.com/questions/13728392/moving-average-or-running-mean
    
def running_mean(l, N):
    # Also works for the(strictly invalid) cases when N is even.
    if (N//2)*2 == N:
        N = N - 1
    front = np.zeros(N//2)
    back = np.zeros(N//2)
    for i in range(1, (N//2)*2, 2):
        front[i//2] = np.convolve(l[:i], np.ones((i,))/i, mode = 'valid')
    for i in range(1, (N//2)*2, 2):
        back[i//2] = np.convolve(l[-i:], np.ones((i,))/i, mode = 'valid')
    return np.concatenate([front, np.convolve(l, np.ones((N,))/N, mode = 'valid'), back[::-1]])

def mkdir(dir):  
    try:os.makedirs(dir) 
    except OSError:pass

mkdir('IndivVAchars');mkdir('ReferenceShot')

SmoothCoefficient=100
HWNapetovyDelic=100 
HWProudovyResistor=200 # Ohm
HFSweepFrequency=1000 #Hz
CurrentScale=1000 # so .. mA

FigSize=10,8


# # Get data

# In[10]:

#ShotNo=string.replace(os.path.basename(os.getcwd()),'#','') # get ShotNumber from dir name where we are
ShotNo='0' #Last shot
#ShotNo='23034' #Last shot
#Plasma parameters
PlasmaStart=int(float(np.loadtxt(urlopen(baseURL+ShotNo+'/plasma_start')))*1e6) # in us
PlasmaEnd=int(float(np.loadtxt(urlopen(baseURL+ShotNo+'/plasma_end')))*1e6) # in us
ShotNumber=int(np.loadtxt(urlopen(baseURL+ShotNo+'/shotno')))
ReferenceShot=ShotNumber # sorry




# In[13]:

#plt.ylim(ylim_min,ylim_max);
os.system('wget '+baseURL + str(ReferenceShot) + '/bpp_current -O ReferenceShot/bpp_current.txt');
os.system('wget '+baseURL + str(ReferenceShot) + '/bpp_voltage -O ReferenceShot/bpp_voltage.txt');

VAchar_voltage=np.loadtxt('ReferenceShot/bpp_voltage.txt')[PlasmaStart:PlasmaEnd]
VAchar_current=np.loadtxt('ReferenceShot/bpp_current.txt')[PlasmaStart:PlasmaEnd]


if is_interactive():plt.figure(figsize=(10, 8), dpi= 80, facecolor='w', edgecolor='k')
plt.plot(VAchar_voltage[:,1],VAchar_current[:,1],'.')
plt.xlabel('bpp_voltage U [V]');plt.ylabel('bpp_current U [V]')
plt.title('BPP raw XY data: bpp_current vs bpp_voltage')
plt.savefig('ReferenceShot/RawData.jpg', bbox_inches='tight')
if is_interactive():plt.show();
plt.close();


# In[9]:

print "Referencni vyboj ..."
os.system('wget http://golem.fjfi.cvut.cz/shots/' + str(ReferenceShot) + '/basicdiagn/Btoroidal.npz -O ReferenceShot/Btoroidal.npz');
os.system('wget http://golem.fjfi.cvut.cz/shots/' + str(ReferenceShot) + '/basicdiagn/Iplasma.npz -O ReferenceShot/Iplasma.npz');
os.system('wget http://golem.fjfi.cvut.cz/shots/' + str(ReferenceShot) + '/basicdiagn/Uloop.npz -O ReferenceShot/Uloop.npz');
os.system('wget http://golem.fjfi.cvut.cz/shots/' + str(ReferenceShot) + '/basicdiagn/graphpres.png -O ReferenceShot/graphres.png');


Btoroidal=np.load('ReferenceShot/Btoroidal.npz');
Iplasma=np.load('ReferenceShot/Iplasma.npz');
Uloop=np.load('ReferenceShot/Uloop.npz');

f,ax = plt.subplots(3,sharex=True);plt.subplots_adjust(hspace=0.001)
ax[0].set_title('#' + str(ReferenceShot))
ax[0].plot(Uloop['data']);ax[0].set_ylabel('$U_l$ [V]')
ax[1].plot(Btoroidal['data']);ax[1].set_ylabel('$B_t$ [T]')
ax[2].plot(Iplasma['data']/1000);ax[2].set_ylabel('$I_p$ [kA]')
plt.savefig('ReferenceShot/ReferenceShot.jpg', bbox_inches='tight')
#plt.show();
plt.close();

f = plt.figure(figsize=(20.0, 5.0))
f,ax = plt.subplots(5,sharex=True);plt.subplots_adjust(hspace=0.001)
f.set_size_inches(FigSize)
ax[0].set_title('#' + str(ReferenceShot))
ax[0].plot(Uloop['data']);ax[0].set_ylabel('$U_l$ [V]')
ax[1].plot(Btoroidal['data']);ax[1].set_ylabel('$B_t$ [T]')
ax[2].plot(Iplasma['data']/1000);ax[2].set_ylabel('$I_p$ [kA]')
ax[3].plot(np.loadtxt('ReferenceShot/bpp_voltage.txt')*100);ax[3].set_ylabel('$U$ [V]')
ax[4].plot(np.loadtxt('ReferenceShot/bpp_current.txt')/HWProudovyResistor*CurrentScale);ax[4].set_ylabel('$I$ [mA]')
ax[4].set_ylim(-5,10)
ax[4].set_xlim(PlasmaStart*8/10,PlasmaEnd*11/10)
plt.savefig('ReferenceShot/ReferenceShotWithBPP.jpg', bbox_inches='tight')
if is_interactive():plt.show();
plt.close();


# # Finding sweeping intervals

# In[85]:

VAchar_voltage_smoothed=running_mean(VAchar_voltage[:,1]*HWNapetovyDelic,SmoothCoefficient)
VAchar_current_smoothed=running_mean(VAchar_current[:,1]/HWProudovyResistor*CurrentScale,SmoothCoefficient)
ylim_min=np.min(VAchar_current_smoothed)
ylim_max=np.max(VAchar_current_smoothed)
xlim_min=np.min(VAchar_voltage_smoothed)
xlim_max=np.max(VAchar_voltage_smoothed)

#ack http://stackoverflow.com/questions/4624970/finding-local-maxima-minima-with-numpy-in-a-1d-numpy-array
x = np.linspace(PlasmaStart,PlasmaEnd,(PlasmaEnd-PlasmaStart))
data=VAchar_voltage_smoothed
extrems_tmp = (np.diff(np.sign(np.diff(data))).nonzero()[0] + 1) # local min+max
minims = ((np.diff(np.sign(np.diff(data))) > 0).nonzero()[0] + 1) # local min
maxims = ((np.diff(np.sign(np.diff(data))) < 0).nonzero()[0] + 1) # local max


# graphical output...
from pylab import *
plt.figure(figsize=(FigSize), dpi= 80, facecolor='w', edgecolor='k')
plt.plot(x,data)
plt.plot(x,VAchar_voltage[:,1]*HWNapetovyDelic)
plt.plot(x[minims], data[minims], "o", label="min")
plt.plot(x[maxims], data[maxims], "o", label="max")
plt.ylabel('bpp_voltage U [V]');plt.xlabel('time [us]')
plt.title('Maxima and minima identification')
plt.legend()
plt.savefig('ReferenceShot/MaximMinims.jpg', bbox_inches='tight')
if is_interactive():plt.show()

# Sometimes double extrems appear, removal:
#print extrems
extrems=[]
for i in range(len(extrems_tmp)-1): 
    if abs(extrems_tmp[i]-extrems_tmp[i+1])>10: 
        extrems.append(extrems_tmp[i])
#print extrems_corr        


# In[104]:

AllVAchars=[]
#for i in range(2):  # tuning purposes
for i in range(len(extrems)-1): 
    fig, ax1 = plt.subplots();
    AllVAchars.append([VAchar_voltage[extrems[i]:extrems[i+1]][:,1]*HWNapetovyDelic,VAchar_current_smoothed[extrems[i]:extrems[i+1]]])
    ax1.plot(AllVAchars[i][0],AllVAchars[i][1],'.')
    plt.grid(True)
    ax1.set_xlim(xlim_min,xlim_max);
    ax1.set_ylim(ylim_min,ylim_max);
    ax1.set_title('VAchar #'+str(ShotNumber)+' t=<'+str(PlasmaStart+extrems[i])+','+str(PlasmaStart+extrems[i+1])+'> us')
    ax1.set_xlabel('$U$ [V]')
    ax1.set_ylabel('$I$ [mA]')
    ax1.axhline(0);
    ax2 = fig.add_axes([0.25, 0.45, 0.3, 0.18]);
    ax2.set_xlim(PlasmaStart,PlasmaEnd);
    ax2.set_ylim(0,6000);
    ax2.set_yticks([j for j in xrange(0,6000,2000)])
    ax2.set_yticklabels([j for j in xrange(0,6,2)])
    ax2.set_xticks([j for j in xrange(PlasmaStart,PlasmaEnd,5000)])
    ax2.set_xticklabels([j for j in xrange(PlasmaStart/1000,PlasmaEnd/1000,2)])
    ax2.plot(Iplasma['data']);
    ax2.axvline(PlasmaStart+(extrems[i]+extrems[i+1])/2);
    ax2.set_ylabel('$I_p$ [kA]')
    ax3 = fig.add_axes([0.25, 0.63, 0.3, 0.18]);
    ax3.set_yticks([(j/10.0) for j in xrange(0,6,2)])
    ax3.set_xticklabels([]);
    ax3.set_ylim(0,0.55);
    ax3.set_xlim(PlasmaStart,PlasmaEnd);
    ax3.plot(Btoroidal['data']);
    ax3.axvline(PlasmaStart+(extrems[i]+extrems[i+1])/2);
    ax3.set_ylabel('$B_t$ [T]')
    ax2.set_xlabel('$t$ [ms]')
    #plt.show(); #just for tuning
    plt.savefig('IndivVAchars/VA'+str(i).zfill(2) +':'+str(extrems[i])+'_'+str(extrems[i+1])+'.jpg', bbox_inches='tight')
    plt.close();


# In[ ]:




# In[34]:

os.system('rm index.html');
fileid = open('index.html','a+')
fileid.write('<html><head><title>Title</title><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><style></style></head><body><center>')
fileid.write('<h2>BPP experiments @ the tokamak GOLEM</h2>')
fileid.write('<h2>Experimental setup</h2>')
fileid.write('<h3>Ball pen probe @ North-East port</h3>')
fileid.write('<center><a href="/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/setup/ExpSetup-BPP.png"><img src="/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/setup/ExpSetup-BPP.png" width="30%"></a><br><a href="/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/setup/">setup</a>')
fileid.write('<h2>The discharge #'+str(ReferenceShot)+'</h2>')
fileid.write('<h3>Basic diagnostics</h3><a href="http://golem.fjfi.cvut.cz/shots/'+str(ReferenceShot)+'/"><img src="ReferenceShot/graphres.png" width="30%"></a><br/><h3>The BPP data vs diagnostics</h3><a href="http://golem.fjfi.cvut.cz/shots/'+str(ReferenceShot)+'/DAS/1011Papouch_St.ON/"><img src="ReferenceShot/ReferenceShotWithBPP.jpg" width="30%"></a><br/><a href="ReferenceShot/">Reference shot data</a><br/><h3>The BPP raw data</h3><a href="http://golem.fjfi.cvut.cz/shots/'+str(ReferenceShot)+'/DAS/1011Papouch_St.ON/"><img src="ReferenceShot/RawData.jpg" width="30%"></a><br/>')
fileid.write('<h2>Data manipulation</h2><h3>Maxima and minima localization in bpp_voltage</h3><img src="ReferenceShot/MaximMinims.jpg" width="30%"><br/>')
fileid.write('<h3>Final movie</h3>')
fileid.write('<img src="finalmovie.gif"><br/><a href="IndivVAchars/">Individual figures</a></br>')
for i in range(len(extrems)-1): 
    fileid.write('<img src="IndivVAchars/VA'+str(i).zfill(2) +':'+str(extrems[i])+'_'+str(extrems[i+1])+'.jpg" width="20%">')
fileid.write('<h2>Experimental Photo</h2>') 
fileid.write('<center><img src="setup/2017-01-19_23-11-12_w.jpg"></center>\n\n')
fileid.write('<h2>Resources</h2>')
fileid.write('<ul><li><a href="VAcharSweeped.ipynb">Jupyter notebook</a><li><li><a href="VAcharSweeped.py">Jupyter notebook python export</a><li><li><a href="VAcharSweeped.html">Jupyter notebook html export</a><li><li><a href="VAcharSweeped.pdf">Jupyter notebook pdf export</a><li></ul>')             
fileid.write('</body></html>')
fileid.close()


# In[ ]: