Private
Public Access
1
0
This commit is contained in:
Sander Roosendaal
2017-06-17 20:27:18 +02:00
parent 87222b88c4
commit b03d88ba64
3 changed files with 66 additions and 38 deletions

View File

@@ -1,6 +1,60 @@
import pandas as pd
import numpy as np
from scipy.interpolate import griddata
from scipy import optimize
def cpfit(powerdf):
# Fit the data to thee parameter CP model
fitfunc = lambda pars,x: abs(pars[0])/(1+(x/abs(pars[2]))) + abs(pars[1])/(1+(x/abs(pars[3])))
errfunc = lambda pars,x,y: fitfunc(pars,x)-y
p0 = [500,350,10,8000]
p1 = p0
thesecs = powerdf['Delta']
theavpower = powerdf['CP']
if len(thesecs)>=4:
p1, success = optimize.leastsq(errfunc, p0[:], args = (thesecs,theavpower))
else:
factor = fitfunc(p0,thesecs.mean())/theavpower.mean()
p1 = [p0[0]/factor,p0[1]/factor,p0[2],p0[3]]
p1 = [abs(p) for p in p1]
fitt = pd.Series(10**(4*np.arange(100)/100.))
fitpower = fitfunc(p1,fitt)
fitpoints = fitfunc(p1,thesecs)
fitpoints0 = fitpoints.copy()
dd = fitpoints-theavpower
ddmin = dd.min()
frac = abs(ddmin)/fitpoints.mean()
while frac>0.0001:
fitpoints = fitpoints*(fitpoints.mean()-ddmin)/(fitpoints.mean())
dd = fitpoints-theavpower
ddmin = dd.min()
frac = abs(ddmin)/fitpoints.mean()
print ddmin,frac
ratio = fitpoints.mean()/fitpoints0.mean()
return p1,fitt,fitpower,ratio
def getlogarr(maxt):
maxlog10 = np.log10(maxt-5)
logarr = np.log10(5.)+np.arange(50)*maxlog10/50.
logarr = [int(10.**(la)) for la in logarr]
logarr = pd.Series(logarr)
logarr.drop_duplicates(keep='first',inplace=True)
logarr = logarr.values
return logarr
def getsinglecp(df):
thesecs = df['TimeStamp (sec)'].max()-df['TimeStamp (sec)'].min()
@@ -9,13 +63,7 @@ def getsinglecp(df):
else:
maxt = 1000.
maxlog10 = np.log10(maxt)
logarr = np.arange(50)*maxlog10/50.
logarr = [int(10.**(la)) for la in logarr]
logarr = pd.Series(logarr)
logarr.drop_duplicates(keep='first',inplace=True)
logarr = logarr.values
logarr = getlogarr(maxt)
dfnew = pd.DataFrame({
@@ -42,7 +90,7 @@ def getcp(dfgrouped,logarr):
tmax = tt.max()
if tmax > 500000:
newlen = int(tmax/5000.)
newlen = int(tmax/2000.)
else:
newlen = len(tt)
if newlen < len(tt):