RP.SVJcall <- function(spot,strike,timetomat,R,V,kappa.V,theta.V,sigma.V,rho,mu.J,sigma.J,lambda, greek=1){

  i<-1i
  tau<-timetomat
  lnspot<-log(spot)
  lndisc<--R*tau
  
  integrand<-function(phi,j){

    iphi<-i*phi

    xi.V<-sqrt((kappa.V - rho*sigma.V*(2-j+iphi))^2-sigma.V^2*iphi*(iphi+3-2*j))
   Constterm<-(theta.V/sigma.V^2)*(2*log(1-((xi.V-kappa.V+rho*sigma.V*(2-j+iphi))*(1-exp(-tau*xi.V)))/(2*xi.V))+tau*(xi.V-kappa.V+rho*sigma.V*(2-j+iphi)))

   Jumpterm<-lambda*tau*(1+mu.J)^(2-j)*((1+mu.J)^(iphi)*exp(0.5*iphi*sigma.J^2*(iphi+3-2*j))-1)-lambda*iphi*mu.J*tau

   Vterm<-(iphi*(iphi+3-2*j)*(1-exp(-xi.V*tau)))/(2*xi.V-(xi.V-kappa.V+rho*sigma.V*(2-j+iphi))*(1-exp(-tau*xi.V)))
   
   fj<-exp(-Constterm+iphi*(lnspot-lndisc)+V*Vterm+Jumpterm)

   integrand<-Re(exp(-iphi*log(strike))*fj/(iphi))
  }

  P1<-0.5+(1/pi)*integrate(integrand,lower=0.000001,upper=10000,subdivisions=10000,j=1)$value
  P2<-0.5+(1/pi)*integrate(integrand,lower=0.000001,upper=10000,subdivisions=10000,j=2)$value
  
  if (greek==2) {RP.SVJcall<-P1}
  else 	{RP.SVJcall<-spot*P1-exp(-R*tau)*strike*P2}
}