/***************************************************************************\
*  This algorithm calculates the integrated flux for a (sub)pixel element.  *
\***************************************************************************/

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

float fqromo(float (*func)(float r, int null1, float *a, float
	*fnull, int *inull, int null2, float cnull), float a, float b, int
	null1, float *ta, float *fnull, int *inull, int null2, float kappa);

float fluxinteg (float rmax, float *a, float (*func)(float r,
    int null1, float *a, float *fnull, int *inull, int null2, float kappa),
    float kappa)
{
    float sersicfunc (float r, int null1, float *a, float *fnull,
	int *inull, int null2, float kappa);

    int *inull;
    float *fnull, flux;

    flux = fqromo (sersicfunc, 0., rmax, -999, a, fnull, inull, 0, kappa);
    return (2 * PI * flux * a[9]);
}

/***************************************************************************\
*  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 as Fourier modes.     *
\***************************************************************************/

float sersicfunc (float r, int null1, float *a, float *fnull, int
    *inull, int null2, float kappa)
{
    float integrand;

    integrand = r * exp (-kappa * ( pow(r/a[4], 1./a[5]) - 1 ));

    return (integrand);   
}