Blog Entry © Friday, October 23, 2026, by James Pate Williams, Jr. Equation of State Calculations

// EquationOfState.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Friday, October 23, 2026
// by James Pate Williams, Jr. BA, BS, MSwE, PhD
// Reference: osti.gov/servlets/purl/4390578/			

#include <math.h>
#include <iostream>
#include <vector>
using namespace std;

const int N = 224;
const int nRows = 14;
const int nCols = 16;
const double d = 1.0 / nRows;

typedef struct tagPoint2d
{
    double x, y;
} Point2d, *PPoint2d;

double d0 = 0.0, AA0 = 0.0, AA2 = 0.0, maxD = 0.0, pi = 0.0;
int numberSteps = 0;
vector<double> Nm(65);
vector<Point2d> particlePositions(N);

static double Distance(double x1, double y1, double x2, double y2)
{
    // 2-D Euclidean norm
    return sqrt(pow(x1 - x2, 2.0) + pow(y1 - y2, 2.0));
}

static double Random01()
{
    // random number in the interval [0, 1]
    double value = static_cast<double>(rand()) / RAND_MAX;
    return value;
}

static double Random11()
{
    // random number in the interval [-1, 1]
    double value = -1.0 + 2.0 * static_cast<double>(rand()) / RAND_MAX;
    return value;
}

static void InitializeData(double K, double nu)
{
    pi = 4.0 * atan(1.0);
    d0 = d * (1.0 - pow(2.0, nu - 8.0));
    AA0 = 1.0 / (sqrt(3.0) * d0 * d0 * N / 2.0);
    maxD = d - d0;
    AA2 = (K * K - 1.0) * pi * d0 * d0 / 64.0;
}

static void InitializeLattice()
{
    // initialize to a regular 2-D Cartesian lattice
    for (size_t r = 0; r < nRows; r++)
    {
        for (size_t c = 0; c < nCols; c++)
        {
            size_t address = r * nRows + c;
            particlePositions[address].x = r / 14.0;
            particlePositions[address].y = c / 16.0;
        }
    }
}

static int Move()
{
    double r = 1.0 / nRows;
    double c = 1.0 / nCols;
    int collisions = 0;

    for (int n = 0; !collisions && n < N; n++)
    {
        double xj = particlePositions[n].x;
        double yj = particlePositions[n].y;
        double xk = 100.0 * r * Random11() + xj;
        double yk = 100.0 * c * Random11() + yj;
        double rr =
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0);

        if (rr < d)
            collisions++;

        if (!collisions)
        {
            particlePositions[n].x = xk;
            particlePositions[n].y = yk;
        }
    }

    return collisions;
}

static double A33()
{
    double a33 = 0.0;
    
    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0, f31 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 3 + 1;
        int k = rand() % 3 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a33 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a33 += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a33 += 1.0;
    }

    a33 /= (pi * pi * pow(d, 4.0));
    a33 /= numberSteps;
    return a33;
}

static double A44()
{
    double a44 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 4 + 1;
        int k = rand() % 4 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a44 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a44 += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a44 += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a44 += 1.0;
    }

    a44 /= (pi * pi * pow(d, 4.0));
    a44 /= numberSteps;
    return a44;
}

static double A45()
{
    double a45 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f13 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 4 + 1;
        int k = rand() % 4 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a45 += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a45 += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a45 += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a45 += 1.0;
        if (j == 1 && k == 3 && rjk < d)
            a45 += 1.0;
    }

    a45 /= (pi * pi * pow(d, 4.0));
    a45 /= numberSteps;
    return a45;
}

static double A46()
{
    double a46 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f13 = 0.0, f24 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 4 + 1;
        int k = rand() % 4 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a46 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a46 += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a46 += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a46 += 1.0;
        else if (j == 1 && k == 3 && rjk < d)
            a46 += 1.0;
        else if (j == 2 && k == 4 && rjk < d)
            a46 += 1.0;
    }

    a46 /= (pi * pi * pow(d, 4.0));
    a46 /= numberSteps;
    return a46;
}

static double A55()
{
    double a55 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f24 = 0.0, f31 = 0.0;
        double f14 = 0.0, f34 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a55 += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a55 += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a55 += 1.0;
        if (j == 4 && 1 == 5 && rjk < d)
            a55 += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a55 += 1.0;
    }

    a55 /= (pi * pi * pow(d, 4.0));
    a55 /= numberSteps;
    return a55;
}

static double A56p()
{
    double a55p = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f24 = 0.0, f31 = 0.0;
        double f14 = 0.0, f34 = 0.0;
        double f53 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a55p += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a55p += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a55p += 1.0;
        if (j == 4 && 1 == 5 && rjk < d)
            a55p += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a55p += 1.0;
    }

    a55p /= (pi * pi * pow(d, 4.0));
    a55p /= numberSteps;
    return a55p;
}

static double A56pp()
{
    double a55pp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f35 = 0.0, f51 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a55pp += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a55pp += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a55pp += 1.0;
        if (j == 3 && 1 == 5 && rjk < d)
            a55pp += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a55pp += 1.0;
    }

    a55pp /= (pi * pi * pow(d, 4.0));
    a55pp /= numberSteps;
    return a55pp;
}

static double A57p()
{
    double a57p = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f35 = 0.0;
        double f45 = 0.0, f51 = 0.0;
        double f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a57p += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a57p += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a57p += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a57p += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a57p += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a57p += 1.0;
    }

    a57p /= (pi * pi * pow(d, 4.0));
    a57p /= numberSteps;
    return a57p;
}

static double A57pp()
{
    double a57pp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a57pp += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a57pp += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a57pp += 1.0;
        if (j == 4 && k == 5 && rjk < d)
            a57pp += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a57pp += 1.0;
        if (j == 3 && k == 1 && rjk < d)
            a57pp += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a57pp += 1.0;
    }

    a57pp /= (pi * pi * pow(d, 4.0));
    a57pp /= numberSteps;
    return a57pp;
}

static double A57ppp()
{
    double a57ppp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f35 = 0.0, f51 = 0.0;
        double f31 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a57ppp += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a57ppp += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a57ppp += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a57ppp += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a57ppp += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a57ppp += 1.0;
        if (j == 3 && k == 1 && rjk < d)
            a57ppp += 1.0;
    }

    a57ppp /= (pi * pi * pow(d, 4.0));
    a57ppp /= numberSteps;
    return a57ppp;
}

static double A58p()
{
    double a58p = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f35 = 0.0, f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a58p += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a58p += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a58p += 1.0;
        if (j == 4 && k == 5 && rjk < d)
            a58p += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a58p += 1.0;
        if (j == 3 && k == 1 && rjk < d)
            a58p += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a58p += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a58p += 1.0;
    }

    a58p /= (pi * pi * pow(d, 4.0));
    a58p /= numberSteps;
    return a58p;
}

static double A58pp()
{
    double a58pp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0, f25 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a58pp += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a58pp += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a58pp += 1.0;
        if (j == 4 && k == 5 && rjk < d)
            a58pp += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a58pp += 1.0;
        if (j == 3 && k == 1 && rjk < d)
            a58pp += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a58pp += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a58pp += 1.0;
        if (j == 2 && k == 5 && rjk < d)
            a58pp += 1.0;
    }

    a58pp /= (pi * pi * pow(d, 4.0));
    a58pp /= numberSteps;
    return a58pp;
}

static double A59()
{
    double a59 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0, f25 = 0.0;
        double f35 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a59 += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a59 += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a59 += 1.0;
        if (j == 4 && k == 5 && rjk < d)
            a59 += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a59 += 1.0;
        if (j == 3 && k == 1 && rjk < d)
            a59 += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a59 += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a59 += 1.0;
        if (j == 2 && k == 5 && rjk < d)
            a59 += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a59 += 1.0;
    }

    a59 /= (pi * pi * pow(d, 4.0));
    a59 /= numberSteps;
    return a59;
}

static double A510()
{
    double a510 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0, f25 = 0.0;
        double f35 = 0.0, f24 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a510 += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a510 += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a510 += 1.0;
        if (j == 4 && k == 5 && rjk < d)
            a510 += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a510 += 1.0;
        if (j == 3 && k == 1 && rjk < d)
            a510 += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a510 += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a510 += 1.0;
        if (j == 2 && k == 5 && rjk < d)
            a510 += 1.0;
        if (j == 2 && k == 4 && rjk < d)
            a510 += 1.0;
    }

    a510 /= (pi * pi * pow(d, 4.0));
    a510 /= numberSteps;
    return a510;
}

int main()
{
    double K = 0.0, nu = 0.0;
    int N12 = 0;

    cout << "K  = ";
    cin >> K;
    cout << "nu = ";
    cin >> nu;
 
    InitializeData(K, nu);

    for (int i = 0; i < 10000; i++)
    {
        if (Move() == 0)
            N12++;
    }
    
    cout << "N(1/2) = ";
    cout << N12 << endl;
    
    double C0 = 64.0 * N12 / (N * N * (K * K - 1.0));
    double A0 = 0.5 * sqrt(3.0) * d0 * d0 * N;
    double a33 = 0.0, a44 = 0.0, a45 = 0.0, a46 = 0.0, a55 = 0.0;
    double a56p = 0.0, a56pp = 0.0, a57p = 0.0, a57pp = 0.0;
    double a57ppp = 0.0, a58p = 0.0, a58pp = 0.0;
    double a59 = 0.0, a510 = 0.0;

    cout << "C0 = " << C0 << endl;
    cout << "A0 = " << A0 << '\t';
    cout << "1 / A0 = " << 1.0 / A0 << endl;

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a33 += A33();     
    
    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a44 += A44();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a45 += A45();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a46 += A46();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a55 += A55();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a56p += A56p();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a56pp += A56pp();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a57p += A57p();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a57pp += A57pp();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a57ppp += A57ppp();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a58p += A58p();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a58pp += A58pp();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a59 += A59();

    for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
        a510 += A510();

    a33 /= 4.0;
    a44 /= 4.0;
    a45 /= 4.0;
    a46 /= 4.0;
    a55 /= 4.0;
    a56p /= 4.0;
    a56pp /= 4.0;
    a57p /= 4.0;
    a57pp /= 4.0;
    a57ppp /= 4.0;
    a58p /= 4.0;
    a58pp /= 4.0;
    a59 /= 4.0;
    a510 /= 4.0;

    cout << "A(3, 3)    = " << a33 << endl;
    cout << "A(4, 4)    = " << a44 << endl;
    cout << "A(4, 5)    = " << a45 << endl;
    cout << "A(4, 6)    = " << a46 << endl;
    cout << "A(5, 5)    = " << a55 << endl;
    cout << "A(5, 6')   = " << a56p << endl;
    cout << "A(5, 6'')  = " << a56pp << endl;
    cout << "A(5, 7')   = " << a57p << endl;
    cout << "A(5, 7'')  = " << a57pp << endl;
    cout << "A(5, 7''') = " << a57ppp << endl;
    cout << "A(5, 8')   = " << a58p << endl;
    cout << "A(5, 8'')  = " << a58pp << endl;
    cout << "A(5, 9)    = " << a59 << endl;
    cout << "A(5, 10)   = " << a510 << endl;

    double C1 = pi / sqrt(3.0);
    double C2 = N12 / ((4.0 * pi * pi * A33()) / 9.0);
    double C3a = (6.0 * A45() - 3.0 * A44() - A46()) / pow(3.0, 1.5);
    double C3b = N12 / (pi * pi * pi * C3a);
    double C40 = 8.0 * pi * pi * pi / 135.0;
    double C4a = 12.0 * a55 - 60.0 * a56p - 10.0 * a56pp;
    double C4b = 30.0 * a57p + 60.0 * a57pp + 10.0 * a57ppp;
    double C4c = -30.0 * a58p - 15.0 * a58pp;
    double C4d = 10.0 * a59 - a510;
    double C4 = C4a + C4b + C4c + C4d;
    C4 = C40 * N12 / C4;

    cout << "C(1) = " << C1 << endl;
    cout << "C(2) = " << 2.0 * C2 << endl;
    cout << "C(3) = " << 9.0 * C3b << endl;
    cout << "C(4) = " << C4 << endl;

    double X3 = C1 * A0 + 2.0 * C2 * A0 * A0 + 9.0 * C3b * A0 * A0 * A0;
    X3 += C4 * pow(A0, 4.0);
    X3 += 1.0;

    cout << "X(3) = " << X3 << endl;
    return 0;
}