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) ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.