int main()
{
	FILE *planets,*results;
	int i;
	struct planet ss[9];

	if((planets=fopen("planets.txt","r"))==NULL)
	{
		printf("Error loading planets.txt\n");
		printf("Please press enter to continue");
		getchar();
		return(0);
	}

	if((results=fopen("results.txt","w"))==NULL)
	{
		printf("Error loading results.txt\n");
		printf("Please press enter to continue");
		getchar();
		return(0);
	}


	for(i=0;i<9;i++)
	{
		fscanf(planets,"%s %lf %lf %lf",&ss[i].name,&ss[i].distance,&ss[i].mass,&ss[i].radius);
	}

	printf("Name       Distance   Mass   Radius\n");
	for(i=0;i<9;i++)
	{
		printf("%10s %g %g %g\n",ss[i].name,ss[i].distance,ss[i].mass,ss[i].radius);
	}

	for(i=0;i<9;i++)
	{
		fprintf(results,"%s %g %g\n",ss[i].name,ylength(&ss[i]),density(&ss[i]));
	}
        printf("Please press enter to continue");
        getchar();
        return(0);
}

doublereal WaterPropsIAPWS::psat(doublereal temperature, int waterState) {
    doublereal densLiq = -1.0, densGas = -1.0, delGRT = 0.0;
    doublereal dp, pcorr;
    if (temperature >= T_c) {
        densGas = density(temperature, P_c, WATER_SUPERCRIT);
        setState_TR(temperature, densGas);
        return P_c;
    }
    doublereal p = psat_est(temperature);
    bool conv = false;
    for (int i = 0; i < 30; i++) {
        if (method == 1) {
            corr(temperature, p, densLiq, densGas, delGRT);
            doublereal delV = M_water * (1.0/densLiq - 1.0/densGas);
            dp = - delGRT * Rgas * temperature / delV;
        } else {
            corr1(temperature, p, densLiq, densGas, pcorr);
            dp = pcorr - p;
        }
        p += dp;

        if ((method == 1) && delGRT < 1.0E-8) {
            conv = true;
            break;
        } else {
            if (fabs(dp/p) < 1.0E-9) {
                conv = true;
                break;
            }
        }
    }
    // Put the fluid in the desired end condition
    if (waterState == WATER_LIQUID) {
        setState_TR(temperature, densLiq);
    } else if (waterState == WATER_GAS) {
        setState_TR(temperature, densGas);
    } else {
        throw Cantera::CanteraError("WaterPropsIAPWS::psat",
                                    "unknown water state input: " + Cantera::int2str(waterState));
    }
    return p;
}

bool AtmosphericMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
        MediumState &state, MediumSample &sample) const
{
    if (state.bounce > _maxBounce)
        return false;

    Vec3f p = (ray.pos() - _center);
    float t0 = p.dot(ray.dir());
    float  h = (p - t0*ray.dir()).length();

    float maxT = ray.farT() + t0;
    if (_absorptionOnly) {
        sample.t = ray.farT();
        sample.weight = std::exp(-_sigmaT*densityIntegral(h, t0, maxT));
        sample.pdf = 1.0f;
        sample.exited = true;
    } else {
        int component = sampler.nextDiscrete(3);
        float sigmaTc = _sigmaT[component];
        float xi = 1.0f - sampler.next1D();

        float t = inverseOpticalDepth(h, t0, sigmaTc, xi);
        sample.t = min(t, maxT);
        sample.weight = std::exp(-_sigmaT*densityIntegral(h, t0, sample.t));
        sample.exited = (t >= maxT);
        if (sample.exited) {
            sample.pdf = sample.weight.avg();
        } else {
            float rho = density(h, sample.t);
            sample.pdf = (rho*_sigmaT*sample.weight).avg();
            sample.weight *= rho*_sigmaS;
        }
        sample.weight /= sample.pdf;
        sample.t -= t0;

        state.advance();
    }
    sample.p = ray.pos() + sample.t*ray.dir();
    sample.phase = _phaseFunction.get();

    return true;
}

                    double l1Distance(const Kernel &rhs, long nr = 1000)
                    {
                        double a,b,x;
                        a = min_ - 4.0*width_;
                        b = rhs.min_ - 4.0*rhs.width_;
                        double dm = (a < b ? a : b);
                        a = max_ + 4.0*width_;
                        b = rhs.max_ + 4.0*rhs.width_;
                        double dM = (a > b ? a : b);
                        double dD=(dM-dm)/(nr-1);
                        const function::Linear<double> lt(0, dm, nr-1, dM);

                        dM = 0.0;
                        for (long i=0; i<nr; ++i)
                        {
                            x = lt(i);
                            dM += fabs(density(x)-rhs.density(x))*dD;
                        }
                        return dM;
                    }

int MixtureFugacityTP::phaseState(bool checkState) const
{
    int state = iState_;
    if (checkState) {
        double t = temperature();
        double tcrit = critTemperature();
        double rhocrit = critDensity();
        if (t >= tcrit) {
            return FLUID_SUPERCRIT;
        }
        double tmid = tcrit - 100.;
        if (tmid < 0.0) {
            tmid = tcrit / 2.0;
        }
        double pp = psatEst(tmid);
        double mmw = meanMolecularWeight();
        double molVolLiqTmid = liquidVolEst(tmid, pp);
        double molVolGasTmid = GasConstant * tmid / (pp);
        double densLiqTmid = mmw / molVolLiqTmid;
        double densGasTmid = mmw / molVolGasTmid;
        double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
        doublereal rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);

        double rho = density();
        int iStateGuess = FLUID_LIQUID_0;
        if (rho < rhoMid) {
            iStateGuess = FLUID_GAS;
        }
        double molarVol = mmw / rho;
        double presCalc;

        double dpdv = dpdVCalc(t, molarVol, presCalc);
        if (dpdv < 0.0) {
            state =  iStateGuess;
        } else {
            state = FLUID_UNSTABLE;
        }

    }
    return state;
}

void XZHydrostatic_Omega<EvalT, Traits>::
evaluateFields(typename Traits::EvalData workset)
{
#ifndef ALBANY_KOKKOS_UNDER_DEVELOPMENT
  for (int cell=0; cell < workset.numCells; ++cell) {
    for (int qp=0; qp < numQPs; ++qp) {
      for (int level=0; level < numLevels; ++level) {
        ScalarT                               sum  = -0.5*divpivelx(cell,qp,level) * E.delta(level);
        for (int j=0; j<level; ++j)           sum -=     divpivelx(cell,qp,j)     * E.delta(j);
        for (int dim=0; dim < numDims; ++dim) sum += Velocity(cell,qp,level,dim)*gradp(cell,qp,level,dim);
        omega(cell,qp,level) = sum/(Cpstar(cell,qp,level)*density(cell,qp,level));
      }
    }
  }

#else
  Kokkos::parallel_for(XZHydrostatic_Omega_Policy(0,workset.numCells),*this);
  cudaCheckError();

#endif
}

        result_type result(Args const &args) const
        {
            if (this->is_dirty)
            {
                this->is_dirty = false;

                std::size_t cnt = count(args);
                range_type histogram = density(args);
                typename range_type::iterator it = histogram.begin();
                while (this->sum < 0.5 * cnt)
                {
                    this->sum += it->second * cnt;
                    ++it;
                }
                --it;
                float_type over = numeric::average(this->sum - 0.5 * cnt, it->second * cnt);
                this->median = it->first * over + (it + 1)->first * (1. - over);
            }

            return this->median;
        }

double N_from_mu(Minimizer *min, Grid *potential, const Grid &constraint, double mu) {
  Functional f = constrain(constraint, OfEffectivePotential(HS + IdealGas()
                                                            + ChemicalPotential(mu)));
  double Nnow = 0;
  min->minimize(f, potential->description());
  for (int i=0;i<numiters && min->improve_energy(false);i++) {
    Grid density(potential->description(), EffectivePotentialToDensity()(1, potential->description(), *potential));
    Nnow = density.sum()*potential->description().dvolume;
    printf("Nnow is %g vs %g\n", Nnow, N);
    fflush(stdout);
    
    density.epsNativeSlice("papers/contact/figs/box.eps", 
                           Cartesian(0,ymax+2,0), Cartesian(0,0,zmax+2), 
                           Cartesian(0,-ymax/2-1,-zmax/2-1));
    density.epsNativeSlice("papers/contact/figs/box-diagonal.eps", 
                           Cartesian(xmax+2,0,zmax+2),  Cartesian(0,ymax+2,0),
                           Cartesian(-xmax/2-1,-ymax/2-1,-zmax/2-1));
    //sleep(3);
  }
  return Nnow;
}

doublereal SingleSpeciesTP::cv_mole() const
{
    /*
     *  For single species, we go directory to the general Cp - Cv relation
     *
     *  Cp = Cv + alpha**2 * V * T / beta
     *
     * where
     *     alpha = volume thermal expansion coefficient
     *     beta  = isothermal compressibility
     */
    doublereal cvbar = cp_mole();
    doublereal alpha = thermalExpansionCoeff();
    doublereal beta = isothermalCompressibility();
    doublereal V = molecularWeight(0)/density();
    doublereal T = temperature();
    if (beta != 0.0) {
        cvbar -= alpha * alpha * V * T / beta;
    }
    return cvbar;
}

void WaterSSTP::getEntropy_R_ref(doublereal* sr) const
{
    doublereal p = pressure();
    double T = temperature();
    double dens = density();
    int waterState = WATER_GAS;
    double rc = m_sub.Rhocrit();
    if (dens > rc) {
        waterState = WATER_LIQUID;
    }
    doublereal dd = m_sub.density(T, OneAtm, waterState, dens);

    if (dd <= 0.0) {
        throw CanteraError("setPressure", "error");
    }
    m_sub.setState_TR(T, dd);

    doublereal s = m_sub.entropy();
    *sr = (s + SW_Offset)/ GasConstant;
    dd = m_sub.density(T, p, waterState, dens);
}

void IonFlow::electricFieldMethod(const double* x, size_t j0, size_t j1)
{
    for (size_t j = j0; j < j1; j++) {
        double wtm = m_wtm[j];
        double rho = density(j);
        double dz = z(j+1) - z(j);

        // mixture-average diffusion
        double sum = 0.0;
        for (size_t k = 0; k < m_nsp; k++) {
            m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
            m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
            sum -= m_flux(k,j);
        }

        // ambipolar diffusion
        double E_ambi = E(x,j);
        for (size_t k : m_kCharge) {
            double Yav = 0.5 * (Y(x,k,j) + Y(x,k,j+1));
            double drift = rho * Yav * E_ambi
                           * m_speciesCharge[k] * m_mobility[k+m_nsp*j];
            m_flux(k,j) += drift;
        }

        // correction flux
        double sum_flux = 0.0;
        for (size_t k = 0; k < m_nsp; k++) {
            sum_flux -= m_flux(k,j); // total net flux
        }
        double sum_ion = 0.0;
        for (size_t k : m_kCharge) {
            sum_ion += Y(x,k,j);
        }
        // The portion of correction for ions is taken off
        for (size_t k : m_kNeutral) {
            m_flux(k,j) += Y(x,k,j) / (1-sum_ion) * sum_flux;
        }
    }
}

// Calculate the modes of the density for the distribution
// and their associated errors
double
vpScale::MeanShift(vpColVector &error)
{

  int n = error.getRows()/dimension;
  vpColVector density(n);
  vpColVector density_gradient(n);
  vpColVector mean_shift(n);

  int increment=1;

  // choose smallest error as start point
  int i=0;
  while(error[i]<0 && error[i]<error[i+1])
    i++;

  // Do mean shift until no shift
  while(increment >= 1 && i<n)
  {
    increment=0;
    density[i] = KernelDensity(error, i);
    density_gradient[i] = KernelDensityGradient(error, i);
    mean_shift[i]=vpMath::sqr(bandwidth)*density_gradient[i]/((dimension+2)*density[i]);

    double tmp_shift = mean_shift[i];

    // Do mean shift
    while(tmp_shift>0 && tmp_shift>error[i]-error[i+1])
    {
      i++;
      increment++;
      tmp_shift-=(error[i]-error[i-1]);
    }
  }

  return error[i];

}

    ShardDefinition::ShardDefinition(std::istream& input)
    {
        CsvTsv::InputColumn<size_t> 
            minPostings("MinPostings",
                        "Minimum number of postings for any document in shard.");
        CsvTsv::InputColumn<double>
            density("Density",
                    "Target density for RowTables in the shard.");

        CsvTsv::CsvTableParser parser(input);
        CsvTsv::TableReader reader(parser);
        reader.DefineColumn(minPostings);
        reader.DefineColumn(density);

        reader.ReadPrologue();
        while (!reader.AtEOF())
        {
            reader.ReadDataRow();

            AddShard(minPostings, density);
        }
        reader.ReadEpilogue();
    }

Real SodProblem::valueExact(Real t, const Point& p, int eq)
{
	switch (eq) {
	case 0:
		return density(t, p);
		break;
	case 1:
		return momentumX(t, p);
		break;
	case 2:
		return momentumY(t, p);
		break;
	case 3:
		return momentumZ(t, p);
	case 4:
		return energyTotal(t, p);
		break;
	default:
		return 0.0;
		mooseError("不可用的分量" << eq);
		break;
	}
}

void XZHydrostatic_GeoPotential<EvalT, Traits>::
evaluateFields(typename Traits::EvalData workset)
{
  const Eta<EvalT> &E = Eta<EvalT>::self();

  ScalarT sum;
  for (int cell=0; cell < workset.numCells; ++cell) {
    for (int node=0; node < numNodes; ++node) {
      for (int level=0; level < numLevels; ++level) {
        
        sum =
        PhiSurf(cell,node) +
        0.5 * Pi(cell,node,level) * E.delta(level) / density(cell,node,level);
        for (int j=level+1; j < numLevels; ++j) sum += Pi(cell,node,j)     * E.delta(j)     / density(cell,node,j);

        Phi(cell,node,level) = sum;
        
        
        //std::cout <<"Inside GeoP, cell, node, PhiSurf(cell,node)="<<cell<<
        //", "<<node<<", "<<PhiSurf(cell,node) <<std::endl;
      }
    }
  }
}

Real
RichardsDensityConstBulk::d2density(Real p) const
{
  return density(p)/_bulk/_bulk;
}

/*
 * Format a summary of the mixture state for output.
 */
void MolalityVPSSTP::reportCSV(std::ofstream& csvFile) const
{


    csvFile.precision(3);
    int tabS = 15;
    int tabM = 30;
    int tabL = 40;
    try {
        if (name() != "") {
            csvFile << "\n"+name()+"\n\n";
        }
        csvFile << setw(tabL) << "temperature (K) =" << setw(tabS) << temperature() << endl;
        csvFile << setw(tabL) << "pressure (Pa) =" << setw(tabS) << pressure() << endl;
        csvFile << setw(tabL) << "density (kg/m^3) =" << setw(tabS) << density() << endl;
        csvFile << setw(tabL) << "mean mol. weight (amu) =" << setw(tabS) << meanMolecularWeight() << endl;
        csvFile << setw(tabL) << "potential (V) =" << setw(tabS) << electricPotential() << endl;
        csvFile << endl;

        csvFile << setw(tabL) << "enthalpy (J/kg) = " << setw(tabS) << enthalpy_mass() << setw(tabL) << "enthalpy (J/kmol) = " << setw(tabS) << enthalpy_mole() << endl;
        csvFile << setw(tabL) << "internal E (J/kg) = " << setw(tabS) << intEnergy_mass() << setw(tabL) << "internal E (J/kmol) = " << setw(tabS) << intEnergy_mole() << endl;
        csvFile << setw(tabL) << "entropy (J/kg) = " << setw(tabS) << entropy_mass() << setw(tabL) << "entropy (J/kmol) = " << setw(tabS) << entropy_mole() << endl;
        csvFile << setw(tabL) << "Gibbs (J/kg) = " << setw(tabS) << gibbs_mass() << setw(tabL) << "Gibbs (J/kmol) = " << setw(tabS) << gibbs_mole() << endl;
        csvFile << setw(tabL) << "heat capacity c_p (J/K/kg) = " << setw(tabS) << cp_mass() << setw(tabL) << "heat capacity c_p (J/K/kmol) = " << setw(tabS) << cp_mole() << endl;
        csvFile << setw(tabL) << "heat capacity c_v (J/K/kg) = " << setw(tabS) << cv_mass() << setw(tabL) << "heat capacity c_v (J/K/kmol) = " << setw(tabS) << cv_mole() << endl;

        csvFile.precision(8);

        vector<std::string> pNames;
        vector<vector_fp> data;
        vector_fp temp(nSpecies());

        getMoleFractions(&temp[0]);
        pNames.push_back("X");
        data.push_back(temp);
        try {
            getMolalities(&temp[0]);
            pNames.push_back("Molal");
            data.push_back(temp);
        } catch (CanteraError& err) {
            err.save();
        }
        try {
            getChemPotentials(&temp[0]);
            pNames.push_back("Chem. Pot. (J/kmol)");
            data.push_back(temp);
        } catch (CanteraError& err) {
            err.save();
        }
        try {
            getStandardChemPotentials(&temp[0]);
            pNames.push_back("Chem. Pot. SS (J/kmol)");
            data.push_back(temp);
        } catch (CanteraError& err) {
            err.save();
        }
        try {
            getMolalityActivityCoefficients(&temp[0]);
            pNames.push_back("Molal Act. Coeff.");
            data.push_back(temp);
        } catch (CanteraError& err) {
            err.save();
        }
        try {
            getActivities(&temp[0]);
            pNames.push_back("Molal Activity");
            data.push_back(temp);
            size_t iHp = speciesIndex("H+");
            if (iHp != npos) {
                double pH = -log(temp[iHp]) / log(10.0);
                csvFile << setw(tabL) << "pH = " << setw(tabS) << pH << endl;
            }
        } catch (CanteraError& err) {
            err.save();
        }
        try {
            getPartialMolarEnthalpies(&temp[0]);
            pNames.push_back("Part. Mol Enthalpy (J/kmol)");
            data.push_back(temp);
        } catch (CanteraError& err) {
            err.save();
        }
        try {
            getPartialMolarEntropies(&temp[0]);
            pNames.push_back("Part. Mol. Entropy (J/K/kmol)");
            data.push_back(temp);
        } catch (CanteraError& err) {
            err.save();
        }
        try {
            getPartialMolarIntEnergies(&temp[0]);
            pNames.push_back("Part. Mol. Energy (J/kmol)");
            data.push_back(temp);
        } catch (CanteraError& err) {
            err.save();
        }
        try {
//.........这里部分代码省略.........

/**
  * Format a summary of the mixture state for output.
  */
std::string MolalityVPSSTP::report(bool show_thermo) const
{


    char p[800];
    string s = "";
    try {
        if (name() != "") {
            sprintf(p, " \n  %s:\n", name().c_str());
            s += p;
        }
        sprintf(p, " \n       temperature    %12.6g  K\n", temperature());
        s += p;
        sprintf(p, "          pressure    %12.6g  Pa\n", pressure());
        s += p;
        sprintf(p, "           density    %12.6g  kg/m^3\n", density());
        s += p;
        sprintf(p, "  mean mol. weight    %12.6g  amu\n", meanMolecularWeight());
        s += p;

        doublereal phi = electricPotential();
        sprintf(p, "         potential    %12.6g  V\n", phi);
        s += p;

        size_t kk = nSpecies();
        vector_fp x(kk);
        vector_fp molal(kk);
        vector_fp mu(kk);
        vector_fp muss(kk);
        vector_fp acMolal(kk);
        vector_fp actMolal(kk);
        getMoleFractions(&x[0]);
        getMolalities(&molal[0]);
        getChemPotentials(&mu[0]);
        getStandardChemPotentials(&muss[0]);
        getMolalityActivityCoefficients(&acMolal[0]);
        getActivities(&actMolal[0]);

        size_t iHp = speciesIndex("H+");
        if (iHp != npos) {
            double pH = -log(actMolal[iHp]) / log(10.0);
            sprintf(p, "                pH    %12.4g  \n", pH);
            s += p;
        }

        if (show_thermo) {
            sprintf(p, " \n");
            s += p;
            sprintf(p, "                          1 kg            1 kmol\n");
            s += p;
            sprintf(p, "                       -----------      ------------\n");
            s += p;
            sprintf(p, "          enthalpy    %12.6g     %12.4g     J\n",
                    enthalpy_mass(), enthalpy_mole());
            s += p;
            sprintf(p, "   internal energy    %12.6g     %12.4g     J\n",
                    intEnergy_mass(), intEnergy_mole());
            s += p;
            sprintf(p, "           entropy    %12.6g     %12.4g     J/K\n",
                    entropy_mass(), entropy_mole());
            s += p;
            sprintf(p, "    Gibbs function    %12.6g     %12.4g     J\n",
                    gibbs_mass(), gibbs_mole());
            s += p;
            sprintf(p, " heat capacity c_p    %12.6g     %12.4g     J/K\n",
                    cp_mass(), cp_mole());
            s += p;
            try {
                sprintf(p, " heat capacity c_v    %12.6g     %12.4g     J/K\n",
                        cv_mass(), cv_mole());
                s += p;
            } catch (CanteraError& err) {
                err.save();
                sprintf(p, " heat capacity c_v    <not implemented>       \n");
                s += p;
            }
        }

        sprintf(p, " \n");
        s += p;
        if (show_thermo) {
            sprintf(p, "                           X        "
                    "   Molalities         Chem.Pot.    ChemPotSS    ActCoeffMolal\n");
            s += p;
            sprintf(p, "                                    "
                    "                      (J/kmol)      (J/kmol)                 \n");
            s += p;
            sprintf(p, "                     -------------  "
                    "  ------------     ------------  ------------    ------------\n");
            s += p;
            for (size_t k = 0; k < kk; k++) {
                if (x[k] > SmallNumber) {
                    sprintf(p, "%18s  %12.6g     %12.6g     %12.6g   %12.6g   %12.6g\n",
                            speciesName(k).c_str(), x[k], molal[k], mu[k], muss[k], acMolal[k]);
                } else {
                    sprintf(p, "%18s  %12.6g     %12.6g          N/A      %12.6g   %12.6g \n",
                            speciesName(k).c_str(), x[k], molal[k], muss[k], acMolal[k]);
//.........这里部分代码省略.........

void SingleSpeciesTP::getStandardVolumes(doublereal* vbar) const
{
    vbar[0] = molecularWeight(0) / density();
}

double run_walls(double eta, const char *name, Functional fhs, double teff) {
  //printf("Filling fraction is %g with functional %s at temperature %g\n", eta, name, teff);
  //fflush(stdout);
  double kT = teff;
  if (kT == 0) kT = 1;

  Functional f = OfEffectivePotential(fhs);

  const double zmax = width + 2*spacing;
  Lattice lat(Cartesian(dw,0,0), Cartesian(0,dw,0), Cartesian(0,0,zmax));
  GridDescription gd(lat, dx);

  Grid constraint(gd);
  constraint.Set(notinwall);
  f = constrain(constraint, f);

  // We reuse the potential, which should give us a better starting
  // guess on each calculation.
  static Grid *potential = 0;
  if (strcmp(name, "hard") == 0) {
    // start over for each potential
    delete potential;
    potential = 0;
  }
  if (!potential) {
    potential = new Grid(gd);
    *potential = (eta*constraint + 1e-4*eta*VectorXd::Ones(gd.NxNyNz))/(4*M_PI/3);
    *potential = -kT*potential->cwise().log();
  }

  // FIXME below I use the HS energy because of issues with the actual
  // functional.
  const double approx_energy = fhs(kT, eta/(4*M_PI/3))*dw*dw*width;
  const double precision = fabs(approx_energy*1e-5);
  printf("\tMinimizing to %g absolute precision from %g from %g...\n", precision, approx_energy, kT);
  fflush(stdout);

  Minimizer min = Precision(precision,
                            PreconditionedConjugateGradient(f, gd, kT,
                                                            potential,
                                                            QuadraticLineMinimizer));
  took("Setting up the variables");
  for (int i=0;min.improve_energy(false) && i<100;i++) {
  }
  took("Doing the minimization");
  min.print_info();

  Grid density(gd, EffectivePotentialToDensity()(kT, gd, *potential));
  //printf("# per area is %g at filling fraction %g\n", density.sum()*gd.dvolume/dw/dw, eta);
  
  char *plotname = (char *)malloc(1024);

  sprintf(plotname, "papers/fuzzy-fmt/figs/walls%s-%06.4f-%04.2f.dat", name, teff, eta);
  z_plot(plotname, Grid(gd, 4*M_PI*density/3));
  free(plotname);

  {
    GridDescription gdj = density.description(); 
    double sep =  gdj.dz*gdj.Lat.a3().norm();
    int div = gdj.Nz;
    int mid = int (div/2.0);
    double Ntot_per_A = 0;
    double mydist = 0;
   
    for (int j=0; j<mid; j++){
      Ntot_per_A += density(0,0,j)*sep;
      mydist += sep;
    }
    
    double Extra_per_A = Ntot_per_A - eta/(4.0/3.0*M_PI)*width/2;
    
    FILE *fout = fopen("papers/fuzzy-fmt/figs/wallsfillingfracInfo.txt", "a");
    fprintf(fout, "walls%s-%04.2f.dat  -  If you want to match the bulk filling fraction of figs/walls%s-%04.2f.dat, than the number of extra spheres per area to add is %04.10f.  So you'll want to multiply %04.2f by your cavity volume and divide by (4/3)pi.  Then add %04.10f times the Area of your cavity to this number\n",
	    name, eta, name, eta, Extra_per_A, eta, Extra_per_A);
    
    int wallslen = 20;
    double Extra_spheres =  (eta*wallslen*wallslen*wallslen/(4*M_PI/3) + Extra_per_A*wallslen*wallslen);  
    fprintf (fout, "For filling fraction %04.02f and walls of length %d you'll want to use %.0f spheres.\n\n", eta, wallslen, Extra_spheres);
    
    fclose(fout); 
  }
  
  {
    //double peak = peak_memory()/1024.0/1024;
    //double current = current_memory()/1024.0/1024;
    //printf("Peak memory use is %g M (current is %g M)\n", peak, current);
  }
  
  took("Plotting stuff");
  printf("density %g gives ff %g for eta = %g and T = %g\n", density(0,0,gd.Nz/2),
         density(0,0,gd.Nz/2)*4*M_PI/3, eta, teff);
  return density(0, 0, gd.Nz/2)*4*M_PI/3; // return bulk filling fraction
}