#!/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