/**********************************************************************\ * This subroutine consolidates all the derivatives having to do with * * the radial coordinate into one place. This makes calculating the * * the derivatives with respect to all the parameters more elegant * * because the calculations here can be shared by all profiles, * * such as Sersic, moffat, Nuker, etc.. What do I mean by this? To * * calculate the derivative of df/dPA, for example, you use the chain * * rule: * * df/dPA = df/dr * (dr/dx * dx/dPA + dr/dy * dr/dPA + * * dr/dtheta*dtheta/dPA) * * It's useful to note that beyond the first term, dfdr, every * * model profile will share the same algorithm that calculates * * everything else in the parenthesis. That's what this subroutine * * provides upon output, and why I named this subroutine dr/dblah.... * \**********************************************************************/ #include #include "nrutil.h" #include "structs.h" #include "const.h" #define BARDEF 0.349066 #define MINR 1e-7 void drdblah (float x, float y, float c, 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) { float or = 1e20, ir = 0., skypa = 0., incl = 0., rotang = 0., gsr = 0., stval = 0., ns = 0., rt = 1., dRdrdist = 0., dRdrt = 0., fabsra; float xd, yd, xd0, yd0, rdist, cosPA, sinPA, xdp, ydp, ydp9, axrotang=0., absxdpc, absydpc, fabsrc, am0, xytheta, foursum, fpa, drcoeff, xdprime, ydprime, xp, yp, cosR, sinR; float dFdtheta, dam_dgamma, dam_ddelta, dam_dtheta, dgamdcosPA, dgamdsinPA, ddeldsinPA, ddeldcosPA, sech2, dgamdxd, dgamdyd, ddeldxd, ddeldyd, datan; float dxd_dxc, dyd_dxc, drdistdxc, dRdxc, dcosRdxc, dsinRdxc, dcosPAdxc, dsinPAdxc, dgamdxc, ddeldxc, dtheta_dxc; float dxd_dyc, dyd_dyc, drdistdyc, dRdyc, dcosRdyc, dsinRdyc, dcosPAdyc, dsinPAdyc, dgamdyc, ddeldyc, dtheta_dyc; float dcosPA_dPA, dsinPA_dPA, dgamdPA, ddeldPA, dtheta_dPA; float sech2S; float dgsr_dir, dstval_dir, dRdir, dcosPAdir, dsinPAdir, dgamdir, ddeldir, dtheta_dir, dRdns, dcosPAdns, dsinPAdns, dgamdns, ddeldns, dtheta_dns; float dgsr_dor, dstval_dor, dRdor, dcosPAdor, dsinPAdor, dgamdor, ddeldor, dtheta_dor; float dgsrdra, dstvaldra, dRdra, dcosRdra, dsinRdra, dgamdra, ddeldra, dtheta_dra, datanh; float drdist_di, dRdi, dcosPAdi, dsinPAdi, dyd_di, dgamdi, ddeldi, dtheta_di; float dxd_dskypa, dyd_dskypa, drdist_dsky, dR_dskypa, dcosPA_dskypa, dsinPA_dskypa, dgam_dskypa, ddel_dskypa, dtheta_dskypa; float dR_dal, dcosRdal, dsinRdal, dgamdal, ddeldal, dtheta_dal; float damdq, dtheta_dq; float rdist_or, rdist_rt, asympl, tanhr, logr, al, dpl_drdist, dydp_dxdp=0., ddeldbn; float cosskypa, sinskypa, cosi; int n, mode, sgn; /*******************************************************\ * Simplify the nomenclature of some PA rotation terms * \*******************************************************/ if (r != NULL) { if (strncmp (r->objtype, "tanh", 4) == 0) { ir = r->a[1]; /* Bar radius in pixels */ or = r->a[2]; /* outer radius in pixels */ } else if (strncmp (r->objtype, "log", 3) == 0) { ns = r->a[1]; /* Transition sharpness from bar to spiral */ rt = r->a[2]; /* Radius of transition */ }; rotang = -r->a[3]/360.*PI; /* Total rotation out to "or" in degrees */ fabsra = fabs (rotang); incl = r->a[9] / 180. * PI; skypa = -r->a[10]/180. * PI; }; /******************************************************\ * Rotate a disk in projection to the apparent sky PA * \******************************************************/ cosskypa = cos(skypa); sinskypa = sin(skypa); cosi = cos (incl); /* Cosine of the disk inclination */ xd0 = x-a[1]; yd0 = y-a[2]; xp = xd0 * cosskypa - yd0 * sinskypa + a[1]; yp = xd0 * sinskypa + yd0 * cosskypa + a[2]; /* Project the model by the inclination angle to line of *\ \* sight. This is effectively the axis ratio. */ yp = (yp - a[2])/cosi + a[2]; /******************************************\ * Calculate the circular center distance * \******************************************/ xd = xp - a[1]; yd = yp - a[2]; if (xd == 0.) xd = MINR; if (yd == 0.) yd = MINR; /***************************************************************\ * For the hyperbolic tangent rotation: * * Rotate the position angle according to central distance. * * First calculate the speed of coordinate rotation as a * * function of central distance. This controls the tightness * * of the spiral arm winding. The coefficient "rotang" * * has units of degrees -- that is to say, from the inner * * radius R = ir out to R = or the PA will have rotated through * * by that many degrees total. So then the average rotation * * *rate* for the tanh function is simply: * * rate = rotang / (or - ir) [degrees / pixel]. * * * * The outer radius (or) is defined to be the location where * * the position angle of the winding spiral arm reaches 4% * * the asymptotic value. * * * * The inner radius (ir) is where the bar is 4% of the final * * asymptotic value, and can be used to model a galaxy bar. * \***************************************************************/ rdist = sqrt(xd*xd+yd*yd); if (r != NULL && strncmp (r->objtype, "tanh", 4) == 0) { /* 20 degrees = 0.349066 radian bar definition. */ al = r->a[4]; gsr = (2 - atanh (2 *BARDEF/(fabsra + BARDEF) - 1.00001)) * or/(or-ir); stval = 2 - gsr; rdist_or = rdist / or; sech2 = pow(cosh (gsr * (rdist_or) + stval), -2); asympl = pow( rdist_or + 1., al) / pow(2, al); dpl_drdist = al / pow(2, al) * pow (rdist_or + 1, al - 1) / or; tanhr = rotang * (tanh (gsr*(rdist_or)+stval) - tanh(stval)); axrotang = tanhr * asympl; dRdrdist = rotang * gsr / or * sech2 + tanhr * dpl_drdist; } else if (r != NULL && strncmp (r->objtype, "log", 3) == 0) { al = r->a[4]; rdist_rt = rdist / rt; dRdrt = -rotang * ns / log(2.) / rdist * pow(rdist_rt, ns+1) / (pow (rdist_rt, ns) + 1); asympl = pow( rdist_rt + 1., r->a[4]) / pow(2, r->a[4]); dpl_drdist = al / pow(2, al) * pow (rdist_or + 1, al - 1) / or; logr = rotang / log(2.) * log (pow(rdist_rt, ns) + 1); axrotang = logr * asympl; dRdrdist = rotang / log(2.) * ns * pow(rdist_rt, ns-1) / rt / (pow (rdist_rt, ns) + 1) + tanhr * dpl_drdist; }; cosR = cos (axrotang); sinR = sin (axrotang); cosPA = (cosR * OcosPA - sinR *OsinPA); sinPA = sinR * OcosPA + cosR * OsinPA; xdp = xd * cosPA - yd * sinPA; absxdpc = FMAX(MINR, pow(fabs(xdp), c)); ydp = xd * sinPA + yd * cosPA; /**********************************************\ * Take care of coordinate bending transform * \**********************************************/ if (b != NULL) { for (n=1; n <= b->npars; n++) { if (!(b->a[n] ==0 && b->ia[n] ==0) && b->ia[n] != -1) ydp += b->a[n] * pow (xdp, n); dydp_dxdp += b->a[n] * pow (xdp, n-1); }; }; /*************************\ * Traditional ellipsoid * \*************************/ ydp9 = ydp/a[9]; /********************************************\ * Diskyness/boxyness coordinate transform. * \********************************************/ absydpc = FMAX(MINR, pow(fabs(ydp9), c)); fabsrc = FMAX(MINR, absxdpc + absydpc); if (xdp == 0.) xdp = MINR; if (ydp == 0.) ydp = MINR; am0 = FMAX(MINR, pow(fabsrc, 1./c)); /**************************************\ * Fourier index coordinate transform * \**************************************/ xytheta = atan2 (ydp9, xdp); foursum = 1.; dFdtheta = 0.; if (f != NULL) { for (n=1; n <= f->npars; n+=2) { if (!(f->a[n] ==0 && f->ia[n] ==0) && f->ia[n] != -1) { mode = n/2+1; fpa = -f->a[n+1] / 180. * PI; foursum += f->a[n]* cos (mode *(xytheta + fpa)); dFdtheta += -f->a[n] * mode * sin (mode *(xytheta + fpa)); }; }; }; *am = am0 * fabs(foursum); sgn = foursum/fabs(foursum); /*****************************\ * For classical parameters, * * note that: * * gamma = xdp * * delta = ydp * \*****************************/ dgamdxd = cosPA; dgamdyd = -sinPA; ddeldxd = sinPA + dydp_dxdp * cosPA; /* The latter comes from bending */ ddeldyd = cosPA - dydp_dxdp * sinPA; datan = 1./ (1+ pow(ydp9/xdp,2)); drcoeff = *am / (absxdpc + absydpc); dam_dgamma = xdp * pow(fabs(xdp), c-2) * drcoeff; dam_ddelta = ydp9 * pow(fabs(ydp9), c-2) * drcoeff / a[9]; dam_dtheta = am0 * dFdtheta * sgn; dgamdcosPA = xd; dgamdsinPA = -yd; ddeldsinPA = xd - dydp_dxdp * yd; ddeldcosPA = yd + dydp_dxdp * xd; if (ia[1] > 0) { dxd_dxc = -cosskypa; dyd_dxc = -sinskypa / cosi; drdistdxc = 1./rdist * (xd * dxd_dxc + yd * dyd_dxc); /* For PA rotation */ dRdxc = dRdrdist * drdistdxc; dcosRdxc = -sinR * dRdxc; dsinRdxc = cosR * dRdxc; dcosPAdxc = OcosPA * dcosRdxc - OsinPA * dsinRdxc; dsinPAdxc = OcosPA * dsinRdxc + OsinPA * dcosRdxc; dgamdxc = dgamdxd * dxd_dxc + dgamdyd * dyd_dxc + dgamdcosPA * dcosPAdxc + dgamdsinPA * dsinPAdxc; ddeldxc = ddeldxd * dxd_dxc + ddeldyd * dyd_dxc + ddeldcosPA * dcosPAdxc + ddeldsinPA * dsinPAdxc; dtheta_dxc = 1./xdp / a[9] * (ddeldxc - ydp/xdp * dgamdxc) * datan; drd[1] = dam_dgamma * dgamdxc + dam_ddelta * ddeldxc + dam_dtheta * dtheta_dxc; }; if (ia[2] > 0) { dxd_dyc = sinskypa; dyd_dyc = -cosskypa / cosi; drdistdyc = 1./rdist * (xd * dxd_dyc + yd * dyd_dyc); /* For PA rotation */ dRdyc = dRdrdist * drdistdyc; dcosRdyc = -sinR * dRdyc; dsinRdyc = cosR * dRdyc; dcosPAdyc = OcosPA * dcosRdyc - OsinPA * dsinRdyc; dsinPAdyc = OcosPA * dsinRdyc + OsinPA * dcosRdyc; dgamdyc = dgamdxd * dxd_dyc + dgamdyd * dyd_dyc + dgamdcosPA * dcosPAdyc + dgamdsinPA * dsinPAdyc; ddeldyc = ddeldxd * dxd_dyc + ddeldyd * dyd_dyc + ddeldcosPA * dcosPAdyc + ddeldsinPA * dsinPAdyc; dtheta_dyc = 1./xdp / a[9] * (ddeldyc - ydp/xdp * dgamdyc) * datan; drd[2] = dam_dgamma * dgamdyc + dam_ddelta * ddeldyc + dam_dtheta * dtheta_dyc; }; if (ia[9] > 0) { damdq = -drcoeff * absydpc/a[9]; dtheta_dq = - ydp9/xdp/a[9] * datan; drd[9] = damdq + dam_dtheta * dtheta_dq; }; if (ia[10] > 0) { dcosPA_dPA = -OsinPA * cosR - OcosPA * sinR; dsinPA_dPA = -OsinPA * sinR + OcosPA * cosR; dgamdPA = dgamdcosPA * dcosPA_dPA + dgamdsinPA * dsinPA_dPA; ddeldPA = ddeldsinPA * dsinPA_dPA + ddeldcosPA * dcosPA_dPA; dtheta_dPA = 1./a[9] / xdp * (ddeldPA - ydp/xdp * dgamdPA) * datan; drd[10] = -PI/180. * (dam_dgamma * dgamdPA + dam_ddelta * ddeldPA + dam_dtheta * dtheta_dPA); }; /*****************\ * Bending terms * \*****************/ if (b != NULL) { for (n=1; n <= b->npars; n++) { if (b->ia[n] > 0) { ddeldbn = pow (xdp, n); bdrd[n] = (dam_ddelta + dam_dtheta * datan / xdp / a[9]) * ddeldbn; }; }; }; /*******************************************\ * "Traditional" diskyness/boxyness term * \*******************************************/ if (c0 != NULL && c0->ia[1] > 0) { xdprime = log(FMAX(MINR, fabs(xdp))) * absxdpc; ydprime = log(FMAX(MINR, fabs(ydp9))) *absydpc; *drdc = *am * (-1./c/c * log(fabsrc) + 1./c/ fabsrc * (xdprime + ydprime)); }; /*****************\ * Fourier terms * \*****************/ if (f != NULL) { for (n=1; n <= f->npars; n+=2) { if (!(f->ia[n]==0 && f->a[n]==0) && f->ia[n] != -1) { mode = n/2+1; fpa = -f->a[n+1] / 180. * PI; /* sgn is for derivative of absolute value */ if (f->ia[n] > 0) fdrd[n] = am0 * cos (mode * (xytheta+fpa)) * sgn; if (f->ia[n+1] > 0) fdrd[n+1] = -mode * am0 * f->a[n] * sin(mode * (xytheta + fpa)) * -PI/180 * sgn; }; }; }; /*********************************\ * Position angle rotation terms * \*********************************/ if (r != NULL) { if (strncmp (r->objtype, "tanh", 4) == 0) { sech2S = pow(cosh (stval), -2.); if (r->ia[1] > 0) { dgsr_dir = gsr / (or - ir); dstval_dir = -dgsr_dir; dRdir = rotang * (sech2 * (rdist_or * dgsr_dir + dstval_dir) - sech2S * dstval_dir); dcosPAdir = (-OcosPA * sinR - OsinPA * cosR) * dRdir; dsinPAdir = ( OcosPA * cosR - OsinPA * sinR) * dRdir; dgamdir = xd * dcosPAdir - yd * dsinPAdir; ddeldir = xd * dsinPAdir + yd * dcosPAdir; dtheta_dir = 1./ xdp / a[9] * (ddeldir - ydp / xdp * dgamdir) * datan; rdrd[1] = dam_dgamma * dgamdir + dam_ddelta * ddeldir + dam_dtheta * dtheta_dir; }; } else if (strncmp (r->objtype, "log", 3) == 0) { if (r->ia[1] > 0) { dRdns = rotang / log(2.) * log(rdist_rt) / (1 + pow(rt / rdist, ns)); dcosPAdns = (-OcosPA * sinR - OsinPA * cosR) * dRdns; dsinPAdns = ( OcosPA * cosR - OsinPA * sinR) * dRdns; dgamdns = xd * dcosPAdns - yd * dsinPAdns; ddeldns = xd * dsinPAdns + yd * dcosPAdns; dtheta_dns = 1./ xdp / a[9] * (ddeldns - ydp / xdp * dgamdns) * datan; rdrd[1] = dam_dgamma * dgamdns + dam_ddelta * ddeldns + dam_dtheta * dtheta_dns; }; }; /*******************************************************************\ * Below are mostly common between tanh and log rotation functions * \*******************************************************************/ if (r->ia[2] > 0) { if (strncmp (r->objtype, "tanh", 4) == 0) { dgsr_dor = gsr * (1./or - 1./(or - ir)); dstval_dor = -dgsr_dor; dRdor = rotang * (sech2 * (rdist_or * dgsr_dor - gsr*rdist_or/or + dstval_dor) - sech2S * dstval_dor) * asympl - al * axrotang * rdist_or / or; } else if (strncmp (r->objtype, "log", 3) == 0) { /* "or" here is "rt", the only radius term */ dRdor = dRdrt - al * axrotang * rdist_or / or; }; dcosPAdor = (-OcosPA * sinR - OsinPA * cosR) * dRdor; dsinPAdor = ( OcosPA * cosR - OsinPA * sinR) * dRdor; dgamdor = xd * dcosPAdor - yd * dsinPAdor; ddeldor = xd * dsinPAdor + yd * dcosPAdor; dtheta_dor = 1. / xdp / a[9] * (ddeldor - ydp/xdp * dgamdor) * datan; rdrd[2] = dam_dgamma * dgamdor + dam_ddelta * ddeldor + dam_dtheta * dtheta_dor; }; if (r->ia[3] > 0) { if (strncmp (r->objtype, "tanh", 4) == 0) { datanh = 1./(1-pow (2*BARDEF /(fabsra+BARDEF)-1.00001, 2.)); dgsrdra = datanh * or / (or - ir) / (fabsra+BARDEF)/ (fabsra+BARDEF) * 2 * BARDEF * rotang / fabsra; dstvaldra = -dgsrdra; dRdra = axrotang / rotang + rotang * (sech2 * (dgsrdra * rdist_or + dstvaldra) - sech2S * dstvaldra); } else if (strncmp (r->objtype, "log", 3) == 0) { /* rotation angle here is just the angle *\ \* at truncation radius */ dRdra = axrotang / rotang; }; dcosRdra = -sinR * dRdra; dsinRdra = cosR * dRdra; dgamdra = xd * (OcosPA * dcosRdra - OsinPA * dsinRdra) - yd * (OcosPA * dsinRdra + OsinPA * dcosRdra); ddeldra = xd * (OcosPA * dsinRdra + OsinPA * dcosRdra) + yd * (OcosPA * dcosRdra - OsinPA * dsinRdra); dtheta_dra = 1. / xdp / a[9] * (ddeldra - ydp / xdp * dgamdra) * datan; rdrd[3] = -PI/360.* (dam_dgamma * dgamdra + dam_ddelta * ddeldra + dam_dtheta * dtheta_dra); }; if (r->ia[4] > 0) { dR_dal = axrotang * (-log(2.) + log (rdist_or + 1)); dcosRdal = -sinR * dR_dal; dsinRdal = cosR * dR_dal; dgamdal = xd * (OcosPA * dcosRdal - OsinPA * dsinRdal) - yd * (OcosPA * dsinRdal + OsinPA * dcosRdal); ddeldal = xd * (OcosPA * dsinRdal + OsinPA * dcosRdal) + yd * (OcosPA * dcosRdal - OsinPA * dsinRdal); dtheta_dal = 1. / xdp / a[9] * (ddeldal - ydp / xdp * dgamdal) * datan; rdrd[4] = (dam_dgamma * dgamdal + dam_ddelta * ddeldal + dam_dtheta * dtheta_dal); }; if (r->ia[9] > 0) { dyd_di = yd * tan (incl); drdist_di = 1./rdist * yd * dyd_di; dRdi = dRdrdist * drdist_di; dcosPAdi = (-OcosPA * sinR - OsinPA * cosR) * dRdi; dsinPAdi = ( OcosPA * cosR - OsinPA * sinR) * dRdi; dgamdi = xd * dcosPAdi - yd * dsinPAdi - dyd_di * sinPA; ddeldi = xd * dsinPAdi + yd * dcosPAdi + dyd_di * cosPA; dtheta_di = 1. / xdp / a[9] * (ddeldi - ydp / xdp * dgamdi) * datan; rdrd[9] = PI / 180. *(dam_dgamma * dgamdi + dam_ddelta * ddeldi + dam_dtheta * dtheta_di); }; if (r->ia[10] > 0) { dxd_dskypa = (-xd0 * sinskypa - yd0 * cosskypa); dyd_dskypa = (xd0 * cosskypa - yd0 * sinskypa)/cosi; drdist_dsky = 1./rdist * (xd * (-xd0 * sinskypa - yd0 * cosskypa) + yd/cosi * (xd0 * cosskypa - yd0 * sinskypa)); dR_dskypa = dRdrdist * drdist_dsky; dcosPA_dskypa = (-OcosPA * sinR - OsinPA * cosR) * dR_dskypa; dsinPA_dskypa = ( OcosPA * cosR - OsinPA * sinR) * dR_dskypa; dgam_dskypa = cosPA * dxd_dskypa - sinPA * dyd_dskypa + xd * dcosPA_dskypa - yd * dsinPA_dskypa; ddel_dskypa = sinPA * dxd_dskypa + cosPA * dyd_dskypa + xd * dsinPA_dskypa + yd * dcosPA_dskypa; dtheta_dskypa = 1. / xdp / a[9] * (ddel_dskypa - ydp / xdp * dgam_dskypa) * datan; rdrd[10] = -PI/180. * (dam_dgamma * dgam_dskypa + dam_ddelta * ddel_dskypa + dam_dtheta * dtheta_dskypa); }; }; }