#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 splint(float xa[], float ya[], float y2a[], int n, float x, float *y); float gammln(float xx); void writefits(char *phead, struct image *img, char *object, int add); void trunc_pars (char objtype[], int ia[], struct fitpars *fpar, struct derivs *df, struct trunc_links *trlinks, struct trunc_val *trval); struct fitpars *trunc_spiral (struct fitpars *r); float fluxinteg (float r, 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); void highpars (struct fitpars *high, struct derivs *df, struct fitpars **bending, struct derivs **bdf, 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 *a, int *ia, float OcosPA, float OsinPA, struct fitpars *b, struct fitpars *c0, struct fitpars *f, struct fitpars *r, float *am, float *drd, float *bdrd, 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 get_trunc_cons (char *objtype, int normtype, struct trunc_links *trlink, struct trunc_val *trval, float a[], float kap, float OcosPA, float OsinPA, float *trunc_norm, float *dfnorm); void trunc_dfdr_calc (struct trunc_val *trval, struct trunc_links *trlink, float flux0, float *PP); void trunc_derivs (int ix, int iy, float am, float flux, float dfdr, struct trunc_val *trval, struct trunc_links *tr); float ratio (float c); 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 free_fitvecs (float **drd, struct fitpars **tilt, struct fitpars *b, float **bdrd, struct fitpars *f, float **fdrd, float **dIedfp, struct fitpars *r, float **rdrd, float **dIedrp, int normtype, struct trunc_links *trlinks, struct trunc_val *trval); void sersic (char objtype[], float a[], int ia[], int normtype, struct fitpars *high, struct image *model, struct derivs *df) { float gamm2nln(float n); float getinterp(float n); extern float Y2[]; extern float KAPPA[]; extern float TWOn[]; extern int ninterp; extern struct inpars input; float sconst=3.60729; float kap=7.66925; unsigned long i, ix, iy, totsamp; int j, n, mode, nsubsamp, os; float xp, yp, am, ta[NPARS+1], xf, yf, xc, yc, rr, h, err, invp5, magfac, kapinvp5, kapfac, dkapdp5 = 0., flux, c = 2., rc, dfdr, step, dc = 0.1, dfdIe, sumdfdIe, drdc, OsinPA, OcosPA, zpt, uarea = 1., flux0; float cosi=1., sini; float area_ratio=1, dratio_dc=0., totcnt, dgammlndn; float dIedn=0., dIedq=0., dIedc=0., dIedre=0.; float *drd, *bdrd, *fdrd, *rdrd, *dIedfp, *dIedrp; float trunc_norm=1., dfnorm=0., normrad, PP; int iatmp[10], *iaptr; float atmp[10], *aptr; struct fitpars *b=NULL, *f=NULL, *c0=NULL, *r=NULL, *tilt=NULL, *tiltptr=NULL; struct derivs *dptr, *bdf=NULL, *fdf=NULL, *c0df=NULL, *rdf=NULL; struct trunc_links trlink = {NULL, 1, 1, 1, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0., NULL, NULL, NULL, NULL, NULL, NULL}, *trlptr; struct trunc_val trval = {0, 0., 0., 0., 0., 0., 0., 0., 0., 0., NULL}, *trvptr; float dfridr (float (*somefunc)(float n), float x, float h, float *err); os = input.sampfac; magfac = log(10.) / -2.5; 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, &b, &bdf, &f, &fdf, &c0, &c0df, &r, &rdf); drd = vector (1, NPARS); if (b != NULL) bdrd = vector (1, b->npars); if (f != NULL) { dIedfp = vector (1, f->npars); fdrd = vector (1, f->npars); for (i=1; i <= f->npars; i++) dIedfp[i] = 0.; }; if (r != NULL) { dIedrp = vector (1, r->npars); rdrd = vector (1, r->npars); for (i=1; i <= r->npars; i++) dIedrp[i] = 0.; }; if (c0 != NULL) c += c0->a[1]; /*----------------------------------------------------------------\ | Figure out whether we're trying to fit a truncation function. | | If so, an outer truncation is signified by trval->innout = 1, | | whereas an inner truncation is by trval->innout = 0. | | | | If the trunc parameters are modifying an inclined spiral | | light component, we have to remember what the inclination and | | sky PA are in "tilt", so we can tilt the truncation in | | the same way. Otherwise, decoupling the two (truncation and | | light profile) makes things rather unintuitive for creating | | rings and such. | \----------------------------------------------------------------*/ trunc_pars (objtype, ia, high, df, &trlink, &trval); if (r != NULL) /* A spiral component. */ tilt = trunc_spiral (r); /*----------------------------------------------------------------\ | Pre-compute some things we will frequently use later on but | | which need to be calculated one time only. | \----------------------------------------------------------------*/ if (ta[5] >= 5.) kap = 1.9992 * ta[5] - 0.3271; else splint (TWOn, KAPPA, Y2, ninterp, 2*ta[5], &kap); if (strncmp(objtype, "sersic1", 7) != 0 && strncmp(objtype, "sersic2", 7) != 0 && strncmp(objtype, "devauc1", 7) != 0 && strncmp(objtype, "devauc2", 7) != 0) sconst = ta[5] * exp(gammln(2.*ta[5])) * exp(kap)* pow(kap, -2.*ta[5]); if (ia[5] > 0) { if (ta[5] > 5.) dkapdp5 = 1.9992; else { h = 0.1; /* Hardwire inital step size for derivative */ dkapdp5 = dfridr (getinterp, ta[5], h, &err); }; if (ta[5] > 3.) h = log(ta[5]) * 0.1; /* Hardwire initial step size for deriv. */ else h = 0.1; if (strncmp(objtype, "sersic1", 7) != 0 && strncmp(objtype, "sersic2", 7) != 0 && strncmp(objtype, "devauc1", 7) != 0 && strncmp(objtype, "devauc2", 7) != 0) dgammlndn = dfridr (gamm2nln, ta[5], h, &err); }; invp5 = 1./ta[5]; kapinvp5 = kap * invp5 / pow(ta[4], invp5+1.); kapfac = kap * invp5 * invp5; /*-------------------------------------------------------------------\ | Differentiate between two scenarios: | | sersic = total magnitude as free parameter. [normtype = 0] | | sersic1 = central surface brightness as free parameter. [=1] | | sersic2 = surface brightness at r_e as free parameter. [=2] | | sersic3 = surf. brightness at r_break as free param.. [=3] | \-------------------------------------------------------------------*/ if (normtype == 0) { /*-------------------------------------------------------------\ | Note that the order of the next 3 "if" statements matters, | | because area_ratio carries through! | \-------------------------------------------------------------*/ if (c0 != NULL && f == NULL) { area_ratio = ratio (c); dratio_dc = dfridr (ratio, c, dc, &err); }; if (f != NULL) { uarea = fourarea (ta, f->a, f->ia, f->npars, c); area_ratio = a[9] * PI / uarea; /* uarea = unit area */ }; if (r != NULL) { cosi = cos (r->a[9] / 180. * PI); sini = sin (r->a[9] / 180. * PI); area_ratio = area_ratio / cosi; }; /*---------------------------------------------------------------\ | a[3] is the galaxy total brightness while ta[3] | | is the count rate at the effective radius. | | | | Notice the last term ta[9] in ta[3]. | | 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. | \---------------------------------------------------------------*/ totcnt = pow (10., (model->magzpt - a[3])/ 2.5); ta[3] = totcnt / (2. * PI * sconst * ta[4] * ta[4] * ta[9] ) * area_ratio; /*---------------------------------------------------------\ | Calculate some of the cross coupling terms between Ie | | with all the other parameters. | \---------------------------------------------------------*/ if (ia[4] > 0) dIedre = -2 * ta[3] / ta[4]; if (ia[5] > 0) dIedn = -ta[3] * (invp5 + dgammlndn + dkapdp5 - 2*(log(kap) + ta[5] / kap * dkapdp5) ); if (ia[9] > 0) dIedq = -ta[3] / ta[9]; if (f != NULL) { for (i=1; i <= f->npars; i+=2) { mode = i/2 + 1; if (f->ia[i] > 0) dIedfp[i] = -ta[3] / uarea * dratio_dfp (ampderiv, mode, ta, f->a, f->ia, f->npars, c); if (f->ia[i+1] > 0) dIedfp[i+1] = -ta[3] / uarea * dratio_dfp (phderiv, mode, ta, f->a, f->ia, f->npars, c); }; if (ia[9] > 0) dIedq = -ta[3] / uarea * dratio_dfp (qderiv, 0, ta, f->a, f->ia, f->npars, c); if (c0 != NULL && c0->ia[1] > 0) dIedc = ta[3] / uarea * dratio_dfp (c0deriv, 0, ta, f->a, f->ia, f->npars, c); } else if (c0 != NULL) dIedc = ta[3] / area_ratio * dratio_dc; if (r != NULL) { if (r->ia[9] > 0) dIedrp[9] = ta[3] / cosi * sini / 180. * PI; }; } else ta[3] = pow (10., (ta[3] - model->muzpt)/-2.5); /*-----------------------------\ | For the azimuthal profiles | \-----------------------------*/ OcosPA = cos (ta[10]); OsinPA = sin (ta[10]); /*----------------------------------------------------------------------\ | Calculate truncation transition parameters and constants for | | derivative calculations. | \----------------------------------------------------------------------*/ get_trunc_cons ("sersic", normtype, &trlink, &trval, ta, kap, OcosPA, OsinPA, &trunc_norm, &dfnorm); /*----------------------------------------------------------------------\ | The normalization radius: For central surface brightness and | | effective radius, r_norm = 0. and a[4], respectively. But, for | | normtype = 3 (truncation), r_norm = break radius of the *first* | | truncation component in the link. | \----------------------------------------------------------------------*/ if (normtype == 1) normrad = 0.; /* Central surface brightness */ if (normtype == 2) normrad = ta[4]; /* Effective radius */ if (normtype == 3) { if (trlink.lrtype == 2 && trlink.innout == 0) normrad = trlink.a[5]; /* Outer break radius */ else normrad = trlink.a[4]; /* Break radius */ }; /*------------------------------------------\ | 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, ta, ia, OcosPA, OsinPA, b, c0, f, r, &am, drd, bdrd, &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; sumdfdIe = 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, ta, ia, OcosPA, OsinPA, b, c0, f, r, &am, drd, bdrd, &drdc, fdrd, rdrd); /*-------------------------------------------------------\ | If the truncation has its own PA and axis ratio, | | calculate the truncation shape separately from | | the light profile model. Note that if the light | | profile is a tilted spiral component, the truncation | | is also tilted in the same way, given by parameters | | stored in "tilt." | \-------------------------------------------------------*/ trlptr = &trlink; trvptr = &trval; while (trlptr != NULL && trlptr->objtype != NULL) { if (trlptr->tilt == 1) tiltptr = tilt; else tiltptr = NULL; /*--------------------------------------------------\ | If the truncation type is "radial2", i.e. inner | | and outer break radius, and if it's an inner | | truncation, then the radius scale paramter is | | the outer break radius. | \--------------------------------------------------*/ if (trlptr->innout == 0 && trlptr->lrtype == 2) { for (i=1; i<=10; i++) { atmp[i] = trlptr->a[i]; iatmp[i] = trlptr->ia[i]; }; atmp[4] = trlptr->a[5]; atmp[5] = trlptr->a[4]; iatmp[4] = trlptr->ia[5]; iatmp[5] = trlptr->ia[4]; aptr = atmp; iaptr = iatmp; } else { aptr = trlptr->a; iaptr = trlptr->ia; }; drdblah (xp, yp, aptr, iaptr, trvptr->cosPA, trvptr->sinPA, trlptr->b, trlptr->c0, trlptr->f, tiltptr, &trvptr->am, trlptr->drd, trlptr->bdrd, &trlptr->drdc, trlptr->fdrd, NULL); trlptr = trlptr->next; trvptr = trvptr->next; }; /*-----------------------\ | Classical parameters | \-----------------------*/ if (am < 1./(nsubsamp*2.)) { dfdIe = fluxinteg (1./nsubsamp/sqrt(PI * ta[9]), ta, sersicfunc, kap); } else { rr = pow(am/ta[4], invp5); dfdIe = exp(-kap*(rr-1.)) / totsamp; }; /* Just the Sersic term without the truncation funcs. */ dfdIe = dfdIe / trunc_norm; flux0 = ta[3] * dfdIe; /*-------------------------------------------------------\ | What does dfdr operate on when there's a truncation? | | | | Keep in mind of the following equation: | /---/ \-----\ | f = Ie exp(-k[(r/re)^(1/n) - 1]) * { P1(r1) * P2(r2) * ... } | \---\ /-----/ | where Pn are the hyperbolic tangent truncation | | functions, and rn's are the radial coordinates when | | they have their own q and PA. Thus, for example, | | the deriv with respect to parameter z is: | /---/ \-----\ | df/dz = df/dr * dr/dz + df/dr1 * dr1/dz + df/dr2 * dr2/dz ... | \---\ /-----/ | The function values df/dr1, df/dr2, etc. are stored | | in trval->dfdam, whereas dr1/dz, dr2/dz, etc. are | | calculated by drdblah and stored into trlink->drd | | structure. | | | | So here are the scenarios: | | | | 1. If trunc has its own q and PA: dfdr only | | operates on the Sersic radial coordinate, am. | | The trunc radial coordinate is stored into | | trval->am. | | | | 2. If Sersic is *not* a spiral AND trunc | | doesn't have own q and PA: "am" is the same | | for both, so dfdr operates on both the Sersic | | and truncation profiles. | | | | 3. If Sersic is a spiral: trunc inherits q and | | PA from the Sersic, but not the higher order | | modes (or else truncation has a weird shape). | | Thus, because of the inheritance, dfdr | | will operate on trlink->am, but only on free | | parameters q and PA, and nothing else. | \-------------------------------------------------------*/ dfdr = flux0 * -kap/ta[5]/ta[4] * pow(am/ta[4],invp5-1); /*--------------------------------------------------------\ | A profile is truncated by multipling by a hyperbolic | | tangent calculated in PP. Note that outer truncation | | and inner truncation are actually two independent | | "components", components and each one has its own | | inner (a[3]) and outer (a[4]) transition radii. | \--------------------------------------------------------*/ if (trlink.objtype != NULL && normtype > 0) { PP = 1; trunc_dfdr_calc (&trval, &trlink, flux0, &PP); dfdIe = dfdIe * PP; dfdr = dfdr * PP; }; flux = ta[3] * dfdIe; /*---------------------------------------------------\ | Note that "trlink.drd" would only exist if | | there is a truncation, and "trlinks.ia" is | | what tells galfit that the truncation has its | | own x and y (ia = 1) or if it is the same as | | same as the Sersic component (ia == 3). | | ---- This is set in "trunc_cent.c" ---- | | If the truncation has its own x and y, then | | df/dx and df/dy for the Sersic component doesn't | | act on it here, but it does in the truncation | | derivs below. | | | | Note that below dfdr_P = df/dP * dP/dr. | \---------------------------------------------------*/ if (ia[1] > 0) { df->dpm[1][iy][ix] += dfdr * drd[1]; trlptr = &trlink; trvptr = &trval; while (trlptr != NULL && trlptr->drd != NULL) { if (trlptr->ia[1] == 3) df->dpm[1][iy][ix] += trvptr->dfdam * trlptr->drd[1]; trlptr = trlptr->next; trvptr = trvptr->next; }; }; if (ia[2] > 0) { df->dpm[2][iy][ix] += dfdr * drd[2]; trlptr = &trlink; trvptr = &trval; while (trlptr != NULL && trlptr->drd != NULL) { if (trlptr->ia[2] == 3) df->dpm[2][iy][ix] += trvptr->dfdam * trlptr->drd[2]; trlptr = trlptr->next; trvptr = trvptr->next; }; }; if (ia[3] > 0) df->dpm[3][iy][ix] += flux * magfac; if (ia[4] > 0) { df->dpm[4][iy][ix] += (flux* kapinvp5 * pow(am,invp5) + dfdIe * dIedre) * os; if (normtype == 3) df->dpm[4][iy][ix] += flux * dfnorm * normrad / ta[4] * os; }; if (ia[5] > 0) { df->dpm[5][iy][ix] += (flux * (kapfac * pow(am/ta[4], invp5) * log(am/ta[4]) - dkapdp5 * (pow(am/ta[4], invp5) - 1.)) + dfdIe * dIedn); if (normtype == 1) df->dpm[5][iy][ix] -= flux * dkapdp5; else if (normtype == 3) df->dpm[5][iy][ix] -= flux * (kapfac * pow(normrad/ta[4], invp5)*log(normrad/ta[4])- dkapdp5 * (pow(normrad/ta[4], invp5) - 1.)); }; if (ia[9] > 0) { df->dpm[9][iy][ix] += (dfdr * drd[9] + dfdIe * dIedq); trlptr = &trlink; trvptr = &trval; while (trlptr != NULL && trlptr->drd != NULL) { if (trlptr->ia[9] == -1) df->dpm[9][iy][ix] += trvptr->dfdam * drd[9]; trlptr = trlptr->next; trvptr = trvptr->next; }; }; if (ia[10] > 0) { df->dpm[10][iy][ix] += dfdr * drd[10]; trlptr = &trlink; trvptr = &trval; while (trlptr != NULL && trlptr->drd != NULL) { if (trlptr->ia[10] == -1) df->dpm[10][iy][ix] += trvptr->dfdam * drd[10]; trlptr = trlptr->next; trvptr = trvptr->next; }; }; /*-------------------------------------------\ | Traditional diskyness/boxyness parameter | \-------------------------------------------*/ if (c0 != NULL && c0->ia[1] > 0) c0df->dpm[1][iy][ix] += (dfdr * drdc + dfdIe * dIedc); /*-------------------------------------------\ | Bending, Fourier and PA rotational terms | \-------------------------------------------*/ if (b != NULL) high_derivs (ix,iy, b, bdf, dfdr, bdrd, 0., NULL); if (f != NULL) high_derivs (ix,iy, f, fdf, dfdr, fdrd, dfdIe, dIedfp); if (r != NULL) high_derivs (ix,iy, r, rdf, dfdr, rdrd, dfdIe, dIedrp); /*-------------------------------------------------------\ | Hybrid or meta-profile terms. Note that trlink | | outer taper and inner taper correspond to two | | different, and independent components so that | | trlink->a for inner and trlink->a for outer are | | mutually exclusive parameters. Note also that a[4] | | a[5] do not "mix" between the two. | \-------------------------------------------------------*/ trlptr = &trlink; trvptr = &trval; while (trlptr != NULL && trlptr->a != NULL) { trunc_derivs (ix,iy, trvptr->am, flux, trvptr->dfdam, trvptr, trlptr); trlptr = trlptr->next; trvptr = trvptr->next; }; sumdfdIe += dfdIe; yc = yc + 1.; }; xc = xc + 1.; }; flux = ta[3] * sumdfdIe; df->dpm[0][iy][ix] = flux; /* ... for convlution */ trlptr = &trlink; while (trlptr != NULL && trlptr->df != NULL) { if (normtype > 0) trlptr->df->dpm[0][iy][ix] = 0.; /* Trunc. are not real components */ trlptr = trlptr->next; }; }; }; free_fitvecs (&drd, &tilt, b, &bdrd, f, &fdrd, &dIedfp, r, &rdrd, &dIedrp, normtype, &trlink, &trval); } /* #if DEBUG sprintf (model->name, "model-0.fits"); Output image writefits ("test.fits", model, "", 0); #endif */ float gamm2nln(float n) { float n2, g2nln; n2 = 2. * n; g2nln = gammln(n2); return (g2nln); } float getinterp(float n) { extern float TWOn[], Y2[], KAPPA[]; extern int ninterp; float n2, y; n2 = 2. * n; splint (TWOn, KAPPA, Y2, ninterp, n2, &y); return (y); }