/* ---------------------------------------------------------------------- 
**  Figure 19.8:  A heat flow simulation.
**  Figure 19.9:  Parameters for the simulation.
**  Figure 19.10: Iteration for the heat-flow simulation, sim_cool().
**  Figure 19.11: Print function for the heat-flow simulation.
** ----------------------------------------------------------------------
*/
#include "tools.h"

#define MAX_STEPS 1000
#define DEPTH_LIMIT .25

void setup( double* p_deep, double* p_init, double* p_last, double* p_surf );
double* sim_cool( int slot, double initial, double goal, double surf, int* actual);  
void print_array( double nodes[], int last );

void main( void )
{   
    double init_temp, last_temp; /* Beginning and ending temperatures. */
    double surf;                 /* Temperature at surface of slab. */
    double depth;                /* Depth at which to monitor temperature. */   
    int slot;                    /* Array slot at given depth. */

    int steps, seconds;          /* Simulation steps, time to reach goal. */
    double* nodes;               /* Pointer to array holding final temperatures. */   
    
    banner();     
    puts("\nSimulation of Heat Conduction in a Slab\n");
    setup( &depth, &init_temp, &last_temp, &surf );
    slot = ceil( (depth - 0.003) / 0.006 );
        
    printf("\nGoal temperature = %g at node %i\n\n", last_temp, slot);
    nodes = sim_cool( slot, init_temp, last_temp, surf, &steps );

    seconds = 30 * steps;
    printf( "\nTemperature of \%g reached at node \%i \n"  
            "\tin \%i seconds (= \%.1f minutes or \%.2f hours).\n",   
            nodes[slot], slot, seconds, seconds/60.0, seconds/3600.0 );
    printf("\nFinal nodal temperatures after %i steps:\n    ", steps);
    print_array( nodes, steps ); 
    free( nodes );
    bye();
}
/* ---------------------------------------------------------------------- */
void 
setup( double* p_deep, double* p_init, double* p_last, double* p_surf )
{
    double deep;                /* Local variable for depth (input) */
    do {
        printf( "Enter depth to monitor (up to %g meters): ", DEPTH_LIMIT );
        scanf( "%lg", &deep );
    } while (deep<0 || deep > DEPTH_LIMIT);
    *p_deep = deep; 

    printf( "Enter initial temperature of the slab: " );
    scanf( "%lg", p_init );
    
    printf( "Enter target temperature for depth %g: ", *p_deep );
    scanf( "%lg", p_last );

    printf( "Enter temperature of cold outer surface: " );
    scanf( "%lg", p_surf );
}

/* ---------------------------------------------------------------------- */
double*                     /* Pointer to array of nodes holding answers. */
sim_cool(  int slot,        /* Number of node being monitored.            */ 
           double initial,  /* Initial temperature of the infinite slab.  */ 
           double goal,     /* Stop when chosen slot <= goal temp.        */ 
           double surf,     /* Temperature of the cold outer surface.     */ 
           int* actual )    /* Actual number of iterations carried out.   */
{   
    int n = 3 * slot;       /* Current allocation length.  */
    int p;                  /* Counter for time steps.     */ 
    int k;                  /* Index for node array.       */ 
    double* swap;           /* To swap old and new arrays. */

    double* old = (double*) malloc( n * sizeof(double) ); 
    double* new = (double*) malloc( n * sizeof(double) );
    
    new[0] = old[0] = surf;   
    for (k = 1; k < n; ++k) new[k] = old[k] = initial;

    for (p = 1; p <= MAX_STEPS; ++p) {   
        if (p == n - 1) {   /* If no neighbor exists on right, make one */
            n *= 2;
            new = realloc (new, n * sizeof(double) );
            old = realloc (old, n * sizeof(double) );
            for (k = p + 1; k < n; ++k) new[k] = old[k] = initial;
        } 
        for (k = 1; k <= p; ++k) new[k] = (old[k-1] + old[k+1]) / 2.0;
        printf( "%3i.", p );  
        print_array( new, slot ); 

        if( new[slot] <= goal ) break;
        swap = old; old = new; new = swap; 
    }
    *actual = p;       /* Actual number of iterations executed. */
    free( old );
    return new;        /* This array holds results when goal is reached. */
}
/* ---------------------------------------------------------------------- */
void 
print_array( double nodes[], int last )
{
    int k;
    for (k = 0; k <= last; ++k) {
        printf( "[%2i] = %6.2f  ", k, nodes[k] ); 
        if (k % 5 == 4) printf( "\n    " ); /* Newline after five items. */
    }
    printf( "\n\n" );
}