Найти площадь поверхности сегмента эллипсоида

#include "stdafx.h"
#define M_PI       3.14159265358979323846
#include<iostream>
#include<math.h>
#include<iomanip>
using namespace std;                         
 
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_ellint.h>
 
float a = 2, b = 1.2, c = 0.7; // ellipsoid radiuses
 
// the function
 
double sq(double x, void* data) 
{
    double d = 1 - c*c / (a*a);
    double eps = 1 - c*c / (b*b);
    float x2 = x*x;
    float dx = 1 - d*x2;
    float m = eps*(1 - x2) / dx;
    return sqrtf(dx)*gsl_sf_ellint_Ecomp(m, GSL_PREC_DOUBLE);
}
 
float Romberg(float a, float b) // calculates integral using Romberg method
 
{
    float R[6][6];
    int i, j, k;
    int n = 5;
 
    float h, sum;
 
    h = b - a;
 
    R[0][0] = (h / 2)*(sq(a, NULL) + sq(b, NULL));
 
    for (i = 1; i <= n; i++)
 
    {
 
        h = h / 2;
 
        sum = 0;
 
        for (k = 1; k <= int(pow(2.0, i) - 1); k += 2)
 
            sum += sq(a + k*h, NULL);
 
        R[i][0] = R[i - 1][0] / 2 + sum*h;
 
        for (j = 1; j <= i; j++)
 
            R[i][j] = R[i][j - 1] + (R[i][j - 1] - R[i - 1][j - 1]) / (pow(4.0, j) - 1);
 
    }
    return (R[n][n]);
}
 
 
int _tmain(int argc, _TCHAR* argv[])
{
    float ac = a*a - c*c;
    float mu = asinf(sqrtf(ac) / a);
    float k = a / b*sqrtf(b*b - c*c) / sqrtf(ac);
    float square = M_PI*(c*c + b / (sqrtf(ac))*
        (c*c*gsl_sf_ellint_F(mu, k, GSL_PREC_DOUBLE) + ac*gsl_sf_ellint_E(mu, k, GSL_PREC_DOUBLE)));
 
    float squareRomberg = 4 * a*b*Romberg(0, 1);
 
    gsl_integration_workspace * w
        = gsl_integration_workspace_alloc(1000);
 
    double result, error;
 
    gsl_function F;
    F.function = &sq;
    F.params = NULL;
 
    gsl_integration_qags(&F, 0, 1, 0, 1e-7, 1000,
        w, &result, &error);
    float squareGauss = result * 4 * a*b;
 
    cout << "Right square = " << square << endl;
    cout << "Square using Romberg = " << squareRomberg << endl;
    cout << "Square using Gauss (GSL) = " << squareGauss << endl;
 
    system("PAUSE");
 
    return 0;
 
}