// This program uses the Newton-Raphson method to find the root of a function
//
// General code for Function class written by Paul M. Eugenio 
// Main and Class  Modified by: Sean Kuvin
// 2-19-08

#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;


const double gkRad2Deg = 180.0 / M_PI;   // rad2pi converter
                                         // M_PI defines Pi in math.h

//  A simple Function class 
//  for functions which 
//  we want to find a root 
//
// The member function f(theta) describes 
// angle between two spatially separated 
// identical springs which are attached 
// to a hanging mass.   

class Function {
  double iMass;              // mass of object in kg
  double iHalfLength;        // half the distance between supports in m
  double iK;                 // spring constant k  in units of [N/kg]
  double iG;                 // acceleration due to gravity in units [m/s^2]
  void setG(double aG)       // private set and get for iG;
               {iG = aG;}
  double getG(){return iG;}

public:
  Function();
  Function(double aMass, double aHalfLength, double aSpringConstant);
  void setMass(double aMass){iMass = aMass;}
  void setSpringConstant(double aK){iK = aK;}
  void setHalfLength(double aHalfLength){iHalfLength = aHalfLength;}
  double getMass(){return iMass;}
  double getSpringConstant(){return iK;}
  double getHalfLength(){ return iHalfLength;}
  double f(double aTheta);
  double fprime(double aTheta);   //added for new method
  void printValue(double aTheta);
};

// Default constructor with defaults parameter values 
Function::Function() {
  setG( 9.81 );
  setMass( 5.0 );
  setHalfLength( 0.3 );
  setSpringConstant( 1000.0 );
}

// Constructor with defining values for Mass, HalfLength, SpringConstant
//
Function::Function(double aMass, double aHalfLength,double aSpringConstant) {
  setG( 9.81 );
  setMass( aMass );
  setHalfLength( aHalfLength );
  setSpringConstant( aSpringConstant );
}

// Print function 
//
void Function::printValue(double aTheta){
  cerr << endl;
  cerr << "\t f(" << aTheta << ")= " << f(aTheta) << endl;
  cerr << "\t The value of theta in degrees is " << aTheta * gkRad2Deg << endl;
}

// The function f(theta) is the function for which we want
// to find its root:
//
// A function for a mass hanging from two equal springs
//
double Function::f(double aTheta){

  double value =  tan( aTheta ) - sin(aTheta) - ( getMass() * getG() ) /
    ( 2.0 * getSpringConstant() * getHalfLength() );

  return value;
}

double Function::fprime(double aTheta){

  double value =( 1/(cos( aTheta)))*(1/(cos( aTheta ))) - cos(aTheta);

  return value;
}


////////////////////////////////////////////


int main() {
  int count = 0;                    // iteration counter
  const double delta = 0.17e-3;     // our desired accuracy
  double fX;
  double fthetaZero;
  double thetaZero;
  double interceptB;
  double slopeM;
  double theta = 0.785;	        // Begins at 45 degrees since we can take 
				// an arbitrary value between 1 and 90
  Function func;

    do{
    count++;
	fX=func.f(theta);
	slopeM=func.fprime(theta);
	interceptB= fX-(slopeM*theta);
	thetaZero=-1*(interceptB/slopeM);
	fthetaZero=func.f(thetaZero);	
	theta=thetaZero;

    cout << setw(2) << count << setw(16) << fthetaZero << endl;  

  }while( fabs(fthetaZero) > delta );            // loop until f(x)=0 to 
                                                 // the desired accuracy

  func.printValue( thetaZero );
}