code: def greenhall_sz(t, F, alpha, d):
      
      if d == 1:
          a = 2*greenhall_sx(t, F, alpha)
          b = greenhall_sx(t-1.0, F, alpha)
          c = greenhall_sx(t+1.0, F, alpha)
          return a-b-c
      elif d == 2:
          a = 6*greenhall_sx(t, F, alpha)
          b = 4*greenhall_sx(t-1.0, F, alpha)
          c = 4*greenhall_sx(t+1.0, F, alpha)
          dd = greenhall_sx(t-2.0, F, alpha)
          e = greenhall_sx(t+2.0, F, alpha)
          return a-b-c+dd+e
      elif d == 3:
          a = 20.0*greenhall_sx(t, F, alpha)
          b = 15.0*greenhall_sx(t-1.0, F, alpha)
          c = 15.0*greenhall_sx(t+1.0, F, alpha)
          dd = 6.0*greenhall_sx(t-2.0, F, alpha)
          e = 6.0*greenhall_sx(t+2.0, F, alpha)
          f = greenhall_sx(t-3.0, F, alpha)
          g = greenhall_sx(t+3.0, F, alpha)
          return a-b-c+dd+e-f-g
  
      assert(0)
docs: Eqn (9) from Greenhall2004