Untitled

function fun(m, ind, n, par)
  integer, intent(in) :: m
  integer, intent(in) :: ind(6), n
  double precision, intent(inout), optional :: par(*)
  integer :: i
  double precision :: b(32)

  b = [2.999314767377287d0, 2.993844170297569d0, 2.982962913144534d0, 2.9667902132486006d0, 2.945503262094184d0, &
       2.9193352839727122d0, 2.8885729807284855d0, 2.853553390593274d0, 2.814660195524919d0, 2.7723195175075137d0, &
       2.7269952498697734d0, 2.6791839747726502d0, 2.6294095225512604d0, 2.5782172325201156d0, 2.526167978121472d0, &
       2.473832021878528d0, 2.4217827674798844d0, 2.3705904774487396d0, 2.3208160252273498d0, 2.2730047501302266d0, &
       2.2276804824924867d0, 2.185339804475081d0, 2.146446609406726d0, 2.1114270192715145d0, 2.080664716027288d0, &
       2.054496737905816d0, 2.033209786751399d0, 2.0170370868554657d0, 2.0061558297024313d0, 2.000685232622713d0]

  fun = exp((b(ind(1)) - b(ind(4)))**2 + (b(ind(2)) - b(ind(5)))**2 + (b(ind(3)) - b(ind(6)))**2)
end function fun