// Funciones para análisis de Decaimiento y Equilibrio Radiactivo
//
// desarrollado por: Marcelo M. David B.
// e-mail: marcelod@cin.edu.uy
//
// destinado para: Curso Básico de Metodología de Radioisótopos
//
// Unidad de Instrumentación Nuclear y Computación
// Centro de Investigaciones Nucleares
// Facultad de Ciencias - Universidad de la República
//
// Copyleft, Licencia GNU http://www.gnu.org
//


function [eq,Teq]=equilibrio(Ta,uTa,Tb,uTb,Aa0,Ab0,uGr)
  // Análisis del equilibrio de un sistema generador de Radionucleidos
  // eq=equilibrio(Ta,uTa,Tb,uTb,Aa0,Ab0,uGr)
  // Ta = T1/2 del padre
  // Tb = T1/2 del hijo
  // uTa = unidad del T1/2 del padre ("s","m","h","d")
  // uTb = unidad del T1/2 del hijo  ("s","m","h","d")
  // Aa0 = Actividad inicial del padre
  // Ab0 = Actividad inicial del hijo
  // uGr = unidad de tiempo para graficar
  
  select uGr
  case "s" then
    select uTa
      case "m" then
        Ta=Ta*60;
      case "h" then
        Ta=Ta*60*60;
      case "d" then
        Ta=Ta*60*60*24;
    end
    select uTb
      case "m" then
        Tb=Tb*60;
      case "h" then
        Tb=Tb*60*60;
      case "d" then
        Tb=Tb*60*60*24;
    end
  case "m" then
    select uTa
      case "s" then
        Ta=Ta/60;
      case "h" then
        Ta=Ta*60;
      case "d" then
        Ta=Ta*60*24;
    end
    select uTb
      case "s" then
        Tb=Tb/60;
      case "h" then
        Tb=Tb*60;
      case "d" then
        Tb=Tb*60*24;
    end
  case "h" then
    select uTa
      case "s" then
        Ta=Ta/(60*60);
      case "m" then
        Ta=Ta/60;
      case "d" then
        Ta=Ta*24;
    end
    select uTb
      case "s" then
        Tb=Tb/(60*60);
      case "m" then
        Tb=Tb/60;
      case "d" then
        Tb=Tb*24;
    end
  case "d" then
    select uTa
      case "s" then
        Ta=Ta/(60*60*24);
      case "m" then
        Ta=Ta/(60*24);
      case "h" then
        Ta=Ta/24;
    end
    select uTb
      case "s" then
        Tb=Tb/(60*60*24);
      case "m" then
        Tb=Tb/(60*24);
      case "h" then
        Tb=Tb/24;
    end
  end

	T=max([Ta,Tb]);
	t=[0:T/100:T*10]';
	la=log(2)/Ta;
	lb=log(2)/Tb;
	Aa=Aa0*exp(-la*t);
	Abi=Ab0*exp(-lb*t);
	if Ta==Tb
	 Ta=Tb+0.01;
	end
	Ab=(Ta/(Ta-Tb))*Aa0*(exp(-la*t)-exp(-lb*t))+Abi;
	At=Aa+Ab;
	Ar=Ab-Aa;
	plot2d(t,[Aa Ab At],leg="Padre@Hijo@Total")
	xtitle("Estudio del equilibrio y actividad del sistema","tiempo("+uGr+")","Actividad(Bq)");
	if (Ta<Tb)
		eq="No hay equilibrio";
		Tieq="";
	elseif (Ta>=100*Tb)
		eq="Equilibrio secular";
	else
		eq="Equilibrio transitorio";
	end
	
	if (Ta>=Tb)
	i=5;

 // El equilibrio se alcanza cuando la desinteg del hijo es igual a la del padre. Al dividir se cancelan las exponenciales siendo el resultado una constante
  while (abs(1-(Ab(i)-Ab(i-3))/(Aa(i)-Aa(i-3)))>0.005)
    i=i+1;
  end
	Teq=i*T/100;
	Tieq=["Tiempo para alcanzar el equilibrio: " + string(Teq)+" "+uGr,string(Teq/Tb)+" veces Tb"];
	end
//	ans=x_message([eq,Tieq]);
endfunction

function A=actividad(Tm,uTm,A0,tt,ut,f)
  // Devuelve la actividad remanente después de un tiempo tt
  
  // Uso:
  // A=actividad(Tm,uTm,A0,tt,ut)
  // Tm = T1/2 del Radionucleido
  // uTm = unidades de Tm
  // A0 = actividad inicial del Radionucleido
  // tt = tiempo transcurrido
  // ut = unidades de tt
  // f = bandera de graficar (f=0 no grafica. f=1 grafica)
  
  bTm=Tm;
  select ut
  case "s" then
    select uTm
      case "m" then
        Tm=Tm*60;
      case "h" then
        Tm=Tm*60*60;
      case "d" then
        Tm=Tm*60*60*24;
    end
  case "m" then
    select uTm
      case "s" then
        Tm=Tm/60;
      case "h" then
        Tm=Tm*60;
      case "d" then
        Tm=Tm*60*24;
    end
  case "h" then
    select uTm
      case "s" then
        Tm=Tm/(60*60);
      case "m" then
        Tm=Tm/60;
      case "d" then
        Tm=Tm*24;
    end
  case "d" then
    select uTm
      case "s" then
        Tm=Tm/(60*60*24);
      case "m" then
        Tm=Tm/(60*24);
      case "h" then
        Tm=Tm/24;
    end
  end
  
  t=[0:tt/100:max([Tm*5 tt+tt/10])];
  if (tt<0)
      t=[tt+tt/50:-tt/100:max([Tm*5 -tt-tt/10])];
  end
  lm=log(2)/Tm;
  Ac=A0*exp(-lm*t);
  A=A0*exp(-lm*tt);
  if (f==1)
    plot2d(t,Ac)
    plot(tt,A,'*')
    xtitle("Actividad del Radionucleido de T1/2 = "+string(bTm)+uTm,"Tiempo("+ut+")","Actividad(Bq)")
  end
endfunction

function [At,t]=mezclar(A0,Ta,uTa,B0,Tb,uTb,uT)
  // Genera la mezcla de dos radionucleidos devolviendo la actividad total del sistema
  
  // Uso:
  // A=mezclar(A0,Ta,uTa,B0,Tb,uTb)
  // A0 = Actividad inicial del RN A
  // Ta = T1/2 de RN A
  // uTa = unidades de Tma
  // B0 = Actividad inicial del RN B
  // Tb = T1/2 de RN B
  // uTb = unidades de Tmb
  //
  // A = Actividad de la mezcla
  // uT = unidades del tiempo de la actividad
  
  bTa=Ta;
  bTb=Tb;
 select uT
  case "s" then
    select uTa
      case "m" then
        Ta=Ta*60;
      case "h" then
        Ta=Ta*60*60;
      case "d" then
        Ta=Ta*60*60*24;
    end
    select uTb
      case "m" then
        Tb=Tb*60;
      case "h" then
        Tb=Tb*60*60;
      case "d" then
        Tb=Tb*60*60*24;
    end
  case "m" then
    select uTa
      case "s" then
        Ta=Ta/60;
      case "h" then
        Ta=Ta*60;
      case "d" then
        Ta=Ta*60*24;
    end
    select uTb
      case "s" then
        Tb=Tb/60;
      case "h" then
        Tb=Tb*60;
      case "d" then
        Tb=Tb*60*24;
    end
  case "h" then
    select uTa
      case "s" then
        Ta=Ta/(60*60);
      case "m" then
        Ta=Ta/60;
      case "d" then
        Ta=Ta*24;
    end
    select uTb
      case "s" then
        Tb=Tb/(60*60);
      case "m" then
        Tb=Tb/60;
      case "d" then
        Tb=Tb*24;
    end
  case "d" then
    select uTa
      case "s" then
        Ta=Ta/(60*60*24);
      case "m" then
        Ta=Ta/(60*24);
      case "h" then
        Ta=Ta/24;
    end
    select uTb
      case "s" then
        Tb=Tb/(60*60*24);
      case "m" then
        Tb=Tb/(60*24);
      case "h" then
        Tb=Tb/24;
    end
  end
  
  T=max([Ta,Tb]);
  Te=min([Ta,Tb]);
	t=[0:T/100:T*5]';
	la=log(2)/Ta;
	lb=log(2)/Tb;
  
  A=A0*exp(-la*t);
  B=B0*exp(-lb*t);
  At=A+B;
//  lt=log((At(int(Te*100/T)))/At(int(10*Te*100/T)))/(t(int(10*Te*100/T))-t(int(Te*100/T)));
  plot2d(t,[At A B],leg="Total@A "+string(bTa) + uTa+"@B "+string(bTb) + uTb)
  xtitle("Mezcla de dos RN de T1/2 : " + string(bTa) + uTa + " y " + string(bTb) + uTb,"Tiempo ("+uT+")","Actividad (Bq)");
endfunction

function f=fmc(p,At)
  // función para ser usada en la minimización cuadrática por la función separar
  // p parámetros A0 lA B0 lB
  // At matriz de dos columnas: Actividad y tiempo
  
  m=max(size(At)); 
  f=zeros(m,1);
  for i=1:m
    f(i)=At(i,1) - (p(1)*exp(-At(i,2)*p(2)) + p(3)*exp(-At(i,2)*p(4)));
  end
endfunction

function [Ao,TmA,Bo,TmB]=separar(A,T,uT,Tms,uTms)
  // Determina, si hay una posible mezcla de dos Radionucleidos, el T1/2 del segundo
  
  // Uso:
  // [Ao,TmA,Bo,TmB]=separar(A,T,uT)
  // A : Vector de actividad del sistema en tiempos T
  // T : Vector de tiempo de medidas de actividad
  // uT : unidad de T
  // Tms : Tiempo medio estimado de un RN de la mezcla (Tms=0 indica desconocer el T1/2 estimado)
  // uTms : unidades de Tms
  // Ao : Actividad inicial de RN A
  // TmA : T1/2 de RN A en uT
  // Bo : Actividad inicial de RN B
  // TmB : T1/2 de RN B en uT
  
if ((length(A)<>length(T))|(length(A)<4))
    disp("Los vectores de actividad y de tiempo deben tener la misma longitud y al menos 4 valores")
else
  
  select uT
  case "s" then
    select uTms
      case "m" then
        Tms=Tms*60;
      case "h" then
        Tms=Tms*60*60;
      case "d" then
        Tms=Tms*60*60*24;
    end
  case "m" then
    select uTms
      case "s" then
        Tms=Tms/60;
      case "h" then
        Tms=Tms*60;
      case "d" then
        Tms=Tms*60*24;
    end
  case "h" then
    select uTms
      case "s" then
        Tms=Tms/(60*60);
      case "m" then
        Tms=Tms/60;
      case "d" then
        Tms=Tms*24;
    end
  case "d" then
    select uTms
      case "s" then
        Tms=Tms/(60*60*24);
      case "m" then
        Tms=Tms/(60*24);
      case "h" then
        Tms=Tms/24;
    end
  end
  
  [r e]=size(A);
  if (e>r)
    A=A';
  end
  [r e]=size(T);
  if (e>r)
    T=T';
  end
  [r e]=min(T)
  T=T-r;
  At=[A,T];
  p(1)=A(e)/2;
  p(2)=0;
  if (Tms>0)
    p(2)= log(2)/Tms;  
  end
  p(3)=A(e)/2;
  p(4)=log(A(1)/A(length(A)))/(T(1)-T(length(A)));
  [A p]=leastsq(list(fmc,At),p);
  Ao=p(1);
  TmA=log(2)/p(2);
  Bo=p(3);
  TmB=log(2)/p(4);
    
end

endfunction

function teff=biologico(Ao,Ta,uTa,Tb,uTb,uteff)
  // Determina el T1/2 efectivo tomando en consideración el T1/2 radiactivo y biológico
  
  // Uso:
  // Ao : Actividad Inicial inyectada
  // Ta : T1/2 radiactivo
  // uTa : unidades de TmA
  // Tb : T1/2 biológico
  // uTb : unidades de Tmb
  // uteff : unidades de teff
  // teff : T1/2 efectivo
 
  bTa=Ta;
  bTb=Tb;
  select uteff
  case "s" then
    select uTa
      case "m" then
        Ta=Ta*60;
      case "h" then
        Ta=Ta*60*60;
      case "d" then
        Ta=Ta*60*60*24;
    end
    select uTb
      case "m" then
        Tb=Tb*60;
      case "h" then
        Tb=Tb*60*60;
      case "d" then
        Tb=Tb*60*60*24;
    end
  case "m" then
    select uTa
      case "s" then
        Ta=Ta/60;
      case "h" then
        Ta=Ta*60;
      case "d" then
        Ta=Ta*60*24;
    end
    select uTb
      case "s" then
        Tb=Tb/60;
      case "h" then
        Tb=Tb*60;
      case "d" then
        Tb=Tb*60*24;
    end
  case "h" then
    select uTa
      case "s" then
        Ta=Ta/(60*60);
      case "m" then
        Ta=Ta/60;
      case "d" then
        Ta=Ta*24;
    end
    select uTb
      case "s" then
        Tb=Tb/(60*60);
      case "m" then
        Tb=Tb/60;
      case "d" then
        Tb=Tb*24;
    end
  case "d" then
    select uTa
      case "s" then
        Ta=Ta/(60*60*24);
      case "m" then
        Ta=Ta/(60*24);
      case "h" then
        Ta=Ta/24;
    end
    select uTb
      case "s" then
        Tb=Tb/(60*60*24);
      case "m" then
        Tb=Tb/(60*24);
      case "h" then
        Tb=Tb/24;
    end
  end


  leff = log(2)/Ta + log(2)/Tb
  T=max([Ta Tb]);
  t=[0:T/100:T*7];
  A=Ao*exp(-leff*t);
  Arn=Ao*exp(-t*log(2)/Ta)
  teff=log(2)/leff;
  plot2d(t,[A;Arn]',leg="Actividad efectiva@Actividad RN")
  xtitle("Estudio de Actividad efectiva de RN de T1/2 : " + string(bTa) + uTa + " y T1/2 biol:" + string(bTb) + uTb,"Tiempo ("+uteff+")","Actividad (Bq)");
endfunction