func EnvSplitTcB(baseEnv *tempAll.Environment, TcFactors, BeFields []float64, epsAbs, epsRel float64) ([]*tempAll.Environment, error) {
	TcEnv := baseEnv.Copy()
	TcEnv.Be_field = 0.0
	TcEnv.Mu_b = 0.0
	_, err := tempCrit.CritTempSolve(TcEnv, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	Tc := 1.0 / TcEnv.Beta
	omegaFit, err := tempCrit.OmegaFit(TcEnv, tempCrit.OmegaPlus)
	if err != nil {
		return nil, err
	}
	TcEnv.A, TcEnv.B = omegaFit[0], omegaFit[2]
	TcEnv.PairCoeffsReady = true

	result := []*tempAll.Environment{}
	for _, TcFactor := range TcFactors {
		env := TcEnv.Copy()
		T := TcFactor * Tc
		env.Beta = 1.0 / T
		env.Temp = T
		// fix (D1, Mu_h) appropriate for Beta
		//_, err := SolveD1Mu_h(env, epsAbs, epsRel)
		_, err := SolveD1Mu_hMu_b(env, epsAbs, epsRel)
		if err != nil {
			return nil, err
		}
		// keep (D1, Mu_h) independent of magnetic field
		BeNum := len(BeFields)
		thisEnv_BeSplit := env.MultiSplit([]string{"Be_field"}, []int{BeNum}, []float64{BeFields[0]}, []float64{BeFields[BeNum-1]})
		result = append(result, thisEnv_BeSplit...)
	}
	return result, nil
}

func innerMu_h(env *tempAll.Environment, k vec.Vector) float64 {
	sxy := math.Sin(k[0]) - math.Sin(k[1])
	numer := sxy * sxy * math.Tanh(env.Beta*env.Xi_h(k)/2.0)
	denom := env.Mu_b + 2.0*env.Xi_h(k)
	//denom := env.Mu_b - 2.0*env.Be_field*env.A + 2.0*env.Xi_h(k)
	return numer / denom
}

func AbsErrorBeta(env *tempAll.Environment, variables []string) solve.Diffable {
	F := func(v vec.Vector) (float64, error) {
		if v.ContainsNaN() {
			fmt.Printf("got NaN in AbsErrorBeta (v=%v)\n", v)
			return 0.0, errors.New("NaN in input")
		}
		env.Set(v, variables)
		if !env.FixedPairCoeffs || !env.PairCoeffsReady {
			// Before we evaluate error in Beta, Mu_h and D1 should have
			// appropriate values.
			eps := 1e-9
			_, err := D1MuSolve(env, eps, eps)
			if err != nil {
				return 0.0, err
			}
		}
		// Beta equation error = x - x1 - x2
		x1 := X1(env)
		x2, err := tempCrit.X2(env)
		if err != nil {
			fmt.Printf("error from X2(): %v\n", err)
			return 0.0, err
		}
		lhs := env.X
		rhs := x1 + x2
		return lhs - rhs, nil
	}
	h := 1e-5
	epsabs := 1e-4
	return solve.SimpleDiffable(F, len(variables), h, epsabs)
}

// Solve the (D1, Mu_h, Beta) system with x and F0 fixed.
func D1MuBetaSolve(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	// our guess for beta should be above beta_c
	if env.A == 0.0 && env.B == 0.0 {
		D1, Mu_h, F0 := env.D1, env.Mu_h, env.F0
		env.F0 = 0.0 // F0 is 0 at T_c
		_, err := tempCrit.CritTempSolve(env, epsAbs, epsRel)
		if err != nil {
			return nil, err
		}
		fmt.Printf("%v; Tc = %f\n", env, 1.0/env.Beta)
		omegaFit, err := tempCrit.OmegaFit(env, tempCrit.OmegaPlus)
		if err != nil {
			return nil, err
		}
		env.A, env.B = omegaFit[0], omegaFit[2]
		env.PairCoeffsReady = true
		env.Beta += 0.1
		// we are at T < T_c; uncache env
		env.D1, env.Mu_h, env.F0 = D1, Mu_h, F0
	}
	//fmt.Printf("%v; Tc = %f\n", env, 1.0 / env.Beta)
	// solve low temp system for reasonable values of D1 and Mu_h first
	_, err := D1MuSolve(env, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	// solve the full low temp system
	system, start := D1MuBetaSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

// Calculate x - (x_1 + x_2) with Mu_h fixed.
func AbsErrorX(env *tempAll.Environment, variables []string) solve.Diffable {
	F := func(v vec.Vector) (float64, error) {
		if v.ContainsNaN() {
			fmt.Printf("got NaN in AbsErrorX (v=%v)\n", v)
			return 0.0, errors.New("NaN in input")
		}
		env.Set(v, variables)
		// Before we evaluate error in X, Mu_b and D1 should have
		// appropriate values.
		system, start := D1Mu_bSystem(env)
		eps := 1e-9
		_, err := solve.MultiDim(system, start, eps, eps)
		if err != nil {
			return 0.0, err
		}
		if env.Mu_b > 0.0 {
			fmt.Println("Warning: got Mu_b > 0 in AbsErrorX")
			env.Mu_b = 0.0
		}
		// evaluate X error
		x1 := tempPair.X1(env)
		x2, err := tempCrit.X2(env)
		if err != nil {
			fmt.Printf("error from X2(): %v\n", err)
			return 0.0, err
		}
		lhs := env.X
		rhs := x1 + x2
		return lhs - rhs, nil
	}
	h := 1e-5
	epsabs := 1e-4
	return solve.SimpleDiffable(F, len(variables), h, epsabs)
}

func AbsErrorMu_h(env *tempAll.Environment, variables []string) solve.Diffable {
	F := func(v vec.Vector) (float64, error) {
		env.Set(v, variables)
		lhs := env.X
		rhs := X1(env)
		return lhs - rhs, nil
	}
	h := 1e-6
	epsabs := 1e-4
	return solve.SimpleDiffable(F, len(variables), h, epsabs)
}

func SolveNoninteracting(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	env.F0 = 0.0
	env.Mu_h = 0.3
	env.Beta = 50.0
	system, start := NoninteractingSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

// Evaluate the anomalous retarded pair Green's function,
// Pi^A(k, omega)_{xx, xy, yy}. k must be a two-dimensional vector.
func PiAnom(env *tempAll.Environment, k vec.Vector, omega float64) vec.Vector {
	piInner := func(q vec.Vector, out *vec.Vector) {
		// Do vector operations on out to avoid allocation:
		//  first case, out = k/2 + q
		(*out)[0] = k[0]/2.0 + q[0]
		(*out)[1] = k[1]/2.0 + q[1]
		Delta1 := env.Delta_h(*out)
		E1 := env.BogoEnergy(*out)
		//  second case, out = k/2 - q
		(*out)[0] = k[0]/2.0 - q[0]
		(*out)[1] = k[1]/2.0 - q[1]
		Delta2 := env.Delta_h(*out)
		E2 := env.BogoEnergy(*out)
		// Get part of result that's the same for all (xx, xy, yy):
		t1 := math.Tanh(env.Beta * E1 / 2.0)
		t2 := math.Tanh(env.Beta * E2 / 2.0)
		common := -Delta1 * Delta2 / (4.0 * E1 * E2) * ((t1+t2)*(1.0/(omega+E1+E2)-1.0/(omega-E1-E2)) + (t1-t2)*(1.0/(omega-E1+E2)-1.0/(omega+E1-E2)))
		// Set out = result:
		sx := math.Sin(q[0])
		sy := math.Sin(q[1])
		(*out)[0] = sx * sx * common
		(*out)[1] = sx * sy * common
		(*out)[2] = sy * sy * common
	}
	return bzone.VectorAvg(env.PointsPerSide, 2, 3, piInner)
}

// Partial derivative of Mu_h with respect to T; x and V held constant.
func dMu_hdT(env *tempAll.Environment) (float64, error) {
	ct := 0
	// F gets Mu_h given Beta
	F := func(Beta float64) (float64, error) {
		ct += 1
		// save the environment state before changing it
		// (don't want one call of F to affect the next)
		oD1, oMu_h, oBeta, oMu_b := env.D1, env.Mu_h, env.Beta, env.Mu_b
		env.Beta = Beta
		// fix free variables
		eps := 1e-9
		_, err := D1MuF0Solve(env, eps, eps)
		if err != nil {
			return 0.0, err
		}
		Mu_h := env.Mu_h
		// restore the environment
		env.D1, env.Mu_h, env.Beta, env.Mu_b = oD1, oMu_h, oBeta, oMu_b
		return Mu_h, nil
	}
	h := 1e-4
	epsAbs := 1e-5
	deriv, err := solve.OneDimDerivative(F, env.Beta, h, epsAbs)
	//fmt.Println("MuT ct", ct)
	return -math.Pow(env.Beta, 2.0) * deriv, err
}

// Partial derivative of F (some function of env) with respect to Mu_h;
// T and V held constant.
func dFdMu_h(env *tempAll.Environment, F envFunc) (float64, error) {
	ct := 0
	// G gets F given Mu_h (allow x to vary; constant Beta)
	G := func(Mu_h float64) (float64, error) {
		ct += 1
		// save the environment state before changing it
		// (don't want one call of F to affect the next)
		oD1, oMu_h, oX, oMu_b := env.D1, env.Mu_h, env.X, env.Mu_b
		env.Mu_h = Mu_h
		// fix free variables2
		eps := 1e-9
		_, err := D1F0XSolve(env, eps, eps)
		if err != nil {
			return 0.0, err
		}
		vF, err := F(env)
		if err != nil {
			return 0.0, err
		}
		// restore the environment
		env.D1, env.Mu_h, env.X, env.Mu_b = oD1, oMu_h, oX, oMu_b
		return vF, nil
	}
	h := 1e-4
	epsAbs := 1e-5
	deriv, err := solve.OneDimDerivative(G, env.Mu_h, h, epsAbs)
	//fmt.Println("dF_dMu ct", ct)
	return deriv, err
}

// Partial derivative of F with respect to T; Mu_h and V held constant.
func dFdT(env *tempAll.Environment, F envFunc) (float64, error) {
	ct := 0
	// G gets F given Beta (allow x to vary; constant Mu_h)
	G := func(Beta float64) (float64, error) {
		ct += 1
		// save the environment state before changing it
		// (don't want one call of F to affect the next)
		oD1, oBeta, oX, oMu_b := env.D1, env.Beta, env.X, env.Mu_b
		env.Beta = Beta
		// fix free variables
		eps := 1e-9
		_, err := SolveD1Mu_bX(env, eps, eps)
		if err != nil {
			return 0.0, err
		}
		vF, err := F(env)
		if err != nil {
			return 0.0, err
		}
		// restore the environment
		env.D1, env.Beta, env.X, env.Mu_b = oD1, oBeta, oX, oMu_b
		return vF, nil
	}
	h := 1e-4
	epsAbs := 1e-5
	deriv, err := solve.OneDimDerivative(G, env.Beta, h, epsAbs)
	fmt.Println("dF_dT ct", ct)
	return -math.Pow(env.Beta, 2.0) * deriv, err
}

// For use with solve.Iterative:
func CritTempStages(env *tempAll.Environment) ([]solve.DiffSystem, []vec.Vector, func([]vec.Vector)) {
	vars0 := []string{"D1", "Mu_h"}
	vars1 := []string{"Beta"}
	diffD1 := tempPair.AbsErrorD1(env, vars0)
	diffMu_h := tempPair.AbsErrorBeta(env, vars0)
	system0 := solve.Combine([]solve.Diffable{diffD1, diffMu_h})
	diffBeta := AbsErrorBeta(env, vars1)
	system1 := solve.Combine([]solve.Diffable{diffBeta})
	stages := []solve.DiffSystem{system0, system1}
	start := []vec.Vector{[]float64{env.D1, env.Mu_h}, []float64{env.Beta}}
	accept := func(x []vec.Vector) {
		env.D1 = x[0][0]
		env.Mu_h = x[0][1]
		env.Beta = x[1][0]
	}
	return stages, start, accept
}

// Calculate Mu_b - (-Omega_+(0))
func AbsErrorMu_b(env *tempAll.Environment, variables []string) solve.Diffable {
	F := func(v vec.Vector) (float64, error) {
		if v.ContainsNaN() {
			fmt.Printf("got NaN in AbsErrorMu_b (v=%v)\n", v)
			return 0.0, errors.New("NaN in input")
		}
		env.Set(v, variables)
		zv := vec.ZeroVector(3)
		omega0, err := tempCrit.OmegaPlus(env, zv)
		if err != nil {
			return 0.0, err
		}
		lhs := env.Mu_b
		rhs := -omega0
		return lhs - rhs, nil
	}
	h := 1e-5
	epsabs := 1e-4
	return solve.SimpleDiffable(F, len(variables), h, epsabs)
}

// Solve the (D1, Mu_h, Beta) system with x and Mu_b fixed.
func FlucTempSolve(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	// fix pair coefficients
	if env.A == 0.0 && env.B == 0.0 && env.FixedPairCoeffs {
		D1, Mu_h, Mu_b, Beta := env.D1, env.Mu_h, env.Mu_b, env.Beta
		env.Mu_b = 0.0 // Mu_b is 0 at T_c
		_, err := tempCrit.CritTempSolve(env, epsAbs, epsRel)
		if err != nil {
			return nil, err
		}
		omegaFit, err := tempCrit.OmegaFit(env, tempCrit.OmegaPlus)
		if err != nil {
			return nil, err
		}
		env.A, env.B = omegaFit[0], omegaFit[2]
		env.PairCoeffsReady = true
		// uncache env
		env.D1, env.Mu_h, env.Mu_b, env.Beta = D1, Mu_h, Mu_b, Beta
	}
	// our guess for beta should be a bit above Beta_p
	pairSystem, pairStart := tempPair.PairTempSystem(env)
	_, err := solve.MultiDim(pairSystem, pairStart, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	env.Beta += 0.1
	// solve fluc temp system for reasonable values of Mu_h and D1 first
	system, start := FlucTempD1MuSystem(env)
	_, err = solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	// solve the full fluc temp system
	system, start = FlucTempFullSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

// Calculate U_{1}/N = 1/N \sum_k \epsilon_h(k) f_h(\xi_h(k))
func HolonEnergy(env *tempAll.Environment) (float64, error) {
	inner := func(k vec.Vector) float64 {
		return env.Epsilon_h(k) * env.Fermi(env.Xi_h(k))
	}
	dim := 2
	avg := bzone.Avg(env.PointsPerSide, dim, inner)
	return avg, nil
}

// Solve the environment under the conditions at T = T_c.
func CritTempSolve(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	// our guess for beta should be a bit above Beta_p
	pairSystem, pairStart := tempPair.PairTempSystem(env)
	_, err := solve.MultiDim(pairSystem, pairStart, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	env.Beta += 0.1
	// solve crit temp system for reasonable values of Mu and D1 first
	system, start := CritTempD1MuSystem(env)
	_, err = solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	// solve the full crit temp system
	system, start = CritTempFullSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

// Solve the (D1, Mu_h, F0) system with x and Beta fixed.
func D1MuF0Solve(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	if env.A == 0.0 && env.B == 0.0 {
		// We must have T < T_c < T_p (Beta > Beta_c > Beta_p).
		// Getting Beta_p is fast, so do that first.
		D1, Mu_h, F0, Beta := env.F0, env.Mu_h, env.F0, env.Beta // cache env
		env.F0 = 0.0                                             // F0 is 0 at T_c and T_p
		_, err := tempPair.PairTempSolve(env, epsAbs, epsRel)
		if err != nil {
			return nil, err
		}
		if Beta < env.Beta {
			return nil, fmt.Errorf("Beta = %f less than Beta_p in env %s", Beta, env.String())
		}
		_, err = tempCrit.CritTempSolve(env, epsAbs, epsRel)
		if err != nil {
			return nil, err
		}
		if Beta < env.Beta {
			return nil, fmt.Errorf("Beta = %f less than Beta_c in env %s", Beta, env.String())
		}
		fmt.Printf("%v; Tc = %f\n", env, 1.0/env.Beta)
		omegaFit, err := tempCrit.OmegaFit(env, tempCrit.OmegaPlus)
		if err != nil {
			return nil, err
		}
		env.A, env.B = omegaFit[0], omegaFit[2]
		env.PairCoeffsReady = true
		// we are at T < T_c; uncache env
		env.D1, env.Mu_h, env.F0, env.Beta = D1, Mu_h, F0, Beta
	}
	// solve low temp system for reasonable values of D1 and Mu_h first
	_, err := D1MuSolve(env, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	// solve the full low temp system
	system, start := D1MuF0System(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

func innerX1(env *tempAll.Environment, k vec.Vector) float64 {
	return env.Fermi(env.Xi_h(k))
}

func innerD1(env *tempAll.Environment, k vec.Vector) float64 {
	return math.Sin(k[0]) * math.Sin(k[1]) * env.Fermi(env.Xi_h(k))
}

func innerBeta(env *tempAll.Environment, k vec.Vector) float64 {
	sxy := math.Sin(k[0]) - math.Sin(k[1])
	return sxy * sxy * math.Tanh(env.Beta*env.Xi_h(k)/2.0) / env.Xi_h(k)
}