(******************************************************************************) (* listModel.math: Code to list Plummer-smoothed gamma and NFW models. *) (******************************************************************************) (* rhoGamma: density profile for gamma models. *) rhoGamma[gamma_, r_, a_:1, M_:1] := (3 - gamma) a M / (4 Pi r^gamma (r + a)^(4 - gamma)) (* massGamma: cumulative mass profile for gamma models. *) massGamma[gamma_, r_, a_:1, M_:1] := M (r / (r + a))^(3 - gamma) (******************************************************************************) (* rhoTaper: density profile in taper region (r > b). Equation is rho = (1 + mu) rhostar (b / r)^2 Exp[-r / rstar] *) rhoTaper[r_, taper_] := (1 + taper[[5]]) taper[[4]] (taper[[1]] / r)^2 Exp[-r / taper[[3]]] (* massTaper: cumulative mass profile in taper region (b to r where r > b). *) massTaper[r_, taper_] := (1 + taper[[5]]) 4 Pi taper[[4]] taper[[1]]^2 taper[[3]] * (Exp[2+taper[[2]]] - Exp[-r / taper[[3]]]) (* taperParamsGamma: build list of parameters for tapered Gamma model. *) taperParamsGamma[gamma_, a_, b_, M_] := Module[{beta, rstar, rhostar, mu}, beta = (b / rhoGamma[gamma, b, a, M]) * (D[rhoGamma[gamma, r, a, M], r] /. r -> b); rstar = - b / (2 + beta); rhostar = rhoGamma[gamma, b, a, M] Exp[- (2 + beta)]; mu = M / (massGamma[gamma, b, a, M] + 4 Pi b^2 rstar rhoGamma[gamma, b, a, M]) - 1; List[b, beta, rstar, rhostar, mu]] (******************************************************************************) rhoGammaTaper[gamma_, r_, a_, M_, taper_] := If[r <= taper[[1]], (1 + taper[[5]]) rhoGamma[gamma, r, a, M], rhoTaper[r, taper]] massGammaTaper[gamma_, r_, a_, M_, taper_] := If[r <= taper[[1]], (1 + taper[[5]]) massGamma[gamma, r, a, M], (1 + taper[[5]]) massGamma[gamma, taper[[1]], a, M] + massTaper[r, taper]] (******************************************************************************) rhoGammaSoft[gamma_, r_, a_, M_, epsilon_] := (3 epsilon^2 / 2) * NIntegrate[rhoGamma[gamma, Sqrt[R^2 + (z - r)^2], a, M] * R (R^2 + z^2 + epsilon^2)^(-5/2), {R, 0, Infinity}, {z, -Infinity, Infinity}, Method -> "DuffyCoordinates", MaxRecursion -> 32] rhoGammaTaperSoft[gamma_, r_, a_, M_, taper_, epsilon_] := Module[{alpha = Max[10, (taper[[1]] + r) / taper[[1]]]}, rhoGammaTaperSoft1[gamma, r, a, M, taper, epsilon] + Re[rhoTaperSoft2[r, a, M, taper, epsilon, alpha]] + rhoTaperSoft3[r, a, M, taper, epsilon, alpha] + rhoTaperSoft4[r, a, M, taper, epsilon, alpha]] rhoGammaTaperSoft1[gamma_, r_, a_, M_, taper_, epsilon_] := (3 epsilon^2 / 2) (1 + taper[[5]]) * NIntegrate[rhoGamma[gamma, Sqrt[R^2 + (z - r)^2], a, M] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - taper[[1]], r + taper[[1]]}, {R, 0, Sqrt[taper[[1]]^2 - (z - r)^2]}, MinRecursion -> 64, MaxRecursion -> 128] rhoTaperSoft2[r_, a_, M_, taper_, epsilon_, alpha_:10] := (3 epsilon^2 / 2) * NIntegrate[rhoTaper[Sqrt[R^2 + (z - r)^2], taper] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - taper[[1]], r, r + taper[[1]]}, {R, Sqrt[taper[[1]]^2 - (z - r)^2], alpha taper[[1]]}, MinRecursion -> 64, MaxRecursion -> 128] rhoTaperSoft3[r_, a_, M_, taper_, epsilon_, alpha_:10] := (3 epsilon^2 / 2) * If[r <= taper[[1]], NIntegrate[rhoTaper[Sqrt[R^2 + (z - r)^2], taper] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - (alpha + 1) taper[[1]], r - taper[[1]]}, {R, 0, alpha taper[[1]]}, Method -> {"GlobalAdaptive", MinRecursion -> 32, MaxRecursion -> 128}], NIntegrate[rhoTaper[Sqrt[R^2 + (z - r)^2], taper] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - (alpha + 1) taper[[1]], 0, r - taper[[1]]}, {R, 0, Sqrt[z^2 + epsilon^2]/2, alpha taper[[1]]}, Method -> {"GlobalAdaptive", MinRecursion -> 32, MaxRecursion -> 128}]] rhoTaperSoft4[r_, a_, M_, taper_, epsilon_, alpha_:10] := (3 epsilon^2 / 2) * NIntegrate[rhoTaper[Sqrt[R^2 + (z - r)^2], taper] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r + taper[[1]], r + (alpha + 1) taper[[1]]}, {R, 0, alpha taper[[1]]}, MinRecursion -> 64, MaxRecursion -> 128] (******************************************************************************) (* listGammaModel: tabulate a smoothed Gamma model. *) listGammaModel[gam_:1, a_:1, b_:0, M_:1, eps_:1/256, file_:"", iter_:{2, -15, 12, 1/2}] := Module[{ostr = If[file != "", OpenWrite[file], \$Output], base = First[iter], taper = If[b > 0, taperParamsGamma[gam, a, b, M], {}], r, rho, mass, rhoEps, rhoApp}, WriteString[ostr, StringForm["# listGammaModel[``, ``, ``, ``, ``, ``, ``]\n", InputForm[gam], InputForm[a], InputForm[b], InputForm[M], InputForm[eps], InputForm[file], InputForm[iter]]]; WriteString[ostr, "# r", " rho(r)", " mass(r)", " rho_eps(r)", " rho(r_eps)\n"]; If[Length[iter] == 5 && iter[[5]], If[b == 0, rhoEps = rhoGammaSoft[gam, 0.0, a, M, eps]; rhoApp = rhoGamma[gam, eps, a, M]]; If[b != 0, rhoEps = rhoGammaTaperSoft[gam, 0.0, a, M, taper, eps]; rhoApp = rhoGammaTaper[gam, eps, a, M, taper]]; WriteString[ostr, " ", NumStr[0], " Infinity", " ", NumStr[0], " ", NumStr[rhoEps], " ", NumStr[rhoApp], "\n"]]; Do[r = base^x; If[b == 0, rho = rhoGamma[gam, r, a, M]; mass = massGamma[gam, r, a, M]; rhoEps = rhoGammaSoft[gam, r, a, M, eps]; rhoApp = rhoGamma[gam, Sqrt[r^2 + eps^2], a, M]]; If[b != 0, rho = rhoGammaTaper[gam, r, a, M, taper]; mass = massGammaTaper[gam, r, a, M, taper]; rhoEps = rhoGammaTaperSoft[gam, r, a, M, taper, eps]; rhoApp = rhoGammaTaper[gam, Sqrt[r^2 + eps^2], a, M, taper]]; WriteString[ostr, " ", NumStr[r], " ", NumStr[rho], " ", NumStr[mass], " ", NumStr[rhoEps], " ", NumStr[rhoApp], "\n"], {x, iter[[2]], iter[[3]], iter[[4]]}]; If[file != "", Close[ostr]]] (******************************************************************************) (* rhoNFW: density profile for NFW models. *) rhoNFW[r_, a_:1, rho0_:1/(2Pi)] := (rho0 a^3)/(r (a + r)^2) (* massNFW: cumulative mass profile for NFW models. *) massNFW[r_, a_:1, rho0_:1/(2Pi)] := 4 Pi rho0 a^3 (Log[(a + r)/a] - r/(a + r)) (* rhoNFWTaper: density profile for NFW models with SW tapering. *) rhoNFWTaper[r_, a_, rho0_, taper_] := If[r <= taper[[1]], rhoNFW[r, a, rho0], rhoNFWTaper1[r, taper[[1]], taper[[2]], taper[[3]], taper[[4]]]] rhoNFWTaper1[r_, b_, eta_, a_, rhostar_] := rhostar (b / r)^eta Exp[- r / a] (* massNFWTaper: cumulative mass profile for NFW models with SW tapering. *) massNFWTaper[r_, a_, rho0_, taper_] := If[r <= taper[[1]], massNFW[r, a, rho0], massNFW[taper[[1]], a, rho0] + massNFWTaper1[r, taper[[1]], taper[[2]], taper[[3]], taper[[4]]]] massNFWTaper1[r_, b_, eta_, a_, rhostar_] := 4 Pi rhostar a^3 (b/a)^eta (Gamma[3-eta, b/a] - Gamma[3-eta, r/a]) (* taperParamsNFW: build list of parameters for tapered NFW model. *) taperParamsNFW[a_, b_, rho0_] := Module[{eta, rhostar}, eta = (a + 3 b) / (a + b) - b/a; rhostar = rhoNFW[b, a, rho0] Exp[b/a]; List[b, eta, a, rhostar]] (******************************************************************************) rhoNFWSoft[r_, a_, rho0_, epsilon_, acc_:12] := (3 epsilon^2 / 2) * NIntegrate[rhoNFW[Sqrt[R^2 + (z - r)^2], a, rho0] * R (R^2 + z^2 + epsilon^2)^(-5/2), {R, 0, Infinity}, {z, -Infinity, Infinity}, Method -> "DuffyCoordinates", MaxRecursion -> 32, AccuracyGoal -> acc] rhoNFWTaperSoft[r_, a_, rho0_, taper_, epsilon_] := Module[{alpha = Max[10, (taper[[1]] + r) / taper[[1]]]}, rhoNFWTaperSoft1[r, a, rho0, taper, epsilon] + Re[rhoNFWTaperSoft2[r, taper, epsilon, alpha]] + rhoNFWTaperSoft3[r, taper, epsilon, alpha] + rhoNFWTaperSoft4[r, taper, epsilon, alpha]] rhoNFWTaperSoft1[r_, a_, rho0_, taper_, epsilon_] := (3 epsilon^2 / 2) * NIntegrate[rhoNFW[Sqrt[R^2 + (z - r)^2], a, rho0] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - taper[[1]], r + taper[[1]]}, {R, 0, Sqrt[taper[[1]]^2 - (z - r)^2]}, MinRecursion -> 64, MaxRecursion -> 128] rhoNFWTaperSoft2[r_, taper_, epsilon_, alpha_:10] := (3 epsilon^2 / 2) * NIntegrate[rhoNFWTaper1[Sqrt[R^2 + (z - r)^2], taper[[1]], taper[[2]], taper[[3]], taper[[4]]] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - taper[[1]], r, r + taper[[1]]}, {R, Sqrt[taper[[1]]^2 - (z - r)^2], alpha taper[[1]]}, MinRecursion -> 64, MaxRecursion -> 128] rhoNFWTaperSoft3[r_, taper_, epsilon_, alpha_:10] := (3 epsilon^2 / 2) * If[r <= taper[[1]], NIntegrate[rhoNFWTaper1[Sqrt[R^2 + (z - r)^2], taper[[1]], taper[[2]], taper[[3]], taper[[4]]] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - (alpha + 1) taper[[1]], r - taper[[1]]}, {R, 0, alpha taper[[1]]}, Method -> {"GlobalAdaptive", MinRecursion -> 32, MaxRecursion -> 128}], NIntegrate[rhoNFWTaper1[Sqrt[R^2 + (z - r)^2], taper[[1]], taper[[2]], taper[[3]], taper[[4]]] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r - (alpha + 1) taper[[1]], 0, r - taper[[1]]}, {R, 0, Sqrt[z^2 + epsilon^2]/2, alpha taper[[1]]}, Method -> {"GlobalAdaptive", MinRecursion -> 32, MaxRecursion -> 128}]] rhoNFWTaperSoft4[r_, taper_, epsilon_, alpha_:10] := (3 epsilon^2 / 2) * NIntegrate[rhoNFWTaper1[Sqrt[R^2 + (z - r)^2], taper[[1]], taper[[2]], taper[[3]], taper[[4]]] * R (R^2 + z^2 + epsilon^2)^(-5/2), {z, r + taper[[1]], r + (alpha + 1) taper[[1]]}, {R, 0, alpha taper[[1]]}, MinRecursion -> 64, MaxRecursion -> 128] (******************************************************************************) (* listNFWModel: tabulate a smoothed NFW model. *) listNFWModel[a_:1, b_:0, rho0_:1/(2Pi), eps_:1/256, file_:"", iter_:{2, -15, 12, 1/2}, acc_:12] := Module[{ostr = If[file != "", OpenWrite[file], \$Output], base = First[iter], taper = If[b > 0, taperParamsNFW[a, b, rho0], {}], r, rho, mass, rhoEps, rhoApp}, WriteString[ostr, StringForm["# listNFWModel[``, ``, ``, ``, ``, ``, ``]\n", InputForm[a], InputForm[b], InputForm[rho0], InputForm[eps], InputForm[file], InputForm[iter], InputForm[acc]]]; WriteString[ostr, "# r", " rho(r)", " mass(r)", " rho_eps(r)", " rho(r_eps)\n"]; If[Length[iter] == 5 && iter[[5]], If[b == 0, rhoEps = rhoNFWSoft[0.0, a, rho0, eps, acc]; rhoApp = rhoNFW[eps, a, rho0]]; If[b != 0, rhoEps = rhoNFWTaperSoft[0.0, a, rho0, taper, eps]; rhoApp = rhoNFWTaper[eps, a, rho0, taper]]; WriteString[ostr, " ", NumStr[0], " Infinity", " ", NumStr[0], " ", NumStr[rhoEps], " ", NumStr[rhoApp], "\n"]]; Do[r = base^x; If[b == 0, rho = rhoNFW[r, a, rho0]; mass = massNFW[r, a, rho0]; rhoEps = rhoNFWSoft[r, a, rho0, eps, acc]; rhoApp = rhoNFW[Sqrt[r^2 + eps^2], a, rho0]]; If[b != 0, rho = rhoNFWTaper[r, a, rho0, taper]; mass = massNFWTaper[r, a, rho0, taper]; rhoEps = rhoNFWTaperSoft[r, a, rho0, taper, eps]; rhoApp = rhoNFWTaper[Sqrt[r^2 + eps^2], a, rho0, taper]]; WriteString[ostr, " ", NumStr[r], " ", NumStr[rho], " ", NumStr[mass], " ", NumStr[rhoEps], " ", NumStr[rhoApp], "\n"], {x, iter[[2]], iter[[3]], iter[[4]]}]; If[file != "", Close[ostr]]] (******************************************************************************) NumStr[x_,prec_:8] := ToString[PaddedForm[N[x], {prec+1,prec}, ExponentFunction -> (#&), NumberFormat -> (SequenceForm[#1, If[StringTake[#3,1] == "-", "e", "e+"], #3] &)]] (******************************************************************************) listStandardModels[iter_:{2, -15, 5, 1/8, True}] := Module[{}, Print[listGammaModel[1, 1, 0, 1, 1/1, "modH_001.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/2, "modH_002.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/4, "modH_004.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/8, "modH_008.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/16, "modH_016.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/32, "modH_032.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/64, "modH_064.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/128, "modH_128.txt", iter]]; Print[listGammaModel[1, 1, 0, 1, 1/256, "modH_256.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/1, "modJ_001.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/2, "modJ_002.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/4, "modJ_004.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/8, "modJ_008.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/16, "modJ_016.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/32, "modJ_032.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/64, "modJ_064.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/128, "modJ_128.txt", iter]]; Print[listGammaModel[2, 1, 0, 1, 1/256, "modJ_256.txt", iter]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/1, "modNFW_001.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/2, "modNFW_002.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/4, "modNFW_004.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/8, "modNFW_008.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/16, "modNFW_016.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/32, "modNFW_032.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/64, "modNFW_064.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/128, "modNFW_128.txt", iter, 16]]; Print[listNFWModel[1, 0, 1/(2Pi), 1/256, "modNFW_256.txt", iter, 16]]] listTaperedModels[iter_:{2, -20, 12, 1/32, True}] := Module[{}, Print[listGammaModel[1, 1, 100, 1, 1/16, "modHt_016.txt", iter]]; Print[listGammaModel[1, 1, 100, 1, 1/64, "modHt_064.txt", iter]]; Print[listGammaModel[1, 1, 100, 1, 1/256, "modHt_256.txt", iter]]; Print[listGammaModel[2, 1, 100, 1, 1/16, "modJt_016.txt", iter]]; Print[listGammaModel[2, 1, 100, 1, 1/64, "modJt_064.txt", iter]]; Print[listGammaModel[2, 1, 100, 1, 1/256, "modJt_256.txt", iter]]] listTaperedNFWModels[iter_:{2, -15, 7, 1/32, True}] := Module[{}, Print[listNFWModel[1, 5, 1/(2Pi), 1/16, "modNFW_016_05.txt", iter]]; Print[listNFWModel[1, 5, 1/(2Pi), 1/64, "modNFW_064_05.txt", iter]]; Print[listNFWModel[1, 5, 1/(2Pi), 1/256, "modNFW_256_05.txt", iter]]; Print[listNFWModel[1, 10, 1/(2Pi), 1/16, "modNFW_016_10.txt", iter]]; Print[listNFWModel[1, 10, 1/(2Pi), 1/64, "modNFW_064_10.txt", iter]]; Print[listNFWModel[1, 10, 1/(2Pi), 1/256, "modNFW_256_10.txt", iter]]; Print[listNFWModel[1, 15, 1/(2Pi), 1/16, "modNFW_016_15.txt", iter]]; Print[listNFWModel[1, 15, 1/(2Pi), 1/64, "modNFW_064_15.txt", iter]]; Print[listNFWModel[1, 15, 1/(2Pi), 1/256, "modNFW_256_15.txt", iter]]] (******************************************************************************)