a = 2.2
r = 1.0
e = 1.0
U(x,e,r,a) = x<r-.39*a ? 81.*e - 563.*(x-r+.39*a): e*(a/(x-r+a) )**9 - e*(a/(x-r+a) )**3 - 1/3.**9 + 1/3.**3;
Up(x,e,r,a) = x<r-.39*a ? 81.*e - 563.*(x-r+.39*a): e*(a/(x-r+a) )**9 + e*(a/(x-r+a) )**3 - 1/3.**9 + 1/3.**3;
f(x,e,r,a) = x<r-.39*a ? 563. : e*9.*(a/(x-r+a) )**10/a - e*3.*(a/(x-r+a) )**4/a;
fp(x,e,r,a) = x<r-.39*a ? 563. : e*9.*(a/(x-r+a) )**10/a + e*3.*(a/(x-r+a) )**4/a;
U1(x,e,r)  = x<r-.9   ? 28.98*e - 100*(x-r+.9): e*(r/x)**9 - e*(r/x)**3;
f1(x,e,r)  = x<r-.9   ? 100.: .9*e*(r/x)**10/r - .3*e*(r/x)**4/r;
E(r,e,a) =  r < a ? 0. : (r-0.4)*(94.6*e+100.*r-40.) - 50.*(r-0.4)**2 \
	- e*a**9/(8.*(2.*r+a)**8) + e*a**9/(8.*(-.4+a)**8) \
	+ e*a**3/(2.*(2.*r+a)**2) - e*a**3/(2.*(-.4+a)**2) \
	+ e*( a**9/(2.*r+a)**9 - a**3/(2.*r+a)**3 ) * (2.*r+0.4);
E2(r,e,a) = - e*a**9/(8.*(2.*r+a)**8) + e*a**9/(8.*(a)**8) \
        + e*a**3/(2.*(2.*r+a)**2) - e*a**3/(2.*(a)**2) \
        + e*( a**9/(2.*r+a)**9 - a**3/(2.*r+a)**3 ) * (2.*r);
E0(r,e) = e*r*( 1./8. - 1./2. + 1./18. - 1./(3**8*8) ) + e*(1./3.**8 - 1./3.**6) * 2.*r;
set term post
set output 'Integrale.ps'
plot [0:3][] E0(x,1.0), E2(x,1.,1.), E2(x,1.,2.2)
set output 'Potenziale.ps'
plot [0:6][-2:6] U(x,1.0,1.0,1.0), U1(x,1.0,1.0)
set term x11
#a=2.4;e=1.0;r=1.0;
#81.0 - 563.0*($t-1.0+0.39*2.4)
#2641.8/($t+1.4)^9 - 13.8/($t+1.4)^3 + 0.0370
#9.0*2641.8/($t+1.4)^10 - 3.0*13.8/($t+1.4)^4
	
