In [1]:
from __future__ import print_function
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
import sys
from pylab import *
from scipy.signal import argrelextrema
import scipy.optimize as spo
from scipy.optimize import fsolve
from numpy import nan

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

def is_interactive(): # are we in jupyter ??
    import __main__ as main
    return 0 


if is_interactive():
    %matplotlib inline

Probe_position='The BPP top from the plasma centre $d=50$ mm'
Base='http://golem.fjfi.cvut.cz/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/Sessions/090217Bscan_HS'
BaseURL = "http://golem.fjfi.cvut.cz/utils/data/" #global

Datum=240217
SmoothCoefficient=100
HWNapetovyDelic=100 
HWProudovyResistor=20 # Ohm
HFSweepFrequency=1000 #Hz
CurrentScale=1000 # so .. mA
VoltageShift=20.0

FigSize=10,8
def FigDef():plt.figure(figsize=(10, 8), dpi= 80, facecolor='w', edgecolor='k')

VoltageTag='bpp_lp_voltage'
CurrentTag='bpp_lp_current'
PotentialTag='bpp_voltage'




 
In [2]:
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
    

# integration of signal with the rid of offset
def integrate(array, dt = 1e-6, aver_line = 4500, coeff = 1):
    aver = np.sum(array[0:aver_line])/aver_line
    return np.cumsum(array - aver) * dt * coeff 


def GetData(ShotNo,VacuumShot):
    mkdir(str(ShotNo))
    mkdir(str(ShotNo)+'/Basics')
    ShotNoStr=str(ShotNo)
    os.system('wget http://golem.fjfi.cvut.cz/shots/' + ShotNoStr + '/basicdiagn/Btoroidal.npz -O '+ShotNoStr+'/Basics/Btoroidal.npz');
    os.system('wget http://golem.fjfi.cvut.cz/shots/' + ShotNoStr + '/basicdiagn/Iplasma.npz -O '+ShotNoStr+'/Basics/Iplasma.npz');
    os.system('wget http://golem.fjfi.cvut.cz/shots/' + ShotNoStr + '/basicdiagn/Uloop.npz -O '+ShotNoStr+'/Basics/Uloop.npz');
    os.system('wget http://golem.fjfi.cvut.cz/shots/' + ShotNoStr + '/basicdiagn/graphpres.png -O '+ShotNoStr+'/Basics/graphres.png');
    #Plasma parameters
    PlasmaStart=int(float(np.loadtxt(urlopen(BaseURL+ShotNoStr+'/plasma_start')))*1e6) # in us
    #PlasmaStart=10000
    PlasmaEnd=int(float(np.loadtxt(urlopen(BaseURL+ShotNoStr+'/plasma_end')))*1e6) # in us
    #PlasmaEnd=25000
    np.savetxt(ShotNoStr+'/Basics/PlasmaStart',[PlasmaStart], delimiter=" ", fmt="%s")
    np.savetxt(ShotNoStr+'/Basics/PlasmaEnd',[PlasmaEnd], delimiter=" ", fmt="%s")

    os.system('wget '+BaseURL + ShotNoStr + '/' + CurrentTag + ' -O '+ShotNoStr+'/Basics/VAchar_current.txt');
    os.system('wget '+BaseURL + ShotNoStr + '/' + VoltageTag + ' -O '+ShotNoStr+'/Basics/VAchar_voltage.txt');
    os.system('wget '+BaseURL + ShotNoStr + '/' + PotentialTag + ' -O '+ShotNoStr+'/Basics/Potential.txt');

    os.system('wget '+BaseURL + ShotNoStr + '/mirnov_5 -O '+ShotNoStr+'/Basics/mc5_p.txt');
    os.system('wget '+BaseURL + str(VacuumShot) + '/mirnov_5 -O '+ShotNoStr+'/Basics/mc5_v.txt');
    os.system('wget '+BaseURL + ShotNoStr + '/mirnov_13 -O '+ShotNoStr+'/Basics/mc13_p.txt');
    os.system('wget '+BaseURL + str(VacuumShot) + '/mirnov_13 -O '+ShotNoStr+'/Basics/mc13_v.txt');
    
def DischargeBasics(ShotNo,VacuumShot,PlasmaStart, PlasmaEnd, Time, Uloop, Btoroidal, Iplasma, \
                    VAchar_voltage, VAchar_current,Potential, Potential_smoothed, dz):
    
    f = plt.figure(figsize=(20.0, 5.0))
    f,ax = plt.subplots(7,sharex=True);plt.subplots_adjust(hspace=0.001)
    f.set_size_inches(FigSize)
    ax[0].set_title('#' + str(ShotNo))
    ax[0].plot(Time/1000,Uloop['data'][PlasmaStart:PlasmaEnd],label='Loop voltage');ax[0].set_ylabel('$U_l$ [V]');ax[0].set_ylim(0,);ax[0].legend(loc=0)
    ax[1].plot(Time/1000,Btoroidal['data'][PlasmaStart:PlasmaEnd],label='Toroidal magnetic field');ax[1].set_ylabel('$B_t$ [T]');ax[1].set_ylim(0,);ax[1].legend(loc=0)
    ax[2].plot(Time/1000,Iplasma['data'][PlasmaStart:PlasmaEnd]/1000,label='Plasma current');ax[2].set_ylabel('$I_p$ [kA]');ax[2].set_ylim(0,);ax[2].legend(loc=0)
    ax[3].plot(Time/1000,VAchar_voltage[0:,1]*HWNapetovyDelic,label='VA char: U');ax[3].set_ylabel('$U$ [V]');ax[3].legend(loc=0)
    ax[4].plot(Time/1000,VAchar_current[0:,1]/HWProudovyResistor*CurrentScale,label='VA char: I');ax[4].set_ylabel('$I$ [mA]');ax[4].legend(loc=0)
    ax[4].plot(Time/1000,VAchar_current[0:,1]*0)
    ax[5].plot(Time/1000,Potential[0:,1]*HWNapetovyDelic,label='Floating potential');ax[5].legend(loc=0)
    ax[5].plot(Time/1000,Potential_smoothed*0)
    ax[5].plot(Time/1000,Potential_smoothed,label='Floating potential smoothed');ax[5].set_ylabel('$<U_{fl}>$ [V]');ax[5].legend(loc=1)
    ax[6].plot(Time/1000,dz[PlasmaStart:PlasmaEnd],label='Plasma position');ax[6].set_ylabel('PlPos [mm]');ax[6].legend(loc=0)
    ax[6].set_xlabel('$t$ [ms]');#ax[6].set_ylim(-200,200)
    ax[6].set_xlim(PlasmaStart/1000,PlasmaEnd/1000)
    plt.savefig(str(ShotNo)+'/Basics/BasicDiagWithBPP.jpg', bbox_inches='tight')
    if is_interactive():plt.show();
    plt.close();
    
    if is_interactive():plt.figure(figsize=(10, 8), dpi= 80, facecolor='w', edgecolor='k')
    plt.plot(VAchar_voltage[:,1],VAchar_current[:,1],'.')
    plt.xlabel('VA char - voltage U [V]');plt.ylabel('VA char - current U [V]')
    plt.title('VA char raw XY data: current vs. voltage')
    plt.savefig(str(ShotNo)+'/Basics/RawData.jpg', bbox_inches='tight')
    if is_interactive():plt.show();
    plt.close();
    
def FindingSweepIntervals(ShotNo,PlasmaStart, PlasmaEnd, Time, VAchar_voltage):
    #ack http://stackoverflow.com/questions/4624970/finding-local-maxima-minima-with-numpy-in-a-1d-numpy-array
       
    data=running_mean(VAchar_voltage[:,1]*HWNapetovyDelic,SmoothCoefficient)
    
    Extrems_tmp = (np.diff(np.sign(np.diff(data))).nonzero()[0] + 1) # local min+max

    # 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])>20 and abs(data[Extrems_tmp[i]])>60 : 
            Extrems.append(Extrems_tmp[i])

    np.savetxt(str(ShotNo)+'/Data/Extrems',Extrems, fmt="%s")

    # graphical output...
    
    plt.figure(figsize=(FigSize), dpi= 80, facecolor='w', edgecolor='k')
    plt.plot(Time/1000,data)
    plt.plot(Time/1000,VAchar_voltage[:,1]*HWNapetovyDelic)
    plt.plot(Time[Extrems]/1000, data[Extrems], "o", label="Extrems")
    plt.ylabel('VAchar_voltage U [V]');plt.xlabel('time [us]')
    plt.title('Maxima and minima identification')
    plt.legend()
    plt.savefig(str(ShotNo)+'/Basics/MaximMinims.jpg', bbox_inches='tight')
    if is_interactive():plt.show();
    plt.close();
    return Extrems

    

    
def Fits(IndvSetup,PlasmaParameters,Extrems,Time,VAchar_voltage,VAchar_current,VAchar_voltage_smoothed,VAchar_current_smoothed, Potential_smoothed,Uloop, Btoroidal, Iplasma, dz, IndxDef):
    ShotNo=IndvSetup[0]
    ShotNoStr=str(ShotNo)
    PlasmaStart=PlasmaParameters[0]
    PlasmaEnd=PlasmaParameters[1]
    
    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)
    
    print("Fits .. ...", end=":")
    print(len(Extrems)-1)
    AllVAchars=[]
    def AzoozFn(Vbias, a1, a2, a3, a4):
        return np.exp(a1 * np.tanh((Vbias+a2)/a3)) - a4

    def ClassicalFn(Vbias, Isat, Vfl, Te,):
        return -Isat*(1-np.exp((Vbias-Vfl)/Te))

    def Classical4Fn(Vbias, Isat, Vfl, Te, C):
        return -Isat*(1-C*(Vbias-Vfl))*(1-np.exp((Vbias-Vfl)/Te))

    
    Results=[]
    #Voltage = np.linspace(int(xlim_min),int(xlim_max),1000);dVoltage = Voltage[1]-Voltage[0]
    #AzoozInitialGuess = [3.8,34,42,2.9]
    AzoozInitialGuess = IndvSetup[3]
    #AzoozInitialGuess = [1.5,8,24,1.15] # Absolutne nejslabsi misto
    AFP0=0;AFP1=0;AFP2=0;AFP3=0
    #ClassicalInitialGuess = [4.6,-30,21] # Absolutne nejslabsi misto
    ClassicalInitialGuess = IndvSetup[4]
    ClassicalFitParms=[0,0,0,0]
    #ClassicalInitialGuess = [0.3,-20,20] # Absolutne nejslabsi misto
    CFP0=0;CFP1=0;CFP2=0;CFP3=0;Pocet=0
    if IndxDef[1]==0:
        IndxDef[1]=len(Extrems)-1 # test purposes
    
    for i in range(IndxDef[0],IndxDef[1]): 
        data=[VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],\
              VAchar_current_smoothed[Extrems[i]:Extrems[i+1]]]
        #data=[VAchar_voltage[Extrems[i]:Extrems[i+1]],\
        #      VAchar_current[Extrems[i]:Extrems[i+1]]]
        #Voltage = np.linspace(int(xlim_min),int(xlim_max),len(data[0]));dVoltage = Voltage[1]-Voltage[0]
        Voltage = np.linspace(int(np.min(data[0])),int(np.max(data[0])),len(data[0]));dVoltage = Voltage[1]-Voltage[0]
        #print(Voltage)
        #Voltage = data[0];dVoltage = Voltage[1]-Voltage[0]
        #print(Voltage)
        def plotsetup():
            plt.legend()
            ax1.plot(data[0],data[1],',',label="rawdata")
            ax1.plot(VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],VAchar_current_smoothed[Extrems[i]:Extrems[i+1]],label="data averaged")
            plt.grid(True)
            ax1.set_xlim(xlim_min,xlim_max);
            ax1.set_ylim(ylim_min,ylim_max*1.25); # Kvuli legende ..
            ax1.set_title('VAchar #'+ShotNoStr+' \
                       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
        fig, ax1 = plt.subplots();
        plotsetup()
        ax1.legend()
        plt.savefig(ShotNoStr+'/IndivVAchars/VA'\
                    +str(i).zfill(2) +':'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'.jpg',\
                    bbox_inches='tight')
        TimeLineDir=ShotNoStr+'/TimeLine/'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])
        mkdir(TimeLineDir)
        plt.savefig(TimeLineDir+'/VAcharSmoothed.jpg',bbox_inches='tight')
        plt.close();
        # Azooz fit
        #if is_interactive():
        print(i, end=",")
        #else:
        #    print(i)
        #print data[0][0]
        AzoozFitParms, po_cov = spo.curve_fit(AzoozFn, data[0], data[1], AzoozInitialGuess)
        #print(AzoozFitParms, po_cov)
        AzoozFit=AzoozFn(Voltage,AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3])
        #AzoozFit=AzoozFn(data[0],AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3])
        dAzoozFit_dVoltage = np.gradient(AzoozFit, dVoltage)
        InflPointLoc=argrelextrema(dAzoozFit_dVoltage, np.greater)[0]
        #VflPointLoc=fsolve(AzoozFit, 0.0)
        #VflPointLoc=fsolve(AzoozFn(Voltage,AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3]), 0.0)
        #FitLimit=VflPointLoc
        #print(len(InflPointLoc))
        if len(InflPointLoc) != 1 \
        or (ShotNo==23191 and i==7)\
        or (ShotNo==23198 and i>22):\
        # not reasonable data
            print('(Azooz fit problem)',end="")
            Results.append([\
            i,\
            (PlasmaStart+Extrems[i])/1000.0,\
            (PlasmaStart+(Extrems[i]+Extrems[i+1])/2)/1000.0,\
            (PlasmaStart+Extrems[i+1])/1000.0,\
            nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan
            ])
            np.savetxt(TimeLineDir+'/index.html',['\
                       <html><center>\
                       <h2><a href="../../index.html">Home</a></h2>\
                       <h1>\
                       <a href="../'+str(PlasmaStart+Extrems[i-1])+'_'+str(PlasmaStart+Extrems[i])+'/index.html">&lt;=</a>\
                       '+str(i)+':Particular time <i>t=</i>'+str(Results[j][2])+' us\
                       <a href="../'+str(PlasmaStart+Extrems[i+1])+'_'+str(PlasmaStart+Extrems[i+2 if i<IndxDef[1]-1 else i])+'/index.html">=&gt;</a>\
                       </h1>\
                       <h2>Raw and smoothed data</h2>\
                       <a href="VAcharSmoothed.jpg"><img src="VAcharSmoothed.jpg"  width="30%"/></a>\
                       <br/>!NO FITS!</center>'],fmt="%s")
            
        else: # reasonable data
            #nearest_idx = np.where(abs(data[0]-float(Voltage[InflPointLoc]))==abs(data[0]-float(Voltage[InflPointLoc])).min())[0][0] # Azooz limit
            #nearest_idx = np.where(abs(data[0]-FitLimit)==abs(data[0]-FitLimit).min())[0][0] # HS limit
            #nearest_idx = np.where(abs(AzoozFit-0.0)==abs(AzoozFit-0.0).min())[0][0] # Azooz V_fl limit
            #print(nearest_idx)
            #print(Voltage[nearest_idx])
            errorflag=0
            Shift=int(float(len(data[0]))/(np.max(data[0])-np.min(data[0]))*VoltageShift)
            DataIndex=0
            if data[0][0]<0: # Are we sweeping up or down??
            #if 1==1:
                #Voltage = np.linspace(int(np.max(data[0])),int(np.min(data[0])),len(data[0]));
                nearest_idx = np.where(abs(AzoozFit-0.0)==abs(AzoozFit-0.0).min())[0][0] # Azooz V_fl limit
                VoltageReduced=Voltage[0:nearest_idx+Shift]
                DataIndex=int((nearest_idx+Shift)*1.0*len(data[0])/len(Voltage))
                #print(nearest_idx,nearest_idx+Shift,DataIndex)
                #print(VoltageReduced)
                #VoltageReduced=VoltageDef[0:nearest_idx+Shift]
                #print(VoltageReduced)
                try:
                    #ClassicalFitParms, po_cov = spo.curve_fit(Classical4Fn,data[0][0:nearest_idx+Shift], data[1][0:nearest_idx+Shift], ClassicalInitialGuess)
                    ClassicalFitParms, po_cov = spo.curve_fit(Classical4Fn,data[0][0:DataIndex], data[1][0:DataIndex], ClassicalInitialGuess)
                except RuntimeError: 
                    errorflag=1
            else:
                #Voltage = np.linspace(int(np.min(data[0])),int(np.max(data[0])),len(data[0]));
                #VoltageReduced=data[0][0:nearest_idx]
                nearest_idx = np.where(abs(AzoozFit-0.0)==abs(AzoozFit-0.0).min())[0][0] # Azooz V_fl limit
                #VoltageReduced=Voltage[0:len(data[0])-nearest_idx-Shift]
                VoltageReduced=Voltage[0:nearest_idx+Shift]
                DataIndex=int((len(data[0])-(nearest_idx+Shift))*1.0*len(data[0])/len(Voltage))
                #print(VoltageReduced)
                #VoltageReduced=VoltageDef[len(data[0])-nearest_idx-Shift:]
                #print(VoltageReduced)
                try:
                    #ClassicalFitParms, po_cov = spo.curve_fit(Classical4Fn,data[0][len(data[0])-nearest_idx-Shift:], data[1][len(data[0])-nearest_idx-Shift:], ClassicalInitialGuess)
                    ClassicalFitParms, po_cov = spo.curve_fit(Classical4Fn,data[0][DataIndex:], data[1][DataIndex:], ClassicalInitialGuess)
                except RuntimeError: 
                    errorflag=1
            if errorflag==1:
                    print('(Classic fit problem)',end="")
                    Results.append([\
                    i,\
                    (PlasmaStart+Extrems[i])/1000.0,\
                    (PlasmaStart+(Extrems[i]+Extrems[i+1])/2)/1000.0,\
                    (PlasmaStart+Extrems[i+1])/1000.0,\
                    nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan
                    ])
                    np.savetxt(TimeLineDir+'/index.html',['\
                       <html><center>\
                       <h2><a href="../../index.html">Home</a></h2>\
                       <h1>\
                       <a href="../'+str(PlasmaStart+Extrems[i-1])+'_'+str(PlasmaStart+Extrems[i])+'/index.html">&lt;=</a>\
                       '+str(i)+':Particular time <i>t=</i>'+str(Results[j][2])+' us\
                       <a href="../'+str(PlasmaStart+Extrems[i+1])+'_'+str(PlasmaStart+Extrems[i+2 if i<IndxDef[1]-1 else i])+'/index.html">=&gt;</a>\
                       </h1>\
                       <h2>Raw and smoothed data</h2>\
                       <a href="VAcharSmoothed.jpg"><img src="VAcharSmoothed.jpg"  width="30%"/></a>\
                       <br/>!NO FITS!</center>'],fmt="%s")
            else:    
                ClassicalFit=Classical4Fn(VoltageReduced,ClassicalFitParms[0],ClassicalFitParms[1],ClassicalFitParms[2],ClassicalFitParms[3])
                fig, ax1 = plt.subplots();
                plotsetup()
                ax1.plot(VoltageReduced, ClassicalFit, label="Classic fit")
                ax1.plot(Voltage, AzoozFit,label="Azooz fit")
                ax1.legend()
                ax1.axvline(float(Voltage[InflPointLoc]))
                ax1.axvline(float(Voltage[nearest_idx]))
                ax1.axvline(float(Voltage[nearest_idx+Shift]))
                plt.savefig(ShotNoStr+'/IndivVAcharsFits/VA'+str(i).zfill(2) +':'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'.jpg', bbox_inches='tight')
                plt.savefig(ShotNoStr+'/TimeLine/'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'/VAcharFit.jpg', bbox_inches='tight')
                np.savetxt(TimeLineDir+'/VoltageAzooz.dat',np.transpose([Voltage,AzoozFit]),fmt="%s")
                np.savetxt(TimeLineDir+'/VoltageReducedClassic4.dat',np.transpose([VoltageReduced,ClassicalFit]),fmt="%s")
                np.savetxt(TimeLineDir+'/SmoothedVachar.dat',np.transpose([VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],VAchar_current_smoothed[Extrems[i]:Extrems[i+1]]]),fmt="%s")
                np.savetxt(TimeLineDir+'/RawVachar.dat',np.transpose([VAchar_voltage[Extrems[i]:Extrems[i+1]],VAchar_current[Extrems[i]:Extrems[i+1]]]),fmt="%s")
                np.savetxt(TimeLineDir+'/DataToBeClass4Fitted.dat',np.transpose([data[0][0:DataIndex], data[1][0:DataIndex]]),fmt="%s")
                np.savetxt(TimeLineDir+'/nearest_idx.dat',[nearest_idx],fmt="%s")
                np.savetxt(TimeLineDir+'/DataIndex.dat',[DataIndex],fmt="%s")
                np.savetxt(TimeLineDir+'/Voltage@nearest_idx.dat',[Voltage[nearest_idx]],fmt="%s")
                np.savetxt(TimeLineDir+'/Shift.dat',[Shift],fmt="%s")
                np.savetxt(TimeLineDir+'/Voltage@nearest_idxPlusShift.dat',[Voltage[nearest_idx+Shift]],fmt="%s")
                plt.close()
                # plot Raw and smoothed data
                plt.plot(VAchar_voltage[Extrems[i]:Extrems[i+1]],VAchar_current[Extrems[i]:Extrems[i+1]],'+', label="raw data")
                plt.plot(VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],VAchar_current_smoothed[Extrems[i]:Extrems[i+1]], label="smoothed")
                plt.xlabel('$U$ [V]');plt.ylabel('$I$ [mA]');plt.legend()
                plt.savefig(TimeLineDir+'/RawAndSmoothed.jpg', bbox_inches='tight')
                plt.close()
                # plot Raw and smoothed data with fits
                plt.plot(VAchar_voltage[Extrems[i]:Extrems[i+1]],VAchar_current[Extrems[i]:Extrems[i+1]],',', label="raw data")
                plt.plot(VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],VAchar_current_smoothed[Extrems[i]:Extrems[i+1]], label="smoothed")
                plt.plot(VoltageReduced, ClassicalFit, label="Classic fit")
                plt.plot(Voltage, AzoozFit,label="Azooz fit")
                plt.axvline(float(Voltage[nearest_idx]))
                plt.axvline(float(Voltage[nearest_idx+Shift]))
                plt.axhline(0)
                #plt.plot(VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],AzoozFn(data[0],AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3]), label="Azooz fit")
                #plt.plot(VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],Classical4Fn(data[0],ClassicalFitParms[0],ClassicalFitParms[1],ClassicalFitParms[2],ClassicalFitParms[3]), label="Classic4p fit")
                plt.xlabel('$U$ [V]');plt.ylabel('$I$ [mA]');plt.legend()
                plt.savefig(TimeLineDir+'/RawAndSmoothedWithFits.jpg', bbox_inches='tight')
                plt.close()
                # plot Raw and smoothed data with fits zoom
                plt.plot(VAchar_voltage[Extrems[i]:Extrems[i+1]],VAchar_current[Extrems[i]:Extrems[i+1]],',', label="raw data")
                plt.plot(VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],VAchar_current_smoothed[Extrems[i]:Extrems[i+1]], label="smoothed")
                plt.plot(VoltageReduced, ClassicalFit, label="Classic 4 fit")
                plt.plot(Voltage, AzoozFit,label="Azooz fit")
                plt.axvline(float(Voltage[nearest_idx]))
                plt.axvline(float(Voltage[nearest_idx+Shift]))
                plt.axhline(0)
                plt.ylim(-10,10);
                plt.xlabel('$U$ [V]');plt.ylabel('$I$ [mA]');plt.legend()
                plt.savefig(TimeLineDir+'/RawAndSmoothedWithFitsZoom.jpg', bbox_inches='tight')
                plt.close()
                #plot Azooz stuff
                plt.plot(VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]],AzoozFn(data[0],AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3]), label="Azooz fit")
                plt.plot(Voltage,dAzoozFit_dVoltage*100,label="Azooz 1st derivative")
                plt.xlabel('$U$ [V]');plt.ylabel('$I$ [mA]');plt.legend()
                plt.savefig(TimeLineDir+'/AzoozStuff.jpg', bbox_inches='tight')
                plt.close()
                DischTime=(PlasmaStart+(Extrems[i]+Extrems[i+1])/2)/1000.0
                V_BPP=sum(Potential_smoothed[(Extrems[i]):(Extrems[i+1])])/\
                len(Potential_smoothed[(Extrems[i]):(Extrems[i+1])])
                Results.append([\
                i,\
                (PlasmaStart+Extrems[i])/1000.0,\
                DischTime,\
                (PlasmaStart+Extrems[i+1])/1000.0,\
                float(Btoroidal['data'][DischTime*1000:DischTime*1000+1]),\
                AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3],\
                #int(InflPointLoc),\
                Voltage[nearest_idx],\
                nearest_idx,\
                float(Voltage[InflPointLoc]),\
                ClassicalFitParms[0],ClassicalFitParms[1],ClassicalFitParms[2],ClassicalFitParms[3],\
                (np.exp(AzoozFitParms[0])+AzoozFitParms[3])/(np.exp(-AzoozFitParms[0])+AzoozFitParms[3]),\
                V_BPP,\
                (V_BPP-ClassicalFitParms[1])/ClassicalFitParms[2],\
                sum(dz[(PlasmaStart+Extrems[i]):(PlasmaStart+Extrems[i+1])])/\
                len(dz[(PlasmaStart+Extrems[i]):(PlasmaStart+Extrems[i+1])])\
                ])
                j=i-IndxDef[0] # Tuning purposes ..
                k=0
                np.savetxt(TimeLineDir+'/Results',Results[j],fmt="%s")
                np.savetxt(TimeLineDir+'/index.html',['\
                       <html>\
                       <script type="text/x-mathjax-config">MathJax.Hub.Config({tex2jax: {inlineMath: [[\'$\',\'$\'], [\'\\(\',\'\\)\']]}});</script>\
                       <script src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>\<body><center>\
                       <h2><a href="../../index.html">Home</a></h2>\
                       <h1>\
                       <a href="../'+str(PlasmaStart+Extrems[i-1])+'_'+str(PlasmaStart+Extrems[i])+'/index.html">&lt;=</a>\
                       '+str(i)+':Particular time <i>t=</i>'+str(Results[j][2])+' us\
                       <a href="../'+str(PlasmaStart+Extrems[i+1])+'_'+str(PlasmaStart+Extrems[i+2 if i<IndxDef[1]-1 else i])+'/index.html">=&gt;</a>\
                       </h1>\
                       <h2>Raw and smoothed data with fits</h2>\
                       <a href="RawAndSmoothed.jpg"><img src="RawAndSmoothed.jpg"  width="30%"/></a>\
                       <!--<h2>Raw and smoothed data with fits</h2>-->\
                       <a href="RawAndSmoothedWithFits.jpg"><img src="RawAndSmoothedWithFits.jpg"  width="30%"/></a>\
                       <a href="RawAndSmoothedWithFitsZoom.jpg"><img src="RawAndSmoothedWithFitsZoom.jpg"  width="30%"/></a><br/>\
                       <h2>Smoothed data with fits</h2>\
                       <a href="VAcharFit.jpg"><img src="VAcharFit.jpg"  width="30%"/></a><br/>\
                       <h2>Results</h2>\
                       <b>i</b>:'+'{:2d}'.format(Results[j][0]) +\
                ', <b>Time:</b>From:'+'{:4.2f}'.format(Results[j][1]) +\
                ',via:'+'{:4.2f}'.format(Results[j][2]) +\
                ',to:'+'{:4.2f}'.format(Results[j][3]) +\
                ',$B_t$:'+'{:4.2f}'.format(Results[j][4]) +\
                '</br><b>Azooz fit:</b> $a_1$='+'{:+5.3e}'.format(Results[j][5]) +\
                ',$a_2$='+'{:4.2f}'.format(Results[j][6]) +\
                ',$a_3$='+'{:+5.3e}'.format(Results[j][7]) +\
                ',$a_4$='+'{:4.2f}'.format(Results[j][8]) +\
                '</br><b>Inflex point:</b>'+'{:+4.0f}'.format(Results[j][9]) +\
                ',Infl'+'{:+4.0f}'.format(Results[j][10]) +\
                ',Infl'+'{:+3.2f}'.format(Results[j][11]) +\
                '</br><b>Class fit:</b> $I_{sat}$='+'{:4.2f}'.format(Results[j][12]) +\
                ',$\\varphi_{pl}$='+'{:4.2f}'.format(Results[j][13]) +\
                ',$T_e$='+'{:4.2f}'.format(Results[j][14]) +\
                ',$C$='+'{:4.2f}'.format(Results[j][15]) +\
                ',$Ratio_{sat}$='+'{:4.2f}'.format(Results[j][16]) +\
                ',$\\varphi_{fl}$='+'{:4.2f}'.format(Results[j][17]) +\
                ',$\\alpha$='+'{:4.2f}'.format(Results[j][18]) +\
                ',Plasma position='+'{:4.2f}'.format(Results[j][19]) +\
                '<h2>Azooz stuff</h2>\
                 <br/><a href="AzoozStuff.jpg"><img src="AzoozStuff.jpg" width="30%"/></a><br/>\
                 <h2>Data</h2>\
                 <a href="RawVachar.dat">Raw VA char data</a><br/>\
                 <a href="SmoothedVachar.dat">Smoothed VA char data</a><br/>\
                 <a href="VoltageAzooz.dat">Azooz fit VA data</a><br/>\
                 <a href="VoltageReducedClassic4.dat">Classic 4 par fit VA data</a><br/>\
                 <a href="Results">Final results</a><br/>\
                 </center></body></html>'],fmt="%s")
                AFP0=AFP0+AzoozFitParms[0];AFP1=AFP1+AzoozFitParms[1];AFP2=AFP2+AzoozFitParms[2];AFP3=AFP3+AzoozFitParms[3]
                CFP0=CFP0+ClassicalFitParms[0];CFP1=CFP1+ClassicalFitParms[1];CFP2=CFP2+ClassicalFitParms[2];CFP3=CFP3+ClassicalFitParms[3];Pocet=Pocet+1
        #print(end=",")
        #np.savetxt(TimeLineDir+'/chars.dat',\
        #           [[VAchar_voltage[Extrems[i]:Extrems[i+1]][j],\
        #             VAchar_current[Extrems[i]:Extrems[i+1]][j],\
        #             VAchar_voltage_smoothed[Extrems[i]:Extrems[i+1]][j],\
        #             VAchar_current_smoothed[Extrems[i]:Extrems[i+1]][j],\
        #              AzoozFn(data[0][j],AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3]), \
        #             dAzoozFit_dVoltage[int(1000*(j*1.0)/len(data[0]))],\
        #             Classical4Fn(data[0][j],ClassicalFitParms[0],ClassicalFitParms[1],ClassicalFitParms[2],ClassicalFitParms[3])] \
        #             for j in xrange(len(data[0]))],delimiter=" ",fmt="%.3f")

    np.savetxt(ShotNoStr+'/Data/Results',Results,fmt="%+8.2f") 
    #print(IndvSetup[3][0])
    np.savetxt(ShotNoStr+'/Data/AzoozFitCoefficients',['Fit parameters generation *estimate*->&lt;final&gt;: No:{:d}. <br/>'.format(Pocet)\
                                                       +'Azooz fit: a1:{:4.2f} -> {:4.2f};'.format(IndvSetup[3][0],AFP0/Pocet)\
                                                       +'a2:{:4.2f} -> {:4.2f};'.format(IndvSetup[3][1],AFP1/Pocet)\
                                                       +'a3:{:4.2f} -> {:4.2f};'.format(IndvSetup[3][2],AFP2/Pocet)\
                                                       +'a4:{:4.2f} -> {:4.2f};<br/>'.format(IndvSetup[3][3],AFP3/Pocet)\
                                                       +'Classic fit: Isat:{:4.2f} -> {:4.2f};'.format(IndvSetup[4][0],CFP0/Pocet)\
                                                       +'varphi_pl:{:4.2f} -> {:4.2f};'.format(IndvSetup[4][1],CFP1/Pocet)\
                                                       +'Te:{:4.2f} -> {:4.2f};'.format(IndvSetup[4][2],CFP2/Pocet)\
                                                       +'C:{:4.3f} -> {:4.3f};'.format(IndvSetup[4][3],CFP3/Pocet)\
                                                      ],fmt="%s")
    return Results
    
def htmlgeneration(ShotNo,UBt,VacuumShot,PlasmaStart, PlasmaEnd, Time, Extrems, IndxDef, Results):
    print("html generation ...")
    ShotNoStr=str(ShotNo)

    os.system('rm '+str(ShotNo)+'/index.html');
    fileid = open(str(ShotNo)+'/index.html','a+')
    fileid.write('<html><head><title>VA char @ GOLEM</title>\
    <style>\
    div.grafy30 {display:inline-block;width: 30%;}\
    div.grafy20 {display:inline-block;width: 20%;}\
    .grafy20 img {width: 100%;}\
    .grafy30 img {width: 100%;}\
    </style>\
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">\
    <style></style>\
    <script type="text/x-mathjax-config">MathJax.Hub.Config({tex2jax: {inlineMath: [[\'$\',\'$\'], [\'\\(\',\'\\)\']]}});</script>\
    <script src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>\
    </head><body><center>')
    fileid.write('<h2><a href="../index.html">Home</a></h2>')
    fileid.write('<h1>BPP experiments @ the tokamak GOLEM</h1>')
    fileid.write('<h3>Essential parameters:</h3><table><tr><td><ul>\
    <li>Current resistor $R=$'+str(HWProudovyResistor)+' $\Omega$</li>\
    <li>Voltage divider: 1:'+str(HWNapetovyDelic)+'</li>\
    <li>The shot: <a href="http://golem.fjfi.cvut.cz/shots/'+ShotNoStr+'">'+ShotNoStr+'</a></li>\
    <li>The DAS used: <a href="http://golem.fjfi.cvut.cz/shots/'+ShotNoStr+'/DAS/1011Papouch_St.ON/">Papouch_St</a></li>\
    <li>$U_{B_t}$='+str(UBt)+'</li>\
    <li>Vacuum shot /plasma position issues/:<a href="http://golem.fjfi.cvut.cz/shots/'+str(VacuumShot)+'">'+str(VacuumShot)+'</a></li>\
    <li>Voltage shift='+str(VoltageShift)+'</li>\
    </ul></td></tr></table>')
    fileid.write('<div class="grafy30"><h2>Experimental setup</h2>')
    fileid.write('<h3>Ball pen probe @ North-East port</h3>')
    fileid.write('<center>\
    <a href="http://golem.fjfi.cvut.cz/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/setup/vI/ExpSetup-BPP.png">\
    <img src="../setup/ExpSetup-BPP.png" width="100%"></a><br>\
    <a href="/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/setup/">setup</a></div>')
    fileid.write('<div class="grafy30"><h2>The discharge #'+ShotNoStr+'</h2>')
    fileid.write('<h3>Basic diagnostics</h3>\
    <a href="http://golem.fjfi.cvut.cz/shots/'+ShotNoStr+'/"><img src="Basics/graphres.png" width="100%"></a><br/>\
    Plasma start:<a href="ShotNo/PlasmaStart">'+str(PlasmaStart)+'</a> us\
    Plasma end:<a href="ShotNo/PlasmaEnd">'+str(PlasmaEnd)+'</a> us</div>\
    <h3>The BPP data vs diagnostics</h3>\
    <a href="Basics/BasicDiagWithBPP.jpg"><img src="Basics/BasicDiagWithBPP.jpg" width="30%"></a><br/>\
    <a href="ShotNo/">Reference shot data</a><br/>\
    <div class="grafy30"><h3>The BPP raw data</h3>\
    <a href="ShotNo/RawData.jpg"><img src="Basics/RawData.jpg"></a><br/></div>')
    fileid.write('\
    <div class="grafy30"><h3>Maxima and minima localization in bpp_voltage</h3>\
    <a href="ShotNo/MaximMinims.jpg"><img src="Basics/MaximMinims.jpg"></a><br/>\
    <a href="Data/Extrems">Extrems</a><br/>\
    <a href="Data/minima">minima</a><br/>\
    <a href="Data/maxima">maxima</a></div>\
    ')
    fileid.write('<h3>Final movie</h3>')
    fileid.write('<img src="finalmovie.gif"><br/><a href="IndivVAchars/">Individual figures</a></br>\
                 .dat format: VAchar_voltage,VAchar_current,VAchar_voltage_smoothed,VAchar_current_smoothed,Azooz fit,dAzoozFit/dVoltage,ClassicFit<br/><div>')
    StartIndx=0
    if IndxDef[1]==0:
        IndxDef[1]=len(Extrems)-1 # test purposes
    for i in range(IndxDef[0]-IndxDef[0],IndxDef[1]-IndxDef[0]): 
        if os.path.isfile(str(ShotNo)+'/IndivVAcharsFits/VA'+str(i).zfill(2) +':'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'.jpg'): # Azooz fit problem
            fileid.write('<div class="grafy20">\
            <center><a href="TimeLine/'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'/index.html">\
            <img src="IndivVAcharsFits/VA'+str(i).zfill(2) +':'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'.jpg" width="100%"></a>\
            </br>\
            <b>i</b>:'+'{:2d}'.format(Results[i][0]) +\
            ', <b>Time:</b>From:'+'{:4.2f}'.format(Results[i][1]) +\
            ',via:'+'{:4.2f}'.format(Results[i][2]) +\
            ',to:'+'{:4.2f}'.format(Results[i][3]) +\
            ',$B_t$'+'{:4.2f}'.format(Results[i][4]) +\
            '</br><b>Azooz fit:</b> $a_1$='+'{:+5.3e}'.format(Results[i][5]) +\
            ',$a_2$='+'{:4.2f}'.format(Results[i][6]) +\
            ',$a_3$='+'{:+5.3e}'.format(Results[i][7]) +\
            ',$a_4$='+'{:4.2f}'.format(Results[i][8]) +\
            '</br><b>Class fit:</b> $I_{sat}$='+'{:4.2f}'.format(Results[i][12]) +\
            ',$\\varphi_{pl}$='+'{:4.2f}'.format(Results[i][13]) +\
            ',$T_e$='+'{:4.2f}'.format(Results[i][14]) +\
            ',$Ratio_{sat}$='+'{:4.2f}'.format(Results[i][15]) +\
            ',$\\varphi_{fl}$='+'{:4.2f}'.format(Results[i][16]) +\
            ',$\\alpha$='+'{:4.2f}'.format(Results[i][17]) +\
            '<br><!--<a href=TimeLine/'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'/chars.dat>data</a>--></center></div>')
        else:
            fileid.write('<div class="grafy20">\
            <center><a href="TimeLine/'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'/index.html">\
            <img src="IndivVAchars/VA'+str(i).zfill(2) +':'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'.jpg" width="100%"></a>\
            </br>\
            <b>i</b>:'+'{:2d}'.format(Results[i][0]) +\
            ', <b>Time:</b>From:'+'{:4.2f}'.format(Results[i][1]) +\
            ',via:'+'{:4.2f}'.format(Results[i][2]) +\
            ',to:'+'{:4.2f}'.format(Results[i][3]) +\
            '<br/><b>!!Azooz fit problem! -> No Fits !</b><br/><br/><br/><br/>' +\
            '<br><a href=Data/VA'+str(i).zfill(2) +':'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'.dat>data</a></center></div>')
    fileid.write('</div>\
    <img src="GlobalResults.jpg"></a><br/>')
    fileid.write('<table border="1">\
    <tr><th>$i$</th><th colspan=3>Time</th><th>$B_t$</th>\
    <th colspan=4>Azooz fit parameters<br/>$I(a,V_{bias})=\exp[a_1 \\tanh\{(V_{bias}+a_2)/a_3\}]+a_4$</th>\
    <th colspan=3>InflexPoint</th>\
    <th colspan=4>Classical fit parameters<br/>$I(V_{bias},I_{sat},V_{fl},T_e,C)=-I_{sat}*(1-C*(V_{bias}-V_{fl}))*(1-\exp((V_{bias}-V_{fl})/T_e)$</th>\
    <th>$Ratio_{sat}$</th><th>$\\varphi_{fl}$</th><th>$\\alpha$</th><th>$Pl_{pos}$ [mm]</th></tr>\
    <tr><th></th><th>Start</th><th>Middle</th><th>End</th>\
    <th></th>\
    <th>$a_1$</th><th>$a_2$</th><th>$a_3$</th><th>$a_4$</th>\
    <th>Data row</th><th>Data row</th><th>Voltage</th>\
    <th>$I_{sat}$</th><th>$\\varphi_{pl}$</th><th>$T_e$</th><th>C</th><th></th><th></th><th></th><th></th></tr>')
    for i in range(IndxDef[0]-IndxDef[0],IndxDef[1]-IndxDef[0]):
        fileid.write('<tr>\
        <td>'+'{:2.0f}'.format(Results[i][0]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][1]) +'</td>\
        <td><a href=TimeLine/'+str(PlasmaStart+Extrems[i])+'_'+str(PlasmaStart+Extrems[i+1])+'/index.html>'+'{:4.2f}'.format(Results[i][2]) +'</a></td>\
        <td>'+'{:4.2f}'.format(Results[i][3]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][4]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][5]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][6]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][7]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][8]) +'</td>\
        <td>'+'{:3.0f}'.format(Results[i][9]) +'</td>\
        <td>'+'{:3.0f}'.format(Results[i][10]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][11]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][12]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][13]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][14]) +'</td>\
        <td>'+'{:4.3f}'.format(Results[i][15]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][16]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][17]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][18]) +'</td>\
        <td>'+'{:4.2f}'.format(Results[i][19]) +'</td>\
        </tr>')
    fileid.write('</table><br/><a href="Data/Results">data</a><br/>')
    fileid.write(np.loadtxt(ShotNoStr+'/Data/AzoozFitCoefficients',dtype='str'))
    fileid.write('</div><h2>Experimental Photo</h2>') 
    fileid.write('<center><img src="http://golem.fjfi.cvut.cz/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/setup/vI/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('<h2>References</h2>')
    fileid.write('<ul>\
    <li><a href="https://www.mathworks.com/matlabcentral/fileexchange/19312-langmuir-probe-data-analysis-code">Aasim Azooz: Langmuir probe data analysis code. Mathworks.</a></li>\
    </ul>')
    fileid.write('</body></html>')
    fileid.close()
    
def GlobalGraphGen(ShotNo,VacuumShot,PlasmaStart, PlasmaEnd, Time, Extrems, Uloop, Btoroidal, Iplasma,\
                   VAchar_current, Potential, Potential_smoothed, PlasmaPosition, dz, Results):
    print("Global graph generation .. ")
    ShotNoStr=str(ShotNo)
    Borders=0,len(Extrems) # Initial and final problems out (in future option)
    ResultsArr=np.array(Results)
    f = plt.figure(figsize=(20.0, 10.0))
    f,ax = plt.subplots(17,sharex=True);plt.subplots_adjust(hspace=0.001)
    f.set_size_inches(10,20)
    ax[0].set_title('#' + ShotNoStr)
    ax[0].plot(Time/1000,Uloop['data'][PlasmaStart:PlasmaEnd]);ax[0].set_ylabel('$U_l$ [V]')
    ax[1].plot(Time/1000,Btoroidal['data'][PlasmaStart:PlasmaEnd]);ax[1].set_ylabel('$B_t$ [T]')
    ax[1].plot(ResultsArr[:,2],ResultsArr[:,4],'+')
    ax[2].plot(Time/1000,Iplasma['data'][PlasmaStart:PlasmaEnd]/1000);ax[2].set_ylabel('$I_p$ [kA]')
    
    ax[3].plot(Time/1000,VAchar_current[0:,1]/HWProudovyResistor*CurrentScale,label='VA char: I');ax[3].set_ylabel('$I$ [mA]');ax[3].legend(loc=0)
    ax[3].plot(Time/1000,VAchar_current[0:,1]*0)
    
    ax[4].plot(Time/1000,dz[PlasmaStart:PlasmaEnd]);ax[4].set_ylabel('$<Pl_{pos}>$ [mm]')
    ax[4].plot(Time/1000,dz[PlasmaStart:PlasmaEnd]*0)
    ax[4].plot(ResultsArr[:,2],ResultsArr[:,19],'+')
    
    ax[5].plot(ResultsArr[:,2],ResultsArr[:,5],'+');ax[5].set_ylabel('$a_1$')
    ax[6].plot(ResultsArr[:,2],ResultsArr[:,6],'+');ax[6].set_ylabel('$a_2$')
    ax[7].plot(ResultsArr[:,2],ResultsArr[:,7],'+');ax[7].set_ylabel('$a_3$')
    ax[8].plot(ResultsArr[:,2],ResultsArr[:,8],'+');ax[8].set_ylabel('$a_4$')
    ax[9].plot(ResultsArr[:,2],ResultsArr[:,9],'+');ax[9].set_ylabel('InflPoint [V]')
    ax[10].plot(ResultsArr[Borders[0]:Borders[1],2],ResultsArr[Borders[0]:Borders[1],12],'+');ax[10].set_ylabel('$I_{sat}$')
    ax[11].plot(ResultsArr[Borders[0]:Borders[1],2],ResultsArr[Borders[0]:Borders[1],13],'+');ax[11].set_ylabel('$\\varphi_{pl}$ [V]')
    ax[12].plot(ResultsArr[Borders[0]:Borders[1],2],ResultsArr[Borders[0]:Borders[1],14],'+');ax[12].set_ylabel('$T_e$ [eV]')
    ax[13].plot(ResultsArr[Borders[0]:Borders[1],2],ResultsArr[Borders[0]:Borders[1],15],'+');ax[13].set_ylabel('$C$ []')
    ax[14].plot(ResultsArr[Borders[0]:Borders[1],2],ResultsArr[Borders[0]:Borders[1],16],'+');ax[14].set_ylabel('$Ratio_{sat}$ []')
    ax[15].plot(ResultsArr[Borders[0]:Borders[1],2],ResultsArr[Borders[0]:Borders[1],18],'+');ax[15].set_ylabel('$\\alpha$ []')
    ax[16].plot(Time/1000,Potential_smoothed);ax[16].set_ylabel('$\\varphi_{fl}$ [V]')
    ax[16].plot(Time/1000,Potential_smoothed*0)
    ax[16].plot(ResultsArr[:,2],ResultsArr[:,17],'+')
    
    ax[16].set_xlim(PlasmaStart/1000*8/10,PlasmaEnd/1000*11/10);ax[16].set_xlabel('$t$ [ms]')

    plt.savefig(ShotNoStr+'/GlobalResults.jpg', bbox_inches='tight')
    if is_interactive():plt.show();
    plt.close();    
    
def Initials(ShotNo):
    ShotNoStr=str(ShotNo)
    mkdir(ShotNoStr+'/IndivVAchars');mkdir(ShotNoStr+'/IndivVAcharsFits');
    mkdir(ShotNoStr+'/Basics');mkdir(ShotNoStr+'/Data');mkdir(ShotNoStr+'/TimeLine');
    
    
def PlasmaParametersDef(ShotNo):
    return [\
    int(np.loadtxt(str(ShotNo)+'/Basics/PlasmaStart')),\
    int(np.loadtxt(str(ShotNo)+'/Basics/PlasmaEnd'))\
            ]
    
def PlasmaPositionFn(ShotNo):
    ShotNoStr=str(ShotNo)
    mc5_p=np.loadtxt(ShotNoStr+'/Basics/mc5_p.txt')
    mc5_v=np.loadtxt(ShotNoStr+'/Basics/mc5_v.txt')
    mc13_p=np.loadtxt(ShotNoStr+'/Basics/mc13_p.txt')
    mc13_v=np.loadtxt(ShotNoStr+'/Basics/mc13_v.txt')
    #upper mirnov coil
    mc5_p  = integrate(mc5_p[0:,1], coeff=1/(3705 * 1e-6))
    mc5_v  = integrate(mc5_v[0:,1], coeff=1/(3705 * 1e-6))
    mc5 = mc5_p - mc5_v # pure signal from plasma
    #bottom mirnov coil
    mc13_p = integrate(mc13_p[0:,1], coeff=1/(3705 * 1e-6))
    mc13_v = integrate(mc13_v[0:,1], coeff=1/(3705 * 1e-6))
    mc13 = mc13_p - mc13_v # pure signal from plasma
    dz= 93 * (mc5 - mc13)/(mc5 + mc13) # calculation of displacement
    for i in range(len(dz)): 
        if dz[i]>50:dz[i]=50
        if dz[i]<-50:dz[i]=-50    
    return dz
    
    
def IndividualShot(IndvSetup,IndxDef):
    ShotNo=IndvSetup[0];print( 'Doing', ShotNo)
    VacuumShot=IndvSetup[1]
    ShotNoStr=str(ShotNo)
    Initials(ShotNo)
    PlasmaParameters=PlasmaParametersDef(ShotNo)
    PlasmaStart=PlasmaParameters[0]
    PlasmaEnd=PlasmaParameters[1]
    
    Time = np.linspace(PlasmaStart,PlasmaEnd,(PlasmaEnd-PlasmaStart))
    VAchar_voltage=np.loadtxt(ShotNoStr+'/Basics/VAchar_voltage.txt')[PlasmaStart:PlasmaEnd]
    VAchar_current=np.loadtxt(ShotNoStr+'/Basics/VAchar_current.txt')[PlasmaStart:PlasmaEnd]
    VAchar_voltage_smoothed=running_mean(VAchar_voltage[:,1]*HWNapetovyDelic,SmoothCoefficient)
    VAchar_current_smoothed=running_mean(VAchar_current[:,1]/HWProudovyResistor*CurrentScale,SmoothCoefficient)
    
    Btoroidal=np.load(ShotNoStr+'/Basics/Btoroidal.npz');
    Iplasma=np.load(ShotNoStr+'/Basics/Iplasma.npz');
    Uloop=np.load(ShotNoStr+'/Basics/Uloop.npz');
    
    Potential=np.loadtxt(ShotNoStr+'/Basics/Potential.txt')[PlasmaStart:PlasmaEnd]
    Potential_smoothed=running_mean(Potential[:,1]*HWNapetovyDelic,SmoothCoefficient)

    PlasmaPosition=PlasmaPositionFn(ShotNo)

    
    DischargeBasics(ShotNo,VacuumShot,PlasmaStart, PlasmaEnd, Time, Uloop, Btoroidal, Iplasma,\
                   VAchar_voltage, VAchar_current,Potential, Potential_smoothed, PlasmaPosition)
    Extrems=FindingSweepIntervals(ShotNo,PlasmaStart, PlasmaEnd, Time, VAchar_voltage)
    Results=Fits(IndvSetup,PlasmaParameters,Extrems,Time,VAchar_voltage[:,1]*HWNapetovyDelic,VAchar_current[:,1]/HWProudovyResistor*CurrentScale,VAchar_voltage_smoothed,VAchar_current_smoothed,\
                 Potential_smoothed, Uloop, Btoroidal, Iplasma, PlasmaPosition, IndxDef)
    htmlgeneration(ShotNo,IndvSetup[2],VacuumShot,PlasmaStart, PlasmaEnd, Time, Extrems, IndxDef, Results)
    GlobalGraphGen(ShotNo,VacuumShot,PlasmaStart, PlasmaEnd, Time, Extrems, Uloop, Btoroidal, Iplasma,\
                   VAchar_current,Potential, Potential_smoothed, PlasmaPosition, PlasmaPosition, Results)


def GlobalFigs():
    mkdir('GlobalFigs')
    Time = np.linspace(0,40000,40000)

    for graph in GraphSetup:
        FigDef()
        for shot in Setup:
            data=np.load(str(shot[0])+'/Basics/'+graph[0]+'.npz')
            plt.plot(Time/1000,data['data'],label='#'+str(shot[0])+' ($U_{B_t}$='+str(shot[2])+')');
            plt.legend(loc=0)
        plt.title('Reproducibility for: '+graph[2])
        plt.xlabel('$t$ [ms]')
        plt.ylabel(graph[1])
        plt.savefig('GlobalFigs/Reproducibility_'+graph[0]+'.jpg', bbox_inches='tight')
        if is_interactive():plt.show();
        plt.close();      

def Dependences():
    # Dependences with respect to time (Results[0:,2])
    for result in ResultsDef:
        FigDef()
        for shot in Setup:
            Results=np.loadtxt(str(shot[0])+'/Data/Results')
            plt.plot(Results[0:,2],Results[0:,result[0]],shot[6],label=str(shot[0])+' ($U_{B_t}$='+str(shot[2])+')')
        plt.title('Crosscheck for: '+result[2])
        plt.ylabel(result[2])
        plt.xlabel('$t$ [ms]')
        plt.legend(loc=0)
        plt.savefig('GlobalFigs/TimeCrossCheck_'+str(result[1])+'.jpg', bbox_inches='tight')
        if is_interactive():plt.show();
        plt.close(); 
        os.system('rm GlobalFigs/TimeCrossCheck_'+str(result[1])+'.dat')
        with open('GlobalFigs/TimeCrossCheck_'+str(result[1])+'.dat', "a") as myfile:
            for i in xrange(len(Results[0:,3])):
                k=1;datarow=str(Results[i][3])
                #print(lst)
                for shot in Setup:
                    datarow=datarow+result[4].format(Results[0:,result[0]][i])
                    k=k+1
                myfile.write(datarow+'\n')
    # and Ufl x Plpot speciality:
    FigDef()
    for shot in Setup:
        Results=np.loadtxt(str(shot[0])+'/Data/Results')
        plt.plot(Results[0:,2],Results[0:,13],shot[6],label=str(shot[0])+' $\\varphi_{fl}$')
        plt.plot(Results[0:,2],Results[0:,17],shot[6],label=str(shot[0])+' $\\varphi_{pl}$')
    plt.title('Crosscheck for: $\\varphi_{fl}$ x  $\\varphi_{pl}$')
    plt.ylabel('$V_{fl}$ [V]')
    plt.xlabel('$t$ [ms]')
    plt.legend(loc=0)
    plt.savefig('GlobalFigs/TimeCrossCheck_FlowPot_PlPot.jpg', bbox_inches='tight')
    if is_interactive():plt.show();
    plt.close(); 
        
        

    # Dependences with respect to Bt (Results[0:,4])
    for result in ResultsDef:
        FigDef()
        for shot in Setup:
            Results=np.loadtxt(str(shot[0])+'/Data/Results')
            plt.plot(Results[0:,4],Results[0:,result[0]],shot[6],label=str(shot[0])+' ($U_{B_t}$='+str(shot[2])+')')
        plt.title('Crosscheck for: '+result[2])
        plt.ylabel(result[2])
        plt.xlabel('$B_t$ [T]')
        plt.legend(loc=0)
        plt.savefig('GlobalFigs/BtCrossCheck_'+str(result[1])+'.jpg', bbox_inches='tight')
        if is_interactive():plt.show();
        plt.close(); 
        os.system('rm GlobalFigs/BtCrossCheck_'+str(result[1])+'.dat')
        with open('GlobalFigs/BtCrossCheck_'+str(result[1])+'.dat', "a") as myfile:
            for i in xrange(len(Results[0:,4])):
                k=1;datarow=str(Results[i][4])
                for shot in Setup:
                    datarow=datarow+result[4].format(Results[0:,result[0]][i])
                    k=k+1
                myfile.write(datarow+'\n')
    # and Ufl x Plpot speciality:
    FigDef()
    for shot in Setup:
        Results=np.loadtxt(str(shot[0])+'/Data/Results')
        plt.plot(Results[0:,4],Results[0:,13],shot[6],label=str(shot[0])+' $\\varphi_{fl}$')
        plt.plot(Results[0:,4],Results[0:,17],shot[6],label=str(shot[0])+' $\\varphi_{pl}$')
    plt.title('Crosscheck for: $\\varphi_{fl}$ x  $\\varphi_{pl}$')
    plt.ylabel('$V_{fl}$ [V]')
    plt.xlabel('$B_t$ [T]')
    plt.legend(loc=0)
    plt.savefig('GlobalFigs/BtCrossCheck_FlowPot_PlPot.jpg', bbox_inches='tight')
    if is_interactive():plt.show();
    plt.close(); 
    

def HtmlGeneration():
    print("html generation ...")

    os.system('rm index.html');
    fileid = open('index.html','a+')
    fileid.write('<html><head><title>VA char @ GOLEM</title>\
    <style>\
    div.grafy30 {display:inline-block;width: 30%;}\
    div.grafy20 {display:inline-block;width: 20%;}\
    .grafy20 img {width: 100%;}\
    .grafy30 img {width: 100%;}\
    </style>\
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">\
    <style></style>\
    <script type="text/x-mathjax-config">MathJax.Hub.Config({tex2jax: {inlineMath: [[\'$\',\'$\'], [\'\\(\',\'\\)\']]}});</script>\
    <script src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>\
    </head><body><center>')
    fileid.write('<h2>BPP experiments @ the tokamak GOLEM</h2>')
    fileid.write('<h3>Essential parameters:</h3><table><tr><td><ul>\
    <li>Date:'+Datum+'</li>\
    <li>BPP used for Plasma Float potential $\\varphi_{fl}$, BPP_LP used in sweped mode</li>\
    <li>Probe position:'+Probe_position+'</li>\
    <li>The VA characteristic resistor:'+str(HWProudovyResistor)+'$\Omega$</li>\
    <li>Sweeping frequency:'+str(HFSweepFrequency)+'</li>\
    </ul></td></tr></table>')
    fileid.write('<div class="grafy30"><h2>Experimental setup</h2>')
    fileid.write('<h3>Ball pen probe @ North-East port</h3>')
    fileid.write('<center>\
    <a href="http://golem.fjfi.cvut.cz/wikiraw/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/setup/vI/ExpSetup-BPP.png">\
    <img src="setup/ExpSetup-BPP.png" width="100%"></a><br>\
    <a href="'+Base+'/setup/">setup</a></div>\
    <h2>Discharges involved</h2>\
    <table border="1"><tr>\
    <th>Discharge #</th>\
    <th>Vacuum #</th>\
    <th>$U_{B_t}$ [V]</th>\
    <th>Analysis</th>\
    <th>Comment</th></tr>')
    for shot in Setup: 
        fileid.write('<tr>\
        <td><a href="http://golem.fjfi.cvut.cz/shots/'+str(shot[0])+'">#'+str(shot[0])+'</td>\
        <td><a href="http://golem.fjfi.cvut.cz/shots/'+str(shot[1])+'">#'+str(shot[1])+'</td>\
        <td>'+str(shot[2])+'</td>\
        <td><a href="'+str(shot[0])+'/index.html">link</td>\
        <td>'+str(shot[5])+'</td>\
        </tr>')
    for shot in Setup: 
        fileid.write('</table><div class="grafy30"><h2>'+str(shot[0])+'</h2>\
        <a href="'+str(shot[0])+'/index.html">\
        <img src="'+str(shot[0])+'/Basics/BasicDiagWithBPP.jpg" width="100%"></a></div>')    
    fileid.write('<h2>Reproducibility</h2>')
    for graph in GraphSetup:
        fileid.write('<div class="grafy30"><h2>'+graph[0]+':</h2>\
        <a href="GlobalFigs/Reproducibility_'+graph[0]+'.jpg">\
        <img src="GlobalFigs/Reproducibility_'+graph[0]+'.jpg" width="100%"></a></div>')
    fileid.write('<h2>Fits definition</h2>\
                 Azooz fit definition:$I(a,V_{bias})=\exp[a_1 \\tanh\{(V_{bias}+a_2)/a_3\}]+a_4$<br/>\
                 Classical 4 par fit definition:$I(V_{bias},I_{sat},V_{fl},T_e,C)=-I_{sat}*(1-C*(V_{bias}-V_{fl}))*(1-\exp((V_{bias}-V_{fl})/T_e)$<br/>\
                 Classical 3 par fit definition:$I(V_{bias},I_{sat},V_{fl},T_e)=-I_{sat}*(1-\exp((V_{bias}-V_{fl})/T_e)$')    
    fileid.write('<h2>Time cross compares</h2>')
    for result in ResultsDef:
        fileid.write('<div class="grafy30"><h2>'+str(result[1])+'</h2>\
        <a href="GlobalFigs/TimeCrossCheck_'+str(result[1])+'.jpg">\
        <img src="GlobalFigs/TimeCrossCheck_'+str(result[1])+'.jpg" width="100%"></a></br>\
        <a href="GlobalFigs/TimeCrossCheck_'+str(result[1])+'.dat">data</a></div>')
    fileid.write('<div class="grafy30"><h2> $\\varphi_{fl}$ x  $\\varphi_{pl}$</h2>\
    <a href="GlobalFigs/TimeCrossCheck_FlowPot_PlPot.jpg">\
    <img src="GlobalFigs/TimeCrossCheck_FlowPot_PlPot.jpg" width="100%"></a></br></div>')
    fileid.write('<h2>$B_t$ cross compares</h2>')
    for result in ResultsDef:
        fileid.write('<div class="grafy30"><h2>'+str(result[1])+'</h2>\
        <a href="GlobalFigs/BtCrossCheck_'+str(result[1])+'.jpg">\
        <img src="GlobalFigs/BtCrossCheck_'+str(result[1])+'.jpg" width="100%"></a></br>\
        <a href="GlobalFigs/BtCrossCheck_'+str(result[1])+'.dat">data</a></div>')
    fileid.write('<div class="grafy30"><h2> $\\varphi_{fl}$ x  $\\varphi_{pl}$</h2>\
    <a href="GlobalFigs/BtCrossCheck_FlowPot_PlPot.jpg">\
    <img src="GlobalFigs/BtCrossCheck_FlowPot_PlPot.jpg" width="100%"></a></br></div>')

    fileid.write('<h2>References</h2>')
    fileid.write('<ul>\
    <li><a href="https://www.mathworks.com/matlabcentral/fileexchange/19312-langmuir-probe-data-analysis-code">Aasim Azooz: Langmuir probe data analysis code. Mathworks.</a></li>\
    </ul>')
    fileid.write('</body></html>')
    fileid.close()
    
def MovieGeneration(ShotNo):
    os.system('convert -delay 50 '+str(ShotNo)+'/IndivVAchars/* '+str(ShotNo)+'/finalmovie.gif')
In [4]:
Setup=[\
       [23297,23196,1100,[4,39.2,46,3.7],[6.3,-34.5,21.3,0.02],'','o'],\
       [23295,23194,1200,[4.3,42.3,44.6,5.2],[8.4,-34,21,0.02],'','x'],\
       [23296,23192,1300,[4,36.7,43.7,3.4],[5.5,-32,21.4,0.02],'','*'],\
      ]

GraphSetup=[['Iplasma','$I_{pl}$ [kA]','$I_p$'],['Uloop',"$U_{l}$ [V]",'$U_l$'],['Btoroidal','$B_t$ [T]','$B_t$']]

ResultsDef=[\
            [2,'Time','t [ms]','',' {:2.0f}'],\
            [4,'$B_t$','$B_t$ [T]','$T$',' {:4.2f}'],\
            [5,'$a_1$','$a_1$ []','',' {:4.2f}'],\
            [6,'$a_2$','$a_2$ []','',' {:4.2f}'],\
            [7,'$a_3$','$a_3$ []','',' {:4.2f}'],\
            [8,'$a_4$','$a_4$ []','',' {:4.2f}'],\
            [12,'$I_{sat}$','$I_{sat}$ [mA]','',' {:4.2f}'],\
            [13,'$V_{pl}$','$\\varphi_{pl}$ [V]','',' {:4.2f}'],\
            [14,'$T_e$','$T_e$ [eV]','',' {:4.2f}'],\
            [15,'$C$','$C$ []','',' {:4.3f}'],\
            [16,'$Ratio_{Sat}$','$Ratio_{sat}$ []','',' {:4.2f}'],\
            [17,'$U_{fl}$','$\\varphi_{fl}$ [V]','',' {:4.2f}'],\
            [18,'alpha','$\\alpha$','',' {:4.2f}'],\
            ]


for j in Setup:
    #GetData(j[0],j[1])
    IndividualShot(j,[0,0]) # mode: which analysis to be made, second parameter: 0 for everything 
    #MovieGeneration(j[0])
GlobalFigs() 
Dependences()
HtmlGeneration()

    
Doing 23297
Fits .. ...:39
0,
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:326: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
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,html generation ...
Global graph generation .. 
Doing 23295
Fits .. ...:39
0,1,2,3,4,5,6,7,8,9,10,11,12,13,(Classic fit problem)14,15,16,17,(Azooz fit problem)18,19,20,(Classic fit problem)21,22,23,24,(Classic fit problem)25,26,27,28,29,30,31,(Classic fit problem)32,33,34,35,36,37,38,html generation ...
Global graph generation .. 
Doing 23296
Fits .. ...:41
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,(Classic fit problem)23,24,25,26,27,28,29,(Azooz fit problem)30,(Classic fit problem)31,32,33,34,35,36,37,38,39,40,
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:133: RuntimeWarning: overflow encountered in exp
html generation ...
Global graph generation .. 
html generation ...
In [ ]: