2d:null 2e:["$","title",null,{"children":"Plotting a black hole | Max Bo | Observable"}] 33:T478,Luminet describes the creation of this image in [his own blog](https://blogs.futura-sciences.com/e-luminet/2018/03/07/45-years-black-hole-imaging-1-early-work-1972-1988/): > However, lacking at the time of an appropriate drawing software, I had to create it by hand. Using numerical data from the computer, I drew directly on negative Canson paper with black India ink, placing dots more densely where the simulation showed more light – a rather painstaking process! Next, I took the negative of my negative to get the positive, the black points becoming white and the white background becoming black. The result converged into a pleasantly organic, asymmetrical form, as visually engaging as it was scientifically revealing. Obviously, I wished to reproduce this image. Fortunately, [Bjorge Meulemeester](https://github.com/bgmeulem) beat me to it with his Python project, the aptly named [Luminet](https://github.com/bgmeulem/Luminet). Using his [reproduced data](https://github.com/bgmeulem/Luminet/tree/master/Points), I decided to render it with the new [`Plot.raster`](https://observablehq.com/@observablehq/plot-raster) mark anyway: 34:T815,fluxLines = { const makeFluxLines = (extraMarks) => Plot.plot({ margin: 2, width: 585, aspectRatio: 1, x: { domain: xDomain, ticks: false, label: null }, y: { domain: yDomain, ticks: false, label: null }, marks: [ Plot.contour(points_incl85, { x: "X", y: "Y", fill: "flux_o", interpolate: "barycentric", interval: 0.0000000000004, fillOpacity: 0, stroke: "black" }), ...extraMarks ] }); // Time to label the flux lines const fluxLines = makeFluxLines([]); // Get all the contour paths // Drop the axis we've accidentally scooped up const paths = _.tail([ ...d3 .select(fluxLines) .select('g[aria-label="contour"]') .selectAll('g[aria-label="contour"] > path') ]); // figure out what each of the contour paths actually corresponds to scale wise const [lo, hi] = fluxLines.scale("color").domain; const step = hi / paths.length; const ticks = d3.range(lo + step, hi + step, step); const labels = paths.flatMap((d, i) => { const mkLabel = () => { // choose some random point along each contour const length = d.getTotalLength() * Math.random(); return { x: fluxLines.scale("x").invert(d.getPointAtLength(length).x), y: fluxLines.scale("y").invert(d.getPointAtLength(length).y), text: Math.round(ticks[i] * 10e12).toString() + "e-12" }; }; const minimumLabelDistance = 200; // create labels the number of times that a label can fit into the path with room return d3 .range(0, Math.round(d.getTotalLength() / minimumLabelDistance) + 1) .map(() => mkLabel()); }); // now we've figured out where the labels should be, rerender the plot with the labels added too return makeFluxLines([ Plot.text(labels, { x: "x", y: "y", fontSize: 7, text: "text", fill: "black", stroke: "white" }) ]); }35:T76c,The following is a simplified version of Meulemeester's [plotting code](https://github.com/bgmeulem/Luminet/blob/27b21ca2d2c1701fa6689c15e0a66e95d49c1cea/black_hole.py#L351): ```py def plot_points(self, power_scale=.9, levels=100): """ Plot the points written out by samplePoints() :param levels: amount of levels in matplotlib contour plot :param power_scale: powers_cale to apply to flux. No power_scale = 1. Anything lower than 1 will make the dim points pop out more. :return: """ def plot_direct_image(_ax, points, _levels, _min_flux, _max_flux, _power_scale): # direct image points.sort_values(by="angle", inplace=True) points_ = points.iloc[[b_ <= self.get_apparent_outer_edge_radius(a_) for b_, a_ in zip(points["impact_parameter"], points["angle"])]] fluxes = [(abs(fl + _min_flux) / (_max_flux + _min_flux)) ** _power_scale for fl in points_['flux_o']] _ax.tricontourf(points_['X'], points_['Y'], fluxes, cmap='Greys_r', levels=_levels, norm=plt.Normalize(0, 1), nchunk=2 ) br = self.calc_apparent_inner_disk_edge() _ax.fill_between(br.X, br.Y, color='black', zorder=1) # to fill Delauney triangulation artefacts with black return _ax _fig, _ax = self.get_figure() points1 = pd.read_csv(f"Points/points_incl={round(self.t * 180 / np.pi)}.csv") points2 = pd.read_csv(f"Points/points_secondary_incl={round(self.t * 180 / np.pi)}.csv") max_flux = max(max(points1['flux_o']), max(points2['flux_o'])) min_flux = 0 _ax = plot_direct_image(_ax, points1, levels, min_flux, max_flux, power_scale) _ax.set_xlim((-40, 40)) _ax.set_ylim((-40, 40)) plt.savefig('SampledPoints_incl={}.png'.format(self.t), dpi=300, facecolor='black') plt.show() return _fig, _ax ```36:[{"id":"b7b7b5a68cc0635a36ddbd0234f85334a96ef919cb3ca82ac85bce6181a7384ace33d4557027aaa77f066004859411ad34fee4f37140c623bb66a02bedf8061e","href":"https://static.observableusercontent.com/files/b7b7b5a68cc0635a36ddbd0234f85334a96ef919cb3ca82ac85bce6181a7384ace33d4557027aaa77f066004859411ad34fee4f37140c623bb66a02bedf8061e","name":"points_incl=85.csv","type":"text/csv","lastModified":1675152244214.948,"size":5875580},{"id":"198b5c1c1f44d34cc96f1c344a052d403ff46640d85225ae55b03f38bcbb45cb945404d2ab7a311f5d4b1e7293d3f5ea77f2e9f14dab47c09ac4b3c67a2d602d","href":"https://static.observableusercontent.com/files/198b5c1c1f44d34cc96f1c344a052d403ff46640d85225ae55b03f38bcbb45cb945404d2ab7a311f5d4b1e7293d3f5ea77f2e9f14dab47c09ac4b3c67a2d602d","name":"points_secondary_incl=85.csv","type":"text/csv","lastModified":1675152244485.013,"size":5858475},{"id":"6ed3c5018b2433133bea1e15ac31710c6b7219a924d6003927b707c12587ebffe178a9b30c32151b857ad9949d6084bf5793e609b71a3869a1ccb5265b6a4fb0","href":"https://static.observableusercontent.com/files/6ed3c5018b2433133bea1e15ac31710c6b7219a924d6003927b707c12587ebffe178a9b30c32151b857ad9949d6084bf5793e609b71a3869a1ccb5265b6a4fb0","name":"image.png","type":"image/png","lastModified":1675155116898.94,"size":336156},{"id":"5ba330fbc2c93a12578a32458ed74b4d6751cd3cef61dbac52264bb5fd27092f50b558818af9324dcebc78f1db316979d29c9c5f7c4b8dd5b636746d2f2d6d1a","href":"https://static.observableusercontent.com/files/5ba330fbc2c93a12578a32458ed74b4d6751cd3cef61dbac52264bb5fd27092f50b558818af9324dcebc78f1db316979d29c9c5f7c4b8dd5b636746d2f2d6d1a","name":"image@1.png","type":"image/png","lastModified":1675163208704.842,"size":84090}] 31:{"id":"10c61d86c3bea569","slug":"plotting-a-black-hole","title":"Plotting a black hole","thumbnail":"edf7022d90f454c779da55f3f4a6e6bd178614ae88033b2b585c5b2b77033091","access_level":"public","create_time":"$D2023-01-31T08:03:28.551Z","update_time":"$D2023-06-30T05:16:06.337Z","view_time":null,"first_public_time":"$D2023-01-31T08:03:28.551Z","fork":false,"stars":8,"starred":false,"owner":{"type":"team","id":"cc48459e4293f942","login":"mjbo","name":"Max Bo","avatar_url":"https://avatars-next.observableusercontent.com/05ae9baba40510a4a22fa7d179f29dd4ffd85fcff859c54a31d42e166ae61dfa","tier":"starter_2024"},"version":1406,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"5bb83a6e9493ea1d","login":"mjbo","name":"Max Bo","avatar_url":"https://avatars-next.observableusercontent.com/05ae9baba40510a4a22fa7d179f29dd4ffd85fcff859c54a31d42e166ae61dfa","tier":null}],"members":[{"type":"individual","id":"5bb83a6e9493ea1d","login":"mjbo","name":"Max Bo","avatar_url":"https://avatars-next.observableusercontent.com/05ae9baba40510a4a22fa7d179f29dd4ffd85fcff859c54a31d42e166ae61dfa","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":1,"collection":{"id":"a5615d262872a60a","slug":"misc","title":"Misc","owner":{"id":"cc48459e4293f942","login":"mjbo","name":"Max Bo"}},"type":"notebook","body":{"title":"Plotting a black hole","cells":[{"id":0,"value":"# Plotting a black hole\n\nMy friend recently came into a possession of a pen plotter. In want of something to plot, I remembered one of my favorite \"hand plotted\" works, Jean-Pierre Luminet's [*Image of a spherical black hole with thin accretion disk* (1979)](https://articles.adsabs.harvard.edu/pdf/1979A%26A....75..228L), the first simulated image of a Schwarzschild black hole:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":324,"value":"![](https://blogs.futura-sciences.com/e-luminet/wp-content/uploads/sites/11/2018/03/TN-JPL1978-72dpi.jpg)\n

\n Image\tof\ta\tSpherical\tBlack\tHole\twith\tThin\tAccretion\tDisk.\tAs\texpected,\t the\tmain\t\ncharacteristic\tof\tthe\timage\tis\tthe\tstrong\tasymmetry\tof\tthe\tdisk’s\tbrightness,\tso\tthat\tone\tside\t\nis\tfar\tbrighter\tand\tthe\tother\tis\tfar\tdimmer.\tAnother\tone\tis\tthat,\talthough\tthe\tupper\tside\tof\tthe\t\ndisk\tis\tcompletely\tvisible,\tonly\ta\tsmall\tpart\tof\tthe\tlower\tside\tis\tobservable.\tThis\tis\tdue\tto\tthe\t\nfact\tthat in\ta\trealistic\tsituation,\tthe\tgaseous\tdisk\tis\topaque;\ttherefore\tit\tabsorbs\tthe\tlight\trays\t\nthat\t it\t intercepts.\t It\t follows\t that\t a\t major\t part\t of\t the\t secondary\t image\t is\t occulted\t by\t the\t\nprimary\timage,\tits\tvisible\tpart\tbeing stuck\taround\tthe\tedge\tof\tthe\tblack\thole,\tlike\ta\tgleaming\nhalo. [Source](https://arxiv.org/pdf/1902.11196.pdf)\n

\n\n\n\n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":340,"value":"$33","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":590,"value":"viewof powerScale = Inputs.range([0, 1], {\n label: \"Power scale\",\n value: 0.6,\n step: 0.1\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":6,"value":"Plot.plot({\n x: {\n domain: xDomain,\n ticks: 0,\n label: null\n },\n y: {\n domain: yDomain,\n ticks: 0,\n label: null\n },\n color: {\n range: [\"black\", \"white\"],\n type: \"linear\"\n },\n style: {\n width,\n background: \"black\",\n color: \"white\"\n },\n marks: [\n // Plot.contour(points_incl85, {\n // x: \"X\",\n // y: \"Y\",\n // interval: 0.000000000000004,\n // interpolate: \"random-walk\",\n // fill: (d) => d.flux_o\n // }),\n Plot.raster(points_incl85, {\n x: \"X\",\n y: \"Y\",\n fill: (d) => normalizeFlux(d.flux_o),\n interpolate: \"random-walk\",\n blur: 5,\n imageRendering: \"pixelated\"\n }),\n Plot.raster(points_incl85, {\n x: \"X\",\n y: \"Y\",\n fill: (d) => normalizeFlux(d.flux_o),\n imageRendering: \"pixelated\"\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":418,"value":"We can reproduce some of Luminet's other diagrams","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":411,"value":"

![Curves of constant apparent flux in units of the maximum intrinsic flux, as seen by an observer at 10° above the disk’s plane. The luminosity is maximum in the regions where high intrinsic luminosity combines with high blueshift.](https://blogs.futura-sciences.com/luminet/wp-content/uploads/sites/2/2015/02/courbes-flux-constant.jpg)

Curves of constant apparent flux in units of the maximum intrinsic flux, as seen by an observer at 10° above the disk’s plane. The luminosity is maximum in the regions where high intrinsic luminosity combines with high blueshift. [Source](https://blogs.futura-sciences.com/e-luminet/2018/03/07/45-years-black-hole-imaging-1-early-work-1972-1988/)

","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":474,"value":"with the new [`Plot.contour`](https://observablehq.com/@observablehq/plot-contour) mark:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1164,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":255,"value":"---\n## Implementation\n\n### Data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":257,"value":"https://github.com/bgmeulem/Luminet/tree/master/Points","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":4,"value":"points_incl85 = await FileAttachment(\"points_incl=85.csv\").csv({typed: \"auto\"})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"points_incl85","database":"$undefined"},{"id":535,"value":"Plot.plot({\n x: {\n domain: xDomain\n },\n y: {\n domain: yDomain\n },\n marks: [\n Plot.vector(points_incl85, {\n x: \"X\",\n y: \"Y\",\n stroke: \"flux_o\",\n rotate: (d) => (d.angle * 180) / Math.PI\n })\n ]\n})\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":5,"value":"points_secondary_incl85 = await FileAttachment(\"points_secondary_incl=85.csv\").csv({typed: \"auto\"})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"points_secondary_incl85","database":"$undefined"},{"id":1371,"value":"I have no idea how to unify this as a meaningful visualization: ","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1368,"value":"Plot.plot({\n x: {\n domain: xDomain\n },\n y: {\n domain: yDomain\n },\n marks: [\n Plot.vector(points_secondary_incl85, {\n x: \"X\",\n y: \"Y\",\n stroke: \"flux_o\",\n rotate: (d) => (d.angle * 180) / Math.PI\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":508,"value":"### Plotting","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":493,"value":"xDomain = [-30, 30]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":495,"value":"yDomain = [-17.5, 17.5]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":524,"value":"$35","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":533,"value":"We will attempt to reproduce some of the logic to improve the rendering. ","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":522,"value":"maxFlux = _.maxBy(points_incl85.map((d) => d.flux_o))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":518,"value":"normalizeFlux = (flux_o) => (flux_o / maxFlux) ** powerScale","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W36","stdlib":"1","titling":"auto"}} 32:"$31:body"