implicit none 	
	real d,a0,b0,xs,ys,alpha,beta,betas,cvdeg2rad,lplane
	real ay,by,x1,y1,x2,y2,a1,b1,a2,b2,cvrad2deg
	integer i ; character stitle*64
	logical zoom ; data zoom /.false./
	call mvsetflags('Palette di colori',7.0)
	call mvsetflags('Logo',2.0)
	d=700.					! Distanza sorgente 
	alpha=60 ; alpha=cvdeg2rad(alpha)	! Angolo sorgente al centro
	xs=d*cos(alpha) ; ys=d*sin(alpha)	! Coordinate sorgente
	a0=0 ; b0=tan(alpha)			! Parametri linea centrale
	beta= 2   ; 				! Angolo piano riflessione
	lplane=45.0				! Larghezza del piano
	if (zoom) then
	   call plottaframe(-20.0,20.0,-5.0,35.0)
	   write(stitle,'(a,i2,a)')'Zoom view - Plane angle =',nint(beta),' deg'
	else
	   call plottaframe(-100.0,650.0,-50.0,700.0)
	   write(stitle,'(a,i2,a)') 'Plane angle=',nint(beta),' deg'
	endif
	beta=cvdeg2rad(beta) 
	call gslwsc(2.0)
	call gsplci(6)
	call linea(0.0,0.0,xs,ys)
	call gsplci(7)
	do i=-45,45,10
	   call linearparam(float(i),0.0,xs,ys,a1,b1)
	   x1=float(i) ; x2=xs ; call linea(x1,a1+b1*x1,x2,a1+b1*x2)
	enddo
	call linearparam( 5.0,0.0,xs,ys,a2,b2)
	call gsplci(9) ; call gslwsc(5.0)
	ay=0 ; by=tan(beta)
	x1=-lplane ; x2=lplane ; call linea(x1,ay+by*x1,x2,ay+by*x2)
	call displayexample('example43','Bob Nerini Light path',stitle)
	call intersection(a2,b2,ay,by,x1,y1)
	call angleestimation(a2,b2,x1,betas)
	write(6,'(10x,a,f10.7)') '     True angle: ',cvrad2deg(beta)
	write(6,'(10x,a,f10.7)') 'Estimated angle: ',cvrad2deg(betas)
	end

	subroutine angleestimation(a,b,x,beta)	
	implicit none
	real a,b,x,beta,y
	y=a+b*x
	beta=atan(y/x) 
	return
	end

	subroutine intersection(a1,b1,a2,b2,x,y)
	implicit none
	real a1,b1,a2,b2,x,y
	x=-(a1-a2)/(b1-b2)
	y=a1+b1*x
	return
	end

	subroutine linearparam(x1,y1,x2,y2,a,b)
	implicit none
	real x1,y1,x2,y2,a,b
	b=(y1-y2)/(x1-x2)
	a=y1-b*x1
	return
	end

	real function cvrad2deg(alpha)
	cvrad2deg=45.*alpha/atan(1.)
	return
	end

	real function cvdeg2rad(alpha)
	cvdeg2rad=4.*atan(1.)*alpha/180.
	return
	end