/*------------------------------------------------------*/
/* NewtRaph.cpp											*/
/* Author: Aaron Pulley									*/
/* (804) 714-9052 - awpulley@cavtel.net					*/
/* Created: 12 Feb 2010									*/
/* Finds root of a function by bisection.				*/
/* History of modifications:							*/
/* -	12 Feb 2011 - Changed to Newton-Raphson while	*/
/*		re-taking class									*/
/* -													*/
/*------------------------------------------------------*/
#include <iostream>
#include <cmath>
#include <fstream>
#include <string>
using namespace std;

//function prototype
double fun(double f, double R);
double dfun(double f, double R);


int main()
{
	//user intro
	cout << "This program will use the Newton-Raphson method to determine friction (f)" << endl
		<< "from Reynolds Number (Re) within an error of 0.001%. A file of iterations" <<endl
		<< "will be generated." << endl<< endl;

	//declare variables
	double R, f, fnew, x, dx, E(.001), e;
	int ok(0);

	//creates output file
	string filename;
	cout << endl << "Please give a file name to store under." << endl << endl;
	cin >> filename;
	cout << endl;
	ofstream output(filename.c_str());

	//get f and Re
	cout << "Please enter a value for Re:  ";
	cin >> R;
	cout << "Please enter an initial guess for f:  ";
	cin >> f;
	cout << endl;

	//NR iteration
	for(int i=0; i<20; i++)
	{
		//function values
		x = fun(f, R);
		dx = dfun(f, R);

		//check deriv
		if(dx<.00001)
		{
			cout << endl << "Derivative became too close to 0. Aborted.";
			exit (0);
		}

		//NR
		fnew = f-(x/dx);

		//error
		e = 100*(fnew-f)/fnew;

		//output
		cout << i << "\t" << R << "\t" << x << "\t" << dx << "\t" << f << "\t" << e << endl;
		output << i << "\t" << R << "\t" << x << "\t" << dx << "\t" << f << "\t" << e << endl;

		//end condition
		if(e<E)
		{
			ok = 1;
			break;
		}
		
		//setup for next iteration
		f = fnew;
	}

	//conclude
	output.close();

	if(ok == 1)
	{
		cout << endl << "f = " << fnew << " +/- " << e << endl << endl;
		system("pause");
	}
	else
	{
		cout << endl << "Exceeded 20 iterations. Last value was f = " << fnew << " +/- " << e << endl << endl;
		system("pause");
	}

	return (0);
}

//functions
double fun(double f, double R)
{
	return (4.0*log10(R*sqrt(f))-.4-(1/sqrt(f)));
}

double dfun(double f, double R)
{
	return (.5*pow(f, -1.5) + (2.0/(f*log(10.0))));
}