function  ob = fcone2(rsX,riX,riY,ha,az1,az2,x0,y0,z0,or,pr,tp)
% function  ob = fcone2(rsX,riX,riY,ha,az1,az2,x0,y0,z0,or,pr,tp)
% rsX = demi-axe supérieur en x
% riX = demi-axe inférieur en x
% riY = demi-axe inférieur en y
% ha = hauteur
% az1 = azimut 1
% az2 = azimut 2
% x0 = position en x
% y0 = position en y
% z0 = position en z
% or = orientation (11=axe X pos;12=axe X nég;21=axe Y pos;22=axe Y nég;
%                   31=axe Z pos;32=axe Z nég)
% pr = precision (standard=50)
% tp = type (1=surf;2=surfl)

[ra,th]=meshgrid(rsX:(riX-rsX)/10:riX,az1:(az2-az1)/pr:az2);
x=ra.*cos(th);
y=ra.*sin(th);
z=ha*(riX-ra)/(riX-rsX);

switch tp
case 1   
   switch or
   case 11
      ob=surf(z+x0,x+y0,riY/riX*y+z0);
   case 12
      ob=surf(-z+x0,x+y0,riY/riX*y+z0);
   case 21
      ob=surf(x+x0,z+y0,riY/riX*y+z0);
   case 22
      ob=surf(x+x0,-z+y0,riY/riX*y+z0);
   case 31
      ob=surf(x+x0,riY/riX*y+y0,z+z0);
   case 32
      ob=surf(x+x0,riY/riX*y+y0,-z+z0);
   end   
case 2
   switch or
   case 11
      ob=surfl(z+x0,x+y0,riY/riX*y+z0);
   case 12
      ob=surfl(-z+x0,x+y0,riY/riX*y+z0);
   case 21
      ob=surfl(x+x0,z+y0,riY/riX*y+z0);
   case 22
      ob=surfl(x+x0,-z+y0,riY/riX*y+z0);
   case 31
      ob=surfl(x+x0,riY/riX*y+y0,z+z0);
   case 32
      ob=surfl(x+x0,riY/riX*y+y0,-z+z0);
   end   
end

%© by Cédric Blaser / 27.11.2001