#include #include #include #include "structs.h" #include "mymath.h" #include "nrutil.h" #include "const.h" #include "debug.h" #define OVERSAMPR 10. #define MINR 1.e-15 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 isophote (float a[], int ia[], struct fitpars *high, struct image *model, struct derivs *df) { extern struct inpars input; unsigned long i, ix, iy, totsamp; int j, nsubsamp, os, n, mode, pa; float am, ta[NPARS+1], OcosPA, OsinPA, sumflux, xf, yf, xp, yp, xc, yc, flux, dfdr, c, step, totcnts, dc=0.1, dratio_dc, err; float *drd, *fdrd, *rdrd, *dIedfp, drdc; struct fitpars *f=NULL, *c0=NULL, *r=NULL, *w=NULL; struct derivs *dptr, *fdf=NULL, *c0df=NULL, *rdf=NULL, *wdf=NULL; os = input.sampfac; for (i=1; i <= NPARS; i++) ta[i] = a[i]; /* a[6] is the gradient at each isophote. */ ta[6] = -a[6]; /**********************************************************\ * 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); /*****************************************************************\ * * * Pre-calculate some things we will frequently use later on but * * which needs to be calculated one time only. * * * \*****************************************************************/ if (c0 != NULL) c = c0->a[1] + 2; else c = 2.; ta[10] = -1. * 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++) { sumflux = 0.; drdblah ((float)ix, (float)iy, c, 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 (am >= a[4] - a[5]/2.- 1 && am <= a[4] + a[5]/2.+ 1) { if (am < OVERSAMPR) { if (ta[4] * ta[9] <= 1 && am <= os) { /* If scalelength is small, oversample by a lot more. */ nsubsamp = IMIN(100, (int)(1./fabs(ta[4] * ta[9]) * 20 / os)); } else if (am <= 3. * os) nsubsamp = 20 / os; else { /* If scalelength is large, bigger step */ nsubsamp = IMIN(100, NINT(20./am/os)); }; } else nsubsamp = 1; nsubsamp = IMAX (1, nsubsamp); totsamp = nsubsamp * nsubsamp; step = 1./nsubsamp; xf = ix - 0.5 + 1./(2 * nsubsamp); yf = iy - 0.5 + 1./(2 * nsubsamp); /****************************\ * Now loop for subsampling * \****************************/ 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, c, ta, ia, OcosPA, OsinPA, c0, f, r, &am, drd, &drdc, fdrd, rdrd); /************************\ * Classical parameters * \************************/ /* this needs to be in the loop because it *\ \* modified near the edges. */ dfdr = ta[6] / totsamp; /* Isophote edges are discontinuties... *\ \* treat carefully */ if (am > a[4] + a[5]/2. - 0.5 ) dfdr = -0.01 * (ta[6] * (am-a[4]) + a[3])/ fabs(ta[6] * (am-a[4]) + a[3]) / totsamp; else if (am < a[4] - a[5]/2. + 0.5) dfdr = 0.01 * (ta[6] * (am-a[4]) + a[3])/ fabs(ta[6] * (am-a[4]) + a[3]) / totsamp; if (ia[1] > 0) df->dpm[1][iy][ix] += dfdr * drd[1]; if (ia[2] > 0) df->dpm[2][iy][ix] += dfdr * drd[2]; if (ia[3] > 0) df->dpm[3][iy][ix] += 1. / totsamp; if (ia[6] > 0) df->dpm[6][iy][ix] += -(am-a[4]) / totsamp; if (ia[9] > 0) df->dpm[9][iy][ix] += dfdr * drd[9]; if (ia[10] > 0) df->dpm[10][iy][ix] += dfdr * drd[10]; /********************************************\ * Traditional diskyness/boxyness parameter * \********************************************/ if (c0 != NULL && c0->ia[1] > 0) c0df->dpm[1][iy][ix] += (dfdr * drdc); /***********************************\ * Fourier and PA rotational terms * \***********************************/ if (f != NULL) high_derivs (ix, iy, f, fdf, dfdr, fdrd, 0., fdrd); /* the second fdrd is dummy */ if (r != NULL) high_derivs (ix, iy, r, rdf, dfdr, rdrd, 0., rdrd); /* the second rdrd is dummy */ sumflux += (ta[6] * (am-a[4]) + a[3]); yc = yc + 1.; }; xc = xc+1.; }; /* Sum within annulus only */ if (am <= a[4] + a[5]/2. || am >= a[4] - a[5]/2.) { sumflux = sumflux / totsamp; df->dpm[0][iy][ix] = sumflux; /* ... for convolution */ }; }; }; }; free_vector (drd, 1, NPARS); if (f != NULL) free_vector (fdrd, 1, f->npars); if (r != NULL) free_vector (rdrd, 1, r->npars); /* #if DEBG sprintf (model->name, "model-0.fits"); Output image writefits ("test.fits", model, 0); #endif */ }