R/calcCirf.R

In TIMP: Fitting Separable Nonlinear Models in Spectroscopy and Microscopy

"calcCirf" <-
function (k, x, irfpar, mirf = FALSE, measured_irf = vector(), 
    convalg = 1, shiftmea = vector(), lamb = 1, reftau = 0, 
    doublegaus = FALSE, multiplegaus = FALSE,
    streak = FALSE, streakT = 0, irffun = "gaus") 
{
   
  if (!mirf) {
    if(irffun == "gaus" || irffun == "doublegaus") {
      mu <- irfpar[1]
      tau <- irfpar[2]
      m <- convGausExp(k, x, mu, tau)
      #m <- rep(0, length(x) * length(k))
      #m <- as.matrix(.C("r_calcCirf", m = as.double(m), as.double(k), 
      #                  as.double(x), as.double(tau), as.double(mu),
      #                  as.integer(length(k)), 
      #                  as.integer(length(x)), PACKAGE="TIMP")$m)
      #dim(m) <- c(length(x), length(k))
      ## for streak synchronscan data get backsweep term 
      if(streak) {
        backsweep <- getStreakBacksweep(streakT, k, mu,x)
        m <- m + backsweep
      }
            if(doublegaus || irffun == "doublegaus") {
              scal <- irfpar[4] #note in multiplegaus, scal comes first
              tau <- irfpar[3]
              m2 <- convGausExp(k, x, mu, tau)
              #m2 <- rep(0, length(x) * length(k))
              #m2 <- as.matrix(.C("r_calcCirf", m = as.double(m2),
              #                   as.double(k), as.double(x), as.double(tau2),
              #                   as.double(mu), as.integer(length(k)),
              #                   as.integer(length(x)), PACKAGE="TIMP")$m)
              #dim(m2) <- c(length(x), length(k))
              if(streak) m2 <- m2 + backsweep
              m <- m + (scal * m2)
            }
    }
    
    if(multiplegaus || irffun == "multiplegaus") {
      mu <- irfpar[1]
      tau <- irfpar[2]
      m <- convGausExp(k, x, mu, tau)
      ## for streak synchronscan data get backsweep term 
      if(streak) {
        backsweep <- getStreakBacksweep(streakT, k, mu,x)
        m <- m + backsweep
      }   
          if(length(irfpar)>4) {
              for( i in 1:((length(irfpar)-2)/3)) {
              # scal shift tau
              scal <- irfpar[1+(i-1)*3+2]
              shift <- irfpar[2+(i-1)*3+2]
              tau <- irfpar[3+(i-1)*3+2]
              m2 <- convGausExp(k, x, mu+shift, tau)              
              if(streak) m2 <- m2 + backsweep
              m <- m + (scal * m2)
              }
           }
            
    }
    
    
    
    if(irffun == "step") {
      m <- convolveUnitStep(k, x, irfpar) 
      
    }
  }
  else {
    m <- matrix(0, nrow = length(x), ncol = length(k))
    xspace <- x[2] - x[1]
    if (length(shiftmea) != 0) {
      if (length(shiftmea) == 1) 
        lamb <- 1
      measured_irf <- .C("r_ShiftCurve", source = as.double(measured_irf), 
                         as.double(measured_irf),
                         as.double(shiftmea[lamb]/xspace), 
                         as.integer(length(x)), PACKAGE="TIMP")$source
    }
    for (i in 1:length(k)) {
      m[, i] <- switch(convalg, .C("r_Conv1", result = as.double(m[, i]),
                                   as.double(measured_irf),
                                   as.integer(length(x)), 
                                   as.double(k[i]), as.double(xspace), 
                                   PACKAGE="TIMP")$result, 
                       .C("r_Conv2", result = as.double(m[, i]),
                          as.double(measured_irf), 
                          as.integer(length(x)), as.double(k[i]),
                          as.double(xspace), 
                          PACKAGE="TIMP")$result, 
                       .C("r_Conv3", result = as.double(m[, 
                                     i]), as.double(measured_irf),
                          as.integer(length(x)), 
                          as.double(k[i]), as.double(xspace),
                          as.double(reftau),  
                          PACKAGE="TIMP")$result)
      
        }
    }
    m
}