In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as spo
from scipy.signal import argrelextrema
import os


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

#if is_interactive():
#    %matplotlib inline

#BasePath='/GW/Experiments/EdgePlasmaPhysics/ParticleFlux/BallPenProbe/Experiments/BottomPosition/Sessions/240217BscanWithStabilization'
BasePath='/home/svoboda/mnt/buon/Aktual/BPPs/020317ReproducibleShots@1300Bt'

def AzoozFn(Vbias, a1, a2, a3, a4):
    return np.exp(a1 * np.tanh((Vbias+a2)/a3)) - a4

def Classical3Fn(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))

def mkdir(dir):  
    try:os.makedirs(dir) 
    except OSError:pass
    
def FigDef():plt.figure(figsize=(10, 8), dpi= 80, facecolor='w', edgecolor='k')    
In [2]:
def MakeFits(DataToBeFitted, AzoozInitialGuess, Classic4InitialGuess):
   
     
    AzoozFitParms, po_cov = spo.curve_fit(AzoozFn, DataToBeFitted[:,0],DataToBeFitted[:,1], AzoozInitialGuess)
        
    return AzoozFitParms
    
In [3]:
def Go(ShotNo,Time,AzoozInitialGuess, Classic4InitialGuess):
    FigDef()
    OutputDir=str(ShotNo)+'/'+Time;mkdir(OutputDir)
    Data=np.loadtxt(BasePath+'/'+str(ShotNo)+'/TimeLine/'+Time+'/SmoothedVachar.dat')
    RawData=np.loadtxt(BasePath+'/'+str(ShotNo)+'/TimeLine/'+Time+'/RawVachar.dat')
    Voltage=Data[:,0];Current=Data[:,1]
    LinVoltageSpace = np.linspace(int(np.min(Voltage)),int(np.max(Voltage)),len(Voltage)) # for diferrential analysis
    SmoothedUfl = np.where(abs(Current-0.0)==abs(Current-0.0).min())[0][0]
    Shift=int(float(len(LinVoltageSpace))/(np.max(LinVoltageSpace)-np.min(LinVoltageSpace))*VoltageShift)  
    plt.plot(RawData[:,0],RawData[:,1],'+')
    plt.plot(Voltage,Current)
    plt.axhline(0)
    plt.axvline(Voltage[SmoothedUfl])
    plt.axvline(Voltage[SmoothedUfl+Shift])
    plt.savefig(OutputDir+'/RawVAchar.jpg', bbox_inches='tight')
    plt.close()
    #1. Azooz Fit
    try:
        AzoozFitParms, po_cov = spo.curve_fit(AzoozFn, Voltage,Current, AzoozInitialGuess)
        #AzoozFitParms=MakeFits(Data,AzoozInitialGuess, Classic4InitialGuess)
        AzoozFit=AzoozFn(LinVoltageSpace,AzoozFitParms[0],AzoozFitParms[1],AzoozFitParms[2],AzoozFitParms[3])
        FigDef()
        plt.plot(Voltage,Current)
        plt.plot(Voltage,AzoozFit)
        plt.plot(RawData[:,0],RawData[:,1],',')
        plt.axhline(0)
        np.savetxt(OutputDir+'/AzoozParm.dat',[AzoozInitialGuess,AzoozFitParms],delimiter=" ",fmt="%s")
        plt.savefig(OutputDir+'/AzoozFit.jpg', bbox_inches='tight')
        plt.close()
        #Let's find Azooz fnction inflex point:

        FigDef()
        dVoltage = LinVoltageSpace[1]-LinVoltageSpace[0]
        dAzoozFit_dVoltage = np.gradient(AzoozFit, dVoltage)
        f = plt.figure()
        f,ax1 = plt.subplots(sharex=True);
        plt.axhline(0)
        ax1.plot(Voltage,AzoozFit, 'b-')
        ax2 = ax1.twinx()
        ax2.plot(LinVoltageSpace,dAzoozFit_dVoltage, 'r-')
        # Where is the inflex point??
        InflPointLoc=argrelextrema(dAzoozFit_dVoltage, np.greater)[0]
        plt.axvline(LinVoltageSpace[InflPointLoc])
        # Where the fit intersect the x axis??
        AzoozUfl = np.where(abs(AzoozFit-0.0)==abs(AzoozFit-0.0).min())[0][0]
        plt.axvline(Voltage[AzoozUfl]) # I do not understand why not LinVoltageSpace ...
        # Let us shift it 1-2 Te ... let's try 20 V (VoltageShift)
        

        plt.savefig(OutputDir+'/AzoozFitWfirstDerv.jpg', bbox_inches='tight')
        plt.close()
    except IOError:
        print('Azooz problem')
    
    print(AzoozUfl)
    
    DataIndex=SmoothedUfl+Shift    
    # Now we have to reduce the data set for the Classsic fit:
    print(SmoothedUfl)
        
    VoltageReduced=Voltage[0:DataIndex]
    CurrentReduced=Current[0:DataIndex]
    RawCurrentReduced=RawData[:,1][0:DataIndex]
    RawVoltageReduced=RawData[:,0][0:DataIndex]
    
    # Classic 3 par fit:
    FigDef()
    Classic3FitParms, po_cov = spo.curve_fit(Classical3Fn,VoltageReduced, CurrentReduced, [Classic4InitialGuess[0],Classic4InitialGuess[1],Classic4InitialGuess[2]])
    Classic3Fit=Classical3Fn(VoltageReduced,Classic3FitParms[0],Classic3FitParms[1],Classic3FitParms[2])
    plt.plot(VoltageReduced,Classic3Fit)
    
    
    plt.title('Classic 3 fit @'+str(ShotNo)) 
    plt.plot(VoltageReduced,CurrentReduced)
    plt.plot(RawVoltageReduced,RawCurrentReduced,',')
    plt.axhline(0);plt.axvline(Voltage[DataIndex]);plt.axvline(Voltage[SmoothedUfl])
    
    plt.savefig(OutputDir+'/Classic3Fit.jpg', bbox_inches='tight')
    plt.close()
    
    # Classic 3 par fit:
    FigDef()
    #Classic4FitParms, po_cov = spo.curve_fit(Classical4Fn,VoltageReduced, CurrentReduced, [Classic4InitialGuess[0],Classic4InitialGuess[1],Classic4InitialGuess[2],Classic4InitialGuess[3]])
    Classic4FitParms, po_cov = spo.curve_fit(Classical4Fn,VoltageReduced, CurrentReduced,bounds=((0,-100,0,0),(10,0,30,0.03)))
    Classic4Fit=Classical4Fn(VoltageReduced,Classic4FitParms[0],Classic4FitParms[1],Classic4FitParms[2],Classic4FitParms[3])
    
    plt.plot(VoltageReduced,Classic4Fit)
    
    plt.title('Classic 4 fit @'+str(ShotNo)) 
    plt.plot(VoltageReduced,CurrentReduced)
    plt.plot(RawVoltageReduced,RawCurrentReduced,',')
    plt.axhline(0);plt.axvline(Voltage[DataIndex]);plt.axvline(Voltage[SmoothedUfl])
    
    plt.savefig(OutputDir+'/Classic4Fit.jpg', bbox_inches='tight')
    plt.close()
    
    # ... and html ...
    np.savetxt(OutputDir+'/index.html',['\
                       <html><style>\
                        div.grafy30 {display:inline-block;width: 30%;}\
                        .grafy30 img {width: 100%;}\
                        </style>\<center>\
                       <div class="grafy30">\
                       <a href="RawVAchar.jpg"><img src="RawVAchar.jpg"  width="30%"/></a></div>\
                       <div class="grafy30">\
                       <a href="AzoozFit.jpg"><img src="AzoozFit.jpg"  width="30%"/></a>\
                       Guess:'+str(AzoozInitialGuess)+'<br/>\
                       Fit:'+str(AzoozFitParms)+'<br/>\
                       </div>\
                       <div class="grafy30">\
                       <a href="AzoozFitWfirstDerv.jpg"><img src="AzoozFitWfirstDerv.jpg"  width="30%"/></a>\
                       InflPointLoc:'+str(InflPointLoc[0])+'<br/>\
                       LinVoltageSpace[InflPointLoc]:'+'{:4.2f} V'.format(LinVoltageSpace[InflPointLoc[0]])+'<br/>\
                       AzoozUfl:'+'{:4.2f} V'.format(Voltage[AzoozUfl])+'<br/>\
                       DataIndex'+'({:4.2f} V):'.format(VoltageShift)+'{:4.2f} V'.format(Voltage[DataIndex])+'\
                       </div>\
                       <div class="grafy30">\
                       <a href="Classic3Fit.jpg"><img src="Classic3Fit.jpg"  width="30%"/></a>\
                       Guess:'+str(Classic4InitialGuess)+'<br/>\
                       Fit:'+str(Classic3FitParms)+'<br/>'+\
                       '$I_{sat}=$'+'{:4.2f}'.format(Classic3FitParms[0])+\
                       '</div>\
                       <div class="grafy30">\
                       <a href="Classic4Fit.jpg"><img src="Classic4Fit.jpg"  width="30%"/></a>\
                       Guess:'+str(Classic4InitialGuess)+'<br/>\
                       Fit:'+str(Classic4FitParms)+'<br/>'+\
                       '$I_{sat}=$'+'{:4.2f}'.format(Classic4FitParms[0])+\
                       '</div>\
                       '],fmt="%s")
            
    
    #f.tight_layout()
    
    #plt.show()
In [6]:
VoltageShift=20 # V
#Go(23297,'20711_21213',[4,39.2,46,3.7],[6.3,-34.5,21.3,0.02])
#Go(23295,'14677_15171',[4,39.2,46,3.7],[6.3,-34.5,21.3,0.02])
#Go(23447,'7209_8208',[4,39.2,46,3.7],[6.3,-34.5,21.3,0.02])
#Go(23447,'8208_9220',[4,39.2,46,3.7],[6.3,-34.5,21.3,0.02])
#Go(23447,'9220_10221',[4,39.2,46,3.7],[3.02,-65,20,0.02])
Go(23447,'10221_11234',[4,39.2,46,3.7],[10.5,-63,102,0.02])
0
264
In [5]:
import scipy
scipy.__version__
Out[5]:
'0.18.1'
In [ ]: