tIs a bit neater to create list objects for each of star1 and star2 and then, since the various functions are not used outside of gravitation stick them right in gravitation. Also we can use distance for len as well by giving it a default second arg of 0 and sqr does not seem to do much so we can eliminate that too.
Making the vec computation a function eliminates some redundant code and using outer eliminates the two loops. star1 <- list(mass = 30, center = c(-0.5, 0)) star2 <- list(mass = 1, center = c(0.5, 0)) gravitation2 <- function(v1, v2) { distance <- function(a, b = 0) sqrt(sum((a-b)^2)) norm <- function(x) x/distance(x) vec <- function(x) with(x, mass * norm(center - v) / distance(v, center)^2) v <- c(v1, v2) distance(vec(star1) + vec(star2)) } xvals2 <- yvals2 <- seq(vmin, vmax, step) b2 <- outer(xvals, yvals, Vectorize(gravitation2)) On Mon, Oct 20, 2008 at 1:41 PM, Johannes Bauer <[EMAIL PROTECTED]> wrote: > 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. > ______________________________________________ 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.