Celestia

1
-- Title: Show Redshifts of Galaxies
2

3
function get_distance(obj_pos, ref_pos)
4
   -- Returns distance Earth-object in Mpc.
5
   distance = ref_pos:distanceto(obj_pos) * km2Mpc
6
   return distance
7
end
8

9
function get_z(distance)
10
   -- Returns redshift z.
11
   -- Hubble constant = 73.2 (km/s)/Mpc, WMAP 3years, best
12
   -- http://map.gsfc.nasa.gov/m_mm/tp_links.html
13
   local H0 = 73.2
14
   local c = 299792.458
15
   z0 = (H0 * distance)/c
16
   z = z0*(1 - z0/4)/(1 - z0/2)^2
17
   -- d <= 2 * c/H0
18
   return z
19
end
20

21
function rgb2hex(r,g,b)
22
   -- Converts color code from RGB to Hex
23
   local hex = string.format("#%.2x%.2x%.2x", math.floor(r), math.floor(g), math.floor(b))
24
   return hex
25
end
26

27
function get_color(z_rel)
28
   -- Returns Hex color code from z_rel
29
   if z_rel > 1 then
30
       z_rel = 1
31
   end
32
   --local green = 255 * (1 - z_rel)/(1 + math.sqrt(z_rel))
33
   local green = 255 * (1 - z_rel)^2
34
   local red   = 255 
35
   local blue  = 255* math.sqrt(green/255)
36
   hex_color   = rgb2hex(red, green, blue)
37
   return hex_color
38
end
39
   
40
function mark_galaxies(sel_pos)
41
   zz = z_max(sel_pos)
42
   for dso in celestia:dsos() do
43
      if dso:type() == "galaxy" then
44
         dso_pos = dso:getposition()
45
         local d = get_distance(dso_pos,sel_pos)
46
         local z_rel = get_z(d)/zz
47
         local hex_color = get_color(z_rel)
48
         dso:mark( hex_color, "disk", 1, 1 )
49
      end
50
   end
51
end
52

53
function z_max(ref_pos)
54
   -- determine maximal redshift in catalog wrto ref_pos
55
   z_old = 0
56
   for dso in celestia:dsos() do
57
      if dso:type() == "galaxy" then
58
         dso_pos = dso:getposition()
59
         local d = get_distance(dso_pos, ref_pos)
60
         local z = get_z(d)
61
         if z > z_old then
62
             z_max = z
63
             z_old = z_max
64
             dsomax = dso
65
         end    
66
      end
67
   end
68
   return z_max
69
end
70

71
----------
72
-- main -- 
73
----------
74
celestia:unmarkall()
75

76
celestia:show("markers")
77
km2Mpc = 1/3.08568025e19
78
MW = celestia:find("Milky Way")
79
MW_pos = MW:getposition()
80
--
81
-- select and specially mark Milky Way
82
--
83
celestia:select(MW)
84
celestia:mark(MW)
85
--
86
-- color encode all other galaxies according to their redshift
87
-- relative to Milky Way
88
--
89
mark_galaxies(MW_pos)
90
--
91
-- move observer to a distance of 1000 Mpc from Milky Way
92
--
93
observer = celestia:getobserver()
94
observer:gotodistance(MW, 1000/km2Mpc,5)
95

96
sel_old = MW
97
   
98
while true do   
99
   sel = celestia:getselection()
100
   --
101
   -- specially mark possible new selection, unmark the old one
102
   --
103
   if sel ~= sel_old then
104
       sel:unmark()
105
       sel_old:unmark()
106
       sel:mark("#00FF00","disk",10,1)
107
       sel_old = sel
108
   end    
109
   
110
   sel_pos = sel:getposition()
111
   
112
   -- obs_pos = observer:getposition()
113
   
114
   if sel:type() == "galaxy" then
115
      local d = get_distance(MW_pos,sel_pos)
116
      local z = get_z(d)
117
      celestia:print(sel:name()..": "..string.format("redshift z = %5.3f,  distance = %5.2f  Mpc", z,d).."\nmax. redshift: "..dsomax:name(),5,-1,-1,0,6)
118
   end
119
   wait(0)
120
end
121