Follow the reluctant adventures in the life of a Welsh astrophysicist sent around the world for some reason, wherein I photograph potatoes and destroy galaxies in the name of science. And don't forget about my website, www.rhysy.net



Showing posts with label CG project. Show all posts
Showing posts with label CG project. Show all posts

Saturday, 8 April 2017

Tutorial : Let's Get Dirty

As you may or may not know, I've jumped on board the VR bandwagon. In order to create 360 3D VR content with Blender you have little choice but to use the Cycles rendering engine, which for me at least is a fairly steep learning curve from the old material editor. Recently I got quite frustrated trying to use Cycles to do the trivial task of using one stencil to act as a stencil for another, so here's a quick, basic tutorial on how to do this from the perspective of someone new to Cycles but used to the original Blender Internal engine and material editor. I don't guarantee that this is the best way to do things, just the one that makes the most sense to me.

A very nice tutorial on some slightly more complex materials can be found here, but it's an hour long video. I won't cover nearly as much as that in this short post, but hopefully I'll give enough of the basics that the Cycles workflow won't seem so scary by the end. It really isn't that different to creating standard materials, it's just that you have to think more about what each button is doing and configure those buttons yourself.

I'll assume a basic working knowledge of Blender (where all the buttons are, how to organise windows etc.) but little or nothing about Cycles materials. I'll cover :

  • How to create standard image and procedural textures
  • How to use bump mapping
  • How to mix textures using different shaders or by mixing their colour/alpha information directly
  • How to use one texture as a stencil for another, with some Cycles-specific quirks to watch out for
  • How to arrange things to make the Cycles node window seem less intimidating.


1) Getting started

First, change the rendering engine to Cycles (top menu bar, the drop-down menu that says "Blender Render" by default). Then set up a window environment suitable for editing the materials. In Cycles this is done almost entirely via nodes, so we're going to need a nice big window for that. I like the following setup :
The left section shows the 3D view of the scene, the middle one is the node editor, and the right is the standard properties window.
Notice that on the left I've set the shading method to "rendered" (the little sphere at the bottom of the left panel, next to "object mode"). This gives a realtime preview of the scene. I've got a single Sun lamp pointing at an angle to the plane, but in more complex scenes you'll see the full effects of shading so this can be a great way to see what you're material will really look like. This setting is also available for the standard Blender engine. You also get a preview of the material in the usual panel, just like with the regular rendering engine.

A word of warning : the realtime preview can be unstable. Save your work often, and if in doubt, revert to the classical technique of rendering after making changes to the material. In the Rendering menu there are various options you can set to speed up rendering time; for this example, try lowering the number of samples in the Sampling panel.

Anyway, select your object and add a material just like normal. In the node editor you'll see the following :


What is this ? The "Diffuse BSDF" thing is scary and not at all an intuitive label. What it is is the way in which the diffuse component of the material is constructed. This is then linked to the material's output, in this case to the surface. I won't look at volumetric materials or displacement settings in this tutorial, just the regular surface shaders. We also don't need to worry about any of the material settings on the right, though it can be helpful to open the "preview" panel, here shown closed.

The first thing to understand about this Diffuse node is that it contains more information than colour. If you were to create a node containing only colour information, such as an image file, and link it to the material output surface socket directly, it wouldn't work. It needs that other information contained in the Diffuse node, even if all you wanted to do was set the colour.

Fortunately if you do just want to set the colour, that's easy. You can click on the colour picker in the Diffuse node and set it there directly. Or, more usefully, you can generate the colour from other nodes (in very complex combinations) and use that as the input to the node. Let's start with an image texture node : Shift+A -> Texture -> Image Texture. Use the file selector in the Image Texture node to choose the image file to use. Then connect the Image Texture node to the Diffuse node by clicking on its colour output socket and dragging the line to the input colour socket of the Diffuse node.


But, confusingly, although the material preview will look correct, the render preview won't. Thing is this setup by itself doesn't work. Just like with regular materials, we have to tell Blender what texture coordinates we use. So now add a Texture Coordinate node : Shift+A -> Input -> Texture Coordinate. Draw a line from its "generated" output socket to the "vector" input socket of the Image Texture node, and both the render and material previews should be what we'd expect.


Hooray ! Now you may be wondering about the other sockets on the Diffuse node. "Roughness" changes the shading style but its effects are subtle so let's not worry about that. "Normal" is used for bump maps. The way you set this up is to add a "Bump" node (from the "Vector" section of the Add (shift+A) menu) :


Just like trying to connect the image texture directly to the material output, if you'd tried to connect it
directly to the Normal socket on the Diffuse node, it wouldn't work - the Bump node converts the data into a format it can process correctly. But note how here I've cheated by using the colour of the material as a measure of its height. This works quite well in this case, and in many other cases too, but you could also use a completely different source for the bump information - another image or a procedural texture, for example.

You might also be wondering how we can alter the texture scaling. For that we need to add a Mapping node from the Vector submenu. A useful trick here - first move the Texture Coordinate node a little way off to the left. Then add the Mapping node and position it over the connecting link between the Texture Coordinate and Image Texture nodes. Notice how the link turns orange. When you left click to accept the new node, the links will be set up automatically. Then you can alter the texture scaling parameters in the usual way.

Now it's starting to look quite a bit like a regular Blender material, except that you can both see and control how each part relates to the rest.
One final extremely useful thing to be aware of : frames. Frames are (sort of but not really) meta-nodes, a way to group nodes together so they can be easily moved in blocks for tidy screen organisation. They also let you label blocks of nodes so you can see at a glance what each part is doing. Add a frame from the Layout option of the add menu. Then you can stick the nodes inside it just like parenting objects together : (shift) select the node(s) you want to add in the frame, then shift-select the frame and press CTRL+P. The frame automatically reshapes to hold the nodes. You can un-parent them in the usual way with ALT+P.

Label the frame with the N menu and edit the "Label" parameter near the top. I got rid of the Mapping node here because it's not necessary for this example.
If you select the frame you can move all its nodes around at once. Obviously this isn't much use for this simple material, but, as we shall see, it becomes essential for more complex materials - especially if you re-use them months later !

Now that we have the very basics out of the way, let's get on with setting up multiple textures in different distributions.


2) Multiple textures

Note how I chose the frame to include the whole shader in the above example. That's because one way of mixing different textures is to generate a new shader and then mix them together, so it makes sense to have each frame contain a complete shader. But we don't have to do it that way. A perhaps more intuitive method is to generate the other texture and have the colours combined before they're input into the Diffuse node. You can have shaders made up of multiple textures mixed with other shaders, so things can get extremely elaborate - and far more powerful than what you can do with Blender's standard material editor.

For example, let's add a second texture on top of our image texture. Let's first create a simple procedural noise texture on its own. For this we need two additional nodes : Noise Texture (from the Texture submenu) and the ColourRamp (from the Converter submenu for some reason - I'm not sure the layout of the add menu is particularly sensible, but never mind). Set them up like this :

You can have multiple Material Output nodes if you like, so you can create this node setup without having to remove the existing one. The output nodes that's used (i.e. the material that's displayed) will be the one that's active.
It should be fairly self-explanatory - the parameters are very similar to Blender's standard internal noise textures.

Now comes the interesting part : overlaying the noise texture on top of the original concrete texture. Let's first try this by using a different diffuse shader for each texture and mixing them. We only need to add one more node to do this - the Mix Shader node from the Shader submenu. The base layer (the concrete in this case) goes into the top socket and the overlaid "dirt" - the white fuzzy noise goes into the top. The really important bit is that the "Fac" socket which controls the mixing should come from the Alpha socket of the dirt colour ramp.

This alpha-based mixing was automatic in the original Blender materials. Now we have to specify it, which is a bit more work but also gives us more control.
You can see that this gives us exactly what we'd expect - a concrete texture with some white splodges, with the dominance of the white controlled exactly by its alpha value (if you need to modify that alpha value, send it through a Math node (Converter submenu) before it goes into the Fac socket). Places which should be totally white are totally white. All is well with the world.

Now let's try the other approach of mixing the colours before we send them into a single Diffuse node. We eliminate the Diffuse node of the dirt texture and (optionally) unparent the Diffuse node from the concrete frame. Then we add a MixRGB node (colour submenu) and, equivalent to what we did before, we put the colour output of the concrete texture into the first slot and the colour of the white patches into the second. Again we use the dirt's alpha channel to control the mixing.


Uh-oh ! This has sort-of worked, but very badly. The dirt is mixed, but now its dominance is not properly determined by its alpha. Parts which should be totally white still clearly show the concrete image. Why is this ? Watch what happens if we eliminate the Bump node :


Ahh, now we see what's going on. It wasn't that the colours of the concrete image were showing through, it's that the normal map was causing bumps that were independent of the dirt texture. If we insist on using the mixing colours rather than shaders method, there is a way we can keep the bump mapping - we use the dirt texture to set the strength of the bumps. We need to use the alpha output of the dirt for this, but since we want the white patches to be smooth we need to invert it (Colour-> Invert).


A bit messy, but it works. The nice thing about this is that by altering the "Fac" slider in the Invert node we can change the strength of the effect, so that you could still see the bumps of the concrete in the white patches a little bit if you wanted to. You could get the same effect with the mixing shaders method too, just by dragging the output of the Bump node to the second Diffuse shader as well as the first.

Of course you don't have to use a simple noise texture - you could use a Musgrave, which has a lot more parameters, or a wave (equivalent to the old "Marble" texture)... or more interestingly, you can link another texture to the input scale and distortion sockets. You can get arbitrarily complicated effects like this. And because of the node setup, you're no longer limited to arranging things in the strictly-linear manner of the original material editor.


3) Controlling where the dirt appears

We might be happy having great white splodges all over the place, but we'll often prefer it if they only appeared in a few places on the material. That is, we want patches controlling where the patches appear. Especially for very small patches of dirt - you might not want a uniform coating of dirt, flecks of paint of whatever over your whole material.


This had me stumped for the longest time, and then I felt very silly when I realised how easy it was. Then I felt a bit less silly when some kindly soul pointed out to me that a small mistake can have big consequences here.

All we need to do is make another texture that defines where the splodges appear and use this to multiply the colours of the dirt. So our distribution texture should be white where we want the dirt to appear and black where we don't. Or we can use its alpha channel, which is slightly simpler that setting the RGB colours, but it doesn't really matter. Then we use a Math node to multiply the two alphas (or colours) together, by setting the method to "Multiply" instead of the default "Mix".


The colour of the dirt is generated by a single noise texture, and it's mixed with the base material in the same way as when we had a single texture. The only difference is that now the alpha mixing is done not by the alpha channel of the dirt itself, but by a multiplied combination of the dirt's alpha and the dirt distribution texture alpha. So dirt distribution section only affects the distribution of the dirt and nothing else.

You can download the .blend file here. Play with the scale parameters in both the the dirt distribution and dirt textures to see what effect they have. If you change the dirt distribution scale, you should change the large patches where the dirt appears, whereas if you change the dirt scale, you should change the size of the smaller sub-patches.

But be careful ! There are two pitfalls here :

1) To see the effect really clearly, there should be a large difference in the scale parameters - the dirt scale should be around 10x larger than the dirt distribution (larger scale means more texture within the same area). But if you go to the other extreme, you can break it - a larger scale of the dirt texture than its distribution scale means you won't see any sub-patches.

Keeping the scale of the distribution texture at 1 but altering the scale of the dirt texture from left to right (10, 25, 50).
Keeping the scale of the dirt texture at 50 but varying the scale of the distribution texture from left to right (1, 25, 50). The distribution scale has the strongest impact at small values - once it gets as high as the dirt scale everything becomes uniform, so varying it just shrinks the size of the individual patches rather than controlling where those patches appear.

2) The position of the colourband sliders matters. Keep the distribution colourband unchanged and try moving the dirt colourband sliders around. If you move them too far to the right, the sub-patches can become too small to see or disappear entirely. Too far to the left and the sub-patches become larger than the patches, so changing the scale doesn't have any affect. A similar problem happens if you move the dirt distribution colourband sliders - it can seriously screw everything up.
A good guideline is to start with both the colourbands having markers in the same positions, then tweak it from there. So if you're making your own material, start simple and gradually build up the complexity - make sure each stage has the effect you think it should have.

Keeping all the texture scales constant but varying the position of the dirt colourband controllers (keeping their relative offset the same but moving them both to the right). In the left image the size of the gaps in the dirt are so small that they're barely visible, whereas in the right they're so large that the dirt is barely visible. This is not what happens with a single texture, where you'd have to move the colourband markers to the extreme ends to get such a strong effect - here you get unexpected behaviour if the markers are just a bit too far from the centre. And in those situations altering the scaling of the textures will give some very strange results. 

Well that about wraps this up. There are lots more aspects to creating materials, obviously, but hopefully that's enough to break the ice and show that Cycle materials aren't so scary. Since non-video tutorials are in short supply, I'll try and do more of these assuming that stuff doesn't keep getting in the way.

Saturday, 16 April 2016

Nope


Sometime late last summer I saw a job advert I was morally obligated to apply for. Astronomy Visualisation Specialist ? I am one already ! Experience with visualisation software and generic programming languages, e.g. Python ? The 11,000 lines of Python code for Blender that I wrote ought to tick that box and then some. New ways of visualising data ? Assisting astronomers cre... look, just visit my website. It's all there. All of it. Never have I seen a job description that felt so precisely tailored to me.


The application deadline was 1st October, though I submitted mine well before that. It's normal in astronomy for responses to take at least one month, sometimes two or even three. A few places - which are downright rude - never bother to respond. Still, by early December I was beginning to suspect that despite being objectively very, very qualified for the job, they must have given it to someone else. Well, it did say, "as soon as possible" on the job description.

They hadn't. I'm not sure if you'd call it an early Christmas present or not but I had a Skype interview on 18th December. Which was the day after I moved out of my flat in Prague (roomate left, couldn't possibly afford the place on my own) which had involved a week of hauling heavy suitcases back and forth to move my stuff to the institute. And it was the day after I got back to Cardiff, just to make things as frantic as possible.

Anyway it went well, but unfortunately it went well for everyone else as well. After a rather nervous Christmas, in early January I got en email saying that they'd go to a second stage round where they'd send us all a data set to visualise. Which they duly did a couple of weeks later.

It was actually quite a fun little project to work on, because the data set wasn't in a format I was familiar with. 3D data sets generally describe the density or temperature or whatever at different locations in space. The location is specified by 3 positions : x, y, z. Nothing very complicated about that.


And that's fine if, as is usually the case, your data set describes something that's roughly box-shaped. And by roughly I mean very roughly indeed, like this :

Simulation of a star-forming filamentary cloud, or something.
But this data set didn't use ordinary "Cartesian" coordinates, it used spherical polar coordinates. These aren't difficult either, but they may be unfamiliar. Instead of specifying 3 linear distances from the origin, they specify one distance and two angles :


Why in the world would you want to use such things ? Surely, it's more intuitive to think in terms of distances, not angles ! No, not always. There's a very simple everyday example that should help you understand : maps. With a street map, you could easily specify a position in Cartesian coordinates. You could say, for example, that Cardiff Castle is about 150m north and 25m west of the Revolution bar, if you thought that breaking into the castle on a Saturday night was somehow a good idea.

You could also specify how high the castle keep is, if it was vitally important to reach a precise level for some reason.
On a scale this small, the fact that the Earth is curved doesn't matter. You could hold out your arm and say, "go 50 metres in that direction" and no-one would have any difficulty. But if you said, "go 5,000 miles in that direction", anyone taking you literally would have ended up in space. Of course, they intuitively understand that you mean "along the surface of the Earth", not really, "in the path followed by a perfectly straight laser beam going in that direction". Unless they're a cat, of course.

This is why you don't give cats directions using laser pointers, because if it involves going into space then damn it that's what they'll do.
North, south, east and west are really just angles. Nothing very complicated about that : if you want to go to Australia you can say it's around 140 degrees east of Great Britain and 20 degrees south of the equator. Or you can give this in miles, it's the same thing.

We don't normally specify the distance from the centre r unless you're a mountaineer, pilot, miner, or deep sea diver. OK, you'd normally give distance from sea level rather than the centre of the Earth, but it's the same thing.
Or is it ? Well, not quite. 50 miles north or south is the same everywhere, unless you're so close to a pole you can't actually go that far. E.g. if you head in a northerly direction when you're 25 miles south of the north pole, after 25 miles you'll find yourself heading south. Much worse is the case of walking east or west if you're near the pole. You can't actually go east or west if you're standing on the pole itself, and if you're just a few steps away from the pole, then walking 25 miles east or west is going to involve walking around in a lot of circles until you get dizzy.

And at the north pole it will also involve discovering Santa's secret hideout or being eaten by a polar bear.
Using angles makes things a lot easier for cartographers. Line of longitude (east or west position) have constant angular separation, even though the physical distance between them varies (i.e. 1 degree involves walking a much larger distance at the equator than near the poles). And the mathematics to convert between the two is easy and precise, so if we have two lat-long positions we can easily compute how far we have to travel to get from one to the other - even if we're at weird positions like the poles.

In numerical simulations, polar coordinates have some other advantages, which are a bit more complicated. Imagine if you will a cat on a record player*. If you wanted to specify any point on the cat, you could give its x,y position. Or you could state its r,φ (pronounced "phi") coordinates instead. There's not really any advantage to either... unless the record player is turned on and it stars spinning.

*Conjecture : there is no scientific concept which cannot be explained with the right cat gif.

If that happens it's very much easier to specify how fast each point is moving in the φ direction. If you wanted to specify its velocity using x,y coordinates, you'd have to give two velocities - which is much less intuitive than saying, "it's spinning at such-and-such a speed". And anyway the velocities in the x,y directions are constantly changing and depend on distance from the centre of the record player, whereas the angular speed in the φ direction remains constant everywhere.


"The cat is spinning at 20 rpm (or 120 degrees per second)", vs, "the cat's x velocity is 1 m/s and its y velocity is 0.5 m/s, no wait now it's 0.6 m/s and 0.2 m/s, no wait it's changing again, aaaaaargh !"
Or to illustrate this slightly more scientifically :



When the green point is at the top or bottom of the circle, it has no velocity in the y-direction at all. Similarly, when it's at the extreme left or right, it has no velocity in the x-direction. But it always has a constant angular velocity.

So polar coordinates are much more useful for describing rotating discs. The problem is that of course for rendering images, pixels generally aren't in polar coordinates : we have to convert back to Cartesian. That's a problem when visualising simulations : just like trying to map the spherical Earth with a flat surface, you can get horrible distortions or lose detail if you're not careful.

Converting between polar and Cartesian coordinates is literally like trying to square the circle.
My first approach was to convert the data into the regular Cartesian system : knowing the r,θ,φ coordinates directly from the data, it was easy to convert to x,y,z. Which gave me this :

It's simulation of a protoplanetary disc.

... at which I make a brief interjection because at about that moment I received the following email, which I will keep anonymous because I'm not a total douchebag :
...and I do my first steps in scientific visualization. I am very interested in astronomy, though my main job till now was connected only with graphic design in university sector. I have some 3D modelling experience (several years ago me and my colleague made a fulldome video). Now I learn Blender and try to write my first scripts in Python. Slowly I become the idea, how everything works. Although the more I read, the more question I get...But it is normal I guess :)
Hopefully soon the quantity of my knowledges will transfer to their quality. About a month ago I got a chance to apply for a position as a specialist for astronomical visualization. My interview was quite successful and now we got a test task, which will probably define the proper candidate. I have already an idea, how to solve it. But it would be nice to find someone, who understand the materia, could evaluate my job and give me some practical advices. So I would like to ask you, if you had time and wish to answer some of my questions.
OK... one of the other candidates is asking me for help without even realising that I'm applying for the same job.


I decided the only safe course of action was to make absolutely no response whatsoever, so that's what I did.

Anyway the first result was not bad, but not great. The problem is that when you convert between coordinate systems there's no guarantee the Cartesian data set will be completely filled - especially at large radii. The polar coordinates tell you the centres of each data cell (pixel) and the density (or temperature or whatever) in that entire cell. But the centre of that cell only corresponds to the position of one particular pixel in Cartesian coordinates - it's not the same as checking every pixel in Cartesian coordinates and finding which polar cell they're in. The upshot is that you end up losing detail in the centre (where the polar cells are closer together than the Cartesian cells) and large blank areas at the edges (where the polar cells are further apart than the Cartesian cells).

To illustrate that, let's take another look at the comparison between polar and Cartesian grids. First in the very centre :

Every large square of the Cartesian grid contains multiple points from the polar grid. So multiple polar cells get reduced to a few Cartesian cells - detail is lost.
And now in the outskirts, shading every Cartesian cell that's intersected by at least one polar cell corner point :

Oh noes ! Not every Cartesian cell is filled ! And this only gets worse at larger radii.
That's not a hopeless problem. One nice feature when converting is that you're free to choose how many Cartesian pixels you want very easily, so you can optimise for a balance of detail in the centre vs. empty regions on the outskirts. In principle, you could then fill in the blanks based on the nearest pixel, or accurately determine the value for each Cartesian cell by working out which polar cell it corresponds to. Doable, but not easy - and certainly not doable in the space of an afternoon, which was the stated scope of the exercise.

There's a more fundamental problem : to you show all the detail in the central regions, you'll need a lot more cells in Cartesian coordinates than if you used polar. Large data sets can easily run into hundreds of millions of cells, which means hundreds of millions of pixels : ouch ! Wouldn't it be better if we could somehow have non-square pixels ? Then we could show the data in its original polar form, with no loss of detail and no need to have a single pixel more than we really needed.

It turns out that we can do just that in Blender. Consider a slice right through the centre of the protoplanetary discs. If we pretend that r and φ are really y and x, we get this :


Of course they're not really x and y at all, which means we're looking at something that's weirdly distorted. But we can correct for this. The method I came up with was to assign each pixel to a face in Blender (UV mapping). Then we can move each vertex (i.e. distort each face) to put it back where it would have been in polar coordinates.


Bingo - we can have our cake (original spherical polar coordinates) and eat it too (no need to convert to square pixels). Of course, the real data isn't just one slice - it's lots of slices, each of a constant angle θ. So what we have is a series of cones :



And if we show all the cones, we get this - which is a pretty convincing way to fake a volumetric render, with the gaps between the cones only becoming apparent at certain angles :



In the above, density controls bother temperature and opacity (transparency). But it doesn't have to be density. It could be, for example, this mysterious Q parameter which is apparently heat transport that I know nothing about except that it looks nice :




In principle, we could fill those gaps by switching to spheres when the viewing angle is through the cones. Actually, I started with spheres because I'd already tried this* - I only had the idea to use cones during a long and boring meeting. They look nice enough on their own, though to really get things perfect we'd need to combine the two.

* Spheres are much easier to do because there's no need for UV mapping - Blender can calculate the distortion to a sphere itself just fine, but not cones.



So having figured out this somewhat complicated process, some considerable time passed before I heard anything back. It felt like forever.


Eventually at the beginning of March - a full five months after the application deadline, I got the news that... I was invited to an on-site interview ! Which led me to an odd mix of glee and frustration.


Many wisecracks did ensue, of course. Maybe, said colleagues, they just wanted more data visualised as a sort of way of getting cheap labour. Maybe they hadn't rejected anyone from the first two rounds. Maybe the number of candidates was actually increasing at each stage.

Nonetheless, I want along to said interview about three weeks later and went at it hell for leather. I brought along both glass cubes, 3D glasses, a copy of the Discover magazine in which I nuked a potato, and I even organized most of my Blender files and scripts from the last 14 years. It would be a five year position with the possibility of a permanent contract at the end, with a salary far more... European than those of the Czech Republic. Worth fighting for, even if the "as soon as possible" phrase had long since rung hollow.

Getting to Heidelberg involved a 7 hour bus trip. It was a very nice bus and I had the whole lower compartment to myself, which was nice. With wi-fi. Heck, it was better than most British trains by a considerable margin. There wasn't much to look at on the trip, but I've always thought that it's far easier to make my own entertainment than it is to make my own legroom.


The only real event that happened was that there was a perilously short connection between the bus to Mannheim and the train to Heidelberg, so I ended up jumping on a plausible-looking train (German trains do not indicate the destination and route very clearly) more out of hope than expectation. Which resulted in a fairly tense 20 minutes until the train stopped at the correct destination. After hauling my really quite surprisingly heavy bags to my hotel, I had little enough time in the evening to do anything except a short walk. I didn't get to see any of the pretty parts of Heidelberg.

Of course I did see the Haus der Astronomie the next morning. My interview went as well as it could have gone. With hindsight I wouldn't have changed a damn thing. There wasn't time to show everything, but I left fully satisfied that I'd done as much as I could possibly have done. I answered all their questions. I showed them the extremely heavy data cubes and 3D movies. I was a little peeved that - bizarrely - they really did want someone to start extremely soon, but I'd have been more than prepared to take the job anyway (even though I'd really, really, really like an extended period back in the UK). So off I went back to Prague while they interviewed the other candidate.

Alas the return trip was not without incident. About one hour in to the seven hour trip, the bus ground to a halt due to an accident up ahead. It didn't move an inch for the next four hours. The bus driver let everyone off to walk around (and even some people on to use the toilet). I talked to some otherwise politically like-minded Germans who were, it must be said, none too fond of the Czechs, but even the nicest bus becomes wearying after 11 hours. I eventually crawled back into my room - still dragging my laptop and heavy glass cubes, of course, at about 2:30am.



Then I proceeded to play the waiting game. Again. For another three weeks. Until finally :
Dear Rhys,
thank you for your patience. It was a difficult decision, but in the end we offered the job to the other remaining candidate. This is in no way a reflection on the skills you demonstrated - we were very much impressed by your visualizations, found that you communicated well with the astronomers who would have been your colleagues here, and you would have been a great addition to our team. In the end, it came down to experience - the candidate to whom we offered, and who has now accepted, the job, is older than you and had put those additional years to good use in the visualization field.
We wish you all the best for your future career - you have an unusual and interesting mix of skills, and I hope you will find a good place to put them to the best possible use, either in astronomy or beyond!
It's a very nice rejection letter, but a rejection all the same. All of that effort had been for absolutely nothing except a free but incredibly exhausting 24 hours in Heidelberg. Was it worth it ? I'll let John Cleese answer that. Skip to 46:28 if it doesn't automatically.


On the one hand, the successful applicant is about 10 years older than me, done four-dimensional relativistic raytracing calculations, and has written a freakin' book. Fair enough, I'd have hired him instead of me. On the other hand, it doesn't take six months to decide who's more experienced. You can do that from the start. Especially if you use the words, "as soon as possible" in the advertisement.

Oh well. As Captain Picard once eloquently put it, "shit happens". In another reality, alternate me is taking a short break in Cardiff before preparing to move to yet another country despite never wanting to leave home in the first place. Actual me is now in the more mundane process of searching for a flat before he gets kicked out of the institute's accommodation. Aaaargh.

Tuesday, 16 February 2016

How To Build A Galaxy

Take one part gas, ten parts stars, a hundred parts dark matter, mix together and simmer lightly for a billion years. Serve warm with a booming godlike voice.

Alas, no-one has any idea what the real recipe for galaxies is. Actually that's not true : lots of people have lots of different ideas, but hardly anyone agrees with anyone else. We're reasonably sure of the ingredients, but as to the cookery method we're definitely still at the pre-school level of making mud pies.


Ordinary spiral galaxies contain lots of stars, a fair bit of gas, some dust, and probably lots and lots of mysterious dark matter. But trying to work out the best way to cook all these together and get a galaxy isn't easy. Or rather, making a spiral galaxy is remarkably easy, but figuring out the method nature actually uses is much more difficult.

Sometimes it's best to avoid trying to figure out where galaxies come from and concentrate on what they're doing. After all, we're a lot more sure of what galaxies are than how they're made. For example, we know of several dozen clouds of hydrogen gas without any stars, seemingly just floating freely in space without any nearby galaxies. Most people wave their hands and say, "tidal debris !" by which they mean the gas must have been ripped out galaxies somehow, even if its properties don't match expectations at all.


A few people think that maybe these clouds are themselves "dark galaxies" - clouds of gas held together by dark matter that have never formed stars. I'll discuss the arguments for and against this in a future post. For now, it's enough to know that I want to test one simple hypothesis : could a dark galaxy survive in the chaotic environment of a cluster of galaxies ?

To do this, of course I need a dark galaxy to start with. That means I need to tell the computer to simulate a rotating disc of gas inside a cloud of dark matter. Doesn't sound so hard, does it ?

Building a galaxy turned out to be a fascinating exercise. With massive naivety, I thought it would be better to start with a "simple" uniform density disc of gas and nothing else. It wasn't. Now, I expected that because gas is rather more complicated than stars, it wouldn't be perfectly stable. I thought it would probably develop a few clumps and spiral arms and whatnot. But I didn't expect it to insist on turning into a ring and then explode.

Particle trajectories. You can see the original disc as the faint circle in the centre.
It turned out that I'd made a very rudimentary error here. There is a well-known and very simple equation that gives the velocity for a particle to maintain a circular orbit. In fact, it's so simple I'm going to show it here :
M is the mass of the object, G is the gravitational constant, and r is how far away you are from the centre. You can use this to calculate the speed of satellites around the Earth or planets around the Sun no problem at all.

But I had forgotten a very important assumption behind this equation. Here, M is the mass within the radius r - that is, all the mass that's below you. The assumption is that if you're inside a sphere, all the mass beneath you (closer to the centre) can be treated as though it was an infinitely small point at the centre, whereas all the mass above you can be ignored. It turns out that this is perfectly fine for spherical masses, but it is not true in general. I assumed that since a uniform density disc is just a thin slice of uniform density sphere, it would still work. But no matter what parameters I tried, I always ended up with a ring.


The problem is the assumption that you can treat a disc as being like a slice of a sphere. It turns out that really doesn't work.

Imagine that you're floating around inside a giant shell with the mass of a planet, and you look around with a telescope. Or, to better illustrate it, you shine a torch at the surface. Actually make that two torches, pointing in opposite directions.


Clearly the area illuminated by the more distant region is bigger than the the part of the shell illuminated that's nearer to you. The two flashlights have the same beam angle, it's just distance that makes the difference. It turns out that the gravitational force pulling you towards the more distant but more massive area is exactly cancelled out by the smaller nearer one. Because it's a sphere, it doesn't matter which way you look - the two opposite areas in any direction will always cancel.

The reason for this, essentially, is that the force of gravity depends on the distance squared and the mass, but the mass itself also depends on distance squared. If the more distant area is twice as far away, then the force it exerts will be four times smaller at the same mass... but it's also actually four times more massive. So it exerts exactly the same force as the nearest area ! Everything is in perfect balance.


Now, if you nest a bunch of shells inside one other, you get a filled sphere. By the shell theorem, you won't feel the force of any shells above you - only below. This is ground into us at undergraduate level so much that I completely forgot to check whether this assumption still holds if you're in a disc or a ring.

It doesn't.

The geometry of a ring is fundamentally different to that of a sphere. If you're inside a hoop and illuminate the two opposite sides with your torches, the mass of each segment now depends only on their distance. Not the distance squared, just the distance. But gravity still depends on distance squared... so the two opposing forces do not cancel. If you're floating inside a giant hoop, you'll be pulled towards the nearest part.

If the segment on the far side is twice as far away, the force it exerts at any given mass is still four times smaller... but the actual mass will only be twice as large. Hence the nearby smaller segment exerts more force. Rings are much more similar to Mr Burn's unstable equilibrium state than spheres.
Or imagine an even simpler example : two identical planets floating in space. Unless you're exactly midway between the two, you'll be pulled towards the closer one. In short, the gravitational field depends very strongly on geometry - you can't use the classical sqrt(GM/r) equation to determine the velocity in every case.

Anyway, my first reaction on finding out that my galaxy had exploded was to call the author of the code, my good friend Rory Smith (who we've previously met destroying a galaxy, albeit last time deliberately) and ask him why the hell it wasn't working. His reaction was not encouraging.
Where to start. Setting up disks is a nightmare ! I spent about 6 months of my PhD on it and they still weren't great. And its not just me - DICE is a very clever parallelised code used for setting up galaxies. It runs on 8 cores, chugging away for half an hour,  carefully measuring the potential gradients everywhere on logarithmic grids. Then you run the initial conditions and the disk collapses, bars, forms rings, etc, etc!
Rory also sent me a link to a paper where the authors had calculated how to set the rotational velocity for a uniform disc, just like what I was trying to do. It was even more dispiriting than Rory's own response. I mean just look at it. LOOK AT IT !!!




It doesn't even end with a nice simple equation I can use. No way am I prepared to slog through that many equations, so instead I came up with another approach.


The simulation code can directly measure the accelerations every single particle experiences. They aren't derived from some clever formula, they're directly measured. Fortunately, the equation for circular motion is very simple (v^2 = a*r) once you know the acceleration a in the radial direction, which of course we do. So by feeding in the particle disc into the simulation we can measure the initial accelerations and use them to determine the velocity needed for circular motion. A little bit of vector analysis is needed to make sure we're measuring the acceleration in the right direction, but it's not too difficult.

And so I did that... and lo ! It didn't work. But at least now it didn't work in a slightly different way.


This time the disc simply fragments and breaks apart without any large-scale motion turning it into a ring. The problem is that the gas is cold. Because the particles are distributed at random, by chance there are some regions which are denser than others and tend to collapse. In reality these would form stars, but that won't help with my nice "simple" pure gas disc. Another solution is to make the gas hotter, so that its random motions will smooth out any condensations that start to form. Trouble is we don't want to make the gas so hot that the thing just turns into a sphere or completely evaporates.

After a great deal of trial and error tweaking, the best I could come up with was this :


Any colder than this and the disc becomes a few large fragments which disperse because the system is no longer symmetrical. Any hotter and the whole thing evaporates. As it is there's just one large fragment in the center - so it remains a disc, just not a uniform density one. I concluded (along with some dodgy mathematics) than uniform density discs are probably impossible.

Fortunately, there is a solution. As we examined last time with Rory, observations show that galaxies are spinning much faster than expected. This implies there's a lot of extra mass holding them together - without "dark matter", they'd just fly apart. Adding dark matter just makes everything better. After a few failed attempts where I made some silly mistakes, pretty soon I had this gobsmackingly stable disc :


OK, sure, it's not perfectly stable. It gets a bit denser in the centre, but not very much. Some rings appear but they more or less disappear by the end. And it definitely does not explode or fragment.

The beautiful thing about dark matter is that it provides lots of extra stability to the disc. The extra gravity means the gas can be much hotter without it evaporating, which means it's much less prone to fragmentation. It also changes the rotation velocity curve (that is, how fast it's rotating at different distances from the centre), making it much more stable against differential rotation.

Without dark matter I had to very very carefully tweak parameters to get something that even vaguely works; with dark matter I found I could set things almost however I wanted and the thing would be incredibly well behaved. I could even send it hurtling through a galaxy cluster and nothing much happened. Unlike Mr Burns, who is in terrible danger from even a light breeze, it really is practically indestructible.

Dark matter shown in pink.

But what about that paper with the scary maths ?

The problem here is that despite the mathematical horrors, it's actually an over-simplification. It considers only velocity and gravity - the thermal motions and fluid properties of the gas are ignored. And as we've seen, these can be extremely important. It was also written by non-astronomers who concluded that because the gravitational field of a disc is different to a sphere, there's no need for dark matter. Which, as I've expounded many times previously, is just not true.

And yet... when I measured the rotation velocities needed to support the uniform pure gas disc, I found something very close to what the authors said it should be. That this didn't actually make a stable disc didn't stop me from being intrigued.



The black line is what you'd expect if you could use the shell theorem and ignore the mass beyond any given radius. The blue line is what my code produced based on the actual accelerations experienced by the particles - there's some scatter in it because there are only a finite number of particles in the simulation (which is another complication that I won't go in to). The green line is the mathematical result from the paper, linearly stretched to fit the blue line as best as possible. It works pretty well except at the very end, where the mathematical solution goes to infinity which is obviously nonsense.

This extremely good agreement perked my interest, and reassured me that my code was working correctly. But could it really be that dark matter was an illusion based on a faulty assumption ? The authors claim that a nice flat rotation curve arises if you use a linearly-decreasing density profile instead of a uniform density as above. I wrote a code that could distribute particles according to any input density profile, so I tested the linear density profile too.


Pretty darn good ! This time the red points are measuring the curve for the gas particles and the blue points use massless "test particles" placed beyond the edge of the disc, just to see how the curve would extend. The green is from the paper again - it looks rather flatter in their paper only because they use a different scaling. It does not, however, resemble the flat (or slightly rising) rotation curve from observations. Like this one for M33, from Wikipedia :


But there's worse. The author's assumed that the density of the disc decreases linearly with radius. This is simply not the case at all for real galaxies, which usually have an exponentially decreasing density profile. So I constructed a rotation curve using a realistic density distribution for M33. Very realistic, because I simply took it directly from observational data. I used the Arecibo data I've described previously for the gas and Digital Sky Survey data for the optical. Google searching for papers gave me reasonable values for the total mass in stars. It was quick and dirty, but far superior to assuming a uniform or linearly-decreasing density profile.

The gas density is uniform in the centre then exponential beyond
the stellar disc, while the stellar disc is always exponential. Hence
the combination produces something a bit weird.
Oh. Well that's not very much like the observations at all, is it ? No, no it isn't. An exponential density profile means there's much less material in the outskirts than in the uniform or linear case. So much less, in fact, that the classical v = sqrt(GM/r) becomes a good approximation after all. Why the authors felt they should use a linear density profile I've no idea, because even a Google image search for "galaxy density profile" reveals that their profiles aren't linear at all. That's more supporting evidence for my theory that mathematical abilities are no guarantee of intelligence, or even simple knowledge in this case.

And it gets worse. In all of the above cases, including the uniform and linear densities, I assumed a total mass equivalent to the combined observable gas and stars. Have a look at the rotation velocity values : maxima of 50-70 km/s. The real M33 rotation velocity exceeds 100 km/s even at the highest distance to which we can measure the rotation curve. Unless you want to modify gravity or invoke another force entirely, there's no way to avoid the need for extra mass. You just can't avoid dark matter that easily.


In summary, there's a heck of a lot of physics in just a rotating disc of gas. It's an under-appreciated point that the gravitational field of a disc is very different to a sphere... but accounting for that isn't enough to keep it stable. The distribution of the material matters, but so too does its temperature and rotation. Ironically, it really does seem to need such a large amount of mass that it brings the rotation back to what you'd expect if it was a sphere after all - even if the disc is uniform, when its velocity profile should be very different if you ignore all the other factors at work.

If I'd just started with a disc inside a dark matter cloud, I'd have been completely unaware just how complicated the situation really is. So the old cliché is true : you learn more from mistakes than successes. That's a nice warm fuzzy feeling to end on, so let's throw in a dose of cynicism to keep it real.