function neu(p,min,max)

%     warning off
% 
%     load('4SimulierteDaten100.mat')
%     
%     figure(1)
%     ecdf( MAX ) ;
    
    %figure(2)
    x = min:5:max ;
    y = F_Mp(x,p) 
    %plot(x,y,'r') ;
    %axis([min max 0 1]) ;
    
end

 
function erg = F_Mp(x,p)                % exakte Verteilungsfunktion
	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
			y1 = 1/( C * PIp ) * dblquad( @(r,phi) ( (pi+atan(-t./(r.^p-(-t).^p).^(1/p)) <= phi) | (phi >= 3*pi/2-atan(-t./(r.^p-(-t).^p).^(1/p))) ) * ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(-2/p) .*r.*exp(-r.^p/p) ,-2^(1/p)*t,Inf,0,2*pi) ;
		else                            % Fall t > 0
			y1 = 1/( C * PIp ) * ( dblquad( @(r,phi) ( (pi-atan(t./(r.^p-t.^p).^(1/p)) <= phi) | (phi <= 3/2*pi+atan(t./(r.^p-t.^p).^(1/p))) ) * ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(-2/p) .*r.*exp(-r.^p/p) ,2^(1/p)*t,Inf,0,2*pi) ...
                + dblquad( @(r,phi) ( (pi-atan(t./(r.^p-t.^p).^(1/p)) <= phi) | (phi <= 3/2*pi+atan(t./(r.^p-t.^p).^(1/p))) ) * ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(-2/p) + ( (pi/2-atan(t./(r.^p-t.^p).^(1/p)) <= phi) | (phi <= atan(t./(r.^p-t.^p).^(1/p))) ) * ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(-2/p) .*r.*exp(-r.^p/p) ,t,2^(1/p)*t,0,2*pi) ...
                + quadgk( @(r) r.*exp(-r.^p/p) ,0,t) ) ;
		end
		erg = [ erg , y1 ] ;
	end
end

% function erg = DAF(p,alpha,beta)		% Durchschnittsanteilsfunktion für p > 1
% 	erg = 1/(quadgk( @(phi) ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(1/p).^(-2) ,0,2*pi)) * quadgk( @(phi) ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(1/p).^(-2) ,alpha,beta) ;
% end
% 
% function erg = N(p,phi)					% Normierung für p-verallgemeinerte Sinus/Kosinus
% 	erg = ( abs(sin(phi)).^p + abs(cos(phi)).^p ).^(1/p) ;
% end
% 
% function erg = PI_p(p)					% p-verallgemeinerte Kreiszahl
%     erg = 1/2*quadgk( @(phi)N(p,phi).^(-2) ,0,2*pi) ;
% end
% 
% function erg = g(x,p)					% Dichtegenerierende Funktion
% 	erg = exp(-x/p) ;
% end
%