#include <math.h>
#include <stdlib.h>
#include "structs.h"
#include "nrutil.h"
#include "mymath.h"
#include "const.h"
#include "debug.h"

#define OVERSAMPR 5.
#define peaksamp 50
#define dr 0.05 / peaksamp
#define MINR 1.e-15

void writefits(char *phead, struct image *img, char *object, int add);
void highpars (struct fitpars *high, struct derivs *df, struct fitpars **four,
    struct derivs **fourdf, struct fitpars **c0, struct derivs **c0df,
    struct fitpars **r, struct derivs **rdf);
void drdblah (float x, float y, float c, float *a, int *ia, float OcosPA, 
    float OsinPA, struct fitpars *c0, struct fitpars *f, struct fitpars *r, 
    float *am, float *drd, float *drdc, float *fdrd, float *rdrd);
void high_derivs (int ix, int iy, struct fitpars *h, struct derivs *hdf,
    float dfdr, float *drd, float dfdIe, float *dIedfp);

void fsinc (float a[], int ia[], struct fitpars *high, struct image *model, 
    struct derivs *df)
{
    extern struct inpars input;

    unsigned long i, ix, iy, totsamp;
    float nukefunc (float ta[], float x, float y);

    float am, af, ta[NPARS+1], OcosPA, OsinPA, flux,
	rta4, drdc, dIdr, c0v;
          
    int n, j, os, mode, nsubsamp;
    float *drd, *fdrd, *rdrd, xf, yf, xc, yc, xp, yp, dfdI, sumdfdI, step;
    struct fitpars *f=NULL, *c0=NULL, *r=NULL, *w=NULL;
    struct derivs *dptr, *fdf=NULL, *c0df=NULL, *rdf=NULL, *wdf=NULL;

/*  Ib=ta[3], rb=ta[4], alpha=ta[5], beta=ta[6], gamma=ta[7]  */

    os = input.sampfac;

    for (i=1; i <= NPARS; i++)
        ta[i] = a[i];

    /**********************************************************\
    *  Set pointers to high order correction structures, then  *
    *  create arrays to hold dr/dpar values.                   *
    \**********************************************************/

    highpars (high, df, &f, &fdf, &c0, &c0df, &r, &rdf);
    drd = vector (1, NPARS);
    if (f != NULL) fdrd = vector (1, f->npars);
    if (r != NULL) rdrd = vector (1, r->npars);

    if (c0 != NULL)
        c0v = c0->a[1] + 2;
    else
        c0v = 2.;

    /******************************\
    *  For the azimuthal profiles  *
    \******************************/

    ta[10] = -(ta[10])/ 180. * PI + PI/2.;  /* "Intuitive" orientation */
    OcosPA = cos (ta[10]);
    OsinPA = sin (ta[10]);

    /*******************************************\
    *  Now fill the image and derivative grids  *
    \*******************************************/

    for (iy=1; iy <= df->naxes[2]; iy++) {
        for (ix=1; ix <= df->naxes[1]; ix++) {

            drdblah ((float)ix, (float)iy, c0v, ta, ia, OcosPA,
                OsinPA, c0, f, r, &am, drd, &drdc, fdrd, rdrd);

            /******************************************************\
            *  If we're near the center of the galaxy we'll        *
            *  subdivide the sampling.                             *
            \******************************************************/

            if (fabs(am) < OVERSAMPR) {
                if (ta[4] * ta[9] <= 1. && fabs(am) <= os) {

                    /* If scalelength is small, oversample by a lot more.  */

                    nsubsamp = IMIN(100, (int)(1./fabs(ta[4]*ta[9])*2*
                               OVERSAMPR/os));
                } else  if (fabs(am) <= 4. * os)
                    nsubsamp = 2 * OVERSAMPR / os;
                else {
                    /* If scalelength is large, bigger step */

                    nsubsamp = IMIN(100, NINT(2*OVERSAMPR/fabs(am)/os));
                };
            } else
                nsubsamp = 1;

            nsubsamp = IMAX (1, nsubsamp);
            totsamp = nsubsamp * nsubsamp;

            xf = ix - 0.5 + 1./(2 * nsubsamp);
            yf = iy - 0.5 + 1./(2 * nsubsamp);
            step = 1./nsubsamp;

	    sumdfdI = 0.;

	    /*************************\
	    *  Now loop for sampling  *
	    \*************************/

            xc = 0.;
            while (xc <= nsubsamp - 1) {
                yc = 0.;
                xp = xf + xc * step;
                while (yc <= nsubsamp - 1) {
                    yp = yf + yc * step;

                    if (nsubsamp > 1)
                        drdblah (xp, yp, c0v, ta, ia, OcosPA, OsinPA, c0,
                            f, r, &am, drd, &drdc, fdrd, rdrd);

            	    /************************\
            	    *  Classical parameters  *
            	    \************************/

            	    rta4 = am/ta[4];

		    if (am == 0.) {
	                flux= ta[3];
		        dfdI = 1. / totsamp;
			dIdr = 0.;
		    } else {
		        dfdI = sin(rta4)/rta4 / totsamp;
	                flux= ta[3] * dfdI;
                        dIdr = ta[3] * (cos(rta4)/am - ta[4]/am/am*sin(rta4)) /
				totsamp;
		    };

            	    if (ia[1] > 0) df->dpm[1][iy][ix] += dIdr * drd[1];
            	    if (ia[2] > 0) df->dpm[2][iy][ix] += dIdr * drd[2];
		    if (ia[3] > 0) {
			if (am == 0.)
			    df->dpm[3][iy][ix] += 1/totsamp;
			else
			    df->dpm[3][iy][ix] += dfdI;
		    };

	    	    if (ia[4] > 0) {
			if (am == 0.)
			    df->dpm[4][iy][ix] += 0.;
			else
			    df->dpm[4][iy][ix] += 
	       		        ta[3] * (-cos(rta4)/ta[4] + sin(rta4)/am)/totsamp;
		    };

            	    if (ia[9] > 0) df->dpm[9][iy][ix] += dIdr * drd[9];
		    if (ia[10] > 0) df->dpm[10][iy][ix] += dIdr * drd[10];

            	    /********************************************\
            	    *  Traditional diskyness/boxyness parameter  *
            	    \********************************************/

            	    if (c0 != NULL && c0->ia[1] > 0)
                	c0df->dpm[1][iy][ix] += (dIdr * drdc);

            	    /*****************\
            	    *  Fourier terms  *
            	    \*****************/

            	    if (f != NULL) high_derivs (ix, iy, f, fdf, dIdr, fdrd,
                	0., fdrd);  	/* the second fdrd is dummy */
            	    if (r != NULL) high_derivs (ix, iy, r, rdf, dIdr, rdrd,
                	0., rdrd);  	/* the second rdrd is dummy */

		    sumdfdI += dfdI;
		    yc = yc + 1;
        	};
		xc = xc + 1;
    	    };

            flux = ta[3] * sumdfdI;
            df->dpm[0][iy][ix] = flux;               /* ... for convlution   */
	};
    };

    free_vector (drd, 1, NPARS);
    if (f != NULL) free_vector (fdrd, 1, f->npars);
    if (r != NULL) free_vector (rdrd, 1, r->npars);

}

#if DEBUG
    sprintf (model->name, "model-0.fits");  /* Output image */
    writefits ("test.fits", model, "", 0);
#endif