C     program main
C     double precision gker,s
C     double precision genKer
C     s =  0.38716867
C     gker = genKer(s,s)
C     write(*,*) gker
C     end
      
C     Note that s1 and t1 should be within [0,1] unit
      subroutine genker(s1,t1,val)
      double precision s1,t1,val
      double precision k1s,k1t,k2s,k2t,st,k1st,k4st

      k1s = s1 - 0.5
      k1t = t1 - 0.5

C     must use 1.0/12 instead of 1/12 to avoid the mistake
      k2s = (k1s**2 - 1.0/12)/2
      k2t = (k1t**2 - 1.0/12)/2
     
C     must use 7.0/240 instead of 7/240
      st = abs(s1-t1)
      k1st = st - 0.5
      k4st = (k1st**4 - k1st**2/2 + 7.0/240)/24

C     genKer = k1s*k1t + k2s*k2t - k4st
C      genker = k2s*k2t - k4st
      val = k2s*k2t - k4st
      end
      




