title: How to use plot3d/surface3d (or another function?) to plot 4d function ("fourth dimension" denoted by color scale)?
body: <p>I would like to plot: 
  <code>production.ts(31, .002, 10,12,125313.93,211,95,x,"2014-02-01","2014-05-14",z,y)</code> as function of x,y,z</p>
  
  <p>As something like this plot from Mathematica, (if possible in R): 
  <a href="http://i.stack.imgur.com/3PRaf.png" rel="nofollow">http://i.stack.imgur.com/3PRaf.png</a></p>
  
  <p>I have a function:</p>
  
  <pre><code>library("lubridate"); library("rgl")
  
  production.ts &lt;- function(a, b, z, c, d, e, 
                      f, g, h, j, r, k) {
    elapsed &lt;- (4-z)*10 + (4-c)
    un.days &lt;- 100 - elapsed
    gone.days &lt;- day(as.Date(h))
    rem.days &lt;- day(as.Date(j))
    r.days &lt;- as.numeric(as.Date(j) - as.Date(h))
  
    m.r &lt;- f/100*d
    inputs &lt;- d * a * (gone.days - 1)/365 + r
    prin &lt;- m.r + inputs 
    costs &lt;- (r.days/365 * r + 1) * prin 
    added.p &lt;- a/100*d + r 
    due &lt;- d * 1-un.days
    tomr.f &lt;- 1- due + k^2
    acct.paid &lt;- (d - due)*tomr.f
    net &lt;- added.p + due + acct.paid
    pv.net &lt;- net/(1+r*(e-30-day(as.Date(j)))/365)
    end &lt;- d - due - acct.paid
    more.add.p &lt;- end*a*(rem.days-1)/365
    rem &lt;- (f-g)/100 * end
    total.fv &lt;- pv.net + rem + more.add.p 
    out &lt;- costs - total.fv
  
    out 
  }
  
  
  x&lt;-seq(-10,10,by=.1)
  y&lt;-seq(0,1000,by=.1)
  z&lt;-seq(0,90,by=.1)
  </code></pre>
  
  <p>I have tried: 
  <code>func.3d&lt;-Vectorize(production.ts(31, .002, 10,12,125313.93,211,95,x,"2014-02-01","2014-05-14",z,y))</code></p>
  
  <pre><code>c &lt;- func.3d; c &lt;- cut(c,breaks=64); cols &lt;- rainbow(64)[as.numeric(c)]
  
  open3d()
  plot3d(x, y, z, col=cols,type="s",size=1)
  </code></pre>
  
  <p>But this plots lines and the colors don't line up with the values the function should output. </p>
  
  <p>Does anyone know how I could do this? Thanks, I really appreciate your time!</p>