#!/usr/bin/env python
# -*- coding: utf-8 -*-
import time
t = time.time()
import os,sys
import urllib2
from CalcIonProjections import GapsFilling,CalcProjections
from IDA import IntegLightAnalyz
from CalcDensTemp import CalcDensTemp,plotDensTemp
from SpectrometerControl import Tokamak, OceanOptics_Spec
from numpy import *
import matplotlib
from scipy.linalg import norm
from scipy.signal import fftconvolve
matplotlib.rcParams['backend'] = 'Agg'
#matplotlib.rc('font', size='10')
#matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!!
#from matplotlib.pyplot import *
import matplotlib.pyplot as plt
from pygolem_lite.modules import load_adv,get_page_paths,saveconst,save_adv
from pygolem_lite import Shot
print 'including time', time.time()-t
def LoadData_web(shot_num):
print 'downloading '+str(shot_num)
url = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/diagnostics/Radiation/1111Spectrometer.ON/data/spectra.txt_Data_.txt'
try:
#print './data/PlasmaStart_'+str(shot_num)
plasma_start = loadtxt('./data/PlasmaStart_'+str(shot_num))
plasma_end= loadtxt('./data/PlasmaEnd_'+str(shot_num))
except:
print 'missing plasma start or end = no plasma?'
plasma_start = 0
plasma_end = inf
#print os.path.isfile('data_'+str(shot_num)+'.txt')
for i in range(10):
try:
if not os.path.isfile('data_'+str(shot_num)+'.txt'):
u = urllib2.urlopen(url)
localFile = open('data_'+str(shot_num)+'.txt', 'w')
localFile.write(u.read())
localFile.close()
break
except:
print 'spectra not found, waiting ...'
time.sleep(1)
print 'data waiting time %ds'%i
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
shot = HR2000.loadData('','Python_Data_file','data_'+str(shot_num))
return HR2000, shot, plasma_start,plasma_end
def LoadData():
Data = Shot()
plasma_start = Data['plasma_start']
plasma_end = Data['plasma_end']
plasma = Data['plasma']
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
file_path = Shot().get_data('spectrometr:data', return_path=True)
nmax = 100
for i in range(nmax):
try:
shot = HR2000.loadData('','Python_Data_file',file_path[:-4])
break
except:
print 'no spectroscopic data '+str(i)
time.sleep(1)
if i == nmax-1:
raise Exception('no data found!')
return HR2000, GOLEM, shot, plasma_start,plasma_end
def LoadData_offline(shot_num):
print 'Loading shot ', shot_num
#try:
##print './data/PlasmaStart_'+str(shot_num)
#plasma_start = loadtxt('./data/PlasmaStart_'+str(shot_num))
#plasma_end= loadtxt('./data/PlasmaEnd_'+str(shot_num))
#except:
#print 'missing plasma start or end = no plasma?'
plasma_start = 0
plasma_end = inf
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
shot = HR2000.loadData('','Python_Data_file','./net_data/data_'+str(shot_num))
return HR2000, shot, plasma_start,plasma_end
def LoadDataIDA(plasma_start,plasma_end,spectrometer ): #TODO není dokončené načítání
#LOAD photodiod ========================================
photoLabel = list()
photoData = list()
photoError = list()
photoTvec = list()
photoFilter = list()
photoSensitivity = list()
photoTvecError = list()
correl_list = list() # correlation of the errors in time vectors
Data = Shot()
photo_list = ('photodiode','photodiode_alpha','photodiode alpha')
for photo_name in photo_list:
try:
photodiod = array(Data[photo_name], copy = False).T
except:
print photo_name, 'was not found'
continue
#correct abnormalilies
photo_tvec = photodiod[:,0]
photodiod = photodiod[:,1]
#photodiod[photodiod < 0] = 0
ind = (photo_tvec < plasma_end) & (photo_tvec > plasma_start)
photodiod[~ind] = 0
med = median(photodiod[ind])
photodiod[photodiod>3*med] = 3*med #opraví se tím "EM disrupce" na konci
SNR = 100
background = 0.005
photoTvec.append(photo_tvec)
photoData.append(photodiod)
photoError.append(hypot(photodiod/SNR,background))
photoLabel.append(photo_name)
photoSensitivity.append(loadtxt('bpy47.txt')) #Photodiod_BPY47 Leybold!
foto_filter = array([[200,1],[1200,1]]) #BUG temporal solution, photodiode properties are not loaded from database
if photo_name in ('photodiode_alpha', 'photodiode alpha'):
foto_filter = loadtxt('./filters/656.txt')
foto_filter[:,1] /= 100
photoFilter.append(foto_filter)
photoTvecError.append(1e-4) #[s]
correl_list.append(1)
for j,camera in enumerate(('camera\color_emiss_1','camera\color_emiss_2')): #BUG !!!
try:
tvec,data = Data[camera]
if size(data,-1) != 3:
continue
except:
continue
tvec/= 1000
tvec += plasma_start
camera_sensitivity = loadtxt('filters/camera_filters.txt') #TODO loadovat to ze serveru
for i,c in enumerate(('R','G','B')):
camera_gaps = isnan(data[:,i])
camera_filter = array(((350,1),(750,1))) #no filter
photoTvec.append(copy(tvec))
photoData.append(data[:,i])
SNR = 10
background = median(data[:,i])/10
photoError.append(hypot(data[:,i]/SNR,background))
photoLabel.append('camera %s%d'%(c,j))
photoSensitivity.append(camera_sensitivity[:,2*i:2*(i+1)])
photoFilter.append(camera_filter)
photoTvecError.append(1e-2) #[s]
correl_list.append(ones((3,3)))
PhotodiodesData = (photoData,photoError,correl_list,photoTvec,photoTvecError,photoFilter,photoSensitivity,photoLabel)
#LOAD spectrometer======================================
tvec,projection = load_adv('./data/projection')
projectionError = projection.data_err
tvecError = projection.tvec_err
SpectrSensitiv = spectrometer.relative_calibration
SpectrometerData = (tvec,tvecError,SpectrSensitiv,projection,projectionError)
return SpectrometerData,PhotodiodesData
def make_graphs():
names = ('HI','OI','OII', 'OIII', 'HeI', 'CII', 'CIII', 'NI','NII', 'NIII','Mystery1', 'CIV+NIV+OIV' )
styles = ('b-','r-','r--', 'r-.' , 'g-' , 'k--', 'k-.' , 'y-','y--', 'y-.' , 'c-' , 'm:' )
tvec,plasma = load_adv('./data/plasma')
if not any(plasma):
print 'no radiation observed'
return
tvec,intensities = load_adv('./data/intensities')
tvec,proj_dict = load_adv('./data/projection')
projection = proj_dict
projectionError = proj_dict.data_err
tvec,totalPower = load_adv('./data/TotalPower')
totalPower_err = totalPower.data_err
totalPower = totalPower
wavelength,components = load_adv('./data/components')
wavelength,overburned = load_adv('./data/overburned')
energy_constants = load('./data/energy_constants.npy')
chi2ion = loadtxt('./results/IonChi2.txt')
n_features = size(components,1)
n_components = size(components,0)
n_measurement = size(intensities,1)
chi2 = load('./results/TotalChi2.npy')
tvec,plasma = load_adv('./data/plasma') #BUG opravit to i jinde!!!
plasma = bool_(plasma)
class MyFormatter(plt.ScalarFormatter):
def __call__(self, x, pos=None):
if pos==0:
return ''
else: return plt.ScalarFormatter.__call__(self, x, pos)
majorLocator = plt.MultipleLocator(50)
majorFormatter = plt.FormatStrFormatter('%d')
minorLocator = plt.MultipleLocator(10)
overburned = all(overburned, 1) #dameged pixels and wavelengths shielded by glass
#fit = mean(dot(projection[plasma,:],components).T, axis=1)*(1-int_(overburned))[:,None]
#print overburned
fit = mean(dot(projection[plasma,:],components).T, axis=1)
mean_data = mean(intensities[:,plasma], axis=1)
resid = mean_data-fit
save_adv('./data/residual_spectra', wavelength,resid)
fig = plt.figure(figsize=(15, 3), dpi=80, facecolor='w', edgecolor='k')
ax = fig.add_subplot(111)
ax.plot(wavelength,resid-10,'c',label = 'residuum',lw = 0.3)
ax.plot(wavelength,mean_data,'r',label = 'data',lw = 0.3)
ax.plot(wavelength,fit,'b', label = 'retrofit',lw = 0.3)
ax.text(0.05,0.8,'$\\chi^2$ = %2.1f'%chi2,horizontalalignment='left',
verticalalignment='bottom',transform=ax.transAxes,backgroundcolor = 'w')
ax.xaxis.set_major_locator(majorLocator)
ax.xaxis.set_major_formatter(majorFormatter)
#for the minor ticks, use no labels; default NullFormatter
ax.xaxis.set_minor_locator(minorLocator)
ax.set_xlabel('$\lambda$ [nm]')
ax.set_ylabel('counts/read noise [-]')
leg = ax.legend(loc='upper right', fancybox=True)
leg.get_frame().set_alpha(0.7)
ax.axis('tight')
ax.set_xlim(220,880)
ax.set_ylim(-20, 20) #BUG dát tam třeba průměr spektra?
ax.set_title('Total retrofit')
fig.savefig('graphs/spectrum_retrofit.png',bbox_inches='tight')
fig.clf()
#plasma = bool_(plasma)
if not any(plasma):
plasma[:] = True
#plot line intensites normalized by their mean value in the database
fig = plt.figure( figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')
ax = fig.add_subplot(111)
for i in range(n_components):
ax.errorbar(tvec[plasma]*1e3,projection[plasma,i],fmt = styles[i],\
yerr=projectionError[plasma,i], label=names[i]) #intensity with errorbars
ax.set_title('ion radiation intensity')
leg = ax.legend(loc='upper left', fancybox=True)
leg.get_frame().set_alpha(0.7)
ax.set_ylabel('intenzity')
ax.set_xlabel('t [ms]')
ax.axis('tight')
ax.set_ylim(0,None)
fig.savefig('graphs/RelativeIntensity.png',bbox_inches='tight')
fig.clf()
ax = fig.add_subplot(111)
projection*= energy_constants
projectionError *= energy_constants
for i in range(n_components):
ax.errorbar(tvec[plasma]*1e3,projection[plasma,i]/1e3,
yerr=projectionError[plasma,i]/1e3,fmt=styles[i], label=names[i]) #intensity with errorbars
ax.set_title('radiation power in range 200-1200 nm')
leg = ax.legend(loc='upper left', fancybox=True)
leg.get_frame().set_alpha(0.7)
ax.set_ylabel('P [kW]')
ax.set_xlabel('t [ms]')
ax.axis('tight')
ax.set_ylim(0,None)
fig.savefig('graphs/RadiatedEnergy.png',bbox_inches='tight')
fig.clf()
#show()
tvec_temp,temp_dict = load_adv('./data/temperature')
tmin = amin(tvec_temp)*1e3
tmax = amax(tvec_temp)*1e3
fig = plt.figure(figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k')
ax = fig.add_subplot(111)
ax.errorbar(tvec*1e3,totalPower/1e3, yerr=totalPower_err/1e3)
ax.set_title('Total radiated power in range 200-1200 nm')
ax.set_ylabel('P [kW]')
ax.set_xlabel('t [ms]')
#plt.axis('tight')
ax.set_xlim(tmin-1, tmax+1)
ax.set_ylim(0,None)
fig.savefig('graphs/TotalPower.png',bbox_inches='tight')
fig.clf()
##MAIN PLOT
#data_max = amax(projection)
#colours = ('b','r','g','c')
#fig = figure(num=None, figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k')
#subplots_adjust(hspace=0, wspace = 0)
##Hydrogen
#ax = fig.add_subplot(6,1,1)
#title('radiation energy in range 200-1200 nm')
#ax.xaxis.set_major_formatter( NullFormatter() )
#ax.yaxis.set_major_formatter( MyFormatter() )
#i = 0
#chi2 = chi2ion[i]/n_measurement
#ax.text(.05,.1,'$\\chi^2$ = %2.2f'% chi2,horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes)
#errorbar(tvec*1e3,projection[:,i], yerr=projectionError[ :,i], label = names[i],fmt=styles[i]) #intensity with errorbars
#ylim(0, data_max)
#ylabel('I normalized')
#leg = legend(loc='upper right', fancybox=True)
#leg.get_frame().set_alpha(0.5)
##Helium
#ax = fig.add_subplot(6,1,2)
#i = 4
#ax.xaxis.set_major_formatter( NullFormatter() )
#ax.yaxis.set_major_formatter( MyFormatter() )
#chi2 = chi2ion[i]/n_measurement
#ax.text(.05,.1,'$\\chi^2$ = %2.2f'% chi2,horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes)
#errorbar(tvec*1e3,projection[:,i], yerr=projectionError[ :,i], label = names[i],fmt=styles[i]) #intensity with errorbars
#ylim(0, data_max)
#ylabel('I normalized')
#leg = legend(loc='upper right', fancybox=True)
#leg.get_frame().set_alpha(0.5)
##carbon
#ax = fig.add_subplot(6,1,3)
#index = [5,6]
#ax.xaxis.set_major_formatter( NullFormatter() )
#ax.yaxis.set_major_formatter( MyFormatter() )
#chi2 = chi2ion[i]/n_measurement
#ax.text(.05,.1,'$\\chi^2$ = %2.2f'% chi2,horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes)
#for i in index:
#errorbar(tvec*1e3,projection[:,i], yerr=projectionError[ :,i], label = names[i],fmt=styles[i]) #intensity with errorbars
#ylim(0, data_max)
#leg = legend(loc='upper right', fancybox=True)
#leg.get_frame().set_alpha(0.5)
#ylabel('I normalized')
##oxygen
#ax = fig.add_subplot(6,1,4)
#index = [1,2,3]
#ax.xaxis.set_major_formatter( NullFormatter() )
#ax.yaxis.set_major_formatter( MyFormatter() )
#chi2 = chi2ion[i]/n_measurement
#ax.text(.05,.1,'$\\chi^2$ = %2.2f'% chi2,horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes)
#for i in index:
#errorbar(tvec*1e3,projection[:,i], yerr=projectionError[ :,i], label = names[i],fmt=styles[i]) #intensity with errorbars
#ylim(0, data_max)
#leg = legend(loc='upper right', fancybox=True)
#leg.get_frame().set_alpha(0.5)
#ylabel('I normalized')
##Nitrogen
#ax = fig.add_subplot(6,1,5)
#ax.xaxis.set_major_formatter( NullFormatter() )
#ax.yaxis.set_major_formatter( MyFormatter() )
#index = [7,8,9]
##index = [7,8]
#chi2 = chi2ion[i]/n_measurement
#ax.text(.05,.1,'$\\chi^2$ = %2.2f'% chi2,horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes)
#for i in index:
#errorbar(tvec*1e3,projection[:,i], yerr=projectionError[ :,i], label = names[i],fmt=styles[i]) #intensity with errorbars
#ylim(0, data_max)
#leg = legend(loc='upper right', fancybox=True)
#leg.get_frame().set_alpha(0.5)
#ylabel('I normalized')
##Mystery
#ax = fig.add_subplot(6,1,6)
#index = [10,11]
##index = [9,10]
#chi2 = chi2ion[i]/n_measurement
#ax.text(.05,.1,'$\\chi^2$ = %2.2f'% chi2,horizontalalignment='left',verticalalignment='bottom',transform=ax.transAxes)
#for i in index:
#errorbar(tvec*1e3,projection[:,i], yerr=projectionError[ :,i], label = names[i],fmt=styles[i]) #intensity with errorbars
#ylim(0, data_max)
#leg = legend(loc='upper right', fancybox=True)
#leg.get_frame().set_alpha(0.5)
#ylabel('I normalized')
#xlabel('t [ms]')
#savefig('./graphs/MultiPlot.png', bbox_inches='tight')
##show()
#clf()
def PlotDetailResiduum(GOLEM):
print 'PlotDetailResiduum'
style_dict = {'H':'b-','OI':'r-','OII':'r--', 'OIII':'r-.','OIV':'r:', 'HeI':'g-','HeII':'g--',\
'CI':'k-','CII':'k--', 'CIII':'k-.','CIV':'k:', 'NI':'y-','NII':'y--', 'NIII':'y-.','NIV':'y:','unknow':'c-'}
names = ('HI','OI','OII', 'OIII', 'HeI', 'CII', 'CIII', 'NI','NII', 'NIII','Mystery1', 'CIV+NIV+OIV' )
styles = ('b-','r-','r--', 'r-.' , 'g-' , 'k--', 'k-.' , 'y-','y--', 'y-.' , 'c-' , 'm:' )
wavelength,components = load_adv('./data/components')
tvec,proj_dict = load_adv('./data/projection')
tvec,plasma = load_adv('./data/plasma') #BUG opravit to i jinde!!!
plasma = bool_(plasma)
projection = (proj_dict)[plasma,:]
projectionError = proj_dict.data_err[plasma,:]
projection = mean(projection, axis=0)
tvec,intensities = load_adv('./data/intensities')
intensities = mean(intensities[:,plasma],axis=1)
norm_components = projection[:,None]*components
fit = mean(dot(projection,components), axis = 0)
fit = sum(norm_components, axis = 0)
resid = intensities- fit
majorLocator = plt.MultipleLocator(10)
majorFormatter = plt.FormatStrFormatter('%d')
minorLocator = plt.MultipleLocator(1)
#fig = plt.figure(figsize=(10,5))
fig = plt.figure(figsize=(20,15))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
shift = 0.7 - (wavelength-206)/183*0.1 #by eye
ax.plot(wavelength-shift, resid, 'b:', linewidth=1,label='fit residuum')
ax.plot(wavelength-shift, intensities, 'k:', linewidth=1, label='data')
for shot_spectrum,n,s in zip(norm_components, names, styles):
ind = shot_spectrum!= 0
ax.plot(wavelength[ind]-shift[ind],shot_spectrum[ind],s,linewidth=1, label=n)
#plot lines from database
for i, element in enumerate(GOLEM.elements):
s = style_dict[element[1]]
ax.axvline(x = 0 , label = element[1], c = s[0], ls = s[1:])
for j,line in enumerate(element[0]):
if line< amin(wavelength) or line> amax(wavelength):
continue
ax.axvline(x=line , c=s[0], ls=s[1:])
ax.text(line , i+j, element[1],color=s[0])
ax.set_xlim(amin(wavelength),amax(wavelength))
ax.set_ylim(-1,None)
ax.set_yscale('symlog')
leg = ax.legend(loc='center right', fancybox=True)
leg.get_frame().set_alpha(0.7)
ax.xaxis.set_major_locator(majorLocator)
ax.xaxis.set_major_formatter(majorFormatter)
#for the minor ticks, use no labels; default NullFormatter
ax.xaxis.set_minor_locator(minorLocator)
print 'saving '
fig.savefig('./graphs/spectra.svgz',bbox_inches='tight')
plt.close()
def main():
#print 'start 1',sys.argv
for path in ['graphs', 'results','data','HistoricalAnalysis']:
if not os.path.exists(path):
os.mkdir(path)
#print 'start 2'
if sys.argv[1] == "analysis":
#print 'load'
spectrometr,tok, actualShot,plasma_start,plasma_end = LoadData()
print "loaded data"
noplasma = CalcProjections(spectrometr,actualShot,plasma_start,plasma_end)
print "CalcProjections done"
if noplasma:
print "no plasma"
return
CalcDensTemp()
#print 'start 3'
if sys.argv[1] == "plots":
make_graphs()
plotDensTemp()
#print 'plotDensTemp'
saveconst('status', 0)
os.system('convert -resize 150x120\! ./graphs/density.png icon.png')
#print 'start 4'
if sys.argv[1] == "postanalysis":
spectrometr,tok, actualShot,plasma_start,plasma_end = LoadData() #BUG loadování je šíleně pomalé? neuložit si to někd předtím??
wavelength, _ = actualShot.getData()
SpectrometerData,PhotodiodesData = LoadDataIDA(plasma_start,plasma_end,spectrometr)
IntegLightAnalyz(wavelength,plasma_start,plasma_end, SpectrometerData,PhotodiodesData)
#PlotDetailResiduum(tok)
##print 'start 5'
if __name__ == "__main__":
main()
#TODO 7307 - co je to tam za divné čáry kolem 600nm?
#7960
#TODO 9298 . 750nm . není to CI??
#TODO 9156 .350-400nm co to je
#TODO 9664,9665,9670,9839,9906,10131,10136 úplně noé čáry
#TODO 10097 590 to fitne mizerně
#TODO 9504 - molekulární spektrum
#problémy - heliové výboje
#neznámý prvek, podivné čáry
#NOTE nejošlklivější výboje všech dob 11527,11528,11529
#514.1 should by CI line!?
#CIV čára 580-581 nm - dodělat proč mi tam chybí??
#CV 494.5 slabá, je vidět??
#velmi teplé výboje - kolem 21545
