Biplots for Functional Principal Component Scores

I'm trying to get the biplots between two F. principal components (or harmonics). I provide an example from fda package doc. to solve the riddle:

library(fda)
#BASIS FUNCTIONS
daybasis65 <- create.fourier.basis(c(0, 365), nbasis=65, period=365)
#DATA
harmaccelLfd <- vec2Lfd(c(0,(2*pi/365)^2,0), c(0, 365))
#SMOOTHING BASIS TO THE DATA
harmfdPar <- fdPar(daybasis65, harmaccelLfd, lambda=1e5)
daytempfd <- smooth.basis(day.5, CanadianWeather$dailyAv[,,"Temperature.C"],
daybasis65, fdnames=list("Day", "Station", "Deg C"))$fd
#FPCA
daytemppcaobj <- pca.fd(daytempfd, nharm=4, harmfdPar)
#RETRIEVING COMPONENTS
f_loadings <- daytemppcaobj $harmonics
#BIPLOT BETWEEN 1ST AND 2ND HARMONIC
biplot(f_loadings[1],f_loadings[2])
+Error in 1L:n : argument of length 0

Please take into account that f_loadings and daytempfd are fd objects while daytemppcaobj is a pca.fd object.

Thank you in advance. All type of help is welcome and rewarded!

Please if you have any idea to implement this with ggplot2, please, don't hesitate to write an answer.

For this case biplot is not the right function as it does not admit fd and pcafd objects. FDA package has plotscores which allows to BIplot Principal Component Scores or in other words, draw the scores of the curves on pairs of principal components

plotscores(daytemppcaobj, scores=c(1, 2))