dat<- data.frame(t=seq(0, 2*pi, by=0.1) ) xhrt <- function(t) 16*sin(t)^3 yhrt <- function(t) 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t) dat$y=yhrt(dat$t) dat$x=xhrt(dat$t) with(dat, polygon(x,y, col="hotpink"))
(R code inspired by/lifted from: DWin on StackOverflow)