Revision 4d3d7a16aa5557ba88de7ce1b78253ad68c7a4e7

TrainingCourses/FTTF/2015-2016/BorLeitl/GD_Spectrometry/spectrometerDataProcVer3.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
# coding: utf-8
import numpy as np
import re
#import matplotlib
# matplotlib.use( 'Agg' )
import matplotlib.pyplot as plt
import datetime as dt
from matplotlib.dates import date2num, DateFormatter
import os
import sys
# from os import remove
import urllib
import codecs
from time import sleep
import shutil

from scipy.optimize import curve_fit

class glowDischarge:
    
    def __init__(self):
        """
        """
#         self.GDType = "" #He or H GD possible on GOLEM
        self.dataFile = ""

    def loadOneGdSpectra(self,filePath, name = "spectra"):
        """
        load single spectra file, return data and time hh:mm:ss of spectra generation
        """

        filePath = filePath + name + ".txt"

        with open(filePath) as spectra:
            lines = spectra.readlines()
            time = (lines[2].split(" ")[-1]).split("\n")[0]
            time = np.array(time.split(":"),dtype=int)

        data = np.zeros(2048)

        for i, line in enumerate(lines[(2048+10):(2048*2+10)]): #TODO vylepsit pro vsechny data, s ruznym mnozstvim spekter
            data[i] = float(line)

        return data, time

    def loadAllGdSpectra(self,filePath, name = "spectra"):
        """
        load all spectra in one data file
        """
        files = sorted(np.array(os.listdir(path = filePath), dtype = int))
        files = np.array(files, dtype = str)
        files_num = len(files)
        data = np.zeros([files_num,2048])
        time = np.zeros([files_num,3])
        

        for i, data_file in enumerate(files):
            data[i,:], time[i,:] = self.loadOneGdSpectra(filePath + "/" + data_file + "/")
        return data, time

def Gauss_function(x, a, x0, sigma,c):
    return a * np.exp(-(x - x0)**2 / (2 * sigma**2)) +c

def Gauss_fit(x,y,baseMean):
    """
    fit given data with gauss function - standart least square method is implicitly used by curve_fit from scipy.optimize
    for test if the fit is appropriate, the use of student test or so is necessary 
    """
    #test if x,y have same dimension
    if len(x) != len(y):
        return print("Gauss_fit problem - wrong dimmensions of fitting data")
    
    mean = sum(x * y) / sum(y)
    sigma = np.sqrt(sum(y * (x - mean)**2) / sum(y))
    popt, pcov = curve_fit(Gauss_function, x, y, p0=[max(y), mean, sigma,baseMean])
    
    return popt, pcov
    
def findpeak(data, minima):
    """
    find all local maxima in data with given limit (minima)
    """
    wavelenght_file = "/home/bobon/SCHOOL/GOLEM/PRAKTIKA2/Skripty/"
    #READ WAVELENGTH FILE
    try: 
        channels = np.loadtxt("%swavelenghts.txt" %wavelenght_file)
    except: print("wavelengths file cannot be opened")
    
    datalen = len(data)
    local_max = np.zeros((2,datalen))
    fit_data = np.zeros((9,2,datalen))
    
    pos = 0
    
    for i in range(datalen):
        if (i!=1 and i!= (datalen-1)):
            if ((data[i] >= data[i-1]) and (data[i] >= data[i+1]) and (data[i] > minima)):
                local_max[1,pos] = data[i]
                local_max[0,pos] = channels[i] #uz rovnou vln delka kanalu odpovidajici otypovanemu maximu
#                 local_max[0,pos] = i #pouze cislo kanalu
                
                try:
#                     fit_data[:,0,pos] = np.arange(i-4,i+5) # x take 9 values around with central - prepocitat rovnou na vln delky?
                    fit_data[:,0,pos] = np.array([channels[chan] for chan in np.arange(i-4,i+5)]) # x take 9 values around with central
                    fit_data[:,1,pos] = data[i-4:i+5] # y intensities
#                     print(np.arange(i-2,i+3))
                except: pass
    
                pos += 1

    return local_max, fit_data

def fitApply(data, index):
    """
    """
    #try:
    baseMean = np.mean(data[index,:])
    maxima, fitin_data = findpeak(data[index,:],baseMean+50)
    l = range(len(data[index,:]))

    #fit
    fit_params = np.zeros((4,len(fitin_data[1,0,:])))
    for i in range(len(fitin_data[1,0,:])):
        if fitin_data[2,0,i] == 0:
            break
        x = fitin_data[:,0,i]
        y = fitin_data[:,1,i]
        xfit = np.linspace(fitin_data[0,0,i], fitin_data[-1,0,i],50)

        try:
            fit_params[:,i], co = Gauss_fit(x,y,baseMean)
        except: print("cannot find fit")
        
        
        eps = np.abs(fit_params[3,i] + fit_params[0,i] - np.max(y))/np.max(y)
        #print("eps = %f, wavelenght = %f, intensity = %f,i = %i" %(eps, fit_params[1,i], (fit_params[3,i] + fit_params[0,i]),i))
        if eps >0.5: #pokud maximum  fit[1][0,:] + fit[1][3,:]
            fit_params[3,i] = 0
            fit_params[0,i] = 0

    return fit_params

def fitAllData(data,time):
    """
    store fit data to DataFitParam as [time,fit_parameters]  fit_parameters - fit parameters for finded peaks
    """
    dataFitParam = []
    for i, t in enumerate(time):
        dataFitParam.append([t,fitApply(data,i)])
        
    return dataFitParam

def plotFittedData(time, Intensities, spectrum):
    """
    for individual ploting of fitted data
    """
    
    #PLOT TIMEBASE
    base = dt.datetime(2017, 5, 5, 0, 0, 0)
    time_plot = [base + dt.timedelta(seconds=((time[i,2]))) + dt.timedelta(minutes=(time[i,1])) + dt.timedelta(hours=(time[i,0])) for i,k in enumerate(time[:,1])]
    time_plot = date2num(time_plot)
    index_lastview = 50
#     plt.figure(1)
    fig, ax = plt.subplots(1,1)
    plt.axvline(ymin=.0, ymax = 1, x=time_plot[index_lastview], linewidth=2, color = 'k')
    
    index = np.arange(len(Intensities[0,:])) #index of actual time
    segment1 = (index <= index_lastview)
    segment2 = (index >= index_lastview)
    
    for line in range(len(Intensities[:,1])):
        if np.mean(Intensities[line,:]) > 100:
            plt.plot_date(time_plot[segment1],Intensities[line,segment1], label = '{:.1f}'.format(spectrum[line]), ms = "4")
#             plt.plot_date(time_plot[segment2],Intensities[line,segment2], label = '{:.1f}'.format(spectrum[line]), ms = "1")
            
            #                     #PLOT FORMAT
            fig.autofmt_xdate(rotation = 30)
            ax.xaxis.set_major_formatter( DateFormatter('%H:%M:%S') )
        #     plt.title("channel = %.2f" %channels[line] + 'wavelength = %.2f' %spektrum[i] + ", Intensity threshold = %i" %min_intensity)
            plt.xlabel("time")
            plt.ylabel("Intensity")

    plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
    
    plt.gca().set_prop_cycle(None) #restart of color palete needed
    #plt.hold(True)
    for line in range(len(Intensities[:,1])):
        if np.mean(Intensities[line,:]) > 100:
            plt.plot_date(time_plot[segment2],Intensities[line,segment2], label = '{:.1f}'.format(spectrum[line]), ms = "1")
            plt.hold(True)
    plt.show()
    plt.close()

def loadOneGdSpectra(file, name = "spectra"):
        """
        load single spectra file, return data and time hh:mm:ss of spectra generation
        """

        file = file + name + ".txt"

        with open(file) as spectra:
            lines = spectra.readlines()
            time = (lines[2].split(" ")[-1]).split("\n")[0]
            time = np.array(time.split(":"),dtype=int)
        
        data = np.zeros(2048)
        
        i=0
        for line in lines:
            if line.find("Wavelengths:"):
               shift = i+1
               print("shift = %d" %shift)
               break
            i+=1
            
            
        
        
        for i, line in enumerate(lines[(2048+shift):(2048*2+shift)]): #shift = 10lines for description 11 for tubes files
            data[i] = float(line)

        return data, time


def loadALLGdSpectra(file, name = "spectra"):
    """
    load all spectra in one data file
    """

    files = np.array(os.listdir(path = file),dtype=int)
    files = sorted(files)
    files_num = len(files)
    data = np.zeros([files_num,2048])
    time = np.zeros([files_num,3])
    #print(files)

    for i, data_file in enumerate(files):
        data[i,:], time[i,:] = loadOneGdSpectra(file + "/" + str(data_file) + "/")
    return data, time

#tento kod se spusti pro všechny data - možno brát i po jednom souboru - jen nutno vhodně zadat index času pro zařazení 
def zaradPeaky(peaks,intensities,peakList,intTable, timeIndex):
    """
    try to find wavelength in the list
    """
    
    wavTolerance = 0.3 #tolerance in nanometers
    
    for peakIndex, wavelength in enumerate(peaks):
        minima = 1000
        #if wavelength == 0: #za platnymi daty jsou nuly, zbytecne je brat - ne tak úplně pozor, někde může být nula i tak a pak nebere ostatní data!
            #break
        diff = np.abs(peakList - wavelength)
        try:
            minima = np.min(diff)
            
        except:
            minima = 1000
            pass
        if minima < wavTolerance:
            minIndex = np.argmin(diff)
            intTable[minIndex,timeIndex] = intensities[peakIndex]
        else:
            size = np.size(intTable[0,:])
            newIntVector = np.zeros((1,size))
            newIntVector[0,timeIndex] = intensities[peakIndex]
            newWave = np.array([wavelength])
            peakList = np.concatenate([peakList,newWave])
            intTable = np.concatenate([intTable,newIntVector])
            
    #print(np.sort(peakList))
    return peakList, intTable

def loadElementsData(elements,listPath): #toto by mohlo být venku
        #LOAD ELEMENT SPECTRUM FILE
        tableSpectrum = np.array((len(elements),100))
        for element in enumerate(elements):
            spectrum = np.loadtxt("%s%s.txt" %(listPath, element))
            for j in range(len(spectrum)):
                tableSpectrum[i,j] = spectrum[j]
        return spectrum
    
#V
def identifyElements(wavelenghts,listPath, elements, lineResolution):
    """
    return vector of names for indentified peaks in wavelenghts vector
    
    wavelenghts - vector of wavelenghts for indentification
    listPath - path to elements peaks directory
    elements - vector of names of elements used for identification
    lineResolution
    """
    
    def loadElementsData(elements,listPath):
        """
        LOAD ELEMENT SPECTRUM FILE
        """
        tableSpectrum = np.zeros((len(elements),500))
        for i,element in enumerate(elements):
            spectrum = np.loadtxt("%s%s.txt" %(listPath, element))
            for j in range(len(spectrum)):
                tableSpectrum[i,j] = spectrum[j]
        return tableSpectrum
    
    tableSpectrum = loadElementsData(elements,listPath)
    
    def findEle(wavPeak, lineResolution, tableSpectrum):
        wavDiff = 100
        wavMin = 1000
        corespElement = ""
        
        for i, element in enumerate(elements):
            for eleWav in tableSpectrum[i,:]:
                if eleWav == 0:
                    break
                
                wavDiff = abs(eleWav - wavPeak)
                
                if wavDiff < wavMin and wavDiff < lineResolution:
                #if wavDiff < lineResolution:
                    print(lineResolution)
                    wavMin = wavDiff
                    corespElement = element #string of correspondent element
        if corespElement == "":
            corespElement = "N/I"
        
        return corespElement
    
    identPeaks = []
    
    for wavPeak in wavelenghts:
        identPeaks.append(findEle(wavPeak, lineResolution, tableSpectrum))
    return identPeaks

def plotOneSpectrum(data,fplots,intTable,argsort,identElements, peakList, chanFile,dateMeasurement, numLines = 5):
    """
    plot spectra data file for given element
    data = intensities table
    modified version, plot from intensity table, given wavelenght vector and list of identified elements
    """
     #number of first spectra according to signal intensity the same number is ploted
    #READ channel WAVELENGTH FILE
    try: 
        channels = np.loadtxt("%swavelenghts.txt" %chanFile)
    except IOError as e:
        print("Unable to read wavelenghts file")

    #rescale Intensities to absoluted values
    maxSignal = np.max(intTable)
    
    sizePeakList = np.size(peakList)
    if numLines > sizePeakList:
        numLines = np.size(peakList)
    #intTable = intTable/maxSignal
    

    ##PLOT TIMEBASE
    #base = dt.datetime(2017, 5, 9, 0, 0, 0) #TODO ziskat datum a cas primo ze souboru s daty
    #time_plot = [base + dt.timedelta(seconds=((time[i,2]))) + dt.timedelta(minutes=(time[i,1])) + dt.timedelta(hours=(time[i,0])) for i,k in enumerate(time[:,1])]
    #time_plot = date2num(time_plot)
    #timeLen = len(time)

    #SET PLOT RANGE
    #index = np.arange(timeLen) #index of actual time
    #segment1 = (index <= index_lastview)
    #segment2 = (index >= index_lastview)

    #PLOT INIT
    legend = [identElements[i] + " - %.1f nm" %(peakList[i]) for i in argsort[0:numLines]]
    fig, ax2 = plt.subplots(1,1, sharex = False, sharey = False)
    #fig.subplots_adjust(hspace=0.60)

    #for label in ax1.get_xmajorticklabels():
        #label.set_rotation(30)
        #label.set_horizontalalignment("right")

    #ax1.set_xlabel("time")
    #ax1.set_ylabel("Intensity [a.u.]")
    #ax1.set_ylabel("Intensity")

    #ax1.xaxis.set_major_formatter( DateFormatter('%H:%M') )
    #ax1.set_title("Main spectral lines for %s GD on %s" %(GDGas, dateMeasurement)) #TODO - přidat plyn, ve kterém výboj probíhá - bude asi na volbě zvenčí při spouštění programu - mohlo
    #by se navázat na spouštění samotného výboje, aby se zároveň zjistily charakteristiky jako tlak a teplota - mohl by být třetí graf - 
    #tam by byl - tlak/teplota/ v závislosti na čase - tlak nutno přepočítat na použitý plyn
    #TODO - připojit datum měření

    ax2.set_xlabel("wavelenght [nm]")
    ax2.set_ylabel("intensity [a.u.]")
    ax2.set_ylabel("Intensity")
    #ax2.set_title("Spectrum at %02i:%02i:%02i" %(time[index_lastview,0],time[index_lastview,1],time[index_lastview,2])) #TODO doplnit čas
    #ax2.set_yticks(ticks = [0,0.2,0.4,0.6,0.8,1])
    #ax2.set_ylim(0.,1.) #BOR

    #PLOT - MAIN SPECTRAL LINES
    #for i in range(numLines):
        #ax1.plot_date(time_plot[segment1],intTable[argsort[i],segment1],"+", ms = "4")
    
    #print(len(time_plot[segment1]))
    #print(len(intTable[argsort[i],segment1]))
    
    #ax1.set_prop_cycle(None) #restart of color cycle needed
    #for i in range(numLines):
#             plt.plot_date(time_plot,intTable[argsort[i],0:len(time)])
        #ax1.plot_date(time_plot[segment2],intTable[argsort[i],segment2],"o", ms = "1")
        #plt.hold(True)
    
    #PLOT2 - SPECTRUM
    spectraA = data[0,:]
    #spectraA = data[index_lastview,:]
    ax2.plot(channels,spectraA,"-")
    #plt.plot(channels,spectraA,"-")

    displayAnotateNum = 0
    ax2Legend = []

    print(peakList)
    for i in range(15):
        if i== len(peakList)-1:
            break
        print(i)
        #ax2.axvline(data[i,0], color='k', linestyle='--')
        lineIntensity = intTable[argsort[i],0]
        lineWavelength = peakList[argsort[i]]
        
        if lineIntensity > 0.:
            print(lineIntensity)
            text = "%.1f" %lineWavelength
            
            ax2.annotate(text,xy=(lineWavelength, lineIntensity-0.05), 
                        xytext=(lineWavelength+(+80)/10, lineIntensity), size = 5, rotation = 30,
                        arrowprops=dict(facecolor='black', arrowstyle="->"))
            ax2Legend.append("%s" %identElements[argsort[i]] + "-" + text + " nm")
            displayAnotateNum +=1
            print(text)
        
    if displayAnotateNum >= 15:
        displayAnotateNum = 15

    #ax1 Legend
    #box = ax1.get_position()
    #ax1.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    #ax1.legend(legend[0:numLines], loc=2, borderaxespad=0.,bbox_to_anchor=(1.01, 1),title = "Main lines")
    #ax1.axvline(ymin=.0, ymax = 1, x=time_plot[index_lastview], linewidth=2, color = 'k')

    #ax2 Legend
    box = ax2.get_position()
    ax2LegendText = "".join(["%s \n" %text for text in ax2Legend[0:displayAnotateNum]])
    ax2.set_position([box.x0, box.y0, box.width*0.8, box.height])
    leg = ax2.legend((ax2LegendText[:],), loc=2, borderaxespad=0.,bbox_to_anchor=(1.01, 1), title = "Main lines", markerscale = 5)
    for item in leg.legendHandles:
        item.set_visible(False)

    ##SAVE FIGURE
    fplots = fplots + "/" + str(1)
    plt.savefig(fplots, dpi = 150)
    plt.close("all")
    
#VI
def plotAlongLines(data,fplots,intTable,argsort,identElements, peakList, time, chanFile, index_lastview, dateMeasurement,GDGas, numLines = 5):
    """
    plot spectra data file for given element
    data = intensities table
    modified version, plot from intensity table, given wavelenght vector and list of identified elements
    """
     #number of first spectra according to signal intensity the same number is ploted
    #READ channel WAVELENGTH FILE
    try: 
        channels = np.loadtxt("%swavelenghts.txt" %chanFile)
    except IOError as e:
        print("Unable to read wavelenghts file")

    #rescale Intensities to absoluted values
    maxSignal = np.max(intTable)
    
    sizePeakList = np.size(peakList)
    if numLines > sizePeakList:
        numLines = np.size(peakList)
    #intTable = intTable/maxSignal
    

    #PLOT TIMEBASE
    base = dt.datetime(2017, 5, 9, 0, 0, 0) #TODO ziskat datum a cas primo ze souboru s daty
    time_plot = [base + dt.timedelta(seconds=((time[i,2]))) + dt.timedelta(minutes=(time[i,1])) + dt.timedelta(hours=(time[i,0])) for i,k in enumerate(time[:,1])]
    time_plot = date2num(time_plot)
    timeLen = len(time)

    #SET PLOT RANGE
    index = np.arange(timeLen) #index of actual time
    segment1 = (index <= index_lastview)
    segment2 = (index >= index_lastview)

    #PLOT INIT
    legend = [identElements[i] + " - %.1f nm" %(peakList[i]) for i in argsort[0:numLines]]
    fig, (ax1,ax2) = plt.subplots(2,1, sharex = False, sharey = False)
    fig.subplots_adjust(hspace=0.60)

    for label in ax1.get_xmajorticklabels():
        label.set_rotation(30)
        label.set_horizontalalignment("right")

    ax1.set_xlabel("time")
    #ax1.set_ylabel("Intensity [a.u.]")
    ax1.set_ylabel("Intensity")

    ax1.xaxis.set_major_formatter( DateFormatter('%H:%M') )
    ax1.set_title("Main spectral lines for %s GD on %s" %(GDGas, dateMeasurement)) #TODO - přidat plyn, ve kterém výboj probíhá - bude asi na volbě zvenčí při spouštění programu - mohlo
    #by se navázat na spouštění samotného výboje, aby se zároveň zjistily charakteristiky jako tlak a teplota - mohl by být třetí graf - 
    #tam by byl - tlak/teplota/ v závislosti na čase - tlak nutno přepočítat na použitý plyn
    #TODO - připojit datum měření

    ax2.set_xlabel("wavelenght [nm]")
    #ax2.set_ylabel("intensity [a.u.]")
    ax2.set_ylabel("Intensity")
    ax2.set_title("Spectrum at %02i:%02i:%02i" %(time[index_lastview,0],time[index_lastview,1],time[index_lastview,2])) #TODO doplnit čas
    #ax2.set_yticks(ticks = [0,0.2,0.4,0.6,0.8,1])
   # ax2.set_ylim(0.,1.) #BOR

    #PLOT - MAIN SPECTRAL LINES
    for i in range(numLines):
        ax1.plot_date(time_plot[segment1],intTable[argsort[i],segment1],"+", ms = "4")
    
    #print(len(time_plot[segment1]))
    #print(len(intTable[argsort[i],segment1]))
    
    ax1.set_prop_cycle(None) #restart of color cycle needed
    for i in range(numLines):
#             plt.plot_date(time_plot,intTable[argsort[i],0:len(time)])
        ax1.plot_date(time_plot[segment2],intTable[argsort[i],segment2],"o", ms = "1")
        #plt.hold(True)
    
    #PLOT2 - SPECTRUM
    spectraA = data[index_lastview,:]
    #spectraA = data[index_lastview,:]
    #ax2.plot(channels,spectraA/maxSignal,"-")
    ax2.plot(channels,spectraA,"-")

    displayAnotateNum = 0
    ax2Legend = []
    for i in range(15):
#     ax.axvline(data[i,0], color='k', linestyle='--')
        lineIntensity = intTable[argsort[i],index_lastview]
        if lineIntensity > 0.:
            text = "%.1f" %peakList[argsort[i]]
            ax2.annotate(text,xy=(peakList[argsort[i]], lineIntensity-0.05), 
                        xytext=(peakList[argsort[i]]+(+80+20*(i%3))/10, lineIntensity+0.1*(i%2)+0.2*(i%3)), size = 5, rotation = 0,
                        arrowprops=dict(facecolor='black', arrowstyle="->"))
            ax2Legend.append("%s" %identElements[argsort[i]] + "-" + text + " nm")
            displayAnotateNum +=1
    if displayAnotateNum >= 15:
        displayAnotateNum = 15

    #ax1 Legend
    box = ax1.get_position()
    ax1.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    ax1.legend(legend[0:numLines], loc=2, borderaxespad=0.,bbox_to_anchor=(1.01, 1),title = "Main lines")
    ax1.axvline(ymin=.0, ymax = 1, x=time_plot[index_lastview], linewidth=2, color = 'k')

    #ax2 Legend
    box = ax2.get_position()
    ax2LegendText = "".join(["%s \n" %text for text in ax2Legend[0:displayAnotateNum]])
    ax2.set_position([box.x0, box.y0, box.width*0.8, box.height])
    leg = ax2.legend((ax2LegendText[:],), loc=2, borderaxespad=0.,bbox_to_anchor=(1.01, 1), title = "Main lines", markerscale = 5)
    for item in leg.legendHandles:
        item.set_visible(False)

    #SAVE FIGURE
    fplots = fplots + "/" + str(index_lastview)
    plt.savefig(fplots, dpi = 150)
    plt.close("all")

def getMeasureDate(filePath):
    fileD = "/1/spectra.txt"
    filePath = filePath + fileD
    with open(filePath) as spectra:
            lines = spectra.readlines()[0:10]
            date = lines[2].split("\t")[1].split(" ")[0]
            date = date.split(".")
            measureDate = date[2] + date[1] + date[0][2::]
            return measureDate

#======================================================================================================================================

def makeSpectra():
    #elements = ["H","HeI", "HeII", "CI", "CII", "OI", "OII", "NI"]
    
    GD = glowDischarge()
    try:
        dataPath = sys.argv[1]
    except:
        print("data path missing")
        
    
    #GD.dataFile = dataPath + "mereni_" + measureDate
    GD.dataFile = dataPath + "data"
    #GD.dataFile = dataPath + "mereni_" + measureDate + "_JVPortH2" #4,8 a 10 cm od horního kraje portu - první byl přesaturovaný
    chanFile = dataPath
    SpectralLines_file = dataPath + "SpectralLines/"
    measureDate = getMeasureDate(GD.dataFile) #TODO lepší z konfiguračního souboru nebo jako druhy vstup do programu nebo z prvniho otevreneho souboru - nejlepší varianta
    fplots = dataPath + "PLOTS"
    GDGas = sys.argv[2]
    
    try:
        shutil.rmtree(fplots)
    except:
        pass
    try:
        os.mkdir(fplots)
    except:
        print("can not to create PLOTS file")
    
    data, time = GD.loadAllGdSpectra(GD.dataFile)
    medianData = np.median(data,axis = 1)
    for i in range(len(data[0,:])):
        data[:,i] = data[:,i] - medianData
    
    #II - data fit
    fit = fitAllData(data,time)

    #III - initialization of intensities and peaklist - each column of intTable is for one wavelength
    peakList = np.zeros(0)
    intTable = np.zeros([1,200])
    
    for t,tim in enumerate(time):
        intensities = fit[t][1][0,:] + fit[t][1][3,:] #intensities from fitted data
        peaks = fit[t][1][1,:] #wavelenghts
        peakList, intTable = zaradPeaky(peaks,intensities,peakList,intTable,t)
        
    identPeaks = identifyElements(peakList,SpectralLines_file, elements, lineResolution = 0.3)

    intSum = np.sum(intTable, axis=1)
    sortedInt = np.sort(np.sum(intTable, axis=0))[::-1]
    argsort = np.argsort(intSum, axis=-1)[::-1]
    for i in range(len(time)):
        plotAlongLines(data, fplots, intTable,argsort,identPeaks, peakList, time,chanFile, i, measureDate, GDGas, numLines = 5)
    
def makeOneSpectra():
    """
    plots only one spectrum
    can be used for gas-filled tubes
    """
    #elements = ["H","HeI", "HeII", "CI", "CII","CIII","CIV","CV", "OI", "OII", "NI", "NII","OIV","OV", "OIII", "FeI", "NEI", "NEII"]
    elements = ["NEI", "NEII"]
    GD = glowDischarge()
    try:
        dataPath = sys.argv[1]
    except:
        print("data path missing")
        
    
    #GD.dataFile = dataPath + "mereni_" + measureDate
    GD.dataFile = dataPath + "data"
    #GD.dataFile = dataPath + "mereni_" + measureDate + "_JVPortH2" #4,8 a 10 cm od horního kraje portu - první byl přesaturovaný
    chanFile = dataPath
    SpectralLines_file = dataPath + "SpectralLines/"
    measureDate = getMeasureDate(GD.dataFile) #TODO lepší z konfiguračního souboru nebo jako druhy vstup do programu nebo z prvniho otevreneho souboru - nejlepší varianta
    fplots = dataPath + "PLOTS"
    GDGas = sys.argv[2]
    
    try:
        shutil.rmtree(fplots)
    except:
        pass
    try:
        os.mkdir(fplots)
    except:
        print("can not to create PLOTS file")
    
    data, time = GD.loadAllGdSpectra(GD.dataFile)
    #shift = 7.7 #3.9nm all spectra
    
    medianData = np.median(data,axis = 1)
    for i in range(len(data[0,:])):
        data[:,i] = data[:,i] - medianData
    
    #II - data fit
    fit = fitAllData(data,time)

    #III - initialization of intensities and peaklist - each column of intTable is for one wavelength
    peakList = np.zeros(0)
    intTable = np.zeros([1,200])
    
    intensities = fit[0][1][0,:] + fit[0][1][3,:] #intensities from fitted data
    peaks = fit[0][1][1,:] #wavelenghts
    peakList, intTable = zaradPeaky(peaks,intensities,peakList,intTable,0)
    #peakList = peakList + shift    
    identPeaks = identifyElements(peakList,SpectralLines_file, elements, lineResolution = 1.0)

    intSum = np.sum(intTable, axis=1)
    intTable = intTable[1::]
    #sortedInt = np.sort(np.sum(intTable, axis=0))[::-1]
    argsort = np.argsort(intSum[1::], axis=-1)[::-1]
    #print(len(argsort))
    #print(argsort[0:15])
    print(intTable[:,0])
    print(peakList) #zdá se, že peaklist neodpovídá intenzitám po použití plot spektrum
    #for i in range(15):
    print(argsort)
        #print(intTable[0][argsort[20]])
    plotOneSpectrum(data, fplots, intTable,argsort,identPeaks, peakList, chanFile, measureDate, numLines = 5)
    
    
#makeSpectra()
makeOneSpectra()

#TODO dořeš třídy a přidej novou věc - program, který bude třídy využívat a připlotí 
#TODO funkce na ukládání důležitých informací jako data, time .... potřebuji data tabulku nebo stačí jen tabulka intenzit? 
#přidáno malé kritérium na fit - lepší by byl studentův test
#9:56 - jak řešit problém s nezobrazováním některých celkem dobře určených peaků jako Halfa 656nm? - určeny jsou dobře v intensities - chyba v zařazení???
#11:00 chyba nalezena - problém byl ve hledání peaků - od wavelength = 0 už se nepokračovalo, ovšem to mohlo nastat i za jiných okolností zaradPeaky() line 331