/*****************************************************
*  rootfinderFP.cc: Find the root of a function using
*                   the false position method.
*
*  written by: Paul M Eugenio
*              Department of Physics
*              Florida State University
*
*  modified by: Volker Crede 
*               Shawn Havery
*
*  Computational Physics Project #7
*  Date: 22 Feb 2008
*
******************************************************/                     
#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;

// Global Constants
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);
  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;
}



////////////////////
//  Program Main  //
////////////////////
int main(){
  // variables
  int count = 0;                    // iteration counter
  double thetaMin = 0.17;           // about 1 deg
  double thetaMax = 1.55;           // about 89 deg
  const double delta = 0.17e-3;     // our desired accuracy
  double thetaZero;                 // theta value between Min & Max
  double fZero;                     // the value of f(thetaZero)
  double slope;                     // the slope between Min & Max
  double intercept;                 // the intercept of the above line

  Function func;

  do{
    count++;
    slope = ( func.f(thetaMax) - func.f(thetaMin) ) / ( thetaMax - thetaMin );
                                                     // find the slope
    intercept = func.f(thetaMin) - slope * thetaMin;
                                                     // find the intercept
    thetaZero = -intercept / slope;                  // find the "midpoint"
    fZero = func.f( thetaZero );                     // get f(midpoint)

    cout << setw(2) << count << setw(16) << fZero << endl;  
                                                     // output for plotting
    if( (func.f(thetaMin) * fZero) > 0 )             // compare relative sign
      thetaMin = thetaZero;
    else
      thetaMax = thetaZero;
  }while( fabs(fZero) > delta );                // loop until f(x)=0 to within 
                                                //      the desired accuracy

  func.printValue( thetaZero );
}