OMaWa: OMaWa

estimateF <- function (T, RH){
# function to estimate F in mortality model
  myNorm <- function (mu, sigma, x){
    (1/(sigma * sqrt(pi*2)))*exp(-.5 *((x-mu)/sigma)^2)
  }
  myRh <- 40:100

  aPar <- c(-10.3380054, -10.4669102, -10.9059755, -11.0231063, -11.4011530, -11.8274586,
  -12.1645400, -12.5839053, -12.9565618, -13.4824782, -14.0807376, -14.6205262,
  -15.0394593, -15.5340079, -16.1052630, -16.5788810, -17.3262633, -18.0527660,
  -18.9383688, -19.8676464, -20.7073485, -21.1171309, -21.5736968, -21.6617334,
  -21.7430835, -21.9869836, -22.3929099, -22.6261775, -23.0285699, -23.0525138,
  -23.2420316, -23.4281215, -23.7672734, -24.3003487, -24.4483432, -24.7585358,
  -24.8734425, -25.6732439, -25.7538884, -26.3994288, -26.4420399, -32.0418536,
  -41.1927391, -73.9114453, -12.9157876, -4.2970298, -1.1196798, 0.5213484,
  1.5160074, 2.1778811, 2.6458577, 2.9908043, 3.2526227, 3.4553026,
  3.6143840, 3.7401128, 3.8398532, 3.9187920, 3.9811222, 4.0300582,
  4.0681841)

  bPar <-  c(1581.8644, 1595.2336, 1657.5313, 1667.9590, 1718.5610, 1777.2548,
  1821.0052, 1878.5489, 1926.2877, 1997.6355, 2081.8548, 2153.7098,
  2207.1146, 2276.2420, 2350.1879, 2408.1392, 2513.5856, 2609.7738,
  2734.2640, 2869.6431, 2982.2564, 3056.7456, 3138.6226, 3158.3735,
  3176.9869, 3223.2298, 3297.5007, 3342.1637, 3416.1198, 3424.8489,
  3462.2053, 3498.9484, 3562.6471, 3660.2774, 3690.5694, 3749.4758,
  3773.9340, 3921.7207, 3940.0727, 4059.8509, 4071.2715, 4944.5390,
  6512.4102, 13473.9898, 1697.1185, 749.0488, 469.7643, 343.3771,
  273.3633, 229.7157, 200.3210, 179.4352, 164.0144, 152.3185,
  143.2735, 136.1977, 130.6213, 126.2225, 122.7512, 120.0199,
  117.8809)

  cPar <- c(40.669287, 40.504866, 40.809235, 40.602486, 40.731921, 40.947283, 41.007085,
  41.195142, 41.240273, 41.447332, 41.776378, 41.924945, 41.976609, 42.211424,
  42.308235, 42.298512, 42.652162, 42.809360, 43.164268, 43.601604, 43.775461,
  44.120354, 44.474893, 44.519069, 44.556174, 44.728993, 45.036354, 45.191420,
  45.481229, 45.466308, 45.584791, 45.695318, 45.932225, 46.305757, 46.379981,
  46.578023, 46.624498, 47.216059, 47.232337, 47.660100, 47.645543, 49.478153,
  52.790242, 65.273419, 31.361257, 22.865273, 18.747942, 16.181704, 14.373231,
  13.001804, 11.909914, 11.010555, 10.251701, 9.600868, 9.036236, 8.543245,
  8.110852, 7.730967, 7.396142, 7.100071, 6.837113)

  whRH <- which.min(abs(myRh-RH))
    out <- floor(aPar[whRH]+ bPar[whRH]*myNorm(15, cPar[whRH] ,T))
    if (any(T >= 40)){
      out[T>=40] <- 2
    }
    if (any(T < 5)){
      out[T<5] <- 2
    }
    if (any(RH < 40)){
      out[RH < 40] <- 2
    }
    out[out<2] <- 2
  out
}