/**********************************************************************\ * 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 #include "nrutil.h" #include "structs.h" #include "const.h" #define BARDEF 0.23 /* bar <= 20 degree rotation (determined empirically) */ #define MINR 1e-7 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) { float or = 1e20, ir = 0., skypa = 0., incl = 0., rotang = 0., gsr = 0., stval = 0., ns = 0., rw = 1., dRdrdist = 0., dRdrw = 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_dp, dcosRdp, dsinRdp, dgamdp, ddeldp, dtheta_dp; float damdq, dtheta_dq; float rdist_or, rdist_rw, asympl, tanhr, al, dpl_drdist, ydpn, dydp_dxdp=0., dydp_da4=0., ddeldbn, or_rw, logrd, logro, flog, dlog_drdist, dtheta_da4, damda4, sgn_ydp; float c = 2.; float cosskypa, sinskypa, cosi; int n, mode, sgn; extern struct inpars input; /********************************\ * Diskyness/boxyness parameter * \********************************/ if (c0 != NULL) c += c0->a[1]; /*******************************************************\ * Simplify the nomenclature of some PA rotation terms * \*******************************************************/ if (r != NULL) { ir = r->a[1]; /* Bar radius in pixels */ or = r->a[2]; /* outer radius in pixels */ 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) { 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); tanhr = rotang * (tanh (gsr*(rdist_or)+stval) + 1); }; if (r != NULL && strncmp (r->objtype, "power", 5) == 0) { al = r->a[4]; asympl = pow( rdist_or + 1., al) / pow(2, al); dpl_drdist = al / pow(2, al) * pow (rdist_or + 1, al - 1) / or; axrotang = tanhr * asympl; dRdrdist = rotang * gsr / or * sech2 + tanhr * dpl_drdist; } if (r != NULL && strncmp (r->objtype, "log", 3) == 0) { rw = r->a[4]; /* Logarithmic winding scale radius spiral */ rdist_rw = rdist / rw; or_rw = or / rw; logrd = log (rdist_rw + 1); logro = log (or_rw + 1); flog = logrd / logro; dlog_drdist = 1. / rw / logro / (rdist_rw + 1); axrotang = tanhr * flog; dRdrdist = rotang * gsr / or * sech2 + tanhr * dlog_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) { ydpn = b->a[n] * pow (xdp / a[4], n); ydp += ydpn; dydp_dxdp += ydpn * n / xdp; dydp_da4 += -ydpn * n / a[4]; }; }; }; /*************************\ * 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 + dydp_dxdp * ddeldxd * dxd_dxc; /* <== bending modes */ dtheta_dxc = 1./xdp / a[9] * (ddeldxc - ydp/xdp * dgamdxc) * datan; drd[1] = (dam_dgamma * dgamdxc + dam_ddelta * ddeldxc + dam_dtheta * dtheta_dxc) * input.sampfac; }; 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) * input.sampfac; }; /* if (ia[4] > 0) { dtheta_da4 = dydp_da4/xdp/a[9] * datan; sgn_ydp = ydp/fabs(ydp); damda4 = drcoeff * absydpc / ydp * sgn_ydp * dydp_da4; drd[4] = damda4 + dam_dtheta * dtheta_da4; }; */ 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 / a[4], n); bdrd[n] = (dam_ddelta + dam_dtheta * datan / xdp / a[9]) * ddeldbn * input.sampfac; }; }; }; /*******************************************\ * "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) { 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)); if (strncmp (r->objtype, "power", 5) == 0) dRdir = dRdir * asympl; if (strncmp (r->objtype, "log", 3) == 0) dRdir = dRdir * flog; 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) * input.sampfac; }; if (r->ia[2] > 0) { dgsr_dor = gsr * (1./or - 1./(or - ir)); dstval_dor = -dgsr_dor; if (strncmp (r->objtype, "power", 5) == 0) dRdor = rotang * ( sech2 * (rdist_or * dgsr_dor - gsr*rdist_or/or + dstval_dor) ) * asympl - al * axrotang * rdist_or / or; if (strncmp (r->objtype, "log", 3) == 0) dRdor = rotang * ( sech2 * (rdist_or * dgsr_dor - gsr*rdist_or/or + dstval_dor) ) * flog - axrotang / logro / (or_rw + 1) / rw; 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) * input.sampfac; }; if (r->ia[3] > 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)); 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) { /********************************************************\ * dRdp is dR_dalpha for powerlaw, but dR_dw for log * * rotation function. * \********************************************************/ if (strncmp (r->objtype, "power", 5) == 0) dR_dp = axrotang * (-log(2.) + log (rdist_or + 1)); if (strncmp (r->objtype, "log", 3) == 0) { dR_dp = -axrotang / rw/rw * ( rdist/logrd / (rdist_rw + 1) - or / logro / (or_rw + 1) ); }; dcosRdp = -sinR * dR_dp; dsinRdp = cosR * dR_dp; dgamdp = xd * (OcosPA * dcosRdp - OsinPA * dsinRdp) - yd * (OcosPA * dsinRdp + OsinPA * dcosRdp); ddeldp = xd * (OcosPA * dsinRdp + OsinPA * dcosRdp) + yd * (OcosPA * dcosRdp - OsinPA * dsinRdp); dtheta_dp = 1. / xdp / a[9] * (ddeldp - ydp / xdp * dgamdp) * datan; rdrd[4] = (dam_dgamma * dgamdp + dam_ddelta * ddeldp + dam_dtheta * dtheta_dp); }; 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); }; }; }