sobolp<-function(svec,pvec,g,d,r=0,epsi=0){
# source("~/Art/sobolp.R")

z<-matrix(0,length(svec),length(pvec))
for (i in 1:length(svec)){
   beta<-svec[i]
   for (j in 1:length(pvec)){
      p<-pvec[j]
      q<-p*(beta-d/p+d/2)

if ( beta<= max( r/p+d/2, 
                (r-g)/p+d*(3/(2*p)-1/2),
                d/2+(p-2)*g+(p-1)*r ,
                d*(3/(2*p)-1/2) + r/p - g/(p*(p-1)),
                d/2+2*r*(p-1)/p+g*(p-2)/p) +epsi  ) 
z[i,j]<-1

else{

#upper bound

pote1<-((p-1)/p)*(g+r+d/2)/(beta+g)
pote2<-(1/2)*(beta-r-d/2)/(beta+g)

musd<-2*pi^(d/2)/gamma(d/2)
kup<-musd
ku2<-musd
a1<-(p*(q-r)/(p-1)-q-d)/(q+2*g)
a2<-(q+p*(2*g+r)/(p-1)+d)/(q+2*g)
b1<-(2*(q+g-r)-d)/(q+2*g)  
b2<-(2*(g+r)+d)/(q+2*g)
I1<-(2*pi)^(d*(1-p/(p-1)))*kup*beta(a1,a2)/(q+2*g)
I2<-(2*pi)^(-d)*ku2*beta(b1,b2)/(q+2*g)
vakio<-((g+r+d/2)/(beta-r-d/2))^((beta-r-d/2)/(beta+g))*(beta+g)/(g+r+d/2)
yla<-I1^pote1*I2^pote2*vakio

#lower bound

pote3<-(-1/p)*(g+r+d/2)/(beta+g)
pote4<-(-1/2)*(beta-r-d/2)/(beta+g)
kl3<-musd
kl4<-musd
kl5<-musd
c1<-(p*(q+g)/(p-1)-q-p*(g+r)-d)/(q+2*g)
c2<-(q+p*g/(p-1)+p*(g+r)+d)/(q+2*g)
d1<-(2*(q+g)/(p-1)-2*r-d)/(q+2*g)
d2<-(2*g/(p-1)+2*r+d)/(q+2*g)
e1<-((q+g*(2-p))/(p-1)-2*r-d)/(q+2*g)
e2<-(p*g/(p-1)+2*r+d)/(q+2*g)
I3<-(2*pi)^(d*(1-p/(p-1)))*kl3*beta(c1,c2)/(q+2*g)
I4<-(2*pi)^(d*(1-2/(p-1)))*kl4*beta(d1,d2)/(q+2*g)
I5<-(2*pi)^(-d/(p-1))*     kl5*beta(e1,e2)/(q+2*g)
ala<-I3^pote3*I4^pote4*I5

    z[i,j]<-yla/ala
}

  }
}
return(x=svec,y=pvec,z=z)
}