cycling/ran_cone.m

   1 function rc = ran_cone(w)
   2 %makes cone structure
   3 close all
   4 f = @(x,y,x0,y0,m) m*sqrt((x-x0).^2+(y-y0).^2)+1.5*cos(4*sqrt((x-x0).^2+(y-y0).^2))-1;
   5
   6 [p,q] = size(w.fxlong);
   7 tx = linspace(-10,10,p);
   8 ty = linspace(-10,10,q);
   9 [x,y] = meshgrid(tx,ty);
  10
  11 x0 = 2*randn;
  12 y0 = 4*randn;
  13 z0 = 0;
  14 z_max = 20;
  15 radius = 9+2*rand;
  16
  17 r = sqrt(x0^2+y0^2);
  18 % m <= (z_max-z0)/(9-r)
  19 %m = (1/2+rand)*(z_max-z0)/(radius-r);
  20 m = (1+2*rand)*(z_max-z0)/(radius-r);
  21 %m = 1;
  22
  23 z = f(x,y,x0,y0,m);
  24 z(z>z_max) = z_max;
  25
  26 n = round(1/2*z_max);
  27
  28 base = rand;
  29 new_x0 = x0;
  30 new_y0 = y0;
  31 theta = 2*pi*rand;
  32 for i = 1:n
  33     %vertext of new cone
  34 %     z0 = i/n;
  35 %     z0 = 1/4*z_max*rand;
  36     z0 = z0+rand;
  37     if mod(i,2)
  38         new_x0 = (new_x0+rand)*cos(theta);
  39         new_y0 = (new_y0+rand)*sin(theta);
  40     else
  41         new_x0 = 2*i/n*randn;
  42         new_y0 = 2*i/n*randn;
  43     end
  44
  45     r = sqrt(new_x0^2+new_y0^2);
  46     new_m =(1+4*rand)*(z_max-z0)/(radius+2*randn-r)+base+rand;
  47     c(i,1) = new_x0;
  48     c(i,2) = new_y0;
  49     z2 = f(x,y,new_x0,new_y0,new_m);
  50
  51 %     noise = 10/z_max*randn(p,q);
  52 %     z2 = z2+noise;
  53
  54     alpha = 3/4;
  55     z2(z2>z_max) = z_max;
  56     %figure,mesh(x,y,z),hold on,mesh(x,y,z2),hold off;
  57     %z(z2>z0) = max(z(z2>z0),z2(z2>z0));
  58     z(z2>z) = alpha*z(z2>z)+(1-alpha)*z2(z2>z);
  59     %z = imgaussfilt(z,1/2);
  60     %z = smooth_up(z);
  61     z(z>z_max) = z_max;
  62
  63     %figure(234),contour(x,y,z,20,'k')
  64
  65 end  %for loop
  66
  67 %z = smooth_up(z);
  68 %figure,mesh(x,y,z);
  69 %quick_mesh(z)
  70 %figure,plot(c(:,1),c(:,2));
  71 %figure,contour(x,y,z,20,'k');
  72
  73
  74 rc = z;
  75
  76 end % function