def _fit_orbit_mlogl(new_vxvv,vxvv,vxvv_err,pot,radec,lb,
                     customsky,lb_to_customsky,pmllpmbb_to_customsky,
                     tmockAA,
                     ro,vo,obs):
    """The log likelihood for fitting an orbit"""
    #Use this _parse_args routine, which does forward and backward integration
    iR,ivR,ivT,iz,ivz,iphi= tmockAA._parse_args(True,False,
                                                new_vxvv[0],
                                                new_vxvv[1],
                                                new_vxvv[2],
                                                new_vxvv[3],
                                                new_vxvv[4],
                                                new_vxvv[5])
    if radec or lb or customsky:
        #Need to transform to (l,b), (ra,dec), or a custom set
        #First transform to X,Y,Z,vX,vY,vZ (Galactic)
        X,Y,Z = coords.galcencyl_to_XYZ(iR.flatten(),iphi.flatten(),
                                        iz.flatten(),
                                        Xsun=obs[0]/ro,
                                        Zsun=obs[2]/ro).T
        vX,vY,vZ = coords.galcencyl_to_vxvyvz(ivR.flatten(),ivT.flatten(),
                                              ivz.flatten(),iphi.flatten(),
                                              vsun=nu.array(\
                obs[3:6])/vo,Xsun=obs[0]/ro,Zsun=obs[2]/ro).T
        bad_indx= (X == 0.)*(Y == 0.)*(Z == 0.)
        if True in bad_indx: X[bad_indx]+= ro/10000.
        lbdvrpmllpmbb= coords.rectgal_to_sphergal(X*ro,Y*ro,Z*ro,
                                                  vX*vo,vY*vo,vZ*vo,
                                                  degree=True)
        if lb:
            orb_vxvv= nu.array([lbdvrpmllpmbb[:,0],
                                lbdvrpmllpmbb[:,1],
                                lbdvrpmllpmbb[:,2],
                                lbdvrpmllpmbb[:,4],
                                lbdvrpmllpmbb[:,5],
                                lbdvrpmllpmbb[:,3]]).T
        elif radec:
            #Further transform to ra,dec,pmra,pmdec
            radec= coords.lb_to_radec(lbdvrpmllpmbb[:,0],
                                      lbdvrpmllpmbb[:,1],degree=True,
                                      epoch=None)
            pmrapmdec= coords.pmllpmbb_to_pmrapmdec(lbdvrpmllpmbb[:,4],
                                                    lbdvrpmllpmbb[:,5],
                                                    lbdvrpmllpmbb[:,0],
                                                    lbdvrpmllpmbb[:,1],
                                                    degree=True,
                                                    epoch=None)
            orb_vxvv= nu.array([radec[:,0],radec[:,1],
                                lbdvrpmllpmbb[:,2],
                                pmrapmdec[:,0],pmrapmdec[:,1],
                                lbdvrpmllpmbb[:,3]]).T
        elif customsky:
            #Further transform to ra,dec,pmra,pmdec
            customradec= lb_to_customsky(lbdvrpmllpmbb[:,0],
                                              lbdvrpmllpmbb[:,1],degree=True)
            custompmrapmdec= pmllpmbb_to_customsky(lbdvrpmllpmbb[:,4],
                                                   lbdvrpmllpmbb[:,5],
                                                   lbdvrpmllpmbb[:,0],
                                                   lbdvrpmllpmbb[:,1],
                                                   degree=True)
            orb_vxvv= nu.array([customradec[:,0],customradec[:,1],
                                lbdvrpmllpmbb[:,2],
                                custompmrapmdec[:,0],custompmrapmdec[:,1],
                                lbdvrpmllpmbb[:,3]]).T
    else:
        #shape=(2tintJ-1,6)
        orb_vxvv= nu.array([iR.flatten(),ivR.flatten(),ivT.flatten(),
                            iz.flatten(),ivz.flatten(),iphi.flatten()]).T 
    out= 0.
    for ii in range(vxvv.shape[0]):
        sub_vxvv= (orb_vxvv-vxvv[ii,:].flatten())**2.
        #print(sub_vxvv[nu.argmin(nu.sum(sub_vxvv,axis=1))])
        if not vxvv_err is None:
            sub_vxvv/= vxvv_err[ii,:]**2.
        else:
            sub_vxvv/= 0.01**2.
        out+= logsumexp(-0.5*nu.sum(sub_vxvv,axis=1))
    return -out