hyperbola.m: Generates hyperbolae

function hy = hyperbola( a,b,minx,maxx,miny,maxy,max_val,mean_val )
% generates hyperbolas
% the parametric function x = a cosh t
% and                     y = b sinh t
% are used, where t is in the interval [-1,1]
% max_val is the value of the highest point in the hyperbola, which
% may also be smaller than mean_val, which is the value of the
% background

if (nargin<1) a=1; end;
if (nargin<2) b=1; end;

mint = -1.6;
maxt = 1.6;

if (nargin<3) minx=a; end;
if (nargin<4) maxx=max(a * cosh( mint ), a * cosh( maxt )); end;
if (nargin<5) miny = b * sinh( mint ); end;
if (nargin<6) maxy = b * sinh( maxt ); end;
if (nargin<7) max_val = 10; end;
if (nargin<8) mean_val = 5; end;

deltax = maxx - minx;
deltay = maxy - miny;
delta_val = max_val - mean_val;

% We want the hyperbolas to be in a frame, that is 10% larger on each side
deltax = deltax * 1.2;
deltay = deltay * 1.2;

% Let's create the matrix containing the image
dimx = 30;
dimy = 40;
steps = 160;
hy = mean_val * ones( dimy, dimx );

for t = (mint:(maxt - mint)/steps:maxt)
  x = a * cosh( t );
  y = b * sinh( t );
  matx = max( 1, abs( round(( x - minx * 0.95 ) / deltax * dimx )));
  maty = round( y / deltay * dimy ) + dimy / 2;
  if (maty <= dimy & maty > 0)
    if (matx <= dimx & matx > 0)
		hy(maty,matx) = mean_val + delta_val * 1.0; 
	end;
    if (matx+1 <= dimx & matx+1 > 0)
		hy(maty,matx+1) = mean_val + delta_val * .8; 
	end;
    if (matx-1 <= dimx & matx-1 > 0)
		hy(maty,matx-1) = mean_val + delta_val * .8; 
	end;
    if (matx+2 <= dimx & matx+2 > 0)
		hy(maty,matx+2) = mean_val + delta_val * .7; 
	end;
    if (matx-2 <= dimx & matx-2 > 0)
		hy(maty,matx-2) = mean_val + delta_val * .7; 
	end;
    if (matx+3 <= dimx & matx+3 > 0)
		hy(maty,matx+3) = mean_val + delta_val * .6; 
	end;
    if (matx-3 <= dimx & matx-3 > 0)
		hy(maty,matx-3) = mean_val + delta_val * .6; 
	end;
    if (matx+4 <= dimx & matx+4 > 0)
		hy(maty,matx+4) = mean_val + delta_val * .5; 
	end;
    if (matx+5 <= dimx & matx+5 > 0)
		hy(maty,matx+5) = mean_val + delta_val * .3; 
	end;
  end;
end;

hy = hy';

% dspim(hy,0)

end;