#Arkaitz Rodas Bilbao
#© Derechos reservados excepto normas de uso y distribución de la UPV-EHU.

#--------------------------------------------------------------------------------------------------------
#En este programa calculamos mediante el metodo de Runge kutta de cuarto orden la evolución temporal de 
#la matriz densidad que describe un sistema fisico de dos niveles, tenemos la formula explicita que determina
#su evolucion temporal.
#--------------------------------------------------------------------------------------------------------

disp('Para comenzar el programa introduzca las condiciones inicales de la matriz densidad');
printf('\n');

printf('Introduzca el punto inicial de muestreo');
printf('\n');
printf('Tinical');
to=input(' = ');

printf('Introduzca el tiempo final del muestreo');
printf('\n');
printf('Tfinal');
tfinal=input(' = ');

printf('Introduzca el número de puntos del muestreo');
printf('\n');
printf('N');
N=input(' = ');

m=1;
while (m<2);
 for i= 1:2
    for j=1:2
	  g=1;
	  while (g<2);
       printf('d(%d',i);
	   printf(',%d)',j);
	   d(i,j)=input(' = ');
	   printf('\n');
	   if (abs(d(i,j))>1);
	    printf('Dato mal introducido');
		printf('\n');
	   else;
	    g=g+1;
	   endif;
	  endwhile;
    endfor;
 endfor;
 
 if (trace(d)!=1);
  printf('Matriz densidad mal introducida');
  printf('\n');
 elseif (d(1,2)!=conj(d(2,1)));
  printf('Matriz densidad mal introducida');
  printf('\n');
 else;
  m=m+1;
 endif;

endwhile;

printf('Ahora introduzca las componentes del hamiltoniano, se considera dividido por (h/2pi) para simplificar los calculos');
printf('\n');

l=1;
while (l<2);

printf('Introduzca si el hamiltoniano es dependiente o independiente del tiempo----1-Dependiente;2-Independiente');
printf('\n');

printf('u');
u=input('=');

 if (u==1);

 printf('Introduzca los valores medidos en tablas hij y sus instantes de medicion en otra tabla tl');
 printf('\n');
 
 tl=dlmread('t.txt');
 h11=dlmread('h11.txt');
 h12=dlmread('h12.txt');
 h22=dlmread('h22.txt');
 
#Conocemos la relacion explicita entre los elementos de la matriz derivada temporal y los iniciales
#se introducen ya en el programa para simplificar sus calculos.

k=0;
t=to;

while (k<N);
 
  H(1,1)= interp1 ( tl, h11, t);
  H(1,2)= interp1 ( tl, h12, t);
  H(2,1)= conj(H(1,2));
  H(2,2)= interp1 ( tl, h22, t);
  
  h=(tfinal-to)/N;

  k1=h*(complex(0,-1))*(H*d-d*H);
  
  t=t+h/2;
  H(1,1)= interp1 ( tl, h11, t);
  H(1,2)= interp1 ( tl, h12, t);
  H(2,1)= conj(H(1,2));
  H(2,2)= interp1 ( tl, h22, t);
  
  k2=h*(complex(0,-1))*(H*(d+k1/2)-(d+k1/2)*H);
  k3=h*(complex(0,-1))*(H*(d+k2/2)-(d+k2/2)*H);
  
  t=t+h/2;
  H(1,1)= interp1 ( tl, h11, t);
  H(1,2)= interp1 ( tl, h12, t);
  H(2,1)= conj(H(1,2));
  H(2,2)= interp1 ( tl, h22, t);
  
  k4=h*(complex(0,-1))*(H*(d+k3)-(d+k3)*H);

  d=d+h*(k1+2*k2+2*k3+k4)/6;
  
  d11(k+1)=d(1,1);
  d12(k+1)=d(1,2);
  d21(k+1)=conj(d12(k+1));
  d22(k+1)=d(2,2);

  k=k+1;

endwhile;

printf('La matriz densidad obtenida para este valor de t es d=');
printf('\n');

disp(d);
 
 l=l+1;
 
 elseif (u==2);

 m=1;
 while (m<2);
  for i = 1:2
     for j= 1:2
      printf('H(%d',i);
	  printf(',%d)',j);
	  H(i,j)=input(' = ');
	  printf('\n');
     endfor;
  endfor;
 
  if (H(1,2)!=conj(H(2,1)));
   printf('Hamiltoniano mal introducido');
   printf('\n');
  else;
   m=m+1;
   endif;

 endwhile;
 l=l+1;
 

k=0;

while (k<N);

  h=(tfinal-to)/N;

  k1=h*(complex(0,-1))*(H*d-d*H);
  k2=h*(complex(0,-1))*(H*(d+k1/2)-(d+k1/2)*H);
  k3=h*(complex(0,-1))*(H*(d+k2/2)-(d+k2/2)*H);
  k4=h*(complex(0,-1))*(H*(d+k3)-(d+k3)*H);

  d=d+h*(k1+2*k2+2*k3+k4)/6;
   
  pu(k+1)=trace(d^2);
  t=to+h;
  tu(k+1)=t;
  
  k=k+1;

endwhile;

printf('La matriz densidad obtenida para este valor de t es d=');
printf('\n');

disp(d);

printf('La evolucion de la pureza es');
printf('\n');

plot(pu,tu),xlabel('t'),ylabel('Pureza');
 
 else;
 printf('Dato mal introducido');
 printf('\n');

 endif;

endwhile; 



disp('Fin del programa')