[R] R Newbie Question

[Johannes Bauer]

Hello list,

I just started R today and tried something quite simple. I wanted to
create a colored plot and eventually after hours of fiddling around got
it working. However, my solution seems very suboptimal and I'd really
appreciate your hints on how to improve. I believe that R already offers
many functions I coded (e.g. distance between two vectors, vector
length, vector normalization and so on). I generally didn't even figure
out how to create a simple vector or how to extract a row from a matrix
so the result is a vector (to get scalars, I use the "sum" functions,
which is an incredibly ugly workaround).

To sum the problem up: one has a binary star system with a large star
(e.g. red giant) and a small star (e.g. white dwarf). Gravitation
between them is directly proportional to the mass and indirectly
proportional to the square of the distance. If correctly plotted, one
should be able to see the inner lagrange point L1 which is the point
where the gravitational potentials of the stars "cancel out", e.g. an
object would not be attracted to any star. Well, enough background
information, here's my rookie code - please feel free to comment on
anything :-)

Kind regards,
Johannes



star1center = vector("numeric", 2)
star1center[1] = -0.5
star1center[2] = 0
star1mass = 30

star2center = vector("numeric", 2)
star2center[1] = 0.5
star2center[2] = 0
star2mass = 1

sqr = function(x) {
	return(x * x)
}

distance = function(a, b) {
	return(sqrt(sqr(a[1] - b[1]) + sqr(a[2] - b[2])))
}

len = function(x) {
	return(sqrt(sqr(x[1]) + sqr(x[2])))
}

norm = function(x) {
	return(x / len(x))
}

gravitation = function(invecx, invecy) {
	invec = vector("numeric", 2)
	invec[1] = invecx
	invec[2] = invecy
	vec1 = star1mass * norm(star1center - invec) / sqr(distance(invec,
star1center))
	vec2 = star2mass * norm(star2center - invec) / sqr(distance(invec,
star2center))
	return(len(vec1 + vec2))
}

vmin = -1
vmax = 1
step = 0.1
vals = ((vmax - vmin) / step) + 1

xvals = seq(vmin, vmax, step)
yvals = seq(vmin, vmax, step)

a = expand.grid(seq(vmin, vmax, step), seq(vmin, vmax, step))
b = matrix(seq(1, vals*vals), vals)

for (x in 1:vals) {
	for (y in 1:vals) {
		b[x, y] = gravitation(sum(a[x,][1]), sum(a[y,][1]))
	}
}
filled.contour(xvals, yvals, z = b, color = heat.colors, ylim = c(-1,
1), xlim = c(-1, 1), zlim = c(0, 100), nlevels = 100)