function [s,T0] = tempy0(t,x,nt,m)
%
% Funcion temperatura. 
% 
% [t,T] = temp(t,x,y,nt,m) 
% 
% Entradas: 
%  --- t:   tiempo para el que se quiere obtener la integral.
%  --- x:   Puntos en x
%  --- y:   Puntos en y
%  --- nt:  numero de puntos para el span de t
%  --- m:   Numero de puntos para la integral en tiempo
%
% Salida, tiempo y campo de temperaturas. 
%
% ------------------------ Programa --------------------------
%
% Metodo de simpson para la integral. 
%

% -------------------- Parametros de la simulacion ---------------------

global q lambda a v_delta delta factor
global tc tp ts 

% ------------------- Fin de parametros  -------------------------------

%              Predimensionalizacin y chequeos
%
if (mod(m,2)~=0) ; m=m+1; end

dimx = length(x);

% Variables y vectores auxiliares temporales.

s     = linspace(1e-7,t,nt);           % Variable temporal 
s_t   = linspace(1e-7,t,(nt-1)*m+1);   % Variable de la integral
nth   = (nt-1)*m+1;

f_ts  = f2t(t-s_t ,tc,tp);  

dd1  = 0.5./sqrt(a*s_t );
dd2  = 14*sqrt(a*s_t);
dd3  = 0.25./(a*s_t);

% Integrales en tiempo 

h = s_t(2)-s_t(1);                
h4 = 4*h/3;
h2 = 2*h/3;

% Predimensionalizacion de T0

T0       = zeros(dimx,nt); 

% --------------------- Integral de T0 --------------------------------%

% evaluamos la integral de T0


oo=h4*ones(1,(nt-1)*m-1);
for i=1:dimx
%   fprintf('%i\n',i)
   f_s0 = eval_fT0(x(i),0,s_t,dd1,dd3); % Ev. del integrando
   for ii=1:nt-1
      coef =[h/3 oo(1:(ii*m-1)) h/3 ];
      coef(3:2:end-2) = h2;
      T0(i,ii+1) = sum(coef.*(f_s0(1:ii*m+1).*...
         f_ts((nt-1-ii)*m+1:end)))'; 
   end
end
T0 = T0*factor;

return

%-----------------------------------------------------------------------%
%                                                                       %
% --------------------------- Subrutinas -------------------------------%
%                                                                       %
% ----------------------------------------------------------------------%

%_______________________________________________________________________%
% Funcion eval_fT0. Evalua la funcion presente en la integral de T_0.   |
%                                                                       |
% Entradas:                                                             |
%         x,y  puntos espaciales                                        |
%         s    puntos para la integral temporal                         |
%         dd1  0.5 ./sqrt(a*s_t);                                       |
%         dd3  0.25./(a*s_t);                                           |
%                                                                       |
% USADA POR:                                                            |
%          T0                                                           |
%_______________________________________________________________________%


function f_s0 = eval_fT0(x,y,s,dd1,dd3);

global q lambda a v_delta delta

% global tc tp ts 


f_s0 = exp(-y.*y.*dd3).*dd1.*...
   (erf((delta-x-v_delta.*s).*dd1)+erf((x+v_delta.*s).*dd1));


%¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ Fin funcion eval_fT0  ¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬

%_____________________________________________________________________%
% Funcion f2t. Funci´on presente en la definicion de funcion de frontera
%  d(x,t)
% Entradas:
%         t        vector Temporal
%         tc,tp    tiempos caracteristicos
%_____________________________________________________________________% 

function y =f2t(t,tc,tp)

n=length(t);
cotsup=ceil(max(t)/tc);
cotinf=max(floor(abs(min(t)-tp)/tc),0);
N = cotinf:cotsup;
y = zeros(1,n);
f = zeros(1,n);
for i=N
   f1 = zeros(1,n);
   f1(i*tc+tp-t>0)=1;
   f2 = zeros(1,n);
   f2(t-i*tc>0) = 1;
   y = y + f1.*f2;
end