Вычисление площади криволинейной трапеции

/*вычислить с точностью Е площадь криволинейной трапеции ограниченной осью , Ох, прямыми х=а х=B , и кривой 
y=f(x)=(1-0,25sin^2tx)^1\2 . параметр "t"- корень уравнения t^3-0,39t^2-10,5t+11=0,
принад. отрезку [c, d] и отрезку с точ. Е
 
исх данные:
a=0 b=П\2
с=1
d=2
E=0,001
*/
#include <iostream>
#include <cmath>
#define _PI 3.1415926535897932384626433832795
using std::cin;
using std::cout;
using std::endl;
 
 
//////////////////////////////////////////////////////////////////////////////////////
//      уравнение трапеции(интегрируем)                                             //
//      вход:                                                                       //
//          X-координата  по оси абцисс                                                //
//          T-параметр является корнем уравнения f0                                  //
//      выход:                                                                     //
//          значение уравнения                                                     //
//////////////////////////////////////////////////////////////////////////////////////
float f1(float x, float t);
 
//////////////////////////////////////////////////////////////////////////////////////
//      уравнение (ищем корень)                                                     //
//      вход:                                                                       //
//          T-переменная                                                          //
//      выход:                                                                     //
//          значение уравнения в точке t                                         //
//////////////////////////////////////////////////////////////////////////////////////
float f0(float t);
 
//////////////////////////////////////////////////////////////////////////////////////
//      ПРЯМОУГОЛЬНИКИ (правые) - интегрирование                                  //
//      вход:                                                                       //
//          a-начало интервала интегрирования                                      //
//          б-конец интервала интегрирования                                       //
//          h-точность интегрирования(совпадает с шагом интегрирования               //
//          t-параметр из решения уравнения                                           //
//      выход:                                                                     //
//          интеграл функции f1  сточностью h на интервале [a,b]                    //
//////////////////////////////////////////////////////////////////////////////////////
float rectangle(float a,float b,float h,float t);
 
//////////////////////////////////////////////////////////////////////////////////////
//      ПОЛОВИННОЕ ДЕЛЕНИЕ - поиск корня                                            //
//      вход:                                                                       //
//          a-начало интервала поиска                                              //
//          б-конец интервала поиска                                               //
//          h-точность поиска                                                     //
//      выход:                                                                     //
//          интеграл функции f0  сточностью h на интервале [a,b]                    //
//          если корня нет то значение не попадающее в интервал поиска             //
//////////////////////////////////////////////////////////////////////////////////////
float separ(float a,float b,float e);
 
//////////////////////////////////////////////////////////////////////////////////////
//      главная функция                                                               //
//      вход:                                                                       //
//      выход:                                                                     //
//          0-нет ошибок                                                           //
//          1-ошибка нет корня уравнения fo                                          //
//////////////////////////////////////////////////////////////////////////////////////
int main(){
    float a,b,c,d,e,tmp;
    //a=0,b=_PI/2.,c=1,d=2,e=0.001;
    //cout<<rectangle(a,b,e,separ(c,d,e));
    cout<<"a: "; cin>>a;
    cout<<"b: "; cin>>b;
    cout<<"c: "; cin>>c;
    cout<<"d: "; cin>>d;
    cout<<"e: "; cin>>e;
    if((tmp=separ(c,d,e))<c)
        return 1;
    cout<<rectangle(a,b,e,tmp)<<endl;
    system("pause");
    return 0;
}
 
//////////////////////////////////////////////////////////////////////////////////////
//      уравнение трапеции(интегрируем)                                             //
//      вход:                                                                       //
//          X-координата  по оси абцисс                                                //
//          T-параметр является корнем уравнения f0                                  //
//      выход:                                                                     //
//          значение уравнения                                                     //
//////////////////////////////////////////////////////////////////////////////////////
float f1(float x, float t){ 
    return sqrt((1-0.25*sin(t*x)*sin(t*x)));
}
 
//////////////////////////////////////////////////////////////////////////////////////
//      уравнение (ищем корень)                                                     //
//      вход:                                                                       //
//          T-переменная                                                          //
//      выход:                                                                     //
//          значение уравнения в точке t                                         //
//////////////////////////////////////////////////////////////////////////////////////
float f0(float t){          
    return t*t*t-0.39*t*t-10.5*t+11;
}
 
//////////////////////////////////////////////////////////////////////////////////////
//      ПРЯМОУГОЛЬНИКИ (правые) - интегрирование                                  //
//      вход:                                                                       //
//          a-начало интервала интегрирования                                      //
//          б-конец интервала интегрирования                                       //
//          h-точность интегрирования(совпадает с шагом интегрирования               //
//          t-параметр из решения уравнения                                           //
//      выход:                                                                     //
//          интеграл функции f1  сточностью h на интервале [a,b]                    //
//////////////////////////////////////////////////////////////////////////////////////
float rectangle(float a,float b,float h,float t){
    float s,x1,x2;
    for(s=0,x1=a,x2=b;x1<x2;x1+=h)  //площадь=0,x1=начало интервала,х2=конец,пока х1 < х2
        s+=f1(x1,t);                //увелививаем х1 на величину шага и к площади прибавляем значение
                                    //функции в точке x1
    return s*h;                     //чтобыфиз суммызначений функции получить площадь умножим на величину шага
}
//////////////////////////////////////////////////////////////////////////////////////
//      ПОЛОВИННОЕ ДЕЛЕНИЕ - поиск корня                                            //
//      вход:                                                                       //
//          a-начало интервала поиска                                              //
//          б-конец интервала поиска                                               //
//          h-точность поиска                                                     //
//      выход:                                                                     //
//          интеграл функции f0  сточностью h на интервале [a,b]                    //
//          если корня нет то значение не попадающее в интервал поиска             //
//////////////////////////////////////////////////////////////////////////////////////
float separ(float a,float b,float e){
    float x1=a;
    float x2=b;
    float x;
    if(f0(a)*f0(b)>0)       //функция не пересекает ось X
        return a-1;         //если корня нет то возвращаем значение не попадающее в интервал поиска
    do{
        x=(x1+x2)/2;        //находим середину интервала поиска
        if(f0(x1)*f0(x)>0)  //если значения начала интервала и середины имеют разные знаки
            x1=x;           //теперь середина интервала будет являться началом интервала поиска
        else                //иначе корень на второй половине
            x2=x;           //теперь середина интервала будет являться концом интервала поиска
    }while(fabs(f0(x))>e);
    return x;
}