/***************************************************************************\
*  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, gaussian, and king radial      *
*  profile, 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);

/************************************************************************\
*  Calculate amplitude derivative of the fourier area                    *
\************************************************************************/

float dratio_dfp (float (*func)(float theta, int mode, float *ca,
    float *fa, int *ifa, int nf, float c), int mode, float *ca, float *fa, int
    *ifa, int nf, float c)
{
    float sum = 0.;

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

/************************************************************************\
*  The area element for amplitude derivative of the fourier mode         *
\************************************************************************/

float ampderiv (float theta, int mode, float *a, float *fa, int *ifa, int nf,
     float c)
{
    int i, m, sgnf=1, sgnx=1, sgny = 1, sgntheta=1;
    float foursum = 1., integrand, fpa, coeff, x, y, dfdam, dxda, dyda;

    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( fabs(cos(theta)/sin(theta)), c) + pow (1./a[9], c),
            1./c);

    if (theta < 0.) {
	sgntheta = -1;
	sgny = -1;
    };
    y = 1./coeff / fabs(foursum) * sgntheta;

    m = mode * 2 - 1;
    fpa = -fa[m + 1] / 180. * PI;
    dfdam = cos (mode * (atan2(sin(theta)/a[9], cos(theta)) + fpa));

    if (foursum < 0.) sgnf = -1;
    dyda = -y/fabs(foursum) * dfdam * sgnf; /* * sgny; */

    x = y/tan (theta);

    if (x < 0.) sgnx = -1;
    dxda = dyda / tan (theta); /* * sgnx; */

    integrand = x * dxda + y * dyda;

    return (integrand); 
}


/************************************************************************\
*  The area element for amplitude derivative of the fourier mode         *
\************************************************************************/

float phderiv (float theta, int mode, float *a, float *fa, int *ifa, int nf,
     float c)
{
    int i, m, sgnf=1, sgnx=1, sgny = 1, sgntheta=1;
    float foursum = 1., integrand, fpa, coeff, x, y, dfdfpa, dxdph, dydph;

    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( fabs(cos(theta)/sin(theta)), c) + pow (1./a[9], c),
            1./c);

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

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

    m = mode * 2 - 1;
    fpa = -fa[m + 1] / 180. * PI;
    dfdfpa = mode * fa[m] * sin (mode * (atan2(sin(theta)/a[9], cos(theta)) + 
	     fpa)) / 180. * PI;

    if (foursum < 0.) sgnf = -1;
    dydph = -y/fabs(foursum) * dfdfpa * sgnf; /* * sgny; */

    x = y/tan (theta);

    if (x < 0.) sgnx = -1;
    dxdph = dydph / tan (theta); /* * sgnx; */

    integrand = x * dxdph + y * dydph;

    return (integrand);
}

/************************************************************************\
*  The area element for axis ratio derivative with fourier modes         *
\************************************************************************/

float qderiv (float theta, int mode, float *a, float *fa, int *ifa, int nf,
     float c)
{
    int i, m, sgnf=1, sgnx=1, sgny = 1, sgntheta=1;
    float foursum = 1., integrand, fpa, coeff, x, y,dxdq, dydq,
	dfdq = 0., rc;

    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));
            dfdq += m * fa[i] * sin (m * (atan2(sin(theta)/a[9], 
	            cos(theta)) + fpa)) / (1+pow(tan(theta)/a[9], 2)) * 
		    tan(theta) / a[9]/a[9];

	};
    };

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

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

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

    if (foursum < 0.) sgnf = -1;
    dydq = (y/rc/pow(a[9],c+1) - y/fabs(foursum) * dfdq * sgnf); /* * sgny; */

    x = y/tan (theta);

    if (x < 0.) sgnx = -1;
    dxdq = dydq / tan (theta); /* * sgnx; */

    integrand = x * dxdq + y * dydq;

    return (integrand); 
}

/************************************************************************\
*  The area element for c0 derivative of the fourier area                *
\************************************************************************/

float c0deriv (float theta, int mode, float *a, float *fa, int *ifa, int nf,
     float c)
{
    int i, m, sgnf=1, sgnx=1, sgny = 1, sgntheta=1;
    float foursum = 1., integrand, fpa, coeff, x, y, dxdc, dydc, rc;

    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));
	};
    };


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

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

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

    dydc = -y/fabs(foursum) * dydc; /* * sgny; */

    x = y/tan (theta);

    if (x < 0.) sgnx = -1;
    dxdc = dydc / tan (theta); /* * sgnx; */


/*
         xdp = log(FMAX(MINR, fabs(xdp))) * absxdpc;
         ydp = log(FMAX(MINR, fabs(ydp9))) *absydpc;
         drdc = am * (-1./c/c * log(fabsrc) + 1./c/ fabsrc *
             (xdp + ydp));
*/

    integrand = x * dxdc + y * dydc;

    return (integrand);   
}