def do_mpfit(x,y,covarin,guess,functname=thepolynomial):
    # check if covariance or error bars were given
    covar=covarin
    if np.size(np.shape(covarin)) == 1:
        err=covarin
        print('err')
        #Prepare arguments for mpfit
        fa={'x':double(x),'y':double(y),'err':double(err),'functname':functname}
        costfct=fdeviates
    else:
        print('covar')
        #Fist do a SVD decomposition
        u,s,v=np.linalg.svd(covarin)
        #Prepare arguments for mpfit
        fa={'x':double(x),'y':double(y),'svdvals':double(s),'v':double(v),'functname':functname}
        costfct=chi2svd

    #Run MPFIT
    mpf = mpfit.mpfit(costfct, guess, functkw=fa, autoderivative=1)
    print('Status of the Fit',mpf.status)
    print('Chi2=',mpf.fnorm)
    print('ndf=',mpf.dof)
    print('Fitted params:',mpf.params)

    return(mpf,mpf.params,mpf.perror,mpf.covar)

def test_quadfit():

    x = N.array([-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
          -2.8394297E-01,1.3416969E+00,1.3757038E+00,
          -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
          8.2065094E-01])
    y = N.array([2.3095947E+01,2.6449392E+01,1.0204468E+01,
          5.40507,1.5787588E+01,1.6520903E+01,
          1.5971818E+01,4.7668524E+00,4.9337711E+00,
          8.7348375E+00])
    ey=0.2*N.ones(y.shape,dtype='float64')
    p0=N.array([1.0,1.0,1.0],dtype='float64')  #initial conditions
    pactual=N.array([4.7, 0.0, 6.2]) #actual values used to make data
    parbase={'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}
    parinfo=[]
    for i in range(len(pactual)):
        parinfo.append(copy.deepcopy(parbase))
        parinfo[i]['value'] = p0[i]
    fa = {'x': x, 'y': y, 'err': ey}
    m = mpfit(myfunctquad, p0, parinfo = parinfo, functkw = fa)
    if (m.status <= 0):
        print 'status = ', m.status
    print m.perror
    assert N.allclose(m.params, N.array([ 4.703829, 0.062586,  6.163087], dtype='float64'))
    assert N.allclose(m.perror, N.array([ 0.097512, 0.054802,  0.054433], dtype='float64'))
    chisq = (myfunctquad(m.params, x=x, y=y, err=ey)[1]**2).sum()
    assert N.allclose(N.array([chisq], dtype='float64'), N.array([5.679323], dtype='float64'))
    assert m.dof == 7
    return

def test_gaussfit_fixed():

    x = N.array([-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
        -2.8394297E-01,1.3416969E+00,1.3757038E+00,
        -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
        8.2065094E-01])
    y = N.array([-4.4494256E-02,8.7324673E-01,7.4443483E-01,
        4.7631559E+00,1.7187297E-01,1.1639182E-01,
        1.5646480E+00,5.2322268E+00,4.2543168E+00,
        6.2792623E-01])
    ey=0.5*N.ones(y.shape,dtype='float64')
    p0=N.array([0.0, 1.0, 0.0, 0.1],dtype='float64')  #initial conditions          
    pactual=N.array([0.0, 4.70, 0.0, 0.5]) #actual values used to make data
    parbase={'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}
    parinfo=[]
    for i in range(len(pactual)):
        parinfo.append(copy.deepcopy(parbase))
        parinfo[i]['value'] = p0[i]
    
    parinfo[0]['fixed'] = 1
    parinfo[2]['fixed'] = 1

    fa = {'x': x, 'y': y, 'err': ey}
    m = mpfit(myfunctgauss, p0, parinfo = parinfo, functkw = fa)
    if (m.status <= 0):
        print 'status = ', m.status
    assert N.allclose(m.params, N.array([ 0., 5.059244,  0.,  0.479746 ], dtype='float64'))
    assert N.allclose(m.perror, N.array([ 0., 0.329307,  0.,  0.053804 ], dtype='float64'))
    chisq = (myfunctgauss(m.params, x=x, y=y, err=ey)[1]**2).sum()
    assert N.allclose(N.array([chisq], dtype='float64'), N.array([15.516134], dtype='float64'))
    assert m.dof == 8
    return

def zandtfit(Inclination, Planetaryradius_SI, Keplerlc, Teff, Logg, Metalicity, sigma, Orbitalperiod_SI, KeplerID, Planetnumberer):
	Planetaryradius_Earth = (Planetaryradius_SI) /  6.3675E6 #convert the planetary radius to earth radii
	p0 = [math.cos(1.570795), Planetaryradius_Earth, Teff, Logg, np.log10(Metalicity)] #inital guess (from kepler mast and headers)
	fa = {'Keplerlc':Keplerlc, 'sigma':sigma, 'Orbitalperiod_SI':Orbitalperiod_SI, 'KeplerID':KeplerID, 'Planetnumberer':Planetnumberer} #additional variables to be past to the function
	parinfo = [{'value':p0[0], 'fixed':0, 'limited':[1,1], 'limits':[-1,1], 'step':0.01}, {'value':p0[1], 'fixed':0, 'limited':[1,0], 'limits':[0.01,100], 'step':0},\
		{'value':p0[2], 'fixed':0, 'limited':[1,1], 'limits':[0,40000], 'step':0}, {'value':p0[3], 'fixed':0, 'limited':[1,1], 'limits':[0,5.0], 'step':0},\
		{'value':p0[4], 'fixed':0, 'limited':[1,1], 'limits':[-5,1], 'step':0} ]
	m = mpfit(minimisefunction, p0, functkw = fa, quiet = 0, maxiter = 35, parinfo = parinfo, ftol=0.0001) #run the mpfit for the best fit
	#print m #testing
	bestfitarray = m.params #extract the results 
	errors = m.perror #these need to be properly stored
	covar = m.covar #these need to be properly stored
	if errors == None: #i.e. if mpfits is a pain and decides to not return the errors or covariance properly
		errors = np.zeros(5) #create an empty errors array
	if covar == None: #if no covariance array is produced, create an array of zeros so that there are no problems later
		covar = np.zeros((5,5))
	#print bestfitarray
	Inclination_fit = math.acos(bestfitarray[0]) #Extract the fit inclination
	Planetaryradius_SI_fit =   bestfitarray[1] * 6.3675E6 #Extract the fit planetary radius and convert it to SI units
	Teff_fit = bestfitarray[2] #Extract the fit Teff
	Logg_fit = bestfitarray[3] #Extract the fit logg of the surface gravity
	Metalicity_fit = (10**bestfitarray[4]) #Extract the metalicity and remove the logg
	Mstar_SI, Rstar_SI, Stellardensity_SI = RMrhocalc(Teff_fit, Logg_fit, Metalicity_fit) #Calculate the stellar mass, radius and density
	semimajoraxis_SI = semimajoraxiscalc(Orbitalperiod_SI, Mstar_SI) #calculate the semi-major axis
	T_dur_SI = Transitdurationcalc(Rstar_SI, Planetaryradius_SI_fit, semimajoraxis_SI, Inclination_fit, Orbitalperiod_SI) #calculate the transit duration
	return Inclination_fit, Planetaryradius_SI_fit, Teff_fit, Logg_fit, Metalicity_fit, Mstar_SI, Rstar_SI, Stellardensity_SI, semimajoraxis_SI, T_dur_SI, errors, covar

    def _run_mpfit(self):
        """Run an MPFIT regression"""
        # Make up the PARINFO list

        parinfo = []
        for i in range(len(self._guess)):
            pdict = {
                "value": self._guess[i].copy(),
                "fixed": 0,
                "limited": self._ilimited[i],
                "limits": self._ilimits[i],
            }
            parinfo.append(pdict)

        quiet = 1
        if self.debug:
            quiet = 0

        m = mpfit.mpfit(self._residualsMPFIT, self._guess, parinfo=parinfo, quiet=quiet, debug=self.debug)

        self._result = m.params
        self._niter = m.niter

        if m.perror is not None:
            self._stdev = m.perror
        else:
            self._stdev = zeros(len(self._guess))

        if m.covar is not None:
            self._covar = m.covar
        else:
            self._covar = zeros((len(self._guess), len(self._guess)))

        self._mpfit = m

def newfit2( infile, outfile, xcol, ycol, p0, \
    constrain=[ {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0} ] ) :
  [x,y] = getdata( infile, xcol, ycol )
  err = numpy.ones( len(x), dtype=float )
  parinfo = [ {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0}, {'fixed':0} ]
  fa = {'x':x, 'y':y, 'err':err}
  m = mpfit.mpfit(myfunc, p0, functkw=fa, parinfo=constrain)

  fout = open( outfile, "w" )
  fout.write("# %s\n" %  outfile )
  fout.write("# fit to spectrum %s, xcol=%d, ycol=%d\n" % ( infile, xcol, ycol ) )
  fout.write("# final fit:\n")
  fout.write("#   yavg= %0.5f   %d \n" %  (m.params[0], constrain[0]['fixed'])  )
  fout.write("#   slope = %0.5f   %d\n" % (m.params[1], constrain[1]['fixed'])  )
  fout.write("#   amp = %0.5f   %d\n" %  (m.params[2], constrain[2]['fixed'])  )
  fout.write("#   v0 = %0.3f   %d\n" %  (m.params[3], constrain[3]['fixed'])  )
  fout.write("#   FWHM = %0.3f   %d\n" % (m.params[4], constrain[4]['fixed'])  )
  fout.write("#\n")

  xavg = 0.5 * (x[0] + x[len(x)-1])
  fit = eval('m.params[0] + m.params[1]*(x-xavg) + \
     m.params[2]*numpy.exp(-2.772 * pow( (x-m.params[3])/m.params[4], 2.))')
     # fit to raw data, including slope
  yminusSlope = eval('y - m.params[1]*(x-xavg)')
     # raw data with linear slope removed
  gfit = eval('m.params[0] + m.params[2]*numpy.exp(-2.772 * pow( (x-m.params[3])/m.params[4], 2.))')
     # gaussian fit to data with slope removed
  yminusSlopeNorm = eval('yminusSlope/m.params[0]')
     # normalized absorption 
  gfitNorm = eval('gfit/m.params[0]')
  for i in range(0,len(x)) :
    fout.write(" %10.3f  %9.4f  %9.4f  %9.4f  %9.4f %9.5f %9.5f\n" % \
	  (x[i],y[i],fit[i],yminusSlope[i],gfit[i],yminusSlopeNorm[i],gfitNorm[i] ) )
  fout.close()

def onedgaborfit(xax, data, err=None,
        params=[0,1,0,1,1,0],fixed=[False,False,False,False,False,False],
        limitedmin=[False,False,False,True,True,True],
        limitedmax=[False,False,False,False,False,True], minpars=[0,0,0,0,0,0],
        maxpars=[0,0,0,0,0,360], quiet=True, shh=True,
        veryverbose=False):
    """
    Inputs:
       xax - x axis
       data - y axis
       err - error corresponding to data

       params - Fit parameters: Height of background, Amplitude, Shift, Width, Wavelength, Phase
       fixed - Is parameter fixed?
       limitedmin/minpars - set lower limits on each parameter (default: width>0)
       limitedmax/maxpars - set upper limits on each parameter
       quiet - should MPFIT output each iteration?
       shh - output final parameters?

    Returns:
       Fit parameters
       Model
       Fit errors
       chi2
    """

    def mpfitfun(x,y,err):
        if err is None:
            def f(p,fjac=None): return [0,(y-onedgabor(x,*p))]
        else:
            def f(p,fjac=None): return [0,(y-onedgabor(x,*p))/err]
        return f

    if xax == None:
        xax = np.arange(len(data))

    parinfo = [ {'n':0,'value':params[0],'limits':[minpars[0],maxpars[0]],'limited':[limitedmin[0],limitedmax[0]],'fixed':fixed[0],'parname':"HEIGHT",'error':0} ,
                {'n':1,'value':params[1],'limits':[minpars[1],maxpars[1]],'limited':[limitedmin[1],limitedmax[1]],'fixed':fixed[1],'parname':"AMPLITUDE",'error':0},
                {'n':2,'value':params[2],'limits':[minpars[2],maxpars[2]],'limited':[limitedmin[2],limitedmax[2]],'fixed':fixed[2],'parname':"SHIFT",'error':0},
                {'n':3,'value':params[3],'limits':[minpars[3],maxpars[3]],'limited':[limitedmin[3],limitedmax[3]],'fixed':fixed[3],'parname':"WIDTH",'error':0},
                {'n':4,'value':params[4],'limits':[minpars[4],maxpars[4]],'limited':[limitedmin[4],limitedmax[4]],'fixed':fixed[4],'parname':"WAVELENGTH",'error':0},
                {'n':5,'value':params[5],'limits':[minpars[5],maxpars[5]],'limited':[limitedmin[5],limitedmax[5]],'fixed':fixed[5],'parname':"PHASE",'error':0}]

    mp = mpfit(mpfitfun(xax,data,err),parinfo=parinfo,quiet=quiet)
    mpp = mp.params
    mpperr = mp.perror
    chi2 = mp.fnorm

    if mp.status == 0:
        raise Exception(mp.errmsg)

    if (not shh) or veryverbose:
        print "Fit status: ",mp.status
        for i,p in enumerate(mpp):
            parinfo[i]['value'] = p
            print parinfo[i]['parname'],p," +/- ",mpperr[i]
        print "Chi2: ",mp.fnorm," Reduced Chi2: ",mp.fnorm/len(data)," DOF:",len(data)-len(mpp)

    return mpp,onedgabor(xax,*mpp),mpperr,chi2

def ab_fit(src,xd,td,lguess,xmin_id,xmax_id,evmin1,evmax1,evmin2,evmax2):
	npar    =  len(lguess)
	guessp  = np.array(lguess, dtype='float64')
	plimd   = [[False,False]]*npar
	plims   = [[0.,0.]]*npar
	parbase = {'value': 0., 'fixed': 0, 'parname': '', 'limited': [0, 0], 'limits': [0., 0.]}
	pname   = ['bg']*npar
	pfix    = [False]*npar

	count   = 1
	for i in range(npar):
		if(i==0):
			pname[i] = 'Tcont'

		if( ((i-1)%3 == 0) and (i>0) ):
			pname[i] = 'tau'+str(count)

		if( ((i-2)%3 == 0) and (i>0) ):
			pname[i] = 'v0'+str(count)

		if( ((i-3)%3 == 0) and (i>0) ):
			pname[i] = 'wid'+str(count)
			count    = count + 1

	parinfo = []
	for i in range(len(guessp)):
		parinfo.append(copy.deepcopy(parbase))

	for i in range(len(guessp)):
		parinfo[i]['value']   = guessp[i]
		parinfo[i]['fixed']   = pfix[i]
		parinfo[i]['parname'] = pname[i]
		parinfo[i]['limited'] = plimd[i]

	##  data and 1-sigma uncertainty ##
	x    = xd[xmin_id:xmax_id]
	y    = td[xmin_id:xmax_id]
	sg,\
	erry = get_tb_sigma(xd, td, evmin1, evmax1, evmin2, evmax2)
	erry = erry[xmin_id:xmax_id]

	x    = x.astype(np.float64)
	y    = y.astype(np.float64)
	erry = erry.astype(np.float64)

	fa   = {'x':x, 'y':y, 'err':erry}
	mp   = mpfit(myabfunc, guessp, parinfo=parinfo, functkw=fa, quiet=True)

	## Fit values of Tau ##
	abp   = mp.params
	abper = mp.perror
	fit   = 0.
	for i in range(1, len(abp), 3):
		fit = fit + abp[i]*np.exp(- ( (x-abp[i+1])/(0.6005612*abp[i+2]))**2)  

	fit = np.exp(-fit)
	# fit = abp[1]*fit

	return x,fit,y/abp[0],abp,abper,npar,parbase,pname,parinfo

def delPuncorr_bestfit(k, delP_over_P): 
    '''
    '''
    fa = {'x': k, 'y': delP_over_P}
    
    param_guess = [-4.0]
    bestfit = mpfit.mpfit(nongauss_poly_mpfit, param_guess, functkw=fa, quiet=True)
    return bestfit.params

def mpfitexpr(func, x, y, err, start_params, check=True, full_output=False,
              imports=None, **kw):
    """Fit the used defined expression to the data
    Input:
    - func: string with the function definition
    - x: x vector
    - y: y vector
    - err: vector with the errors of y
    - start_params: the starting parameters for the fit
    Output:
    - The tuple (params, yfit) with best-fit params and the values of func evaluated at x
    Keywords:
    - check: boolean parameter. If true(default) the function will be checked for sanity
    - full_output: boolean parameter. If True(default is False) then instead of best-fit parameters the mpfit object is returned
    - imports: list of strings, of optional modules to be imported, required to evaluate the function
    Example:
    params,yfit=mpfitexpr('p[0]+p[2]*(x-p[1])',x,y,err,[0,10,1])

    If you need to use numpy and scipy functions in your function, then
            you must to use the full names of these functions, e.g.:
            numpy.sin, numpy.cos etc.

    This function is motivated by mpfitexpr() from wonderful MPFIT IDL package
            written by Craig Markwardt

    """

    hash = {}
    hash['numpy'] = numpy
    hash['scipy'] = scipy

    if imports is not None:
        for i in imports:
            #exec '%s=__import__("%s")'%(a,b) in globals(),locals()
            hash[i] = __import__(i)

    def myfunc(p, fjac=None, x=None, y=None, err=None):
        return [0, eval('(y-(%s))/err' % func, hash, locals())]

    myre = "[^a-zA-Z]p\[(\d+)\]"
    r = re.compile(myre)
    maxp = -1
    for m in re.finditer(r, func):
        curp = int(m.group(1))
        maxp = curp if curp > maxp else maxp
    if check:
        if maxp == -1:
            raise Exception("wrong function format")
        if maxp + 1 != len(start_params):
            raise Exception(
                "the length of the start_params != the length of the parameter verctor of the function")
    fa = {'x': x, 'y': y, 'err': err}
    res = mpfit.mpfit(myfunc, start_params, functkw=fa, **kw)
    yfit = eval(func, globals(), {'x': x, 'p': res.params})
    if full_output:
        return (res, yfit)
    else:
        return (res.params, yfit)

def main():

    # Generate a noisy polynomial
    #     [off,  x1,  x2,  x3,  y1,  y2,  y3]
    pIn = [2.0, 1.5, 0.1, 0.3, 1.0, 2.0, 0.05]
    pIn = [1.0, 0.2, 0.0, 0.0, 0.1, 0.0, 0.0]
    shape = (200, 200)
    X, Y, Z, xyData = genpolydata(pIn, shape, 300, 10.2)

    print X.shape
    print Y.shape
    print Z.shape
    print xyData.shape
    # raw_input()

    # Define an function to evaluate the residual
    def errFn(p, fjac=None):
        status = 0
        # poly_surface' returns the 'rfunc' function and the X,Y data is
        # inserted via argument unpacking.
        return status, poly_surface(p)(*[Y, X]) - Z
    
    # Fit the data starting from an initial guess
    mp = mpfit(errFn, parinfo=inParms, quiet=False)
    print 
    for i in range(len(inParms)):
        print "%s = %f +/- %f" % (inParms[i]['parname'],
                                  mp.params[i],
                                  mp.perror[i])
    p1 = mp.params

    #-------------------------------------------------------------------------#
  
    # Plot the original, fit & residual
    fig = pl.figure(figsize=(18,4.3))
    
    ax1 = fig.add_subplot(1,3,1)
    cax1 = ax1.imshow(xyData, origin='lower',cmap=mpl.cm.jet)
    cbar1=fig.colorbar(cax1, pad=0.0)
    ax1.scatter(X, Y, c=Z, s=40, cmap=mpl.cm.jet)
    ax1.set_title("Sampled Data")
    ax1.set_xlim(0, shape[-1]-1)
    ax1.set_ylim(0, shape[-2]-1)
    ax1.set_aspect('equal')
    
    ax2 = fig.add_subplot(1,3,2)
    xyDataFit = poly_surface(p1, shape)
    cax2 = ax2.imshow(xyDataFit, origin='lower', cmap=mpl.cm.jet)
    cbar2=fig.colorbar(cax2, pad=0.0)
    ax2.set_title("Model Fit")
    
    ax3 = fig.add_subplot(1,3,3)
    xyDataRes = xyData - xyDataFit
    cax3 = ax3.imshow(xyDataRes, origin='lower', cmap=mpl.cm.jet)
    cbar2=fig.colorbar(cax3, pad=0.0)
    ax3.set_title("Residual")
    
    pl.show()

def disk_bulge_AIM(Iobs,err,psf,guessP,quiet=0,ftol=1.e-10,return_image=False):

	#Dictionary to pass parameters to function `handleFunc`
	fa = {'y':Iobs,'x':psf,'err':err}

	#Possible constraints for each parameter
	parinfo= [{'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]} for i in range(len(guessP))]
	#Pass initial values to parinfo
	for i in range(len(guessP)):
		parinfo[i]['value'] = guessP[i]
	
	L = guessP[0]

	#Set parameter limits and whether fixed
	parinfo[0]['fixed'] = 1				#xmax
	parinfo[1]['fixed'] = 1				#N
	
	parinfo[2]['limited'] = [1,1]		#disk flux
	parinfo[2]['limits'] = [0.1,1.e10] #
	parinfo[3]['limited'] = [1,1]		#bulge flux
	parinfo[3]['limits'] = [0.1,1.e10] #
	parinfo[4]['limited'] = [1,1]		#q (q=10 corresponds to e = 0.82)
	parinfo[4]['limits'] = [1.0,20.0]
	parinfo[5]['limited'] = [1,1]		#phi (there is 2pi degeneracy, but keep to avoid boundaries)
	parinfo[5]['limits'] = [0.0,6.283185]
	parinfo[6]['limited'] = [1,1]		#disk re
	parinfo[6]['limits'] = [0.05,3.0*L] #
	parinfo[7]['limited'] = [1,1]		#bulge re
	parinfo[7]['limits'] = [0.05,3.0*L] #
	parinfo[8]['limited'] = [1,1]
	parinfo[8]['limits'] = [0.5,4.0]  #n_b
	

	parinfo[9]['fixed'] = 1				#kappa
	parinfo[10]['fixed'] = 1				#Shear 1
	parinfo[11]['fixed'] = 1				#Shear 2
	parinfo[12]['limited'] = [1,1]		#f1
	parinfo[12]['limits'] = [-1.0,1.0]
	parinfo[13]['limited'] = [1,1]		#f2
	parinfo[13]['limits'] = [-1.0,1.0]	
	parinfo[14]['limited'] = [1,1]		#g1
	parinfo[14]['limits'] = [-1.0,1.0]	
	parinfo[15]['limited'] = [1,1]		#g2
	parinfo[15]['limits'] = [-1.0,1.0]
	parinfo[16]['limited'] = [1,1]		#cx
	parinfo[16]['limits'] = [-L,L]
	parinfo[17]['limited'] = [1,1]		#cy
	parinfo[17]['limits'] = [-L,L]


	m = mpfit.mpfit(disk_bulge_func,guessP,parinfo=parinfo,functkw = fa,quiet=quiet,ftol=ftol)

	if return_image == True:
		fit_image = disk_bulge_model(m.params,psf)
		return m,fit_image
	else:
		return m

    def __init__(self,f, xdata, ydata, p0=None, sigma=None, fixed=None,limits=None, maxiter=200,quiet=1):
        args, varargs, varkw, defaults = inspect.getargspec(f)
       
        self.quiet = quiet
        if len(args) < 2:
            msg = "Unable to determine number of fit parameters."
            raise ValueError(msg)
       
        def general_function(params, fjac=None, xdata=None, ydata=None):
            return [0, f(xdata, *params) - ydata]

        def weighted_general_function(params,fjac=None, xdata=None, ydata=None, weights=None):
            return [0, weights * (f(xdata, *params) - ydata)]


        parinfo = [{'value':1., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}  for i in range(len(args)-1)] 
        

        if p0 != None:
            for x,y in zip(parinfo,p0):
                x['value']=y
        if fixed != None:
            for x,y in zip(parinfo,fixed):
                x['fixed']=y


        if limits != None:
            for x,y in zip(parinfo,limits):
                x['limited']=y[0]
                x['limits']=y[1]
            



    #args = (xdata, ydata, f)
        if sigma is None:
            func = general_function
           # print 'here'
          #  my_functkw = {'xdata':xdata,'ydata':ydata,'function':f}
            my_functkw = {'xdata':xdata,'ydata':ydata}
        else:
            func = weighted_general_function
            weights= 1.0/asarray(sigma)
           # my_functkw = {'xdata':xdata,'ydata':ydata,'function':f,'weights':weights}
            my_functkw = {'xdata':xdata,'ydata':ydata,'weights':weights}


        
       
        result = mpfit(func, functkw=my_functkw, parinfo = parinfo,quiet=self.quiet,maxiter=maxiter)
      
        self.params = result.params
        self.errors = result.perror
        self.dof = result.dof
        self.chi2 = result.fnorm
        self.covar = result.covar

def fit_pixel(inputs):
	wl,q_array,u_array,sigma_q,sigma_u=inputs
	pos=[1,1,10]
	funcargs= {'x':wl, 'y':[q_array,u_array],'err':[sigma_q,sigma_u]}
	fit=mpfit.mpfit(alpha_function,pos,functkw=funcargs,autoderivative=1,quiet=1)
	parms=fit.params
	parms[:2]*=1e-7
	dparms=fit.perror
	dparms[:2]*=1e-7
	return np.array([parms,dparms,fit.fnorm])

def dlosenv_peakfit_sigma_env_fit_test(cat_corr, n_NN=3, fit='gauss', **kwargs): 
    """ Linear fit of the best-fit sigma as a function of environment. 
    Bestfit sigma values are computed for the dLOS distribution in bins of 
    galaxy environment. Linear fit is done using MPFit.

    """

    env_low, env_high, fpeaks, fpeak_errs, sigmas, sigma_errs, amps, nbins =  \
            dlos_envbin_peakfit(
                    cat_corr, 
                    n_NN = n_NN,
                    fit = fit, 
                    **kwargs
                    )
        
    env_mid = np.array( 0.5 * (env_low + env_high) )
    p0 = [ -0.03, 4.0 ] # guess
    fa = {'x': env_mid, 'y': sigmas, 'err': sigma_errs}
    
    fit_param = mpfit.mpfit(mpfit_linear, p0, functkw=fa, quiet=1)
        
    best_slope = fit_param.params[0]
    best_yint = fit_param.params[1]

    print 'Entire range'
    print best_slope, best_yint
    
    for range_test in [20.0, 30.0, 40.0, 50.0]: 

        test_range = np.where(
                env_mid < range_test 
                )
        fa = {'x': env_mid[test_range], 'y': sigmas[test_range], 'err': sigma_errs[test_range]}
        
        fit_param = mpfit.mpfit(mpfit_linear, p0, functkw=fa, quiet=1)
            
        best_slope = fit_param.params[0]
        best_yint = fit_param.params[1]
        
        print 'dNN < ', str(range_test)
        print best_slope, best_yint

    return best_slope, best_yint

def Flat_Field( spec_list, flat ):
    # This Function divides each spectrum in spec_list by the flat and writes
    # The new images as fits files. The output is a list of file names of 
    # the flat fielded images. 
    print "\n====================\n" 
    print 'Flat Fielding Images by Dividing by %s\n' % (flat) 
        
    np.seterr(divide= 'warn')
    flat_data = fits.getdata(flat)
    #If flat is a blue spectrum, find the Littro ghost and add those pixels to the header
    if 'blue' in flat.lower():
        fit_data = np.median(flat_data[0][75:85],axis=0)
        low_index = 1210. #Lowest pixel to search within
        high_index = 1650. #highest pixel to search within
        fit_data1 = fit_data[low_index:high_index]
        fit_pix1 = np.linspace(low_index,low_index+len(fit_data1),num=len(fit_data1))
        max_pixel = np.argmax(fit_data1)
        fit_data2 = fit_data1[max_pixel-30:max_pixel+30]
        guess1 = np.zeros(5)
        guess1[0] = np.mean(fit_data2)
        guess1[1] = (fit_data2[-1]-fit_data2[0])/len(fit_data2)
        guess1[2] = np.amax(fit_data2)
        guess1[3] = np.argmax(fit_data2)
        guess1[4] = 4.
        error_fit1 = np.ones(len(fit_data2))
        xes1 = np.linspace(0,len(fit_data2)-1,num=len(fit_data2))
        fa1 = {'x':xes1,'y':fit_data2,'err':error_fit1}
        fitparams1 = mpfit.mpfit(fitgaussslope,guess1,functkw=fa1,quiet=True)
        center_pixel = low_index+max_pixel-30.+fitparams1.params[3]
        littrow_ghost = [np.rint(center_pixel-9.),np.rint(center_pixel+9.)]
        diagnostic[0:len(fit_pix1),17] = fit_pix1
        diagnostic[0:len(fit_data1),18] = fit_data1
        diagnostic[0:len(xes1+low_index+max_pixel-30.),19] = xes1+low_index+max_pixel-30.
        diagnostic[0:len(gaussslope(xes1,fitparams1.params)),20] = gaussslope(xes1,fitparams1.params)
        diagnostic[0,21] = littrow_ghost[0]
        diagnostic[1,21] = littrow_ghost[1]
    else:
        littrow_ghost = 'None'
    f_spec_list = []
    for spec in spec_list:
        spec_data = fits.getdata(spec)
        f_spec_data = np.divide(spec_data, flat_data)
        f_spec_data[ np.isnan(f_spec_data) ] = 0
        print "f"+"%s Mean: %.3f StDev: %.3f" % (spec, np.mean(f_spec_data), np.std(f_spec_data) ) 
        hdu = fits.getheader(spec)
        Fix_Header(hdu)
        hdu.set('DATEFLAT', datetime.datetime.now().strftime("%Y-%m-%d"), 'Date of Flat Fielding')
        hdu.set('LITTROW',str(littrow_ghost),'Littrow Ghost location in Flat')
        hdu.append( ('FLATFLD', flat,'Image used to Flat Field.'), 
               useblanks= True, bottom= True )    
        NewHdu = fits.PrimaryHDU(data= f_spec_data, header= hdu)
        new_file_name= check_file_exist('f'+spec)
        NewHdu.writeto(new_file_name, output_verify='warn', clobber= True)
        f_spec_list.append(new_file_name)
    return f_spec_list

def dlosenv_peakfit_sigma_env_fit(cat_corr, n_NN=3, fit='gauss', **kwargs): 
    """ Linear fit of the best-fit sigma as a function of environment. 
    Bestfit sigma values are computed for the dLOS distribution in bins of 
    galaxy environment. Linear fit is done using MPFit.

    """

    env_low, env_high, fpeaks, fpeak_errs, sigmas, sigma_errs, amps, nbins =  \
            dlos_envbin_peakfit(
                    cat_corr, 
                    n_NN = n_NN,
                    fit = fit, 
                    **kwargs
                    )

    env_mid = np.array( 0.5 * (env_low + env_high) )

    if n_NN == 5: 
        fit_range = np.where(
                env_mid < 40.0
                )
    else: 
        raise NotImplementedError()
        
    p0 = [ -0.03, 4.0 ] # guess
    fa = {'x': env_mid[fit_range], 'y': sigmas[fit_range], 'err': sigma_errs[fit_range]}
    
    fit_param = mpfit.mpfit(mpfit_linear, p0, functkw=fa, quiet=1)
        
    best_slope = fit_param.params[0]
    best_yint = fit_param.params[1]
    
    try: 
        if kwargs['writeout']: 
            dlos_fpeak_envbin_fit_file = ''.join([
                direc('data', cat_corr), 
                'DLOS_sigma_env_d', str(n_NN), 'NN_bin_bestfit.dat'
                ])
        
            data_hdr = "Columns : bestfit_slope, bestfit_yint"
            data_fmt = ['%10.5f', '%10.5f']

            np.savetxt(
                    dlos_fpeak_envbin_fit_file, 
                    np.c_[best_slope, best_yint],
                    fmt = data_fmt, 
                    delimiter = '\t', 
                    header = data_hdr
                    ) 

    except KeyError: 
        pass

    return best_slope, best_yint

def onedmoffatfit(xax, data, err = None,
        params=[0,1,0,1,3],fixed=[False,False,False,False,False],
        limitedmin=[False,False,False,True,True],
        limitedmax=[False,False,False,False,False], minpars=[0,0,0,0,0],
        maxpars=[0,0,0,0,0], quiet=True, shh=True,
        veryverbose=False,
        vheight=True, negamp=False,
        usemoments=False):

    def mpfitfun(x,y,err):
        if err is None:
            def f(p,fjac=None): return [0,(y-onedmoffat(x,*p))]
        else:
            def f(p,fjac=None): return [0,(y-onedmoffat(x,*p))/err]
        return f

    if xax == []:
        xax = numpy.arange(len(data))

    parinfo = [ {'n':0,'value':params[0],'limits':[minpars[0],maxpars[0]],
                     'limited':[limitedmin[0],limitedmax[0]],'fixed':fixed[0],
                     'parname':"HEIGHT",'error':0} ,
                {'n':1,'value':params[1],'limits':[minpars[1],maxpars[1]],
                     'limited':[limitedmin[1],limitedmax[1]],'fixed':fixed[1],
                     'parname':"AMPLITUDE",'error':0},
                {'n':2,'value':params[2],'limits':[minpars[2],maxpars[2]],
                     'limited':[limitedmin[2],limitedmax[2]],'fixed':fixed[2],
                     'parname':"SHIFT",'error':0},
                {'n':3,'value':params[3],'limits':[minpars[3],maxpars[3]],
                     'limited':[limitedmin[3],limitedmax[3]],'fixed':fixed[3],
                     'parname':"ALPHA",'error':0},
                 {'n':4,'value':params[4],'limits':[minpars[4],maxpars[4]],
                     'limited':[limitedmin[4],limitedmax[4]],'fixed':fixed[4],
                     'parname':"BETA",'error':0}]

    mp = mpfit(mpfitfun(xax,data,err), parinfo=parinfo, quiet=quiet)
    mpp = mp.params
    mpperr = mp.perror
    chi2 = mp.fnorm
    try:
        dof = len(data) - len(mpp)
    except TypeError:
        dof = len(data)
    if mp.status == 0:
        raise Exception(mp.errmsg)
    if (not shh) or veryverbose:
        print "Fit status: ",mp.status
        for i,p in enumerate(mpp):
            parinfo[i]['value'] = p
            print parinfo[i]['parname'],p," +/- ",mpperr[i]
        print "Chi2: ",chi2," Reduced Chi2: ",chi2/dof," DOF:",dof
    retax = numpy.arange(xax[0], xax[-1], 0.1)
    return mpp, onedmoffat(retax,*mpp), mpperr, [chi2, dof], retax

def onedgaussfit(xax,data,err=None,params=[0,1,0,1],fixed=[False,False,False,False],limitedmin=[False,False,False,True],
        limitedmax=[False,False,False,False],minpars=[0,0,0,0],maxpars=[0,0,0,0],
        quiet=True,shh=True):
    # """
    # Inputs:
       # xax - x axis
       # data - y axis
       # err - error corresponding to data

       # params - Fit parameters: Height of background, Amplitude, Shift, Width
       # fixed - Is parameter fixed?
       # limitedmin/minpars - set lower limits on each parameter (default: width>0)
       # limitedmax/maxpars - set upper limits on each parameter
       # quiet - should MPFIT output each iteration?
       # shh - output final parameters?

    # Returns:
       # Fit parameters
       # Model
       # Fit errors
       # chi2
    # """

    def mpfitfun(x,y,err):
        if err == None:
            def f(p,fjac=None): return [0,(y-onedgaussian(x,*p))]
        else:
            def f(p,fjac=None): return [0,(y-onedgaussian(x,*p))/err]
        return f

    if xax == None:
        xax = numpy.arange(len(data))

    parinfo = [ {'n':0,'value':params[0],'limits':[minpars[0],maxpars[0]],'limited':[limitedmin[0],limitedmax[0]],'fixed':fixed[0],'parname':"HEIGHT",'error':0} ,
                {'n':1,'value':params[1],'limits':[minpars[1],maxpars[1]],'limited':[limitedmin[1],limitedmax[1]],'fixed':fixed[1],'parname':"AMPLITUDE",'error':0},
                {'n':2,'value':params[2],'limits':[minpars[2],maxpars[2]],'limited':[limitedmin[2],limitedmax[2]],'fixed':fixed[2],'parname':"SHIFT",'error':0},
                {'n':3,'value':params[3],'limits':[minpars[3],maxpars[3]],'limited':[limitedmin[3],limitedmax[3]],'fixed':fixed[3],'parname':"WIDTH",'error':0}]

    mp = mpfit(mpfitfun(xax,data,err),parinfo=parinfo,quiet=quiet)
    mpp = mp.params
    mpperr = mp.perror
    chi2 = mp.fnorm

    if mp.status == 0:
        raise Exception(mp.errmsg)

    if not shh:
        for i,p in enumerate(mpp):
            parinfo[i]['value'] = p
            print parinfo[i]['parname'],p," +/- ",mpperr[i]
        print "Chi2: ",mp.fnorm," Reduced Chi2: ",mp.fnorm/len(data)," DOF:",len(data)-len(mpp)

    return mpp,onedgaussian(xax,*mpp),mpperr,chi2

def fitGauss(dataList,nBins=201):

    hist,histBinEdges = np.histogram(dataList,bins=nBins)
    histBinCenters = histBinEdges[0:-1]+np.diff(histBinEdges)/2.

    amplitude = 0.95*np.max(hist)
    x_offset = histBinCenters[np.argmax(hist)]
    sigma = np.std(dataList)*.8
    y_offset = 1.e-8

    params=[sigma, x_offset, amplitude, y_offset]  # First guess at fit params

    #errs = np.sqrt(hist) #for poisson distributed data
    #errs[np.where(errs == 0.)] = 1.

    nListVals = 1.*np.sum(hist)
    pvals = hist/nListVals
    #assume errors are the multinomial distribution
    errs = np.sqrt(nListVals*pvals*(1-pvals))
    errs[np.where(errs == 0.)] = 1.


    quiet = True

    parinfo = [ {'n':0,'value':params[0],'limits':[sigma/10., 10*sigma], 'limited':[True,True],'fixed':False,'parname':"Sigma",'error':0},
       {'n':1,'value':params[1],'limits':[x_offset-sigma*2, x_offset+sigma*2],'limited':[True,True],'fixed':False,'parname':"x offset",'error':0},
       {'n':2,'value':params[2],'limits':[0.5*amplitude, 2.*amplitude],'limited':[True,True],'fixed':False,'parname':"Amplitude",'error':0},
       {'n':3,'value':params[3],'limited':[False,False],'fixed':True,'parname':"y_offset",'error':0}]

    fa = {'x':histBinCenters,'y':hist,'err':errs}

    m = mpfit.mpfit(gaussian, functkw=fa, parinfo=parinfo, maxiter=1000, quiet=quiet)
    if m.status <= 0:
        print m.status, m.errmsg

    mpp = m.params                                #The fit params
    mpperr = m.perror

    for k,p in enumerate(mpp):
        parinfo[k]['value'] = p
        parinfo[k]['error'] = mpperr[k]
        if k==0: sigma = p
        if k==1: x_offset = p
        if k==2: amplitude = p
        if k==3: y_offset = p

    def gaussFitFunc(x):
        return y_offset + amplitude * np.exp( - (( x - x_offset)**2) / ( 2. * (sigma**2)))
    gaussFit = gaussFitFunc(histBinCenters)
    resolution = np.abs(x_offset/(2.355*sigma))

    return {'gaussFit':gaussFit,'resolution':resolution,'sigma':sigma,'x_offset':x_offset,'amplitude':amplitude,'y_offset':y_offset,'hist':hist,'histBinEdges':histBinEdges,'gaussFitFunc':gaussFitFunc,'histBinCenters':histBinCenters,'parinfo':parinfo}