+ return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None, qstd
+
+ def efjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293.0, return_errors=False):
+ '''
+ Extended Freely-jointed chain function
+ ref: F Oesterhelt, M Rief and H E Gaub, New Journal of Physics 1 (1999) 6.1–6.11
+ Please note that 2-parameter fitting (length and kl) usually does not converge, use fixed kl
+ '''
+ '''
+ clicked_points[0] is the contact point (calculated or hand-clicked)
+ clicked_points[1] and [2] are edges of chunk
+
+ '''
+ #Fixed parameters from reference
+ Kb=(1.38065e-2) #in pN.nm
+ Lp=0.358 #planar, nm
+ Lh=0.280 #helical, nm
+ Ks=150e3 #pN/nm
+
+
+ #STEP 1: Prepare the vectors to apply the fit.
+
+ #indexes of the selected chunk
+ first_index=min(clicked_points[1].index, clicked_points[2].index)
+ last_index=max(clicked_points[1].index, clicked_points[2].index)
+
+ #getting the chunk and reverting it
+ xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
+ xchunk.reverse()
+ ychunk.reverse()
+ #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
+ xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
+ ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
+
+
+ #make them arrays
+ xchunk_corr_up=scipy.array(xchunk_corr_up)
+ ychunk_corr_up=scipy.array(ychunk_corr_up)
+
+ xchunk_corr_up_nm=xchunk_corr_up*1e9
+ ychunk_corr_up_pn=ychunk_corr_up*1e12
+
+
+ #STEP 2: actually do the fit
+
+ #Find furthest point of chunk and add it a bit; the fit must converge
+ #from an excess!
+ xchunk_high=max(xchunk_corr_up_nm)
+ xchunk_high+=(xchunk_high/10.0)
+
+ #Here are the linearized start parameters for the WLC.
+ #[Ns , pii=1/P]
+ #max number of monomers (all helical)for a contour length xchunk_high
+ excessNs=xchunk_high/(Lp)
+ p0=[excessNs,(1.0/(0.7))]
+ p0_plfix=[(excessNs)]
+
+ def dist(px,py,linex,liney):
+ distancesx=scipy.array([(px-x)**2 for x in linex])
+ minindex=np.argmin(distancesx)
+ return (py-liney[minindex])**2
+
+ def deltaG(f):
+ dG0=12.4242 #3kt in pN.nm
+ dL=0.078 #planar-helical
+ return dG0-f*dL
+
+ def Lfactor(f,T=T):
+ Lp=0.358 #planar, nm
+ Lh=0.280 #helical, nm
+ Kb=(1.38065e-2)
+ therm=Kb*T
+ dG=deltaG(f)
+
+ return Lp/(np.exp(dG/therm)+1)+Lh/(np.exp(-dG/therm)+1)
+
+ def coth(z):
+ '''
+ hyperbolic cotangent
+ '''
+ return 1.0/np.tanh(z)
+
+ def x_efjc(params,f,T=T,Ks=Ks):
+ '''
+ efjc function for ODR fitting
+ '''
+
+ Ns=params[0]
+ invkl=params[1]
+ Kb=(1.38065e-2)
+ therm=Kb*T
+ beta=(f/therm)/invkl
+
+ x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
+ return x
+
+ def x_efjc_plfix(params,f,kl_value=pl_value,T=T,Ks=Ks):
+ '''
+ efjc function for ODR fitting
+ '''
+
+ Ns=params
+ invkl=1.0/kl_value
+ Kb=(1.38065e-2)
+ therm=Kb*T
+ beta=(f/therm)/invkl
+
+ x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
+ return x
+
+ #make the ODR fit
+ realdata=scipy.odr.RealData(ychunk_corr_up_pn,xchunk_corr_up_nm)
+ if pl_value:
+ model=scipy.odr.Model(x_efjc_plfix)
+ o = scipy.odr.ODR(realdata, model, p0_plfix)
+ else:
+ print 'WARNING eFJC fit with free pl sometimes does not converge'
+ model=scipy.odr.Model(x_efjc)
+ o = scipy.odr.ODR(realdata, model, p0)
+
+ o.set_job(fit_type=2)
+ out=o.run()
+
+
+ Ns=out.beta[0]
+ Lc=Ns*Lp*1e-9
+ if len(out.beta)>1:
+ kfit=1e-9/out.beta[1]
+ kfitnm=1/out.beta[1]
+ else:
+ kfit=1e-9*pl_value
+ kfitnm=pl_value
+
+ fit_out=[Lc, kfit]
+
+ #Calculate fit errors from output standard deviations.
+ fit_errors=[]
+ fit_errors.append(out.sd_beta[0]*Lp*1e-9)
+ if len(out.beta)>1:
+ fit_errors.append(1e9*out.sd_beta[1]*kfit**2)
+
+
+
+ def efjc_eval(y,params,pl_value,T=T,Lfactor=Lfactor,Ks=Ks):
+ '''
+ Evaluates the eFJC function
+ '''
+ if not pl_value:
+ Ns, invkl = params
+ else:
+ Ns = params
+
+ if pl_value:
+ invkl=1.0/pl_value
+
+ Kb=(1.38065e-2) #boltzmann constant
+ therm=Kb*T #so we have thermal energy
+ beta=(y/therm)/invkl