#include <cassert>
#include <cmath>
#include <cstddef>
#include <iostream>

// We will simulate heat diffusion through a material.
//
// Imagine a bar sampled at equally spaced points.
// An array stores the temperature at each of these
// points.
//
// Given the current temperatures, we will estimate
// the temperatures one second into the future, and
// then repeat this process.
//
// The function update_temperatures(...) will do the
// following:
//
//   1. Create a new array next_temperature[].
//
//   2. Set the first and last entries to the fixed
//      temperatures
//
//          temp_left
//
//      and
//
//          temp_right,
//
//      respectively.
//
//   3. Set each remaining entry using
//
//          next_temperature[k]
//              = temperature[k]
//                + 0.5*(
//                    temperature[k - 1]
//                    - 2*temperature[k]
//                    + temperature[k + 1]
//                );
//
//   4. Once all the new temperatures have been
//      calculated, copy them back into temperature[].
//
// Why do we need a second array? Why can't we simply
// update the entries of temperature[] one at a time?
//
// You can experiment with replacing 0.5 by smaller
// positive values. What happens?
//
// If you replace 0.5 by a value greater than 0.5,
// the simulation may become unstable.
//
// Question:
//
//   Does the change in temperature behave as you
//   might expect?

// Function declarations
void update_temperatures(
    double            temperature[],
    std::size_t const capacity,
    double      const temp_left,
    double      const temp_right
);

int main();

// Function definitions
int main() {
    std::size_t const N{ 20 };

    // We will start with a bar having a uniform
    // temperature of 50 degrees Celsius.
    //
    // This array represents the bar sampled at
    // 20 equally spaced points along its length.
    double temperature[N]{
        50, 50, 50, 50, 50, 50, 50, 50, 50, 50,
        50, 50, 50, 50, 50, 50, 50, 50, 50, 50
    };

    for ( std::size_t time{ 0 }; time <= 200; ++time ) {
        std::cout << time << "\t";

        for ( std::size_t k{ 0 }; k < N; ++k ) {
            if ( temperature[k] < 9.5 ) {
                std::cout << " "
                          << std::round( temperature[k] )
                          << " ";
            } else {
                std::cout << std::round( temperature[k] )
                          << " ";
            }
        }

        std::cout << std::endl;

        // The left-hand end of the bar is touching
        // ice and is therefore maintained at 0 degrees
        // Celsius. The right-hand end is maintained
        // at 23 degrees Celsius.
        update_temperatures( temperature, N, 0, 23 );
    }

    return 0;
}

void update_temperatures(
    double            temperature[],
    std::size_t const capacity,
    double      const temp_left,
    double      const temp_right
) {
    assert( capacity >= 2 );

    // Enter your implementation here.
}












































// Why must we calculate all the new temperatures
// before changing any entries in temperature[]?
//
// Consider what would happen if we immediately changed
// temperature[1]. When calculating temperature[2], we
// would then be using the new value of temperature[1]
// together with the old value of temperature[2] and
// temperature[3].
//
// We want every new temperature to be calculated using
// temperatures from the same instant in time.
//
// Consequently, we first calculate all the new values
// and only then copy them back.












































// There is another interesting observation.
//
// Our formula is
//
//     next_temperature[k]
//         = temperature[k]
//           + 0.5*(
//               temperature[k - 1]
//               - 2*temperature[k]
//               + temperature[k + 1]
//           );
//
// If you simplify this expression, what happens to
// temperature[k]?
//
// You should find that
//
//     next_temperature[k]
//         = 0.5*(
//             temperature[k - 1]
//             + temperature[k + 1]
//         );
//
// Thus, with the coefficient 0.5, the new temperature
// at each interior point is simply the average of the
// temperatures at its two neighboring points.












































// One solution is here:

/*
void update_temperatures(
    double            temperature[],
    std::size_t const capacity,
    double      const temp_left,
    double      const temp_right
) {
    assert( capacity >= 2 );

    std::size_t const MAXIMUM_CAPACITY{ 20 };
    assert( capacity <= MAXIMUM_CAPACITY );

    double next_temperature[MAXIMUM_CAPACITY];

    next_temperature[0] = temp_left;
    next_temperature[capacity - 1] = temp_right;

    for ( std::size_t k{ 1 }; k < capacity - 1; ++k ) {
        next_temperature[k] = temperature[k] + 0.5*(
            temperature[k - 1]
            - 2*temperature[k]
            + temperature[k + 1]
        );
    }

    for ( std::size_t k{ 0 }; k < capacity; ++k ) {
        temperature[k] = next_temperature[k];
    }
}
*/