#!/usr/bin/env python
# -*- coding: utf-8 -*-
from numpy import *
import time
import matplotlib
matplotlib.rcParams['backend'] = 'Agg'
matplotlib.rc('font', size='10')
matplotlib.rc('text', usetex=True)
from matplotlib.pyplot import *
import os
from pygolem_lite.modules import list2array,save_adv,saveconst,load_adv
from pygolem_lite import Shot
LineCount=40000;
TriggerLine=3500;
DeltaTs=1e-6;#[s]
Stab_S_Eff=0.08601;#[m^2]
Mirnov_S_Eff=0.003705;#[m^2]
MirnovPosition=93;#[mm]
MirnovPosition_m=0.093;#[m]
R0=0.4;#m
mu0=0.000001256;
a_limiter=0.085;
a_L=0.085;
def narovnani(field):
aver2 = mean(field[32000:32000+TriggerLine])
field[5000:32000]-= (arange(5000,32000)-5000)*aver2/(32000-5000)
field[32000:]-= aver2
return field
def integrace_mc(field):
#print "====", field, TriggerLine
aver = mean(field[:TriggerLine,:],axis = 0)
integr = cumsum(field-aver, axis = 0, out = field)*DeltaTs/Mirnov_S_Eff
return integr
def integrace_bz(field):
aver = mean(field[:TriggerLine,...],axis = 0)
integr = cumsum(field-aver,axis = 0, out = field)*DeltaTs/Stab_S_Eff
return integr
def load_etalons(Ub_limit,Ust_limit):
path = 'Ub'+str(Ub_limit)+'_Ust'+str(Ust_limit)+'_no-plasma'
try:
tvec,etalon = load_adv('./etalon/'+path)
except:
print 'cant open '+path+' etalon'
saveconst('no_etalon', 0)
exit()
return etalon
def load_pygolem():
Data = Shot()
time.sleep(2) # stupid bug fix with loading of nonsaved file
gd = Data.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)
mc = list2array( Data[das, [m1, m5, m9, m13]][1] )
Ip = array(Data['plasma_current'], copy = False)[1,:].T
try:
Bz = list2array( Data[das, 5][1] )
except:
Bz = zeros(LineCount)
Ub_limit = int(Data['ub'])
Ust_limit = int(Data['ust'])
plasma = Data['plasma']
#Ub_limit = 600#BUG
#Ust_limit = 0
return mc,Bz,Ip,Ust_limit,Ub_limit
def plot_graphs():
mc,Bz,Ip,Stlim,Btlim = load_pygolem()
#print mc,Bz,Ip,Stlim,Btlim
mc = integrace_mc(mc)
Bz = integrace_bz(Bz)
mc = mc[:LineCount,:]
Bz = Bz[:LineCount]
#remove crosstalks
for i in range(4):
mc[:,i] = narovnani(mc[:,i])
etalon = load_etalons(Btlim,Stlim)
mc-= etalon[:,:4] #time mc1 mc5 mc9 mc13 bz
Bz-= etalon[:,-1]
#make calculations
disZ = MirnovPosition*(mc[:,1]-mc[:,3])/(abs(mc[:,1]+mc[:,3]))
disR = MirnovPosition*(mc[:,0]-mc[:,2])/(abs(mc[:,0]+mc[:,2]))
B0 = (mu0*Ip)/(2*pi*MirnovPosition_m)
Lambda = ((mc[:,2]-mc[:,0])/2 - Bz)*(R0/(B0*MirnovPosition_m))-log(MirnovPosition_m/a_L)+1
disR_upgrade=(mc[:,2]-mc[:,0])*MirnovPosition_m/(2*B0)-0.5*(log(MirnovPosition_m/a_L)-1+(Lambda-0.5)*(1+(MirnovPosition_m/a_L)**2))*MirnovPosition_m**2/R0*1000
disR_upgrade[Ip < 50] = 0
disR[Ip < 50] = 0
disZ[Ip < 50] = 0
#plot graphs
tvec = arange(LineCount)*DeltaTs*1000
class MyFormatter(ScalarFormatter):
def __call__(self, x, pos=None):
if pos==0:
return ''
else: return ScalarFormatter.__call__(self, x, pos)
fig = figure(num=None, figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k')
subplots_adjust(hspace=0, wspace = 0)
for i in arange(4):
ax = fig.add_subplot(5,1,i+1)
ax.xaxis.set_major_formatter( NullFormatter() )
ax.yaxis.set_major_formatter( MyFormatter() )
plot(tvec,1000*mc[:,i], label = 'clear signal')
plot(tvec,1000*etalon[:,i+1], label = 'toroidal field')
plot(tvec,1000*(mc[:,i]+etalon[:,i+1]), label = 'signal with toroidal field')
axhline(y = 0, linestyle = '--')
leg = legend(loc='upper left', fancybox=True)
leg.get_frame().set_alpha(0.8)
ylabel('mc'+str(i*4+1)+' B [mT]')
ax.xaxis.set_major_formatter( ScalarFormatter() )
ax.yaxis.set_major_formatter( ScalarFormatter() )
xlabel('t [ms]')
savefig('./mc.png',bbox_inches='tight')
close()
fig = figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k')
plot(tvec,-disZ, label = 'disZ')
plot(tvec,-disR, label = 'disR')
plot(tvec,-disR_upgrade, label = 'disR upgrade')
axhline(y = a_L*1000,linestyle = '--',color = 'g')
axhline(y = -a_L*1000,linestyle = '--',color = 'g')
ylim(-100, 100)
xlim(6, 23)
ylabel('R,Z [mm]')
xlabel('t [ms]')
leg = legend( loc='upper right', fancybox=True)
leg.get_frame().set_alpha(0.8)
savefig('./position.png',bbox_inches='tight')
close()
os.system('convert -resize 150x120\! position.png icon.png')
def main():
plot_graphs()
saveconst('status', 0)
if __name__ == "__main__":
main()
