C     Note that s1 and t1 should be within [0,1] unit
      double precision function genKer(s1,t1)
      double precision s1,t1
      double precision k1s,k1t,k2s,k2t,st,k1st,k4st
      
      k1s = s1 - (0.5d+0)
      k1t = t1 - (0.5d+0)

C     must use 1.0/12 instead of 1/12 to avoid the mistake
      k2s = (k1s**2 - (1.0d+0)/(12.0d+0))/(2.0d+0)
      k2t = (k1t**2 - (1.0d+0)/(12.0d+0))/(2.0d+0)
     

C     must use 7.0/240 instead of 7/240
      st = abs(s1-t1)
      k1st = st - (0.5d+0)
      k4st = (k1st**4 - k1st**2/(2.0d+0)+(7.0d+0)/(240.0d+0))/(24.0d+0)
 
C     penalize linear term
      genKer = k1s*k1t + k2s*k2t - k4st
C     not penalize linear term
C     genKer = k2s*k2t - k4st
      end
      




