/***************************************************************************\
*  This algorithm calculates the total area for a generalized azimuthal     *
*  shape that has a unit radius.  Why do we do this?  For a perfect         *
*  ellipsoid with unit radius=1 and axis ratio q, the area is just:         *
*  area = PI * q.  We also know how to calculate the total flux of an       *
*  ellipsoid with a Sersic, expdisk, moffat, and gaussian radial profiles   *
*  which have simple analytic solutions.  Then, the amazing fact            *
*  is that the ratio of the volume flux between the generalized shape       *
*  and the ellipsoid is just the ratio of their *azimuthal areas*, if all   *
*  the *radial* profile parameters are the same.  This means that we don't  *
*  have to numerically integrate out to inifinity to get the volume flux,   *
*  which can take a long time.  We just have to compare the areas of the    *
*  azimuthal shape with the ellipsoid, and scale the flux accordingly.      *
\***************************************************************************/

#include <math.h>
#include "const.h"

float fqromo(float (*func)(float theta, int mode, float *ca, float *fa, int
	*ifa, int nf, float c), float a, float b, int mode, float *ca, float
	*fa, int *ifa, int nf, float c);

float fourarea(float *ca, float *fa, int *ifa, int nf, float c)
{
    float farea (float theta, int mode, float *a, float *fa, int *ifa, 
	int nf, float c);

    float sum = 0.;

    sum += fqromo (farea, 0., 0.5 * PI, -999, ca, fa, ifa, nf, c);
    sum += fqromo (farea, 0.5 * PI, PI, -999, ca, fa, ifa, nf, c);
    sum += fqromo (farea, -0.5 * PI, 0., -999, ca, fa, ifa, nf, c);
    sum += fqromo (farea, -PI, -0.5 * PI, -999, ca, fa, ifa, nf, c);
    return (sum);
}

/***************************************************************************\
*  Calculate the differential area element of the azimuthal profile.        *
*  "mode" is a dummy variable here so that the derivative integrals can     *
*  use the same qromb algorithm to do the integration.                      *
\***************************************************************************/

float farea (float theta, int mode, float *a, float *fa, int *ifa, int nf, 
    float c)
{
    int i, m;
    float foursum = 1., integrand, fpa, coeff, x, y, sgntheta=1;

    for (i=1; i <= nf; i+=2) {
        if (!(fa[i] ==0 && ifa[i] ==0) && ifa[i] != -1) {
	    m = i/2 + 1;
            fpa = -fa[i+1] / 180. * PI;

	    foursum += fa[i] * cos ( m * (atan2(sin(theta)/a[9], cos(theta)) +
		       fpa));
	};
    };

    coeff = pow( pow( cos(theta)/sin(theta), c) + pow (1./a[9], c), 1./c);

    if (theta < 0.) sgntheta = -1;

    y = 1./coeff / fabs(foursum) * sgntheta;
    x = y/tan (theta);

    integrand = 0.5 * (x*x + y*y);

    return (integrand);   
}