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


#-------------------------------------------------------------------------------------------------------------
# En este programa se calcula la evolucion de dos sistemas y sus fidelidades relativas en el paso del tiempo.
# El calculo de la evolucion es similar al programa rungekutta.m .
# Suponemos que evolucionan bajo el mismo hamiltoniano.
#La evolucion se considera Markoviana, el hamiltoniano de interacción introducido debe cumplir esta condición.
#El sistema de interés sera un qubit, el entorno se considerara sin evolucion.
#-------------------------------------------------------------------------------------------------------------

disp('Guarde el hamiltoniano de interaccion del sistema global');
printf('\n');

cofr;
delta;
interp;
part;
bloq;
sid;

 y=1;
while (y<2)
 

printf('Introduzca si el hamiltoniano es independiente o dependiente en la imagen de Schrodinger: 1---Dependiente; 2---Independiente');
printf('\n');

 u=input (' = ');

 if(u==1);
 
  printf('Introduzca la dimension del "entorno"');
  printf('\n');

  n=input(' = ');
 
  printf('Debe tener los hamiltonianos medidos en el archivo dispuestos en una sola columna');
  printf('\n');
  
  HT=dlmread ("hamiltoniano.txt");
  
  printf('Introduzca el numero de puntos de muestreo y los tiempos de mustreo en la tabla tl');
  printf('\n');
  
  nt=input(' = ');
  
  tl=dlmread("tl.txt");
  
#---------------------------------------------------------------------
#Comprobamos que el conjunto de operadores es hermitico.
#---------------------------------------------------------------------
  
  for r=1: nt
   for i=1: 2*n
    for j=1: 2*n
     if (HT(2*n*(r-1)+i,j)!=conj(HT(2*n*(r-1)+j,i)));
      printf('Dato mal introducido');
	  printf('\n');
	  stop
     else;
   
     endif;
   endfor; 
  endfor;
 endfor;
 
#------------------------------------------------------------------------------------------
#Generamos un cell array para simplificar los calculos de interpolacion del hamiltoniano.
#------------------------------------------------------------------------------------------
 
 for i=1: 2*n
  for j=1: 2*n
   for r=1: nt
    hx(r)=HT(2*n*(r-1)+i,j);
    ht{i,j}=hx;
   endfor;
  endfor;
 endfor;
  
 
 
 t=0;

 printf('Debe tener guardados los archivos con las matrices densidad iniciales del qubit y el entorno en "matriz1.txt" y "matriz2.txt"');
 printf('\n');
 
 d=dlmread("matriz1.txt");
 
 db=dlmread("matriz2.txt");
 
 d2=f(d,db);
 
 d2f=d2;
  
 printf('Introduzca una nueva matriz densidad  en "mat.txt" con la que obtener la evolucion de la fidelidad');
 printf('\n');
 
 du=dlmread("mat.txt");
 
 du2=f(du,db);
 
 du2f=du2;
  
 #------------------------------------------------------------
 #Ahora empleamos el metodo runge-kutta aplicado al problema.
 #------------------------------------------------------------
 
 printf('Introduzca el instante tf donde se desea obtener la matriz densidad');
 printf('\n');
 
 tf=input(' = ');
 
 printf('Introduzca el numero de puntos de la interaccion');
 printf('\n');
 
 N=input(' = ');
 
 printf('Introduzca el instante inicial');
 printf('\n');
 
 to=input(' = ');
 
 t=to;
 
 h=(tf-to)/N;
 
  e=0;
  yu=0;
  yu2=0;
  er=0;
  er2=0;
  eru=0;
  eru2=0;
  yi=0;
  yi2=0;
  tyu=0;
 
 while (e<N);

#--------------------------------------------------
#Interpolamos el hamiltoniano con una subrutina.
#--------------------------------------------------
 
 H1=fi(ht,tl,t,n,nt);
 Ro=H1*d2f-d2f*H1;
 Rou=H1*du2f-du2f*H1;
 Re=H1*d2-d2*H1;
 Reu=H1*du2-du2*H1;
 Q=fu(Ro,n);
 Qu=fu(Rou,n);
 
#--------------------------------
#Coeficiente de runge-kutta.
#--------------------------------
 
p1=-complex(0,1)*Re;
pu1=-complex(0,1)*Reu;
k1=-fc(ht,H1,d,db,tl,t,to,n,nt)-complex(0,1)*Q;
ku1=-fc(ht,H1,du,db,tl,t,to,n,nt)-complex(0,1)*Qu;

 t=t+h/2;
 
#--------------------------------------------------
#Interpolamos el hamiltoniano con una subrutina.
#--------------------------------------------------
 
 H1=fi(ht,tl,t,n,nt);
 Re=H1*(d2+p1/2)-(d2+p1/2)*H1;
 Reu=H1*(du2+pu1/2)-(du2+pu1/2)*H1;
 
#--------------------------------
#Coeficiente de runge-kutta.
#--------------------------------
 
p2=-complex(0,1)*Re;
pu2=-complex(0,1)*Reu;
k2=-fc(ht,H1,d+k1/2,db,tl,t,to,n,nt)-complex(0,1)*Q;
ku2=-fc(ht,H1,du+ku1/2,db,tl,t,to,n,nt)-complex(0,1)*Qu;
 
Re=H1*(d2+p2/2)-(d2+p2/2)*H1; 
Reu=H1*(du2+pu2/2)-(du2+pu2/2)*H1;

#--------------------------------
#Coeficiente de runge-kutta.
#--------------------------------
 
p3=-complex(0,1)*Re;
pu3=-complex(0,1)*Reu;
k3=-fc(ht,H1,d+k2/2,db,tl,t,to,n,nt)-complex(0,1)*Q;
ku3=-fc(ht,H1,du+ku2/2,db,tl,t,to,n,nt)-complex(0,1)*Qu;
 
 
 t=t+h/2;
 
#--------------------------------------------------
#Interpolamos el hamiltoniano con una subrutina.
#--------------------------------------------------
 
 H1=fi(ht,tl,t,n,nt);
 Re=H1*(d2+p3)-(d2+p3)*H1;
 Reu=H1*(du2+pu3)-(du2+pu3)*H1;
 
#--------------------------------
#Coeficiente de runge-kutta.
#--------------------------------
 
p4=-complex(0,1)*Re;
pu4=-complex(0,1)*Reu;
k4=-fc(ht,H1,d+k3,db,tl,t,to,n,nt)-complex(0,1)*Q;
ku4=-fc(ht,H1,du+ku3,db,tl,t,to,n,nt)-complex(0,1)*Qu;
 
 d=d+h*(k1+2*k2+2*k3+k4)/6;
 d2=d2+h*(p1+2*p2+2*p3+p4)/6;
 
 du=du+h*(ku1+2*ku2+2*ku3+ku4)/6;
 du2=du2+h*(pu1+2*pu2+2*pu3+pu4)/6;
 
 ds2=fu(d2,n);
 dsu2=fu(du2,n);
 
#------------------------------------------------------------------------------------------------------------
#Calculamos en columnas los terminos necesarios para obtener la fidelidad y la pureza en funcion del tiempo.
#------------------------------------------------------------------------------------------------------------

 e=e+1;

 yu(e)=trace(d^2);
 yu2(e)=trace(ds2^2);
 er=sqrtm(d);
 er2=sqrtm(ds2);
 eru=er*du*er;
 eru2=er2*dsu2*er2;
 yi(e)=trace(sqrtm(eru));
 yi2(e)=trace(sqrtm(eru2));
 tyu(e)=t;
 

 
 endwhile;
 
 y=y+1;
 
 printf('Mostramos ahora la evolucion de la fidelidad entre los dos qubits');
 printf('\n');
 
 plot( tyu, yi, tyu, yi2),xlabel('t'),title('Medidas de fidelidad');

printf('Las matrices densidad son');
printf('\n');

disp(d);
printf('\n');

disp(ds2);
 
 elseif (u==2);
 
 printf('Empleamos la aproximacion "rotating wave aproximation"');
 printf('\n');
 
 
 printf('Debe tener la matriz densidad del sistema guardada en "sis.txt"');
 printf('\n');
 
 d=dlmread("sis.txt");

 printf('Introduzca una nueva matriz densidad  en "mat.txt" con la que obtener la evolución de la fidelidad');
 
 du=dlmread("mat.txt");

 printf('El hamiltoniano de la interaccion se descompone en el producto tensorial de matrices en los espacios de estados, se guerdan en "matrixA.txt" ordenadas de menor a mayor w para cada i de 1 hasta el numero total de matrices');
 printf('\n');
 
 printf('Introduzca el numero de matrices A');
 printf('\n');
 
 na=input(' = ');
 
 printf('Introduzca el numero de w');
 printf('\n');
 
 nw=input(' = ');
 
 A=dlmread("matrixA.txt");
 
 disp('Se guardan las matrices que componen tau en "matrixG.txt" y "matrixS.txt" con w de menor a mayor, la matriz G debe ser positiva');
 printf('\n');
 
 G=dlmread("matrixG.txt");
 
 S=dlmread("matrixS.txt");
 
 printf('Introduzca el instante tf donde se desea obtener la matriz densidad');
 printf('\n');
 
 tf=input(' = ');
 
 printf('Introduzca el instante inicial de la interaccion');
 printf('\n');
 
 to=input(' = ');
 
 printf('Introduzca el numero de puntos de la integracion');
 printf('\n');
 
 N=input(' = ');
 
 h=(tf-to)/N;
 
#-----------------------------------------------------------------------------
#Con las matrices calculamos bloques con los que poder operar mas facil.
#-----------------------------------------------------------------------------

#---------------------------------------------------------
#Se calcula un cell array para operar de forma matricial.
#---------------------------------------------------------

  for n=1: na
   for k=1: nw
    for i=1:2
     for j=1:2 
     y(i,j,k,n)=A(i+2*(k-1)+2*(n-1)*nw,j);
	 endfor;
    endfor;
	a{k,n}=y(:,:,k,n);
   endfor;
  endfor;
  
 g=fb(G,nw);
 
 s=fb(S,nw);
 
 
  
 
#---------------------------------
#Obtenemos el hamiltoniano H_LS.
#---------------------------------

HLS=0;

for k=1:nw
 for i=1:na
  for j=1:na
   HLS=HLS+s(i,j,k)*transpose(conj(a{k,i}))*a{k,j};
  endfor;
 endfor;
endfor;
 
 
 t=to;
 
e=0;

yu=0;
er=0;
eru=0;
yi=0;
tyu=0;

while (e<N);
 


#----------------------------------------------
#Obtenemos el disipador que es variable con d.
#----------------------------------------------

D=fD(a,d,g,nw,na);
D2=fD(a,du,g,nw,na);
 
#--------------------------------------------------
#Calculamos ahora los coeficientes de runge-kutta.
#--------------------------------------------------
 
R=HLS*d-d*HLS;
R2=HLS*du-du*HLS;

k1=-complex(0,1)*R+D;
ku1=-complex(0,1)*R2+D2;

t=t+h/2;

#----------------------------------------------------
#En los siguientes pasos solo debemos recalcular D.
#----------------------------------------------------
 
 D=fD(a,d+k1/2,g,nw,na);
 D2=fD(a,du+ku1/2,g,nw,na);
 
#--------------------------------------------------
#Calculamos ahora los coeficientes de runge-kutta.
#--------------------------------------------------
 
R=HLS*(d+k1/2)-(d+k1/2)*HLS;
R2=HLS*(du+ku1/2)-(du+ku1/2)*HLS;

k2=-complex(0,1)*R+D;
ku2=-complex(0,1)*R2+D2;


D=fD(a,d+k2/2,g,nw,na);
D2=fD(a,du+ku2/2,g,nw,na);

#--------------------------------------------------
#Calculamos ahora los coeficientes de runge-kutta.
#--------------------------------------------------
 
R=HLS*(d+k2/2)-(d+k2/2)*HLS;
R2=HLS*(du+ku2/2)-(du+ku2/2)*HLS;

k3=-complex(0,1)*R+D;
ku3=-complex(0,1)*R2+D2;

t=t+h/2;

D=fD(a,d+k3,g,nw,na);
D2=fD(a,du+ku3,g,nw,na);

#--------------------------------------------------
#Calculamos ahora los coeficientes de runge-kutta.
#--------------------------------------------------
 
R=HLS*(d+k3)-(d+k3)*HLS;
R2=HLS*(du+ku3)-(du+ku3)*HLS;


k4=-complex(0,1)*R+D;
ku4=-complex(0,1)*R2+D2;
 
#----------------------------------------------------
#Obtenemos la matriz densidad con los coeficientes.
#----------------------------------------------------

 d=d+h*(k1+2*k2+2*k3+k4)/6;
 du=du+h*(ku1+2*ku2+2*ku3+ku4)/6;
 
  e=e+1;
 
 yu(e)=trace(d^2);
 er=sqrtm(d);
 eru=er*du*er;
 yi(e)=trace(sqrtm(eru));
 tyu(e)=t;
 

 
endwhile;
 
 
 printf('Mostramos ahora la evolución de la fidelidad y la pureza del estado del qubit');
 printf('\n');
 
 plot( tyu, yi, tyu, yu),xlabel('t'),title('Medidas de fidelidad y pureza');
 
 printf('Las matrices densidad tras la evolucion son');
 printf('\n');
 
 disp(d)
 printf('\n');
 
 disp(du)
  
  y=y+1;
  
 else;
 
  printf('Dato mal introducido');
  printf('\n');
 
 
 endif;

endwhile;


 


disp('Fin del programa');