// Simple 1D Transient Heat Transfer (parameters a function of temperature)
// Type:	Conduction only
// Shape:	Slab
// BC:		Fixed Temperature (infinite heat sinks, "Dirichlet")
// Params:	Variable! Thermal conductivity, density, & heat capacity all can be a function of temperature, time, and position!!!
// Method:	Explicit (small memory version)
// By:		Wayne Pafko
// Created:	June 13th, 2001
// Updated:	June 25th, 2001

#include <stdio.h>			// The standard C library
#include <time.h>			// How fast is this program? 
#define SUCCESS			1
#define FAILURE			0
#define SMALL			1e-6	// A small number
#define ASCII_SPACE		32		// ASCII decimal value for a space
#define ASCII_TAB		9		// ASCII decimal value for a tab

// Allows the use of "excel" style formulas to define K, CP, & DENSITY (just cut and paste from excel to formula's below)
#define IF(A,B,C)		((A)?(B):(C))

// Solution specific definitions (change these & recompile)
#define NODES	101		// How many nodes do we have (no decimal point please)
#define STEPS	50000	// How many time steps do we have (no decimal point please)
#define	TIME	1000.0	// (sec) Total time
#define WIDTH	10.0	// (cm) Total width of simulation

// Functional relationships (change these & recompile) Excel friendly (paste back and forth)!!!
// In general, each parameter is a function of temperature - T_C (C), time - t_s (s), and position - x_cm (cm)
// Temperature range:	unknown (don't know until we get the solution)
// Time range:			0 s < t < "TIME" s
// Position range:		0 cm < x < "WIDTH" cm.  (Center node is at "WIDTH"/2 * cm)

#define K(T_C,t_s,x_cm)	 (IF((T_C)>=0,IF((T_C)<=100,0.5811+9.888E-4*(T_C),0.467),0.569)+.001*(t_s)+.1*(WIDTH-(x_cm)))	
// (W/m/C) Thermal conductivity

#define CP(T_C,t_s,x_cm)  (IF((T_C)>=0, IF((T_C)<=100, 4.194, 1.827+5.507E-04*(T_C)), 2.098 + 0.007129*(T_C))+.001*(t_s)+.1*(WIDTH-(x_cm)))
// (kJ/kg/C) Heat Capacity

#define DENSITY(T_C,t_s,x_cm)  (IF((T_C)>=0, IF((T_C)<=100, 1.001 - 6.441E-5*(T_C) - 3.622E-6*(T_C)*(T_C), 7.494E-04-1.739E-06*(T_C)+1.700E-09*(T_C)*(T_C)), .8)+.001*(t_s)+.1*(WIDTH-(x_cm)))
// (g/cm3) Density
// note: We use macros instead of "functions" because they are much faster!
// Important! Leave extra parenthesis around the individual variables, otherwise complex expressions won't be inserted into equation correctly!!!

// Problem specific definitions (change these & recompile)
#define INITIAL_TEMP	125.0		// (C) Initial temperature across slab
#define BC_LEFT			200.0			// (C) Left BC
#define BC_RIGHT		50.0		// (C) Right BC
#define TREF			25.0		// (C) Reference temperature (used to calculate enthalpy)
double Tmin_predict	=	0.0;		// (C) Predicted minimum temperature.  Should solution converge?  Also effects parameter output range.
double Tmax_predict	=	200.0;		// (C) Predicted maximum temperature.  Should solution converge?  Also effects parameter output range.

// Unit conversion factors
#define FOURIER_FACTOR				1		// [dimless: (cm2/s) * (s) / (cm2)] Conversion constant for equation k/cp/density*dt/dx/dx
#define THERMAL_DIFFUSIVITY_FACTOR	0.01	// [cm2/s: (W/m/C) / (kJ/kg/C) / (g/cm3)] Conversion constant for equation k/cp/density
#define ENTHALPY_FACTOR				1.0		// [kJ/kg: (kJ/kg/C) * (C)] Conversion constant for equation cp*(T-Tref)
#define M_OVER_CM					0.01	// [dimless: (1 m) / (100 cm)] Conversion constant from meters to cm

// Parameters that effect how the results are saved (but not the solution method)
#define FILE_NAME		"data.txt"	// Save results in this file
#define PRINTSTEP		500			// How many time steps before a result is printed (if "1" see data at every time step, if "2" see every other step)
#define PRINTNODE		10			// How many nodes before a result is printed (if "1" see all the nodes, if "2" see half)
#define PRINTNODESTART	1			// What node should I start printing at?  (usually 1, always >=1)

// Derived values based on above definitions
double dt		= (double)TIME/STEPS;	// (sec) Time step
double dx		= (double)WIDTH/NODES;	// (cm) Distance between nodes
double dtdx2	= (double)TIME/STEPS * NODES/WIDTH * NODES/WIDTH;  // So don't have to keep recalculating (dt/dx/dx) at each time step

// Memory stuff
double storage_a[NODES+1];	// Storage for algorithm; "NODES+1" for safety reasons (don't want to fly past alotted space accidently)
double storage_b[NODES+1];	// More storage for algoritm

// Other useful numbers that will be calculated (leave these alone) 
//double Temp[STEPS][NODES];	// (C) The Solution (Temperature over time for each node, can get very big!)  
double * Temp_now;			// (C) Current Temperature (this time step) (will point to storage_a or storage_b)
double * Temp_past;			// (C) Past Temperature (one time step ago) (will point to storage_a or storage_b)
double aveTemp[STEPS];		// (C) Average temperature of the slab
double fourierDtDx2;		// (dimless) Fourier number for a single step & node (must be small for convergence)
//double fourier[STEPS];		// (dimless) Fourier number (dimensionless time)
double aveH[STEPS];			// (kJ/kg) Average Enthalpy
double TBC_left[STEPS];		// (C) Left BC (may change with time)
double TBC_right[STEPS];	// (C) Right BC (may change with time)
double Tmin;				// (C) Minimum temperature found during solution (Same as TMIN_PREDICT?)
double Tmax;				// (C) Maximum temperature found during solution (Same as TMAX_PREDICT?)

// Time the algorithm
time_t	startTime, endTime;	// "time_t" structure defined in time.h

// Sets initial conditions and boundary conditions
int initial_conditions(void)
{
	long i, t; // Loop variable

	// Set up pointers (to storage space)
	// Establish where Temp_now & Temp_past point to...
	Temp_now = & storage_a[0];
	Temp_past = & storage_b[0];

	// Check expected memory & file usage
	printf("Expected sizes:  memory %.2f mb, file %.2f mb\n",1.1*(sizeof(storage_a)+sizeof(storage_b))/1E6, NODES/PRINTNODE*STEPS/PRINTSTEP*8*8*2/1E6);
	// Initial Conditions for interior nodes
	for(i=0;i<NODES;i++)
	{
		Temp_past[i]=(double)INITIAL_TEMP;  // Really changing values in storage_b
	}

	// Store the highest and lowest temperatures (save in Tmax or Tmin)		
	Tmin=Temp_past[0]; // Give some initial value
	Tmax=Temp_past[0]; // Give some initial value
	for(i=0;i<NODES;i++)  // Loop over all nodes
	{
		if(Temp_past[i]<Tmin) Tmin=Temp_now[i]; // New minimum temperature found, save it
		if(Temp_past[i]>Tmax) Tmax=Temp_now[i]; // New maximum temperature found, save it
	}


	// Set boundary conditions (BC)
	for(t=0;t<STEPS;t++)
	{
		TBC_left[t]=(double)BC_LEFT;
		TBC_right[t]=(double)BC_RIGHT;
	}

	printf("Initial conditions set...\n");

	return(SUCCESS);
};

// Code to solve the problem and give some screen output
int solve_and_save_problem(void)
{
	long i,t;					// Loop variables 
	//double thermalDiffusivity;	// Store thermal diffusivity calculations (more readible code)
	double stable=0.5, converge=0.25, minError=0.166;		// Stability, convergence, and minimum truncation error rules of thumb
	double thermalDiffusivity;	// (cm2/s) Thermal Diffusivity (k/cp/density)
	double k_a, k_b, k_c;		// (W/m/C) Thermal Conductivity at nodes "a", "b", & "c" ("b" is the current node, "a" is one to the left, "c" is one to the right
	double k_ba, k_bc;			// (W/m/C) Thermal Conductivity between nodes "b-a," "b-c," & at a Boundary Condition (BC)
	double T_a, T_b, T_c;		// (C) Temperatures at nodes "a," "b," & "c."
	double x_a, x_b, x_c;		// (cm) Locations of nodes "a," "b," & "c."
	double cp_b;				// (kJ/kg/C) Heat Capacity of node "b"
	double density_b;			// (g/cm3) Density of node "b"
	double q_ba, q_bc, q_net;	// (W/cm2) Heat flux between nodes "b-a," "b-c," & net into node "b"
	double deltaT;				// (C) Temperature increase for node during this timestep
	double Tloop, xloop;		// Loop variable
	FILE *dataFile;				// File input/output

	// PC clock time at start of this function
	time(&startTime);

	// Calculate a few important numbers
	// Estimate this in the middle (ave predicted temp, center node, half full simulation time)
	// Note "\" just means the expression is continued on the next line...
	thermalDiffusivity = K((Tmax_predict+Tmin_predict)/2.0,(float)TIME * dt/2.0,WIDTH/2.0) \
		/ CP((Tmax_predict+Tmin_predict)/2.0,(float)TIME * dt/2.0,WIDTH/2.0) \
		/ DENSITY((Tmax_predict+Tmin_predict)/2.0,(float)TIME * dt/2.0,WIDTH/2.0) * THERMAL_DIFFUSIVITY_FACTOR;
	fourierDtDx2 =  thermalDiffusivity * dtdx2 * (double)FOURIER_FACTOR;
	
	// Print solution properties
	printf("\n---Nodal Representation---\n");
	printf("  Nodes:\t%d\n",(long)NODES);
	printf("  Slab Width:\t%.4f cm\n",(double)WIDTH);
	//printf("  Spacing:\t%.6f cm\n",dx);
	printf("  Steps:\t%d\n",(long)STEPS);
	printf("  Total Time:\t%.2f s\n",(double)TIME);
	//printf("  Time Step:\t%.6f s\n",dt);
	printf("  Fourier #:\t%.3f [must be less than %.3f]\n",fourierDtDx2,stable);
	if(fourierDtDx2<stable) 
	{
		printf("\t\tOK! The solution should converge...\n");
		if(fourierDtDx2>converge) printf("\t\tbut it may oscillate if above %.3f!\n",converge);
		printf("\t\tnote: Truncation error minimized below %.3f\n",minError);

	}
	else 
	{
		printf("\t\tERROR! The solution will not converge!\n\t\tIncrease # of time steps (STEPS) and recompile...\n");
		return(FAILURE);
	}
	
	// Print material properties
	printf("\n---Material Properties---\n");
	printf("  Thermal Conductivity:\t%.4f W/m/C\n",(double)K(Tmax_predict,(float)TIME * dt,0.0));
	printf("  Heat Capacity:\t%.4f kJ/kg/C\n",(double)CP(Tmax_predict,(float)TIME * dt,0.0));
	printf("  Density:\t\t%.4f g/cm3\n",(double)DENSITY(Tmax_predict,(float)TIME * dt,0.0));
	//printf("  Thermal Diffusivity:\t%.6f cm2/s\n",thermalDiffusivity);

	// Print some stuff to the screen
	//printf("\nDoing all the math...\t\t");

//SAVE TO FILE STUFF/////
	printf("\nDoing math and saving in %s...\t    ",FILE_NAME);

	// Open the file (check for failure)
	if( (dataFile = fopen(FILE_NAME,"w"))==FAILURE )
	{
		printf("[failure]\nIs %s open in another program?\n",FILE_NAME);
		return(FAILURE);
	}

	// Start printing to the file
	fprintf(dataFile, "1D Heat Transfer\nby Wayne Pafko\n");

	//Print solution properties
	fprintf(dataFile, "\nNodal Representation\n");
	fprintf(dataFile, "Nodes:\t%d\n",(long)NODES);
	fprintf(dataFile, "Slab Width:\t%.4g\tcm\n",(double)WIDTH);
	fprintf(dataFile, "Spacing:\t%.6g\tcm\n",dx);
	fprintf(dataFile, "Steps:\t%d\n",(long)STEPS);
	fprintf(dataFile, "Total Time:\t%.2f\ts\n",(double)TIME);
	fprintf(dataFile, "Time Step:\t%.6g\ts\n",dt);

	// Calculate & print material properties
	fprintf(dataFile, "\nMaterial Properties\n");

	// Print table of thermal conductivity
	fprintf(dataFile, "Thermal Conductivity (W/m/C)\n");
	fprintf(dataFile, "(C)\t%.4f\t%.4f\t%.4f\t:(cm)\n",0.5*dx,(float)WIDTH/2.0,(float)WIDTH-0.5*dx);
	for(Tloop=Tmin_predict;Tloop<=Tmax_predict+SMALL;Tloop+=(Tmax_predict-Tmin_predict)/4.0)
	{
		fprintf(dataFile, "%.2f\t", Tloop);
		for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
			fprintf(dataFile, "%.4f\t", K(Tloop,0.0,xloop) );
		fprintf(dataFile, "\n");
	}
	fprintf(dataFile, "%.2f\t", Tmax_predict);
	for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
		fprintf(dataFile, "%.4f\t", K(Tmax_predict,(double)TIME,xloop) );
	fprintf(dataFile, "%.2f (s)\n", (double)TIME);

	// Print table of heat capacity
	fprintf(dataFile, "Heat Capacity (kJ/kg/C)\n");
	fprintf(dataFile, "(C)\t%.4f\t%.4f\t%.4f\t:(cm)\n",0.5*dx,(float)WIDTH/2.0,(float)WIDTH-0.5*dx);
	for(Tloop=Tmin_predict;Tloop<=Tmax_predict+SMALL;Tloop+=(Tmax_predict-Tmin_predict)/4.0)
	{
		fprintf(dataFile, "%.2f\t", Tloop);
		for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
			fprintf(dataFile, "%.4f\t", CP(Tloop,0.0,xloop) );
		fprintf(dataFile, "\n");
	}
	fprintf(dataFile, "%.2f\t", Tmax_predict);
	for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
		fprintf(dataFile, "%.4f\t", CP(Tmax_predict,(double)TIME,xloop) );
	fprintf(dataFile, "%.2f (s)\n", (double)TIME);

	// Print table of density
	fprintf(dataFile, "Density (g/cm3)\n");
	fprintf(dataFile, "(C)\t%.4f\t%.4f\t%.4f\t:(cm)\n",0.5*dx,(float)WIDTH/2.0,(float)WIDTH-0.5*dx);
	for(Tloop=Tmin_predict;Tloop<=Tmax_predict+SMALL;Tloop+=(Tmax_predict-Tmin_predict)/4.0)
	{
		fprintf(dataFile, "%.2f\t", Tloop);
		for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
			fprintf(dataFile, "%.4f\t", DENSITY(Tloop,0.0,xloop) );
		fprintf(dataFile, "\n");
	}
	fprintf(dataFile, "%.2f\t", Tmax_predict);
	for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
		fprintf(dataFile, "%.4f\t", DENSITY(Tmax_predict,(double)TIME,xloop) );
	fprintf(dataFile, "%.2f (s)\n", (double)TIME);

	// Print table of thermal diffusivity
	fprintf(dataFile, "Thermal Diffusivity (cm2/s)\n");
	fprintf(dataFile, "(C)\t%.4f\t%.4f\t%.4f\t:(cm)\n",0.5*dx,(float)WIDTH/2.0,(float)WIDTH-0.5*dx);
	for(Tloop=Tmin_predict;Tloop<=Tmax_predict+SMALL;Tloop+=(Tmax_predict-Tmin_predict)/4.0)
	{
		fprintf(dataFile, "%.2f\t", Tloop);
		for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
		{
			thermalDiffusivity=K(Tloop,0.0,xloop)/CP(Tloop,0.0,xloop)/DENSITY(Tloop,0.0,xloop) * THERMAL_DIFFUSIVITY_FACTOR;
			fprintf(dataFile, "%.6f\t", thermalDiffusivity );
		}
		fprintf(dataFile, "\n");
	}
	fprintf(dataFile, "%.2f\t", Tmax_predict);
	for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx)
		{
			thermalDiffusivity=K(Tmax_predict,(double)TIME,xloop)/CP(Tmax_predict,(double)TIME,xloop)/DENSITY(Tmax_predict,(double)TIME,xloop) * THERMAL_DIFFUSIVITY_FACTOR;
			fprintf(dataFile, "%.6f\t", thermalDiffusivity );
		}
	fprintf(dataFile, "%.2f (s)\n", (double)TIME);
	

	// Print headings for main temperature table
	fprintf(dataFile, "\nTemperature Results (C)\n");
	fprintf(dataFile, "\tTref (C):\t%.2f\tx (cm):\t< %.4f\t",(double)TREF,0.0);
	for(i=PRINTNODESTART-1;i<NODES;i+=PRINTNODE) fprintf(dataFile, "%.4f\t",((double)i+0.5)*dx);
	fprintf(dataFile, "> %.4f\n", (double)WIDTH);
	fprintf(dataFile, "step (#)\ttime (s)\tH (kJ/kg)\tAve T (C)\tBC\t");
	for(i=PRINTNODESTART-1;i<NODES;i+=PRINTNODE) fprintf(dataFile, "Node %d\t",i+1);
	fprintf(dataFile, "BC\n");
//DONE SAVING TO FILE FOR NOW/////

	// Now solve the problem
	for(t=0;t<STEPS;t++) // Increment time until end
	{
		// Explicitly calculate the temperature at each node based upon its neighbors one time step ago
		// Special calcualation for edge nodes (use BC)
		// Left node (so "a" is actual outside the nodes...a Boundary Condition)
		T_a=TBC_left[t];	// (C) Temperature at node "a"
		T_b=Temp_past[0];	// (C) Temperature at node "b"
		T_c=Temp_past[0+1];	// (C) Temperature at node "c"

		x_b=(0+0.5)*dx;		// (cm) Location of node "b"
		x_c=(1+0.5)*dx;		// (cm) Location of node "c"

		k_b=(double)K(T_b,(double)t*dt,x_b);	// Current node
		k_c=(double)K(T_c,(double)t*dt,x_c);	// One node to the right (i+1)

		k_ba=k_b;						// (W/m/C) Conduction in series, but no node to the left (take k to be that for node "b")
		k_bc=2*k_b*k_c/(k_b+k_c);		// (W/m/C) Conduction in series, 2/k_bc=1/k_b+1/k_c

		q_ba=k_ba*(T_a-T_b)/dx*(double)M_OVER_CM;	// (W/cm2) Heat flux from "a" to "b." Positive means heat flows into "b"
		q_bc=k_bc*(T_c-T_b)/dx*(double)M_OVER_CM;	// (W/cm2) Heat flux from "b" to "c." Positive means heat flows into "b"
		q_net=q_ba+q_bc;

		cp_b=(double)CP(T_b,(double)t*dt,x_b);				// (kJ/kg/C) Heat Capacity at node "b"
		density_b=(double)DENSITY(T_b,(double)t*dt,x_b);	// (g/cm3) Density at node "b"

		deltaT = q_net/cp_b/density_b*dt/dx;	// (C) Temperature increase for this time step
		Temp_now[0] = T_b + deltaT;				// (C) Temperature for next time step

		//thermalDiffusivity = K(Temp[t][0],(float)t*dt,(0+0.5)*dx)/CP(Temp[t][0],(float)t*dt,(0+0.5)*dx)/DENSITY(Temp[t][0],(float)t*dt,(0+0.5)*dx) * (double)THERMAL_DIFFUSIVITY_FACTOR;
		//Temp[t+1][0] = Temp[t][0] + thermalDiffusivity * dtdx2 * FOURIER_FACTOR * ( Temp[t][1] - 2 * Temp[t][0] + TBC_left[t]);

		// Right node (so "c" is actually outside the nodes...a Boundary Condition) 
		T_a=Temp_past[NODES-1-1];	// (C) Temperature at node "a"
		T_b=Temp_past[NODES-1];	// (C) Temperature at node "b"
		T_c=TBC_right[t];		// (C) Temperature at node "c"

		x_a=((double)(NODES-1-1)+0.5)*dx;		// (cm) Location of node "a"
		x_b=((double)(NODES-1)+0.5)*dx;		// (cm) Location of node "b"
	
		k_a=(double)K(T_a,(double)t*dt,x_a);	// One node to the left (i-1)
		k_b=(double)K(T_b,(double)t*dt,x_b);	// Current node

		k_ba=2*k_b*k_a/(k_b+k_a);		// (W/m/C) Conduction in series, 2/k_ba=1/k_a+1/k_b
		k_bc=k_b;						// (W/m/C) Conduction in series, but no node to the right (take k to be that for node "b")

		q_ba=k_ba*(T_a-T_b)/dx*(double)M_OVER_CM;	// (W/cm2) Heat flux from "a" to "b." Positive means heat flows into "b"
		q_bc=k_bc*(T_c-T_b)/dx*(double)M_OVER_CM;	// (W/cm2) Heat flux from "b" to "c." Positive means heat flows into "b"
		q_net=q_ba+q_bc;

		cp_b=(double)CP(T_b,(double)t*dt,x_b);				// (kJ/kg/C) Heat Capacity at node "b"
		density_b=(double)DENSITY(T_b,(double)t*dt,x_b);	// (g/cm3) Density at node "b"

		deltaT = q_net/cp_b/density_b*dt/dx;	// (C) Temperature increase for this time step
		Temp_now[NODES-1] = T_b + deltaT;			// (C) Temperature for next time step
		//thermalDiffusivity = K(Temp[t][NODES-1],(float)t*dt,((float)NODES-0.5)*dx)/CP(Temp[t][NODES-1],(float)t*dt,((float)NODES-0.5)*dx)/DENSITY(Temp[t][NODES-1],(float)t*dt,((float)NODES-0.5)*dx) * (double)THERMAL_DIFFUSIVITY_FACTOR;
		//Temp[t+1][NODES-1] = Temp[t][NODES-1] + thermalDiffusivity * dtdx2 * (double)FOURIER_FACTOR * ( TBC_right[t] - 2 * Temp[t][NODES-1] + Temp[t][NODES-2]);

		// Calculate interior nodes
		for(i=1;i<NODES-1;i++)  // Loop over internal nodes only (don't write over BC nodes)
		{
			// Note: "b" is the current node, "a" is one to the left, "c" is one to the right

			T_a=Temp_past[i-1];	// (C) Temperature at node "a"
			T_b=Temp_past[i];		// (C) Temperature at node "b"
			T_c=Temp_past[i+1];	// (C) Temperature at node "c"

			x_a=((double)(i-1)+0.5)*dx;		// (cm) Location of node "a"
			x_b=((double)(i)+0.5)*dx;		// (cm) Location of node "b"
			x_c=((double)(i+1)+0.5)*dx;		// (cm) Location of node "c"

			k_a=(double)K(T_a,(double)t*dt,x_a);	// One node to the left (i-1)
			k_b=(double)K(T_b,(double)t*dt,x_b);	// Current node
			k_c=(double)K(T_c,(double)t*dt,x_c);	// One node to the right (i+1)

			k_ba=2*k_b*k_a/(k_b+k_a);		// (W/m/C) Conduction in series, 2/k_ba=1/k_a+1/k_b
			k_bc=2*k_b*k_c/(k_b+k_c);		// (W/m/C) Conduction in series, 2/k_bc=1/k_b+1/k_c

			q_ba=k_ba*(T_a-T_b)/dx*(double)M_OVER_CM;	// (W/cm2) Heat flux from "a" to "b." Positive means heat flows into "b"
			q_bc=k_bc*(T_c-T_b)/dx*(double)M_OVER_CM;	// (W/cm2) Heat flux from "b" to "c." Positive means heat flows into "b"
			q_net=q_ba+q_bc;

			cp_b=(double)CP(T_b,(double)t*dt,x_b);				// (kJ/kg/C) Heat Capacity at node "b"
			density_b=(double)DENSITY(T_b,(double)t*dt,x_b);	// (g/cm3) Density at node "b"

			deltaT = q_net/cp_b/density_b*dt/dx;	// (C) Temperature increase for this time step
			Temp_now[i] = T_b + deltaT;			// (C) Temperature for next time step
			// Code faster if don't have so many intermediate variables (but this is much easier to understand)
			//thermalDiffusivity = K(Temp[t][i],(float)t * dt,((float)i+0.5)*dx)/CP(Temp[t][i],(float)t*dt,((float)i+0.5)*dx)/DENSITY(Temp[t][i],(float)t*dt,((float)i+0.5)*dx) * (double)THERMAL_DIFFUSIVITY_FACTOR;
			//Temp[t+1][i] = Temp[t][i] + thermalDiffusivity * dtdx2 * (double)FOURIER_FACTOR	* ( Temp[t][i+1] - 2 * Temp[t][i] + Temp[t][i-1]);
			// Note to self: For more speed, consider using pointers (Temp[a][b] has to perform multiplication to find address)
		}

		// See if we have a very high or low temperature (save in Tmax or Tmin)		
		for(i=0;i<NODES;i++)  // Loop over all nodes
		{
			if(Temp_now[i]<Tmin) Tmin=Temp_now[i]; // New minimum temperature found, save it
			if(Temp_now[i]>Tmax) Tmax=Temp_now[i]; // New maximum temperature found, save it
		}

//MORE SAVE TO FILE STUFF/////
		// Print the results for this timestep
		if(t % PRINTSTEP==0)  // No remainder, evenly divisible, then save to file
		{
			// Calculate average nodal temperature at this timestep
			for(i=0,aveTemp[t]=0;i<NODES;i++)  // Loop over INTERIOR nodes (not BC nodes)
			{
				aveTemp[t] += Temp_past[i];  // Sum temperatures of the nodes
			}
			aveTemp[t] /= NODES;

			// Calculate average enthalpy at this timestep
			// WRONG!!! FIX THIS!!!!  Need to integrate across nodes!!!! (cp not constant)
			aveH[t]= (double)CP(Tmax_predict,(float)TIME*dt,0.0) * (aveTemp[t] - (double)TREF) * ENTHALPY_FACTOR;

			// Now save it in the file
			fprintf(dataFile, "%6d\t%8.4f\t%8.3f\t%6.2f\t%6.2f\t",t,(double)t * dt, aveH[t], aveTemp[t], TBC_left[t]);
			for(i=PRINTNODESTART-1;i<NODES;i+=PRINTNODE) fprintf(dataFile, "%6.2f\t",Temp_past[i]);
			fprintf(dataFile, "%6.2f\n", TBC_right[t]);

			printf("\b\b\b\b%3d%c",t*100/STEPS,'%');
		}
//DONE SAVING TO FILE STUFF/////

		// Swap Storage Space for next timestep
		if((Temp_now == & storage_a[0]) && (Temp_past == & storage_b[0]))
		{
			Temp_now = & storage_b[0];
			Temp_past = & storage_a[0];
		}
		else if((Temp_now == & storage_b[0]) && (Temp_past == & storage_a[0]))
		{
			Temp_now = & storage_a[0];
			Temp_past = & storage_b[0];
		}
		else
		{
			printf("\nMEMORY ERROR: Wandering pointers...\n");
			return(FAILURE);
		}
	}

	fprintf(dataFile, "\n\nMin Temp (C):\t%.2f\nMax Temp (C):\t%.2f\n",Tmin,Tmax);
	// Close the file
	fclose(dataFile);

	// PC clock time at end of this function
	time(&endTime);

	// Let the user know we're done, and how long it took...
	printf("\b\b\b\b[complete in %d sec]\n", endTime-startTime);

	return(SUCCESS);
};


// Appends the original c source file (if it can find it) to output file
int append_source_file(void)
{
	int letter;						// Store the character (must be int to find EOF) 
	FILE *dataFile, *sourceFile;	// File input/output

	//printf("Data being saved in %s...\t",FILE_NAME);

	// Try to open the source file (contains the original c code)
	// Note: __FILE__ is a preprocessor define that is equal to the source code file name
	if( (sourceFile = fopen(__FILE__,"r"))==FAILURE )
	{
		printf("\nnote: %s missing, couldn't append to %s.\n\n",__FILE__, FILE_NAME);
		return(FAILURE);
	}

	// Open the file (check for failure)
	if( (dataFile = fopen(FILE_NAME,"a"))==FAILURE ) return(FAILURE);

	// Write out some info to the file
	fprintf(dataFile,"\n\n// Source code from %s ",__FILE__);
	fprintf(dataFile,"// Originally compiled on %s %s.\n// Recompile and execute to generate the above data\n\n", __DATE__, __TIME__); 

	// Read and write each character from "sourceFile" to "dataFile"
	while ( (letter=fgetc(sourceFile)) != EOF )
	{
		if(letter==ASCII_TAB) putc(ASCII_SPACE, dataFile);  // Replace a tab with 4 spaces (more Excel friendly)
		else putc(letter, dataFile);  // If it's not a tab, write the letter
	}

	fclose(dataFile);	// Close the data file
	fclose(sourceFile);	// Close the source file
	return(SUCCESS);

}


// Make the user press enter, prevents the window from closing too fast...
void press_enter(void)
{
	char presskey;

	printf("Press <enter> to continue...");
	do{
		presskey=getchar();
		//putchar(presskey);
	} while (presskey!='\n');

	//scanf("%c",&presskey);  //Another way of doing the same thing

}

// The "main" program (start here)
int main (void)
{
	printf("1D Heat Transfer\nby Wayne Pafko\n\n");

	// Run appropriate subroutines (if they fail, don't continue)
	if( initial_conditions() !=SUCCESS) {press_enter(); return(FAILURE);}
	if ( solve_and_save_problem() !=SUCCESS) {press_enter(); return(FAILURE);}
	//if( save_results() !=SUCCESS) {press_enter(); return(FAILURE);}
	append_source_file(); // No big deal if it fails (couldn't find the source file)

	// All done, ready to leave
	press_enter();
	return(SUCCESS);
}

// Created and compiled with Microsoft Visual C++ 4.0