#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 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); 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 dfdI0, float *dI0dfp); 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 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 **dI0dfp, struct fitpars *r, float **rdrd, float **dI0drp, int normtype, struct trunc_links *trlinks, struct trunc_val *trval); void exponential (char objtype[], float a[], int ia[], int normtype, 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, trunc=0; float xp, yp, am, ta[NPARS+1], OcosPA, OsinPA, sumdfdI0, xf, yf, xc, yc, flux, gexp, dfdI0, dfdr, step, err, fpa, foursum, c=2., magfac, drdc, dc=0.1, sumdfdI, uarea=1., flux0; float cosi=1., sini; float area_ratio = 1, dratio_dc = 0., totcnt; float dI0dc=0., dI0drs=0., dI0dq=0.; float *drd, *bdrd, *fdrd, *rdrd, *dI0dfp, *dI0drp; 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; 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) { dI0dfp = vector (1, f->npars); fdrd = vector (1, f->npars); for (i=1; i <= f->npars; i++) dI0dfp[i] = 0.; }; if (r != NULL) { dI0drp = vector (1, r->npars); rdrd = vector (1, r->npars); for (i=1; i <= r->npars; i++) dI0drp[i] = 0.; } if (c0 != NULL) c += c0->a[1]; /*----------------------------------------------------------------\ | Figure out whether we're trying to fit a trunc. If so, an | | outer truncation is innout = 1, whereas if there's an inner | | truncation, it's a innout = 0. If there's both, then it's | | Type 3 ttype. Also set pointers to trunc structures. | | | | 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); /*----------------------------------------------------------------\ | Compute area corrections due to high order terms. Note that | | the order of "if" statements matters. | \----------------------------------------------------------------*/ if (normtype == 0) { 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; }; if (r != NULL) { cosi = cos (fabs(r->a[9]) / 180. * PI); sini = sin (fabs(r->a[9]) / 180. * PI); area_ratio = area_ratio / cosi; }; /*----------------------------------------------------------------\ | 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. | \----------------------------------------------------------------*/ 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. | \---------------------------------------------------------------*/ dI0drs = -2 * ta[3] / ta[4]; dI0dq = -ta[3] / ta[9]; if (f != NULL) { 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] / uarea * dratio_dfp (c0deriv, 0, ta, f->a, f->ia, f->npars, c); } else if (c0 != NULL) dI0dc = ta[3] / area_ratio * dratio_dc; if (r != NULL) { if (r->ia[9] > 0) dI0drp[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 ("exponential", normtype, &trlink, &trval, ta, 0., OcosPA, OsinPA, &trunc_norm, &dfnorm); if (normtype == 1) normrad = 0.; /* Central surface brightness */ if (normtype == 2) normrad = ta[4]; /* Scale 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 && am <= os) { /* If scalelength is small, oversample by a lot more. */ nsubsamp = IMIN(100, (int)(1./fabs(ta[4] * ta[9]) * 20 *OVERSAMPR/ os)); } else if (fabs(am) <= 3. * os) nsubsamp = 20 * 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; 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, 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 | \-----------------------*/ gexp = -am / ta[4]; dfdI0 = exp(gexp) / totsamp / trunc_norm; flux0 = ta[3] * dfdI0; dfdr = flux0 / -ta[4]; /*--------------------------------------------------------\ | 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); dfdI0 = dfdI0 * PP; dfdr = dfdr * PP; }; /*--------------------------------------------------------\ | Now put it all together with the normal parameters. | | Truncation P included here | \--------------------------------------------------------*/ flux = ta[3] * dfdI0; /*--------------------------------------------------------\ | What does dfdr operate on when there's a truncation? | | See "sersic.c" for more details. | \--------------------------------------------------------*/ 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*am/ta[4]/ta[4] + dfdI0*dI0drs) * os; if (normtype == 3) df->dpm[4][iy][ix] += flux * dfnorm * normrad / ta[4] * os; }; if (ia[9] > 0) { df->dpm[9][iy][ix] += (dfdr * drd[9] + dfdI0 * dI0dq); 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 + dfdI0 * dI0dc); /********************************************\ * 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, dfdI0, dI0dfp); if (r != NULL) high_derivs (ix, iy, r, rdf, dfdr, rdrd, dfdI0, dI0drp); /*-------------------------------------------------------\ | 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; }; sumdfdI0 += dfdI0; yc = yc + 1.; }; xc = xc+1.; }; flux = ta[3] * sumdfdI0; 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, &dI0dfp, r, &rdrd, &dI0drp, normtype, &trlink, &trval); } /* #if DEBG sprintf (model->name, "model-0.fits"); Output image writefits ("test.fits", model, 0); #endif */