Play 3d graph with level curves in R software

Question

I tried to add level curves in my Plot as in the second figure below, but I could not, someone knows how I should proceed?

 #Function density probability
    library(pbivnorm)
    bsb <- function(t1,t2){
    a1 <- sqrt(phi1/2)*(sqrt(((phi1+1)*t1)/(phi1*mu1))-sqrt(((phi1*mu1)/((phi1+1)*t1))))
      a2 <- sqrt(phi2/2)*(sqrt(((phi2+1)*t2)/(phi2*mu2))-sqrt(((phi2*mu2)/((phi2+1)*t2))))
      Phi2 <- pbivnorm(a1, a2, rho, recycle = TRUE)
      b1 <- ((phi1+1)/(2*phi1*mu1))*sqrt(phi1/2)*(((phi1*mu1)/((phi1+1)*t1))^(1/2)+((phi1*mu1)/((phi1+1)*t1))^(3/2))
      b2 <- ((phi2+1)/(2*phi2*mu2))*sqrt(phi2/2)*(((phi2*mu2)/((phi2+1)*t2))^(1/2)+((phi2*mu2)/((phi2+1)*t2))^(3/2))
      fdp <- Phi2*b1*b2
      return(fdp)
    }
    t1 <- seq(0.001,5,length=100)
    t2 <- seq(0.001,5,length=100)
    #Parameters
    mu1=5
    phi1=2
    mu2=5
    phi2=2
    rho=0.9

    z<-outer(t1,t2,bsb) # calculate density values

    persp(t1, t2, z, # 3-D plot
          main="Bivariate Birnbaum-Saunders",
          col="lightgray",
          theta=40, phi=10,
          r=10,
          d=0.9,
          expand=0.5,
          ltheta=90, lphi=80,
          shade=0.9,
          ticktype="detailed",
          nticks=5) 

Answer 1

I managed to solve:

source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/panel.3d.contour.R")
library(lattice)
library(manipulate)
library(colorRamps)

#Function density probability
library(pbivnorm)
bsb <- function(t1,t2){
  a1 <- sqrt(phi1/2)*(sqrt(((phi1+1)*t1)/(phi1*mu1))-sqrt(((phi1*mu1)/((phi1+1)*t1))))
  a2 <- sqrt(phi2/2)*(sqrt(((phi2+1)*t2)/(phi2*mu2))-sqrt(((phi2*mu2)/((phi2+1)*t2))))
  Phi2 <- pbivnorm(a1, a2, rho, recycle = TRUE)
  b1 <- ((phi1+1)/(2*phi1*mu1))*sqrt(phi1/2)*(((phi1*mu1)/((phi1+1)*t1))^(1/2)+((phi1*mu1)/((phi1+1)*t1))^(3/2))
  b2 <- ((phi2+1)/(2*phi2*mu2))*sqrt(phi2/2)*(((phi2*mu2)/((phi2+1)*t2))^(1/2)+((phi2*mu2)/((phi2+1)*t2))^(3/2))
  fdp <- Phi2*b1*b2
  return(fdp)
}
#Parameters
mu1=5
phi1=2
mu2=5
phi2=2
rho=0.9
grid <- expand.grid(t1 = seq(0.001,8, by = 0.1),
                    t2 = seq(0.001,8, by = 0.1))
grid$z <- bsb(grid$t1,grid$t2)

manipulate({
  ## Makes the three-dimensional chart
  colr <- colorRampPalette(c(c1, c2, c3), space="rgb")
  arrows <- arr
  wireframe(z ~ t1 + t2,
            data = grid,
            scales = list(arrows = FALSE),
            zlim = extendrange(grid$z, f = 0.25),
            panel.3d.wireframe = "panel.3d.contour",
            nlevels = 8,
            col = "gray40",
            type = c("bottom"),
            col.regions = colr(101),
            drape = TRUE, colorkey=FALSE,
            screen=list(z=z.angle, x=x.angle),
            axis.line = list(col = "transparent"),
            clip = list(panel = "off"),
            par.settings = list(box.3d = list(col=c(1,NA,NA,1,1,NA,NA,NA,NA))))
},
## Controls the value of angles and colors
z.angle=slider(0, 360, step=10, initial=40),
x.angle=slider(-180, 0, step=5, initial=-80),
arr=checkbox(FALSE, "show.arrows"),
c1=picker("transparent","black","red","yellow","orange","green","blue","pink","violet"),
c2=picker("transparent","black","red","yellow","orange","green","blue","pink","violet"),
c3=picker("transparent","black","red","yellow","orange","green","blue","pink","violet")
)