#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 float ratio (float c); float dfridr (float (*ratio)(float n), float x, float h, float *err); 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 *dI0dfp); float fourarea(float *ca, float *fa, int *ifa, int nf, float c); 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 phderiv (float theta, int mode, float *a, float *fa, int *ifa, int nf, float c); float ampderiv (float theta, int mode, float *a, float *fa, int *ifa, int nf, float c); float qderiv (float theta, int mode, float *a, float *fa, int *ifa, int nf, float c); float c0deriv (float theta, int mode, float *a, float *fa, int *ifa, int nf, float c); void exponential (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, mode, n; float xp, yp, am, ta[NPARS+1], OcosPA, OsinPA, sumdfdI0, xf, yf, xc, yc, flux, gexp, dfdI0, dfdr, step, dc=0.1, err, fpa, foursum, c, drdc, dfdIe; float dI0dmag=0., area_ratio = 1, dratio_dc = 0., totcnt, uarea; float dI0dc, dI0drs, dI0dq; float *drd, *fdrd, *rdrd, *dI0dfp; struct fitpars *f=NULL, *c0=NULL, *r=NULL; struct derivs *dptr, *fdf=NULL, *c0df=NULL, *rdf = NULL; os = input.sampfac; for (i=1; i <= NPARS; i++) ta[i] = a[i]; ta[10] = -ta[10]/ 180. * PI + PI/2.; /* "Intuitive" orientation */ /**********************************************************\ * 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 (r != NULL) rdrd = vector (1, r->npars); if (c0 != NULL) c = c0->a[1] + 2; else c = 2.; if (f != NULL) { fdrd = vector (1, f->npars); uarea = fourarea (ta, f->a, f->ia, f->npars, c); area_ratio = a[9] * PI / uarea; }; /*****************************************************************\ * * * Pre-calculate some things we will frequently use later on but * * which needs to be calculated one time only. * * * \*****************************************************************/ /* a[3] is the galaxy total brightness while ta[3] */ /* is the count rate at the disk scalelength. */ /* Notice the last term ta[9] in the first equation here. */ /* That's because the distance from the center is defined */ /* as r = (|x|^c + |y/q|^c)^(1/c), where q is the axis ratio. */ /* This would make the magnitude come out right. */ /* The area_ratio is the ratio of a pure ellipse area */ /* to that of a disky or boxy area. Since we can calculate */ /* the flux of an elliptical model, the flux of a disky/boxy */ /* model is simply the ratio of the two areas. */ if (c0 != NULL && f == NULL) { area_ratio = ratio(c); dratio_dc = dfridr (ratio, c, dc, &err); }; totcnt = pow (10., (model->magzpt - a[3])/2.5); ta[3] = totcnt / (2. * PI * ta[4] * ta[4] * ta[9]) * area_ratio; /**********************************************************\ * Calculate some of the cross coupling terms between I0 * * with all the other parameters. * \**********************************************************/ dI0dmag = ta[3] * log(10.) / -2.5; dI0drs = -2 * ta[3] / ta[4]; dI0dq = -ta[3] / ta[9]; if (f != NULL) { dI0dfp = vector (1, f->npars); for (i=1; i <= f->npars; i+=2) { mode = i/2 + 1; if (f->ia[i] > 0) dI0dfp[i] = -ta[3] / uarea * dratio_dfp (ampderiv, mode, ta, f->a, f->ia, f->npars, c); if (f->ia[i+1] > 0) dI0dfp[i+1] = -ta[3] / uarea * dratio_dfp (phderiv, mode, ta, f->a, f->ia, f->npars, c); }; if (ia[9] > 0) dI0dq = -ta[3] / uarea * dratio_dfp (qderiv, 0, ta, f->a, f->ia, f->npars, c); if (c0 != NULL && c0->ia[1] > 0) dI0dc = ta[3] / area_ratio * dratio_dfp (c0deriv, 0, ta, f->a, f->ia, f->npars, c); } else dI0dc = ta[3] / area_ratio * dratio_dc; /******************************\ * For the azimuthal profiles * \******************************/ 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, 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 < 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; xf = ix - 0.5 + 1./(2 * nsubsamp); yf = iy - 0.5 + 1./(2 * nsubsamp); step = 1./nsubsamp; sumdfdI0 = 0.; /****************************\ * 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 * \************************/ gexp = -am / ta[4]; dfdI0 = exp(gexp) / totsamp; flux = ta[3] * dfdI0; dfdr = flux / -ta[4]; 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] += dfdI0* dI0dmag; if (ia[4] > 0) df->dpm[4][iy][ix] += (flux*am/ta[4]/ta[4] + dfdI0*dI0drs); if (ia[9] > 0) df->dpm[9][iy][ix] += (dfdr * drd[9] + dfdI0 * dI0dq); 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 + dfdI0 * dI0dc); /***********************************\ * Fourier and PA rotational terms * \***********************************/ if (f != NULL) high_derivs (ix, iy, f, fdf, dfdr, fdrd, dfdIe, dI0dfp); if (r != NULL) high_derivs (ix, iy, r, rdf, dfdr, rdrd, 0., rdrd); /* the second rdrd is dummy */ sumdfdI0 += dfdI0; yc = yc + 1.; }; xc = xc+1.; }; flux = ta[3] * sumdfdI0; df->dpm[0][iy][ix] = flux; /* ... for convlution */ }; }; free_vector (drd, 1, NPARS); if (f != NULL) { free_vector (fdrd, 1, f->npars); free_vector (dI0dfp, 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 */ }