' set of functions to calculate age-group etc. WRs
' Howard Grubb 1999-2000
' http://www.rdg.ac.uk/~snsgrubb/athletics/home.html
' see paper:
' Grubb, H.J. (1998) Models for comparing athletic performances.
' The Statistician, 47, 509-521.

' key ones are:
'   =PWRSpeed(dist,time,sex)        % (fraction) WR-speed
'   =PAWRSpeed(dist,time,age,sex)   % (fraction) age-group WR-speed
'   =PWWRSpeed(dist,time,age,sex)   % (fraction) WAVA age-standard
'   =dist_recs(dist,sex)            WR speed (m/s)
'   =dist_arecs(dist,age,sex)       age-group WR speed (m/s)
'   =dist_wrecs(dist,age,sex)       WAVA age-standard (m/s)
' Others:
'   =AgeWRPerf(dist,age,sex)        age-group WR performance (% WR speed)
'   =SAgeWRPerf(dist,age,sex)       as above, for swimming
'   =WAVAFactor(event,age,sex)      WAVA age-factor - all track and field events - event codes
'   =WAVAStandard(event,sex)        WAVA event standard - all track and field events
'   =IWAVAFactor(dist,age,sex)      WAVA age-factor - running only (interpolates distance)
'   =IWAVAStandard(dist,sex)        WAVA event standard - running only (interpolates distances)

' Note that the WAVA functions require the spreadsheet with the WAVA factors
' the other calculations are simple parametric functions, self-contained in the code

Option Explicit

Function AgeWRPerf(dist As Single, age As Integer, sex As Integer) As Single
' uses simple model fitted to age WRs to calculate age-group % WR speed
' for a given distance (km), age (years), sex=0 female, sex=1 male
' returns % (overall) WR-speed for that dist/age/sex world record

  Dim a As Single
  Dim b As Single
  Dim p As Single
  Dim apwr As Single
  
' below are parameters of function fitted to age-group WRs
  a = 1.0064954 ' gradient with age^2
  b = -0.000007905554
  p = 2.5 ' power of age
  If (sex = 0) Then ' female
    a = 1.0056745
    b = -0.0000096131749
    p = 2.5
  End If
  apwr = 1
  If (age > 90) Then age = 90
  If (age >= 35) Then apwr = a + b * (age ^ p)
  AgeWRPerf = apwr ' factor to scale down WR speed by for that age
End Function

Function SAgeWRPerf(dist As Single, age As Integer, sex As Integer) As Single
' as above, but for swimming world records
  Dim a As Single
  Dim b As Single
  Dim p As Single
  Dim apwr As Single
  
  a = 1.0189311
  b = -0.0000075960805
  p = 2.5
  If (sex = 0) Then ' female
    a = 1.0275257
    b = -0.0000088615036
    p = 2.5
  End If
  apwr = 1
  If (age > 90) Then age = 90
  If (age >= 35) Then apwr = a + b * (age ^ p)
  SAgeWRPerf = apwr
End Function

Function WAVAStandard(tevent As String, sex As Integer) As Single
' lookup WAVA standards for all events
' sex=0 female sex=1 male
 
  Dim std As Single
  Dim ifac As Integer
  Dim sexSheet As String
 
  sexSheet = "Men"
  If (sex = 0) Then sexSheet = "Women"
 
' find event
  ifac = -1 ' can we always find event??
  ifac = Application.Match(tevent, Worksheets(sexSheet).Range("Event"), 0)

  std = 0
  If (ifac > 0) Then std = Worksheets(sexSheet).Range("Standard")(ifac).value
 
  WAVAStandard = std
End Function

Function WAVAFactor(tevent As String, age As Integer, sex As Integer) As Single
' lookup WAVA factors for all events
' age in years, sex=0 female sex=1 male
 
  Dim fac As Single
  Dim fac2 As Single
  Dim rage As Single
  Dim rage1 As Single
  Dim prage As Single
  Dim iage As Integer
  Dim iage1 As Integer
  Dim ifac As Integer
  Dim nages As Integer
  Dim sexSheet As String
 
  sexSheet = "Men"
  If (sex = 0) Then sexSheet = "Women"

  nages = Worksheets(sexSheet).Range("Age").Count
' find age
  If (age < 8) Then ' ignore very young
    iage = 1
    iage1 = 1
  Else
    iage = -1
    iage = Application.Match(age, Worksheets(sexSheet).Range("Age"), 1) ' first age below?
    iage1 = iage
    If (Worksheets(sexSheet).Range("Age")(iage) <> age) Then ' interpolate
      If (iage = nages) Then
        iage1 = nages ' don't extrapolate
      Else
        iage1 = iage + 1
      End If
    End If
  End If

' find event
  ifac = -1 ' can we always find event??
  ifac = Application.Match(tevent, Worksheets(sexSheet).Range("Event"), 0)
  
  fac = 1 ' default
  If (iage >= 0) Then
    fac = Worksheets(sexSheet).Range("Factors")((ifac - 1) * nages + iage + 0).value
    If (iage1 <> iage) Then ' interpolate
      rage = Worksheets(sexSheet).Range("Age")(iage)
      rage1 = Worksheets(sexSheet).Range("Age")(iage1)
      prage = (age - rage) / (rage1 - rage)
      fac2 = Worksheets(sexSheet).Range("Factors")((ifac - 1) * nages + iage1 + 0).value
      fac = (1 - prage) * fac + prage * fac2
    End If
  End If
 
  WAVAFactor = fac
End Function

Function IWAVAStandard(dist As Single, sex As Integer) As Single
' interpolate WAVA standards for running events
' sex=0 female sex=1 male
 
  Dim std As Single
  Dim pfac As Single
  Dim ifac As Integer
  Dim ifac1 As Integer
  Dim i As Integer
  Dim nfac As Integer
  Dim Where As Range
  Dim sexSheet As String
 
  sexSheet = "Men"
  If (sex = 0) Then sexSheet = "Women"
 
' find event
  ifac = -1 ' can we always find event??
  i = 12 ' first event with dist(km)
  nfac = 42 ' last event with dist(km)
  Set Where = Worksheets(sexSheet).Range("Dist")
  Do While (i < nfac And Where(i).value < dist)
    i = i + 1
  Loop
  
  If (i = 12) Then ' less than shortest event
    ifac = i
    ifac1 = i + 1
  End If
  If (i >= nfac - 1) Then ' longer than longest event
    ifac1 = nfac - 1
    ifac = nfac - 2
  End If
  If (i > 12 And i < nfac - 1) Then ' within array
    ifac = i - 1
    ifac1 = i
  End If
 
  pfac = 0
  If (Where(ifac).value <> dist) Then '/ need interpolation/extrapolation
    pfac = (dist - Where(ifac).value) / (Where(ifac1).value - Where(ifac).value)
  End If
  
  std = 0
  If (ifac > 0) Then
      std = pfac * Worksheets(sexSheet).Range("Standard")(ifac1).value + _
          (1 - pfac) * Worksheets(sexSheet).Range("Standard")(ifac).value
  End If
 
  IWAVAStandard = std
End Function

Function IWAVAFactor(dist As Single, age As Integer, sex As Integer) As Single
' interpolate WAVA factors for running events
' dist in km, age in years, sex=0 female sex=1 male
 
  Dim fac As Single
  Dim fac2 As Single
  Dim pfac As Single
  Dim rage As Single
  Dim rage1 As Single
  Dim prage As Single
  Dim i As Integer
  Dim iage As Integer
  Dim iage1 As Integer
  Dim ifac As Integer
  Dim ifac1 As Integer
  Dim nfac As Integer
  Dim nages As Integer
  Dim sexSheet As String
  Dim Where As Range
 
  sexSheet = "Men"
  If (sex = 0) Then sexSheet = "Women"
 
  nages = Worksheets(sexSheet).Range("Age").Count
' find age
  If (age < 8) Then ' ignore very young
    iage = 1
    iage1 = 1
  Else
    iage = -1
    iage = Application.Match(age, Worksheets(sexSheet).Range("Age"), 1) ' first age below?
    iage1 = iage
    If (Worksheets(sexSheet).Range("Age")(iage) <> age) Then ' interpolate
      If (iage = nages) Then
        iage1 = nages ' don't extrapolate
      Else
        iage1 = iage + 1
      End If
    End If
  End If
  
' find event
  ifac = -1 ' can we always find event??
  i = 12 ' first event with dist(km)
  nfac = 42 ' last event with dist(km)
  Set Where = Worksheets(sexSheet).Range("Dist")
  Do While (i < nfac And Where(i).value < dist)
    i = i + 1
  Loop
  
  If (i = 12) Then ' less than shortest event
    ifac = i
    ifac1 = i + 1 ' extrapolate
  End If
  If (i >= nfac - 1) Then ' longer than longest event
    ifac1 = nfac - 1
    ifac = nfac - 2 ' extrapolate
  End If
  If (i > 12 And i < nfac - 1) Then ' within array
    ifac = i - 1 ' interpolate
    ifac1 = i
  End If
 
  pfac = 0
  If (Where(ifac).value <> dist) Then '/ need interpolation/extrapolation
    pfac = (dist - Where(ifac).value) / (Where(ifac1).value - Where(ifac).value)
  End If
'  MsgBox ("pfac: " & pfac & " ifac: " & ifac & " dist: " & Where(ifac).Value & " ifac1: " & ifac1 & " dist1: " & Where(ifac1).Value)
  
  fac = 1 ' default
  Set Where = Worksheets(sexSheet).Range("Factors")
  If (iage >= 0) Then ' interpolate distance
    fac = (1 - pfac) * Where((ifac - 1) * nages + iage + 0).value + _
                pfac * Where((ifac1 - 1) * nages + iage + 0)
'    MsgBox ("fac: " & fac)
    If (iage1 <> iage) Then ' interpolate age
      rage = Worksheets(sexSheet).Range("Age")(iage)
      rage1 = Worksheets(sexSheet).Range("Age")(iage1)
      prage = (age - rage) / (rage1 - rage)
'      MsgBox ("prage: " & prage)
      fac2 = (1 - pfac) * Where((ifac - 1) * nages + iage1 + 0).value + _
                   pfac * Where((ifac1 - 1) * nages + iage1 + 0)
'      MsgBox ("fac2: " & fac2)
      fac = (1 - prage) * fac + prage * fac2
    End If
  End If
            
  IWAVAFactor = fac
End Function

Function dist_arecs(dist As Single, age As Integer, sex As Integer) As Single
' given distance (km), age (years), sex=0 female, sex=1 male
' returns speed (m/s) for that dist/age/sex world record
  dist_arecs = dist_recs(dist, sex) * AgeWRPerf(dist, age, sex)
End Function
  
Function dist_wrecs(dist As Single, age As Integer, sex As Integer) As Single
' given distance (km), age (years), sex=0 female, sex=1 male
' returns speed (m/s) for that dist/age/sex WAVA factor
  dist_wrecs = dist * 1000 / IWAVAStandard(dist, sex) * IWAVAFactor(dist, age, sex)
End Function

Function interp_array(array1, array2, value) ' interpolate results

  Dim i As Integer
  Dim indx As Integer
  Dim indx1 As Integer
  Dim prop As Single
  
  indx = -1
  prop = 0
  i = LBound(array1) ' find the first event longer than the chosen distance
  Do While (i < UBound(array1) And array1(i) < value)
    i = i + 1
  Loop

  If (i = LBound(array1)) Then ' smaller than smallest
    indx = 0
    indx1 = 1
  End If
  If (i >= UBound(array1)) Then ' larger than largest
    indx1 = UBound(array1)
    indx = indx1 - 1
  End If
  If (i > LBound(array1) And i < UBound(array1)) Then ' within array
    indx = i - 1
    indx1 = i
  End If
  
  If (array1(indx) <> value) Then ' need interpolation/extrapolation
    prop = (Application.Ln(value) - Application.Ln(array1(indx))) / _
         (Application.Ln(array1(indx1)) - Application.Ln(array1(indx)))
'    prop=(value-array1(indx))/(array1(indx1)-array1(indx))
  End If

  interp_array = prop * array2(indx1) + (1 - prop) * array2(indx)
  
End Function
  
Function dist_recs(dist As Single, sex As Integer) As Single ' compute record for distance
' given distance (km), sex=0 female, sex=1 male
' returns speed (m/s) for that dist/sex world record (smooth model)
' Updated to use model for top 200 times at each distance (September 2000)

  Dim bestdists As Variant
  Dim mbestspeeds As Variant
  Dim wbestspeeds As Variant

  bestdists = Array(100, 200, 400, 800, 1500, 3000, 5000, 10000, 21100, 42200) ' record distances
  mbestspeeds = Array(10.2919, 10.2607, 9.27439, 7.93579, 7.29478, 6.80194, 6.57603, 6.25136, 5.94435, 5.61122)
  wbestspeeds = Array(9.52688, 9.43498, 8.4094, 7.11094, 6.51896, 6.07955, 5.82178, 5.5302, 5.29561, 5.02135)
  If (sex = 0) Then ' female
    dist_recs = interp_array(bestdists, wbestspeeds, (dist * 1000))
  Else
    dist_recs = interp_array(bestdists, mbestspeeds, (dist * 1000))
  End If
End Function
  
' old models
'  mod5 = 1
'
'  If (mod5 = 1) Then ' model (5) from paper
'    a = 10.43
'    b = 3.91
'    c = 4.19
'    If (sex = 0) Then ' female
'      a = 7.99
'      b = 4.21
'      c = 3.9
'    End If
'    recs = a / (Application.Ln(dist * 1000) - b) + c
'  Else  ' model (6) from paper
'    a = 17.29
'    b = 250
'    c = 4.69
'    l = -0.267
'    If (sex = 0) Then ' female
'      a = 16.15
'      c = 4.12
'    End If
'    recs = a * (dist * 1000 - b) ^ l + c
'  End If
'  dist_recs = recs ' WR speed (m/s)
'End Function

Function PWRSpeed(dist As Single, time As Single, sex As Integer) As Single
' dist (km), time(seconds), sex=0 female, sex=1 male
' return percentage (actually fraction) of WR speed
  PWRSpeed = (dist * 1000 / time) / dist_recs(dist, sex)
End Function

Function PAWRSpeed(dist As Single, time As Single, age As Integer, sex As Integer) As Single
' dist(km), time (seconds), sex=0 female, sex=1 male
' return percentage (fraction) of AGE-WR speed
  PAWRSpeed = (dist * 1000 / time) / dist_arecs(dist, age, sex)
End Function

Function PWWRSpeed(dist As Single, time As Single, age As Integer, sex As Integer) As Single
' dist(km), time (seconds), sex=0 female, sex=1 male
' return percentage (fraction) of WAVA AGE-WR speed
  PWWRSpeed = (dist * 1000 / time) / dist_wrecs(dist, age, sex)
End Function

Function DiffSecs(time1 As Variant, time2 As Variant)
  DiffSecs = (Minute(time1) - Minute(time2)) * 60 + _
              Second(time1) - Second(time2)
End Function