function [t,y,dy] = analysis(log_polar,fovea,t,y0,f,dmin_fovea,a,b,dist,distances,plot_window)
%
global h_mixedmirrors_fig;
%
%Instituto Superior Técnico - Instituto de Sistemas e Robótica - 			             09/2001
%OMNIVIEWS PROJECT    -   SVAVISCA sensor									[ccgd,jasv@isr.ist.utl.pt]
%
%
% This program computes the variation between the theorethical and numerical 
% derivative function h(rho) for both cartesian and log polar cameras.
%
% log_polar = 0/1: conventional/log polar cameras
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dist = ipheight for fovea 						%
% dist = d for retina (radius of cylinder)   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


deriv_method = 0;   	% 0 : diff;
							% 1 : centered derivative
 							% 2 : analytic
                      
                      
eval_axis_conv = [0 0.3 0.97 1.03];
eval_axis_log =  [0 0.4 0.87 1.03];

dt = t(2)-t(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%																											 %
% 									CARTESIAN IMAGE													 %
%																											 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

									%				%
									%	FOVEA 	%
									%  ------	%

if log_polar == 0;										%%conventional camera
   
   if fovea == 1;  %%foveal region - calculate the points of the fovea mirror
      
      % Computation of the fovea mirror points 
      [t,y] = ode45 ('foveapoints', t, y0, [], log_polar, f, dmin_fovea, a, b, dist);
      
      %Computation of the function of the sensor
      [h,d_floor,rho] = init_sensor(t,y,f,a,b,dmin_fovea,log_polar);
      
      %Call the derivation function
      if deriv_method == 2,
         n = length(y);
         dy_fovea = ((t.*(d_floor-t)+ y.*(y-dist-f)) ./ (y.*t + (t-d_floor).*(y-dist-f))) +  sqrt(((t.*(d_floor-t)+ y.*(y-dist-f)) ./ (y.*t + (t-d_floor).*(y-dist-f))).^2 + 1);
			min_dy_fovea = 1;
         max_dy_fovea = n;
      else
         [dy_fovea, min_dy_fovea, max_dy_fovea] = derivation(y, dt, deriv_method);
      end
      dy = dy_fovea; %setting output variable
      
      %Calculus of the derivate of the funcion d_floor = a*rho+dmin_fovea
      [dd_floor, min_dd_floor, max_dd_floor]	= derivation(d_floor, dt, 0);
      [drho, min_drho, max_drho]					= derivation(rho,dt, 0);
      
      
      %%%%%%%Computation of the numerical derivations
      % eq. for the numerical derivation of the equation d_floor = a*rho+dmin_fovea
	   fact = t./(y-f);
   	g = (1-dy_fovea.^2);
	  	d_floorn = t + (y+dist).*(fact.*g + 2*dy_fovea)./(g - 2*dy_fovea.* fact);

		%d_floorn = t + y.*(t.*g + 2*dy_fovea.*(y-dist-f))./((y-dist-f).*g - 2*t.*dy_fovea);

      [dd_floorn, min_dd_floorn, max_dd_floorn] = derivation(d_floorn, dt, 0);
      min_dd_floorn = min_dd_floorn + min_dy_fovea;
      max_dd_floorn = max_dd_floorn - ( length(dy_fovea) - max_dy_fovea );
      
      
      %plot
      h_figure_upright = findobj(h_mixedmirrors_fig, 'Tag', plot_window);
      axes(h_figure_upright);
		cla;
      rdfloorn = (dd_floorn(min_dd_floorn: max_dd_floorn)./drho(min_dd_floorn: max_dd_floorn))/a;
      rdfloort = (dd_floor(min_dd_floor: max_dd_floor)./drho(min_dd_floor: max_dd_floor))/a; 
      rate = rdfloorn./rdfloort((min_dd_floorn: max_dd_floorn));
      axis(eval_axis_conv);
		plot(rho(min_dd_floorn: max_dd_floorn), rate, 'bx-.');
      grid on;
      hold on;
      set(h_figure_upright, 'Tag', plot_window);

      
   else %%fovea == 0 - retinal region - calculate the points of the retina mirror

									%				%
									%	RETINA 	%
									%  ------	%
      
      % Computation of the retina mirror points 
      [t,y] = ode45 ('retinapoints', t, y0, [], log_polar, f, dmin_fovea, a, b, dist);
      
      %Computation of the function ht(rho)= a*rho+b
      [ht,d_floor,rho] = init_sensor(t,y,f,a,b,dmin_fovea,log_polar);
      
      if deriv_method == 2,
         n = length(y);
         dy_retina = - ((t.*(dist-t)+(y-f).*(y-ht)) ./ ((y-f).*(dist-t)-t.*(y-ht))) + sqrt(((t.*(dist-t)+(y-f).*(y-ht)) ./ ((y-f).*(dist-t)-t.*(y-ht))).^2 + 1);     
         min_dy_retina = 1;
         max_dy_retina = n;
      else
         [dy_retina, min_dy_retina, max_dy_retina] = derivation(y, dt, deriv_method);
      end
      dy = dy_retina; %setting output variable
      
      
      %Calculus of the derivate of the funcion ht(rho) = a*rho+b
      [dht, min_dht, max_dht]			= derivation(ht, dt, 0);
      [drho, min_drho, max_drho] 	= derivation(rho,dt, 0);
      
      
      % eq. for the numerical derivation (eq.6) of the equation ht(rho)
      k = (2*t.*dy_retina) - ( (y-f).*(1-dy_retina.^2) ) ;
      q =  2*(y-f).*dy_retina + t.*(1-dy_retina.^2) ;
      
      %plot
      nn = length(rho);
      h_figure_upright = findobj(h_mixedmirrors_fig, 'Tag', plot_window);
      axes(h_figure_upright);
      %cla;

      for i=1:length(distances),
         di = distances(i); 
         
         hn = y + ( k./q ).*(di-t) ;
         [dhn, min_dhn, max_dhn] = derivation(hn, dt, 0);
         min_dhn = min_dhn + min_dy_retina;
         max_dhn = max_dhn - ( length(dy_retina) - max_dy_retina );
      
      	% rate r - calculus of rn and rt
	      style = ['r+-'; 'go-'; 'k*-'; 'c+-'; 'md-' ];
      
   	   rn = (dhn(min_dhn: max_dhn)./drho(min_dhn: max_dhn))/a*(200/di);
         rt = (dht(min_dht: max_dht)./drho(min_dht: max_dht))/a; 
         rate = rn./rt((min_dhn: max_dhn));
	   	switch i
   		   case 1,
            plot(rho(min_dhn: max_dhn), rate, 'ro-');
            text(rho(nn-1), rate(max_dhn-min_dhn), sprintf('%d cm', di) );
            grid on;
            axis(eval_axis_conv);
            hold on;
      	   plot(rho(min_dht: max_dht), rt, 'b*-.');
         	text(rho(nn-1), rt(max_dhn-min_dhn), 'Rt' );
	      otherwise
   	      plot(rho(min_dhn: max_dhn), rate, style(mod(i-1,5)+1,:)  );
      	   text(rho(nn-1), rate(max_dhn-min_dhn), sprintf('%d cm', di) );
         end  
   	end
     hold on;
     set(h_figure_upright, 'Tag', plot_window);
     
  
   end %%%%(if fovea == 1) % end the fovea and retina mirror for log_polar = 0 %%%%%
		 %%%%%%%%%%%%%%%%%%% start the fovea and retina mirror for log_polar = 1 %%%%	
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%																											 %
% 								LOG-POLAR IMAGE														 %
%																											 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

else %%%log_polar == 1 - log polar camera
   
									%				%
									%	FOVEA 	%
									%  ------	%
   
   if fovea == 1;  %%foveal region - calculate the points of the fovea mirror

   % Computation of the fovea mirror points 
   [t,y] = ode45 ('foveapoints', t, y0, [], log_polar, f, dmin_fovea, a, b, dist);
   
   %Computation of the function of the sensor
   [h,d_floor,rho] = init_sensor(t,y,f,a,b,dmin_fovea,log_polar);
   
   %Call the derivation function
   dt = t(2)-t(1);
   
   if deriv_method == 2,
      n = length(y);
      dy_fovea = ((t.*(d_floor-t)+ y.*(y-dist-f)) ./ (y.*t + (t-d_floor).*(y-dist-f))) +  sqrt(((t.*(d_floor-t)+ y.*(y-dist-f)) ./ (y.*t + (t-d_floor).*(y-dist-f))).^2 + 1);
      min_dy_fovea = 1;
      max_dy_fovea = n;
   else
      [dy_fovea, min_dy_fovea, max_dy_fovea] = derivation(y, dt, deriv_method);
   end
   dy = dy_fovea; %setting output variable
   
   %Calculus of the derivate of the funcion d_floor = a*log_k(rho/rhof)+dmin_fovea
   [dd_floor, min_dd_floor, max_dd_floor]		= derivation(d_floor, dt, 0);
   [drho, min_drho, max_drho] 					= derivation(rho,dt, 0);
   
   %%%%%%%Computation of the numerical derivations
   % eq. for the numerical derivation of the equation d_floor
   fact = t./(y-f);
  	g = (1-dy_fovea.^2);
  	d_floorn = t + (y+dist).*(fact.*g + 2*dy_fovea)./(g - 2*dy_fovea.* fact);

   [dd_floorn, min_dd_floorn, max_dd_floorn] = derivation(d_floorn, dt, 0);
   min_dd_floorn = min_dd_floorn + min_dy_fovea;
   max_dd_floorn = max_dd_floorn - ( length(dy_fovea) - max_dy_fovea );
   
   
   %plot
   h_figure_downright = findobj(h_mixedmirrors_fig, 'Tag', plot_window);
   axes(h_figure_downright);
   %cla;
   rdfloorn = (dd_floorn(min_dd_floorn: max_dd_floorn)./drho(min_dd_floorn: max_dd_floorn))/a;
   rdfloort = (dd_floor(min_dd_floor: max_dd_floor)./drho(min_dd_floor: max_dd_floor))/a; 
   rate = rdfloorn./rdfloort((min_dd_floorn: max_dd_floorn));
   axis(eval_axis_log);
	plot(rho(min_dd_floorn: max_dd_floorn), rate, 'g+-');
   grid on;
   hold on;
   set(h_figure_downright, 'Tag', plot_window);


   else %%fovea == 0 - retinal region - calculate the points of the retina mirror

									%				%
									%	RETINA 	%
                           %  ------	%
   
   % Computation of the retina mirror points 
   [t,y] = ode45 ('retinapoints', t, y0, [], log_polar, f, dmin_fovea, a, b, dist);
   
   %Computation of the function ht(rho)= a*log(rho/rhoc)+b
   [ht,d_floor,rho] = init_sensor(t,y,f,a,b,dmin_fovea,log_polar);
   
   if deriv_method == 2,
      n = length(y);
      dy_retina = - ((t.*(dist-t)+(y-f).*(y-ht)) ./ ((y-f).*(dist-t)-t.*(y-ht))) + sqrt(((t.*(dist-t)+(y-f).*(y-ht)) ./ ((y-f).*(dist-t)-t.*(y-ht))).^2 + 1);     
      min_dy_retina = 1;
      max_dy_retina = n;
   else
      [dy_retina, min_dy_retina, max_dy_retina] = derivation(y, dt, deriv_method);
   end
   dy = dy_retina; %setting output variable
   
   %Calculus of the derivate of the funcion ht(rho) = a*log(rho/rhoc)+b
   [dht, min_dht, max_dht]			= derivation(ht, dt, 0);
   [drho, min_drho, max_drho] 	= derivation(rho,dt, 0);
   
   
   % eq. for the numerical derivation (eq.6) of the equation ht(rho)
   k = (2*t.*dy_retina) - ( (y-f).*(1-dy_retina.^2) ) ;
   q =  2*(y-f).*dy_retina + t.*(1-dy_retina.^2) ;
   
   %plot
   nn = length(rho);
 	h_figure_downright = findobj(h_mixedmirrors_fig, 'Tag', plot_window);
   axes(h_figure_downright);
   %cla;
   
   for i=1:length(distances),
      di = distances(i); 
      
      hn = y + ( k./q ).*(di-t) ;
      [dhn, min_dhn, max_dhn] = derivation(hn, dt, 0);
      min_dhn = min_dhn + min_dy_retina;
      max_dhn = max_dhn - ( length(dy_retina) - max_dy_retina );
   
      % rate r - calculus of rn and rt
      style = ['r+-'; 'go-'; 'k*-'; 'c+-'; 'md-' ];
      
   	rn = (dhn(min_dhn: max_dhn)./drho(min_dhn: max_dhn))/a;
	  	rt = (dht(min_dht: max_dht)./drho(min_dht: max_dht))/a;
	   rate = rn./rt(min_dhn: max_dhn)*(200/di);
      switch i
      	case 1,
         plot(rho(min_dhn: max_dhn), rate, 'ro-');
         text(rho(nn-1), rate(max_dhn-min_dhn), sprintf('%d cm', di) );
         grid on;
         axis(eval_axis_log);
         hold on;
      otherwise
         plot(rho(min_dhn: max_dhn), rate, style(mod(i-1,5)+1,:)  );
         text(rho(nn-1), rate(max_dhn-min_dhn), sprintf('%d cm', di) );
      end
      hold on;
      set(h_figure_downright, 'Tag', plot_window);
   end  %%%%(if fovea == 1)
end

end %%(log_polar)
