function VF1(p,min,max)

    x = min:.1:max ;
    y = F_Mp(x,p) ;
    plot(x,y,'r') ;
    axis([min max 0 1]) ;
    
end

 
function erg = F_Mp(x,p) 
	erg = [] ;                  
	for i=1:length(x)
		t = x(i) ;
        PIp = quadgk( @(phi) ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(-2/p) ,0,2*pi) ;
        C = quadgk( @(r) r.*exp(-r.^p/p) ,0,Inf) ;
		if( (t==0)||(t<0) )            % Fall t <= 0
            alpha = @(r) atan(-t./(r.^p-(-t).^p).^(1/p)) ;
            integrand = @(r,phi) r.*exp(-r.^p/p)/( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(2/p) ;
			y1 = 1/( C * PIp ) * dblquad( integrand , -2^(1/p)*t , Inf , pi+alpha , 3*pi/2-alpha ) ;
		else                            % Fall t > 0
            beta = @(r) atan(t./(r.^p-t.^p).^(1/p)) ;
            integrand = @(r,phi) r.*exp(-r.^p/p)/( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(2/p) ;
			y1 = 1/( C * PIp ) * ( PIp*quadgk( @(r) r.*exp(-r.^p/p) , 0 , t ) ...
                                 + dblquad( integrand , t , 2^(1/p)*t , pi-beta , 3*pi/2+beta ) + dblquad( integrand , t , 2^(1/p)*t , pi/2-beta , beta ) ...
                                 + dblquad( integrand , 2^(1/p)*t , Inf , pi-beta , 3*pi/2+beta ) ) ;
		end
		erg = [ erg , y1 ] ;
	end
end