function fig411

%  plots ray surfaces

clear *
clf

% get(gcf)
set(gcf,'Position', [1889 957 560 420]);

hold on

% set parameters
nx=20;
nz=nx;

a=1;
b0=1;
y0=-1;
z0=log(4+y0);
b1=-1/(4+y0);

x1=linspace(-1,1,nx);
z1=linspace(-1,1.5,nx);
for ix=1:nx
    for iz=1:nz
	y1(ix,iz) = 10 - (y0 + a*x1(ix)^2 + b0*(z1(iz) - z0)^2 + b1*(z1(iz)-z0));
    end;
end;
xm = 0;  zm = 0.25;  ym = 10 - (y0 + a*xm^2 + b0*(zm - z0)^2 + b1*(zm-z0));
c1 = [xm, ym, zm];
n1 = [2*a*xm, 1, 2*b1*(zm-z0)+b1];

a=3;
b0=1;
y0=4;
z0=log(4+y0);
b1=-1/(4+y0);
x2=linspace(-1,1,nx);
z2=linspace(0,2,nx);
for ix=1:nx
    for iz=1:nz
	y2(ix,iz) = 10 - (y0 + a*x2(ix)^2 + b0*(z2(iz) - z0)^2 + b1*(z2(iz)-z0));
    end;
end;
xm = 0;  zm = 1;  ym = 10 - (y0 + a*xm^2 + b0*(zm - z0)^2 + b1*(zm-z0));
c2 = [xm, ym, zm];
n2 = [2*a*xm, 1, 2*b1*(zm-z0)+b1];

nt=40;
t=linspace(0,2.7,nt);
for it=1:nt
	x3(it) = 0;
	z3(it) = (t(it)-1)*c1(3) - (t(it)-2)*c2(3);
	y3(it) = (t(it)-1)*c1(2) - (t(it)-2)*c2(2);
end;

plot3(x3,z3,y3,'k','Linewidth',1)

%grid on
%box on

text(-0.4,2,-2,'\theta({\bfx})=c','FontSize',20)
text(-0.8,0.5,11.1,'\theta({\bfx})=c+t','FontSize',20)

surf(x1,z1,y1)
surf(x2,z2,y2)
view(-160, 16);
colormap gray
axis off

%axis([-1 1 -1 1 -3 0])
%axes('position',[.13  .43  .8  .53])

i1=18; i2=i1+3;
pb=[x3(i1); y3(i1);  z3(i1)];
pa=[x3(i2); y3(i2);  z3(i2)];
[xc,yc,zc,xc0,yc0,zc0]=cone(pa,pb);
surf(xc,zc,yc,'FaceColor',[0 0 0],'LineStyle','none')
surf(xc0,zc0,yc0,'FaceColor',[0 0 0.65],'LineStyle','none')

i1=35; i2=i1+3;
pb=[x3(i1); y3(i1);  z3(i1)];
pa=[x3(i2); y3(i2);  z3(i2)];
[xc,yc,zc,xc0,yc0,zc0]=cone(pa,pb);
surf(xc,zc,yc,'FaceColor',[0 0 0],'LineStyle','none')
surf(xc0,zc0,yc0,'FaceColor',[0 0 0.65],'LineStyle','none')

function [x,y,z,x0,y0,z0]=cone(pa,pb);
p=pa-pb;
h=norm( p );
r=h/30;
theta=acos(p(3)/h);
phi=acos(p(2)/sqrt(p(1)^2+p(2)^2));
R1=[1 0 0; 0 cos(theta) sin(theta); 0 -sin(theta) cos(theta)];
R2=[cos(phi) sin(phi) 0; -sin(phi) cos(phi) 0; 0 0 1];

% set parameters
nr=20;
nt=nr;
s=linspace(0,r,nr);
t=linspace(0,2*pi,nt);

for ir=1:nr
	for it=1:nt
		xq=s(ir)*cos(t(it));
		yq=s(ir)*sin(t(it));
		zq=h*(1-s(ir)/r);
		zq0=0;

		v=[xq;yq;zq];
		v0=[xq;yq;zq0];

		vv=R2*R1*v;
		vv0=R2*R1*v0;

		x(ir,it)=vv(1)+pb(1);
		y(ir,it)=vv(2)+pb(2);
		z(ir,it)=vv(3)+pb(3);

		x0(ir,it)=vv0(1)+pb(1);
		y0(ir,it)=vv0(2)+pb(2);
		z0(ir,it)=vv0(3)+pb(3);

	end;
end;